基于无人机影像的阿拉善东南部荒漠灌丛植被空间格局研究

程天亮,王新平,马雄忠,潘颜霞

1 中国科学院西北生态环境资源研究院, 兰州 730000 2 中国科学院大学, 北京 100049 3 南宁师范大学地理科学与规划学院, 南宁 530001

干旱荒漠区的低盖度植被常呈现出斑块状空间格局[1],其形成、结构和演替研究,对揭示区域生态系统变化的关键过程有重要意义[2]。灌丛斑块是荒漠中较为稳定的植被存在形式,在防风固沙和沙漠化逆转中发挥着不可替代的作用[3],但气候条件变化和人为活动会导致其格局改变,影响土壤水分、养分等资源的空间分配[4],造成荒漠植被生态系统稳定性的扰动。国外有大量基于空间格局指示生态系统退化程度的研究[5—8],通过引入无标度格局对植被系统的整体性和复杂性进行了良好描述[9],并指出斑块大小空间分布越偏离幂律,则标志着植被系统越接近退化状态[10],但国内的相关研究很少。植物间局部关系对植被斑块的形成有着重要影响[11],关于种群相互作用的研究虽然较为丰富[12—14],但种群空间分布及结构变化与植被空间异质性的关系还不清楚。植被的空间格局是与所处环境长期互馈中形成的,是生态过程和水文过程相互作用的结果,深入探究植被斑块的相关机理,并将其空间配置应用于受损植被恢复中,有着至关重要的生态学意义[15]。

阿拉善荒漠区的极端干旱和风蚀环境对植物的生长和繁殖造成极大压力,植被群落类型少,空间异质性强[16]。内蒙古自治区政府已在当地建立多处自然保护区,并将其纳入生态红线。因此,采用恰当的指标评估该地区生态系统的状态,并以稳定的天然植被格局为参考,恢复和重建退化植被,对维护生态安全屏障意义重大。本文通过无人机影像获取覆沙荒漠灌丛植被的空间结构信息,统计分析斑块大小的空间分布来评估植被系统的状态,结合土壤条件和下垫面粗糙度佐证估计结果,用点格局分析法分析灌木种群的内在关系,重点关注以下问题：(1)用幂律拟合斑块大小空间分布评估生态系统状态的准确性；(2)不同灌木种在灌丛斑块形成中发挥了什么样的作用？以期为生态脆弱区植被的可持续管理提供科学依据。

1 材料与方法

1.1 研究区概况

研究区域位于阿拉善左旗巴润别立镇,腾格里沙漠东南部(图1),属典型大陆干燥性气候,日照时间长,昼夜温差大,风大沙多,干旱少雨,年均降水量约80—150mm,年均气温约8.2℃,年均蒸发量2258.8mm,年均风速4.2m/s。景观类型属覆沙荒漠,土壤为风沙土,植被以半灌木和沙生灌木为优势种,包括柠条锦鸡儿(Caraganakorshinskii)、沙冬青(Ammopiptanthusmongolicus)、蒙古扁桃(Prunusmongolica)、霸王(Zygophyllumxanthoxylon),猫头刺(Oxytropisaciphylla)、驼绒藜(Ceratoideslatens)等[17]。

1.2 研究方法

1.2.1无人机航拍样方植被

2020年5月中旬开展植被调查,用样线拉设30m×30m样方4个,使用大疆Phantom 4 Pro航拍获取样方正射影像,航高20m,航向和旁向重叠80%,样方影像经Photoscan软件拼接后分辨率约0.1cm[18],导入ArcMap定义投影坐标系,通过目视解译和实地验看,识别灌木种类并绘制灌丛斑块,算出各斑块面积,在野外实测灌丛冠幅的长短轴,按椭圆近似计算灌丛斑块面积,与ArcMap中计算的灌丛面积进行配对t检验差异不显著[19](P>0.05) 。计算灌丛植被盖度=总冠幅/样方面积×100%,灌木种密度=某灌木种的丛数/样方面积、灌木种分盖度=某灌木种的总冠幅/样方面积×100%、灌木种频度=某灌木种的丛数/样方内总灌木丛数×100%。

图1 研究区域Fig.1 Study area 图中1—4代表样方

1.2.2灌丛斑块大小分布

根据Kéfi对干旱区植被系统状态的划分[20],用最大似然法拟合以下先验分布：

p(x)=C×e-βx(指数分布)

式中C是常数,α、β为指数,xmin为尾部阈值,xmax为截尾值,μ为平均值,σ为标准差。参考Meloni等[21]的方法,首先估算幂律分布的参数,xmin值由Kolmogorov-Smirnov统计量的最小值确定[22],结合Berdugo等[23]引入的幂律区间(Power Law Rang, PLR)来判断植被斑块拟合纯幂律的优度,幂律区间的计算方法为：

式中xmin为尾部阈值,xsmallest为最小斑块,xlargest为最大斑块。若PLR大于0.57,则原假设为幂律分布为最优拟合,而其他分布为备择假设,此时所有备择假设都要在幂律的最优区间上进行对比,即都以xmin为阈值,通过赤池信息准则(Akaike information criterion, AIC)确定最优拟合[24]。若PLR小于0.57,则认为幂律分布并非最优拟合,将对数正态分布作为原假设,此时其他先验分布不需要设定阈值,而在整个样本区间上进行对比。当幂律分布或截尾幂律为最优拟合时,则说明灌丛植被由大量小型斑块和逐渐减少的大型斑块组成,空间格局无标度或在较大范围内无标度,植被的稳健性较强,处于良好状态；若对数正态分布或指数分布为最优拟合,则意味着灌丛斑块存在特征长度,前者表明有少量大型斑块,但不存在极大斑块,后者表明以小型斑块为主,几乎没有大型斑块,以上两种分布寓意植被的稳健性较差,处于退化状态。采用R语言poweRlaw包进行分析[25]。

1.2.3土壤水分、机械组成测定

选取样方1和样方2在植被样方内以3m为间隔,按照方形网格取土样,每个取样点用100cm3环刀取5—10cm沙土,105℃烘干法测定土壤水分。在环刀取样点周围取3次0—10cm沙土混匀后过2mm筛,采用吸管吸湿法测定土壤机械组成,土壤粒径分级标准采用美国制,计算土壤分形维数：

式中R是土壤机械组成各分组的平均半径,RL是最大组平均半径,等式左边是半径小于某分组的质量百分数,D为分形维度[26]；每个样方取土样162个。

1.2.4表观粗糙度长度测定

在植被样方内用风杯风速计测定不同高度的风速,风杯高度设为1m、2m、3m,以6m间距按方形网格测定12个点[27],每个点测定30分钟[28]。用梯度法计算各点空气动力学粗糙度：

式中z0为空气动力学粗糙度,z1和z2为为两个风杯的高度,v1为z1处测定的风速,v2为z2处测定的风速。为了获得较为稳定的空气动力学粗糙度值,按照不同高度风速比值标准差最小来确定最佳高度组合[29],z1取2m,z2取3m。空气动力学粗糙度表征着在地表以上风速降低为零的高度,在研究范围内取平均值称为表观粗糙度长度[30],能够反映植被与地表组成的下垫面减弱风沙运动的作用[31]。

1.2.5扩展点格局分析

由于灌丛植被盖度低,斑块大小与研究尺度相近,故采用Wiegand等提出的扩展点格局分析法分析灌丛植被间的相互关系[32—33]。扩展点格局基于栅格数据对尺寸有限、形状不规则的植物进行格局分析,用分类的栅格来表示植物,在配对相关函数g(r)的基础上提出圆环二阶统计量O(r)：

O(r)=λ×g(r)

式中λ为考察植被的密度,r为研究尺度,O(r)表示以任意一点为圆心、半径为r、宽为dr的圆环内考察植物点的期望值[14]。

不区分种以灌丛斑块为单变量的格局分析零模型采用完全空间随机模型(complete spatial randomness, CSR),单个灌木种的格局分析除猫头刺外也都采用CSR,而猫头刺由于空间格局具有明显的异质性,因此采用异质泊松过程(heterogeneous poisson process, HPP)；灌木种之间的格局关系即双变量格局分析的零模型采用环形转换(toroidal shift, TS)。用蒙特卡洛法模拟植被格局 199次,将模拟结果和观察结果合并,排序后分别剔除最大和最小的5次构成零模型95%的置信区间。当实际值在置信区间内时,灌丛空间格局为随机分布,大于上限时,O(r)>λ表明在特定尺度上,研究对象的出现频率更高,为聚集分布,小于下限时,O(r)<λ表示植物间存在竞争关系,呈均匀分布[34—35]。通过把ArcMap中勾勒的灌丛矢量文件转换为栅格格式,导入Programita 2018进行分析[36]。

2 结果与分析

2.1 灌丛斑块空间格局及状态评估

植被样方取样点见图1,样方内地形较为平缓,灌丛盖度在5%—20%以内(见表1),各样方中仅零星草本植物,驼绒藜、木蓼和矮锦鸡儿的生长状况较差,样方1和样方2都以柠条锦鸡儿+猫头刺为主要植被类型,样方2中沙冬青在风沙掩埋作用下,分枝向外侧伸展,冠幅较大且能与柠条锦鸡儿相互镶嵌形成大斑块,样方3以霸王+蒙古扁桃为主,且两者常重叠生长,样方4基本为纯的猫头刺,灌丛密度最高。样方灌丛斑块分布情况见图2。

各样方中灌丛斑块分别为93、125、107、372丛,从箱线图可看出(图3),斑块大小分布都严重负偏,中位数与下四分位数非常接近,且尾部异常值较多,属于重尾分布。表2幂律分布拟合结果中样方2的参数值α为1.378与经典文献相近[7](α=1.23),其次是样方3。样方2的阈值最小,尾部样本量多于其他样方,说明样方2斑块分布拟合幂律的结果较优,而样方1和样方4在阈值影响下,样本损失量在80%以上。有且仅有样方2的PLR大于0.57,可将幂律为其斑块大小分布的原假设,而其他样方的幂律区间都小于0.57。基于以上结果,样方2各备择假设都以0.021为阈值进行对比,而样方1、3和4则是以全部样本量拟合对数正态分布、幂律分布、指数截尾的幂律及指数分布。

表1 阿拉善左旗样方概况

图2 样方灌丛斑块分布图Fig.2 Patchiness of shrub in quadrats

图3 灌丛斑块大小箱线图Fig.3 Boxplot of shrub patch size

表2 灌丛斑块大小拟合幂律分布参数表

由表3模型拟合的结果可知,样方2中AIC最小的是指数截尾幂律,而其他三个样方AIC值最小的都是对数正态分布。样方2最优拟合的幂指数为1.078,阈值为0.021m2,截尾值为10.053m2,说明该嵌套函数在区间[0.021,10.053]上服从幂律分布,斑块大小分布无标度。而样方1、3和4的最优拟合为对数生态分布,表示其斑块大小具有特征长度,根据对数生态分布期望的计算公式得出,各样方斑块平均值分别为2.918、0.654、0.361m2,且方差依次减小。图4展示出四种模型对数据的描绘方式各有不同,直观地看,样方1中截尾幂律对尾部数据的描绘最为贴合,然而赤池信息准则不支持其为最优拟合,这是因为对数正态分布对占据较大样本量的小斑块典型特征描述更为精确,同时也能兼顾大斑块的存在(如样方3和4所示)。样方2舍弃部分样本后,在大斑块急速减少前,双对数轴上斑块大小及其频率基本呈线性关系。以上结果表明,样方2灌丛斑块空间格局优于其他三个样方,植被系统稳健性较强,处于良好的状态。

表3 灌丛斑块大小分布模型拟合对比

2.2 不同空间格局灌丛植被生境条件的对比

由于生长季前降水稀少,两样方表层土壤含水量极低,如所示,经独立t检验样方1土壤表层含水量极显著(P<0.01)的低于样方2,说明在长期干旱条件下,样方2的初始土壤水分高于样方1。两样方的土壤中细沙(0.25—0.05mm)含量均为最高(见表4),极粗沙(2.0—1.0mm)含量都低于1%(未列出),样方1中粘粉粒含量极显著(P<0.01)的低于样方2,样方2的土壤分形维度极显著(P<0.01)的高于样方1,说明样方2的土壤质地更细腻,样方1中16份土壤样本分形维度拟合系数的t值不显著,说明土壤微粒不具自相似性,发育情况较差[37]。样方2的表观粗糙度显著(P<0.05)大于样方1,植被与裸地构成的下垫面更为粗糙,削弱风速的能力更强。

图4 灌丛斑块大小分布模型拟合Fig.4 Model fitting of shrub patch size distribution

由图5可知,样方1中10cm土壤水分与土壤沙粒(粗沙加细沙)、粘粉粒都不相关,只与分形维数显著正相关,相关系数为0.26,而样方2中表层土壤水分与土壤沙粒极显著的负相关,与粘粉粒和分形维数都极显著的正相关,土壤水分与其他指标的相关系数在0.42—0.44之间。从两样方的生境条件调查结果来看,样方2抗风蚀能力强,土壤结构更好,更具保水性。这与空间格局对植被系统状态的评估结果一致,说明较优的植被空间格局有利于缓解外界环境压力。

图5 土壤机械组成与含水量相关图Fig.5 Correlogram of Particle size distribution and water content Sand:沙粒;Slit&clay:粘粉粒;D:分形维度;Water:10cm土壤含水量;暖色代表负相关,冷色代表正相关,椭圆上*表示P<0.05相关关系显著,**表示P<0.01相关关系显著

图6 灌丛斑块点格局分析Fig.6 Point pattern analysis of shrub patch 图中红色线条为观测的O(r)值,灰色条带为模拟的95%置信区间;猫头刺空间尺度只关注0—5m

2.3 灌木种内及种间的相互关系

以灌木种为单变量的点格局分析结果表明,样方1中柠条在0—2.5m、4.3—5.4m、5.5—5.6m尺度上处于均匀分布,样方2中在0—1.9 m尺度上为均匀分布(见图6)；样方2和3中沙冬青分别在0—2.8m和0—2.5m尺度上呈均匀分布；样方3中蒙古扁桃在0—1 m尺度上为均匀分布,在1.7m、2.9—3.3m上呈聚集分布,霸王在0—1.5m、8.3—8.7m和9—9.1m上呈均匀分布；在≤10m的其他范围内,以上灌木种都呈随机分布。图2可看到猫头刺具有异质分布的特点,其单变量点格局分析采用核密度估计参数设置为5m 的HPP零模型,主要研究在5m范围内猫头刺的空间格局,其中O(r)值在0—0.2m(样方1)、0—0.3m(样方2)和0—0.5m(样方4)尺度上小于置信区间,样方2中在0.7—1.6m及2.1m尺度上大于置信区间,样方4中2.8—4.8m尺度上大于置信区间,除此之外,在其他≤5m的范围内三个样方中猫头刺的空间格局都符合零模型,说明由于自身的扩繁机制,猫头刺在小尺度上(≤5m)容易形成高密度的种群格局,这种聚集趋势在较纯的群落中更为突出,但在微小尺度上依然以竞争关系为主,导致其不能形成连片的群体。驼绒藜、木蓼和矮锦鸡儿由于个体少、盖度小,未做进一步分析。

双变量的点格局分析结果表明,柠条锦鸡儿和猫头刺在样方1中0—2.3m尺度上存在竞争关系,5.6—6.4m尺度上存在互利关系,样方2中在0—1.5m尺度上存在竞争关系(见图6)。沙冬青和柠条锦鸡儿在4.4—6.3m尺度上存在竞争关系,样方3中沙冬青和蒙古扁桃在0.6—1.8m尺度上存在互利关系,沙冬青和霸王在0.5—1.3m尺度上为互利关系。蒙古扁桃和霸王在0.2—1.6m尺度上存在互利关系。各灌木种在≤10m的其他尺度上关系不明确。

3 讨论

3.1 灌丛斑块空间格局对植被生态系统状态的指示

2007年Kéfi等[7]和Scanlon等[9]分别在地中海和非洲中南部干旱区通过样线横截面实地测量和IKNOS影像(4m分辨率)非监督分类的方法,证明了干旱区植被斑块大小空间分布符合幂律,从而将无标度概念引入到植被空间格局中。根据幂律分布的性质,当幂指数1<α≤2时,其平均值和方差发散,不能算出具体的特征长度,Kéfi等[7]拟合的幂指数都在区间(1,2]上,Scanlon等[9]拟合的6个站点中有5个在该区间内,本文中样方2拟合纯幂律指数为1.374,截尾幂律的指数为1.048,都表征着无标度的空间格局。近年来应用遥感影像进行大尺度植被空间格局的研究中,国外大部分采用谷歌地球[21, 23, 38],因其空间分辨率较高,能识别小型植被斑块及其间隙(如2017年Berdugo等[23]影像分辨率<0.3m)。随着无人机测绘技术的成熟,本文通过合成高重叠率的低空航拍影像,分辨率<0.01m,在小尺度上能够提取更多的斑块信息。在选取遥感与实地验证的尺度时,Xu等[38]影像研究区域为250m×250m,实地调查60个0.5m×0.5m样方,Meloni等[21]影像研究区域最大为1000 m2,田间设80m2样圆5个,Berdugo等[23]田间调查采用30m×30m样方,为了增大样本量、提高统计模型的可信度而将影像面积扩大到50m×50m。本文影像和实地均为30m×30m,样方中灌丛斑块数量最少的接近100丛,符合Clauset等[22]模型对比的样本量要求。Kéfi等[20]在2011年基于渗透理论以无标度空间格局作为生态系统稳健性的指示信号,当外界环境压力增大导致斑块格局破碎化时,植被斑块大小空间分布逐步偏离幂律,因此该指示信号最初还用于估计生态系统到退化临界点的距离。2019年Sankaran等[8]探讨了植被生态系统渐变或突变的内在驱动力后,认为幂律及截尾幂律能够指示植被处于良好状态,而不能用于评价其抵抗力稳定性的强弱。Meloni等[21]指出对数正态分布是除幂律分布外能对干旱区植被斑块进行良好描述的空间模型之一,符合该分布的植被系统其形成受局地随机因子影响,稳健性较差。本文的四个样方中仅样方2符合截尾幂律,而其他样方都符合对数正态分布,说明样方2植被系统的状态较好,灌丛生境条件的对比分析也佐证了空间格局的评估结果。万炜等通过Landsat等影像资料评判该区域沙漠化过程[39],认为在长时间尺度上有先恶化后逆转的趋势。但本文结果则表明在沙漠—荒漠草原过渡的局部地区,灌丛植被系统状态较差,灌丛沙漠化的程度还较为严重。

3.2 灌木种内和种间关系对灌丛空间格局的影响

小尺度植被的局部相互关系是大尺度空间格局形成的重要驱动力[9]。从本文点格局分析的结果看,灌丛斑块的空间格局同时受到种内和种间关系的影响。所有灌木种种群内部在小尺度上都存在显著的竞争关系,表现为自疏过程[35],这可能是导致该地区灌丛斑块破碎化的主要原因。猫头刺在小尺度上符合异质泊松分布,可能是因为其株高一般约15cm,种子在风力作用下传播距离有限；同时猫头刺又在0—5m范围内出现了尺度分离现象[40],即随着尺度的增加有聚集趋势,在样方中有大量连在一起的“豆荚”形斑块(见图2),能形成多样化的空间格局。而蒙古扁桃种子个体大,只能靠重力传播,在平缓地区一定尺度上呈聚集分布。Sheffer等[41]认为景观的异质性会随着植物种类发生变化,小型物种(鳞茎早熟禾Poabulbosa)能够形成多样化的空间格局,而大型物种(刺地榆Sarcopoteriumspinosum)则倾向于表现为类似干旱区的低异质性空间格局；Xu等[38]发现植物间的互利关系有利于形成重尾分布的空间格局。本文样方3中霸王在小尺度上与沙冬青和蒙古扁桃都存在互利关系,因其种子外部包裹着一层轻盈的翅,有利于凭借风媒找到“护理”植物。样方2中沙冬青与柠条锦鸡儿在小尺度上相互关系不明确,或者说二者之间互利和竞争关系共存,彼此镶嵌形成大型斑块,图3样方2中10 m2以上的3个离群点都是由这两个灌木种重叠形成的,反而增加了灌丛斑块的空间异质性,可见在平缓沙地上,大型物种间的互利作用更有利于形成符合幂律分布的空间格局(样方2和3的PLR分别等于0.745和0.405)。

植物通过生态相互作用产生正负反馈,对其所处环境进行调节,根据图灵提出的激活抑制原则,这些尺度依赖的正负反馈使植被由无序的初始状态发展为有序的空间格局[42],不同植物种由于生物遗传特性的差异,调节环境的能力也不同。本文发现在阿拉善高原腾格里沙漠-荒漠草原过渡区,大灌木种更有利于维持植被生态系统的稳定性,一是因其独特的旱生特征对环境的适应性强,如在冬春季风沙作用下,常绿灌木沙冬青能形成较大沙堆[43],进而在沙堆保水作用下萌生新枝,缓解了环境压力；柠条锦鸡儿的树冠呈倒锥形,有利于形成树干径流,促进水分入渗[44]。二是大灌木种树冠对微气候的调节能为中小型物种提供护理植物效应,如本文中蒙古扁桃和沙冬青对霸王都存在互利作用。因此,保护和培育大灌木种对该地区植被的恢复和重建有重要价值。在干旱风蚀环境中,物种丰富度高的灌丛植被种内、种间相互关系更为复杂,空间异质性更强,呈现出点状和条带状斑块共存的双稳态植被格局(如样方2所示),灌丛斑块大小的空间分布在一定范围内表现为无标度。

3.3 灌丛斑块空间格局差异的内在机理

本文中四个样方的气候条件相似,而样方2浅层土壤的平均含水量大于样方1,土壤质地也更为细腻。虽然土壤细粒物质对水分的吸附性更强,但两样方土壤机械组成与水分的相关关系并不一致,这与马生花等在荒漠草原上发现的结果相似[45],说明荒漠区土壤水分不仅与土壤质地有关,还受植被的影响,如局部的遮荫效应会改变土壤水分分布[46]。同时,灌丛是荒漠区下垫面粗糙元的主要组成部分,其密度[47]和水平空间配置[48]会显著影响对风力的阻碍作用。在水资源短缺的制约下,优化的空间格局能在较大尺度上通过增大地表粗糙度,减弱风速,富集精细土壤微粒[49],减缓蒸发[50],以保存可利用土壤水分,减轻环境压力。综上,样方2中的灌丛斑块与地表大气和土壤间产生了较好的正反馈,促使其形成了更为复杂的空间格局。但本文空间梯度少,不足以建立下垫面粗糙度与灌丛斑块格局的数量关系,因此我们将定量分析灌丛格局特征与防风固沙的生态学机制[51]作为下一步研究目标。

4 结论

灌丛斑块的空间格局能够为覆沙荒漠生态系统状态提供准确的指示信息。研究区域内灌丛斑块大小空间分布符合截尾幂律的植被样方,生态系统状态优于符合对数正态分布的植被样方。不同灌木种种内和种间的相互关系对空间格局有显著影响,小尺度上竞争关系为主导是斑块破碎化的首要原因。小灌木(如猫头刺)种内的互利作用能促使其形成多样化的斑块形态,而大灌木(如沙冬青和蒙古扁桃)种间的互利关系,能够增大灌丛斑块的异质性,形成更为复杂的空间格局。通过评估荒漠区生态系统状态,分析保持系统稳态性的内在机制,对受损植被的恢复和重建有重要的参考价值。