沙漠腹地绿洲植被覆盖度提取及植被指数优选

2022-02-19宁周明通魏宣郭玉川

水土保持通报 2022年6期

王宁周明通魏宣郭玉川

〔1.新疆大学地理与遥感科学学院,新疆乌鲁木齐830017;2.绿洲生态教育部重点实验室(新疆大学),新疆乌鲁木齐830017〕

天然植被是区域内生态与环境因素的综合反映,能够维系生态环境各项功能的持续发挥[1]。在极端干旱区,对天然植被生长状况开展野外调查工作量大且相对困难,基于遥感数据获取的植被指数来反演地表植被覆盖状况,是开展地面植被研究的重要方法之一。植被指数(Ⅵ)能够通过特定的表达方式定量描述植被生长活力,而针对不同的应用环境,植被指数又衍生出了多种表现形式[2]。因此,如何选择植被指数研究不同环境下植被覆盖问题受到广泛关注。近年来,国内外学者对稀疏植被覆盖状况反演开展了一些研究。叶静芸等[3]利用地面调查数据,结合Quick Bird影像,提取样方地上植被生物量与植被指数进行回归分析,认为RVI估算乌兰布和沙漠东北缘的荒漠—绿洲过渡带地上生物量效果较好。崔万新等[4]利用无人机可见光波段荒漠植被的光谱特性,反演多种可见光植被指数,结果表明可见光波段差异植被指数、过绿减过红指数、归一化绿蓝差异指数提取荒漠灌木植被精度较好。部分学者[5-6]对不同植被指数模型与实测数据进行线性拟合,选取了当地最适宜的植被指数。唐亮等[7]通过无人机低空遥感获得样方植被覆盖度,与像元二分模型的反演结果(FVCM)进行验证,结果表明FVCM的精确度达到了83.1%。杨红艳等[8]针对内蒙古草原分类方法问题,使用无人机高光谱数据源对草原植被进行机器学习算法分类,随机森林分类精度较高。郝梦宇等[9]利用无人机影像对植被样方进行识别与分类,快速获取植被信息数据。部分学者[10-12]在各自研究区利用卫星、无人机传感器数据反演的植被指数进行分析,以期更好地指示植被生物量的时空变化及反演精度。国外学者对地面植被覆盖状况估算开展了众多研究。Kattenborn等[13]将无人机航拍影像与Sentinel-1和Sentinel-2数据进行空间升尺度匹配,以建立植被覆盖度的估算模型。无人机遥感平台集成多类型传感器,使得小尺度范围内地面植被监测成为可能。有学者[14-15]利用无人机影像获得的多光谱植被指数、可见光植被指数[16]、植被结构及布局[17],在各自研究区内对地面植被进行监测与估算,估测准确度较高。

综上所述,众多学者在不同尺度和不同场景下优选植被指数反演覆盖度方面进行深入研究,基于无人机航拍获取的覆盖度,在大幅减少样地调查工作量的同时可以尽可能降低植被样方的抽样误差,以无人机航拍影像提取的植被覆盖度为基准可以更接近实际植被覆盖情况。将无人机遥感作为数据源提取地面植被信息,相关研究主要集中在植被生长状况较好的区域,沙漠戈壁地区稀疏植被有关报道较少[18]。为此,本研究依托无人机数据提取的植被覆盖度,筛选极端干旱区指示植被覆盖变化的最优植被指数,旨在提高反演植被覆盖变化的精度,为研究维护绿洲生态安全,揭示该地区荒漠—绿洲共生关系提供参考。

1 研究区概况

克里雅河由南向北经过于田绿洲流入塔克拉玛干沙漠腹地,在尾闾处形成达里雅布依天然绿洲(图1),绿洲位于38°16′—38°33′N,81°44′—82°5′E,地势较为平坦,呈叶状分布,进入绿洲后干流由西南向东北和西北方向分支。根据于田县气象站资料,研究区多年平均气温11.76℃,多年平均降水量50.23 mm,多浮尘天气。绿洲内植被群落类型单一,主要以乔木胡杨(Populuseuphratica)、灌木柽柳(Tamarix ramosissima)以及草本植物芦苇(Phragmites australis)为主[19]。

2 材料与方法

2.1 研究思路

本研究选取塔克拉玛干沙漠腹地的克里雅河尾闾绿洲达里雅布依作为研究区,以Image J解译出来的无人机航拍影像样方的植被覆盖度(FVC)为基础,选择Sentinel-2B多光谱影像波段反演7种干旱区常用植被指数,分别与FVC进行Spearman相关分析,根据显著性筛选结果,在R 4.1.1环境中建立逐步多元回归、岭回归、套索回归模型。在最佳回归模型内,利用标准回归系数、确定系数评价指标,确定反演研究区植被覆盖状况的最优植被指数。利用不同年份的航拍数据对本研究结果进行验证。

2.2 数据获取与处理

2.2.1 无人机航拍影像在天气晴朗无风条件下,采用大疆创新科技有限公司(DJI)生产的P4_Multispectral四旋翼无人机平台,搭载了用于可见光(RGB)成像通道的CMOS影像传感器,针对研究区植被覆盖状况开展低空拍摄,航拍时间为2021年8月23—26日,植被样地分布如图1所示。样地航拍位置选择遵循生态学野外调查方法与实地情况,在研究区按上中下游梯度均匀分布,每幅航拍样地范围是1 km×0.3 km,基本覆盖了河床、沙地、灌木丛、芦苇地、林地等地类,能够较好地代表当地植被覆盖实际状况。飞行作业时,飞行器航高设置为100 m,航向重叠率为80%,旁向重叠率为70%,地面分辨率约5 cm,航拍样地总面积约4.5 km2。借助大疆智图软件对无人机多幅可见光图像进行辐射校正及拼接处理,生成15幅样地高清正射影像(digital orthophoto map,DOM)。

图1 研究区遥感影像及各样地位置

2.2.2 植被覆盖度提取本研究以Sentinel-2B影像10 m像元尺度为基准,将无人机航拍样地以10 m×10 m的正方形为单位布设若干样方,为避免样方间距离过近而产生同质性,各样方间隔20 m,以SP5样地为例绘制了样方分布示意图(图2a,2b),Sentinel-2B影像的像元尺度与无人机样方范围保持高度一致(图2c,2d)。

图2 样方分布示意图(以SP5样地为例)

植被样方提取过程在ArcGIS 10.8中利用创建渔网、创建要素类、裁剪等工具完成,在15块样地内共提取到了5 000个植被样方。

样方布设时,遵循与Sentinel-2B卫星像元高度重合的原则,以达到卫星像元植被指数值即为对应样方植被指数值的目的(如图3a所示)。根据已有研究,利用无人机航拍影像在荒漠区划分植被与非植被[20],代替人工调查植被覆盖度是可行的[21],但主要以颜色指数划分阈值来判断植被,增加了区分植被的复杂性。Image J(Fiji)是基于JAVA语言的开源图像处理软件(http:∥imagej.nih.gov/ij/),早期常用于细胞荧光面积、伤口划痕面积计算等生物医学领域,支持图像批量处理。研究区内优势种以乔木胡杨和灌木柽柳为主,俯视形态上植被轮廓较为清晰,与细胞、划痕边缘相似,故本研究利用Image J测量各样方内植被覆盖面积。计算植被覆盖面积占样方总面积的百分比,从而得到各样方的FVC,根据野外实地植被调查,结合生态环境质量评价标准,将研究区植被覆盖度进行划分[22]:劣覆盖度(0%～20%)、低覆盖度(20%～30%)、中覆盖度(30%～60%)、高覆盖度(60%～100%),不同覆盖状况下的解译结果如图3b—3e所示。

图3 样方布设及不同覆盖状况下分类结果

2.2.3 Sentinel-2B卫星遥感影像综合对比常用遥感数据的时空分辨率、应用广泛性以及数据源获取等方面,选择Sentinel-2数据进行常用植被指数的反演,其时间分辨率为5 d,空间分辨率可达10 m。为保证与无人机航拍时间空间的同步性和可比性,下载数据的卫星过顶时间为当地2021年8月19日11:16,数据级别为Level-1C级影像数据,下载自欧洲航空局哥白尼数据中心(https:∥scihub.copernicus.eu/dhus/#/home)。该期幅影像云量少,成像质量高,与无人机航拍时间接近。Level-1C级数据已经过辐射校正与几何精校正处理,只需利用SentiNel Application Platform(SNAP)软件中Sen2Cor插件进行大气校正处理即可。

2.2.4 植被指数选取通过文献梳理及研究区实际,本研究共选用归一化植被指数(NDVI)、差值植被指数(DVI)、比值植被指数(RVI)、土壤调整植被指数(SAVI)、改进型土壤调节植被指数(MSAVI)、全球环境监测指数(GEMI)、大气阻抗植被指数(ARVI)等7种干旱区常用植被指数,其计算公式详见表1。各样方对应的植被指数使用SentiNel Application Platform(SNAP)软件中Vegetation Radiometric Indices功能进行计算。

表1 主要植被指数计算公式

2.3 研究方法

2.3.1 植被覆盖度解译精度评价为客观地掌握Image J对植被覆盖度的提取精度,本研究采用总体分类精度(OA)评价指标对提取结果进行精度评价。由于本研究所使用的无人机航拍影像分辨率是cm级的,所以将人工目视解译识别结果作为真实感兴趣区,与Image J提取结果计算混淆矩阵,得到总体分类精度OA,计算公式为:

式中:Q为植被样方总数;Qij为矩阵中第i行第j列上的频数;k为类别数,仅有植被与裸地之分,k=2。

2.3.2 Spearman相关分析法在进行相关分析前需判断变量的正态性分布。经Shapiro-Wilk检验法(公式2)检验,顺序统计量(W)虽接近1,但各变量参照系数(P)均小于0.05(表2),所以原始数据中变量并不服从正态分布[30]。较Pearson相关分析法而言,Spearman相关分析法不需要各变量服从正态分布,其适用范围更广[31]。本文选用Spearman相关分析法获取样方的植被覆盖度与各植被指数的线性关系。

表2 Shapiro-Wilk正态分布检验结果

式中:Rw为统计数值;ai为常量;Zi为第i个样本;z为样本均值。

相关系数r的大小可以衡量变量间的相关程度,相关系数r的取值范围是-1～1,正值则为正相关,相反则为负相关,相关系数r越大,即认为变量间的关系越紧密。其计算公式为:

式中:n为样本数量;Xi变量为各植被指数对应值;Yi变量为无人机航拍样方提取的植被覆盖度;x,y分别对应变量的平均值。

通过Spearman方法选择显著变量,同时判断自变量间是否存在多重共线性,为植被覆盖度与不同植被指数间的回归分析模型选择提供依据。

2.3.3 植被覆盖度估算模型及验证目前,利用植被指数反演天然植被覆盖度模型主要有线性回归和非线性回归两类[32]。由于自变量间存在较强的共线性,故在R 4.1.1环境中将5 000组数据中随机选取的4 000组不同植被指数值作为自变量,与对应FVC引入逐步多元回归分析模型(MSR)[33]、岭回归分析模型(RR)[34]以及套索回归分析模型(LR)[35],削弱自变量间的多重共线性。

将其余1 000组数据用于模型精度验证,计算决定系数(R2)、希尔不等系数(TIC)、均方根误差(RMSE)、平均绝对误差(MAE),检验回归模型反演精度[36],其计算公式为:

式中:N为参与模型验证的样本数;Ym为第m个样方的实际植被覆盖度值;ÿm为回归模型反演的第m个样方植被覆盖度预测值;y为地上实际植被覆盖度的平均值。TIC值越小,R2值越大,表明模型拟合效果越好。RMSE与MAE能够很好评估数据变化,在相同模拟集和预测集条件下,RMSE与MAE的值越小,表示构建的回归分析模型精度越高。

2.3.4 植被指数优选评价指标研究采用标准回归系数与确定系数两种评价指标来确定最优植被指数。首先,标准化回归系数在一定程度上可以衡量自变量的相对重要程度,其绝对值大小可以说明多元回归模型中自变量的相对重要性[37]。其次,运用LMG方法[38],将决定系数R2自然分解为每个自变量对因变量的贡献程度,并且将自变量的独自贡献量和与其他自变量共同贡献量考虑在内,可以衡量自变量的相对重要程度[39],评价指标计算在R 4.1.1环境中完成。

2.3.5 结果验证本研究采用不同年份航拍数据对结果进行验证。航拍数据为2018年10月8日使用大疆Phantom 4Pro飞行器在研究区内进行植被样地航拍的无人机影像(图4,图内方框为布设的植被样方);卫星数据为2018年10月4日的Sentinel-2B遥感影像,云量少,成图质量较高。数据处理和分析过程与上文步骤一致,在样地内共提取到了660个植被样方。

图4 验证样地及样方分布示意图

3 结果与分析

3.1 植被覆盖度提取精度

本研究采用总体精度OA作为达里雅布依天然绿洲植被覆盖度提取的评价指标。在R 4.1.1环境中,使用caret包的createDataPartition函数进行分段抽样,从每个植被样地中随机抽选100个植被样方进行精度验证,结果详见表3。15个样地平均总体精度为88.67%,各样地分类精度分布在81%～96%,表明Image J可以较高精度提取样方植被覆盖度,能够用于后续无人机航拍样方提取植被覆盖度的处理工作。

表3 植被覆盖度提取精度

将Image J部分提取植被覆盖度有误的样方参照目视解译结果进行修正,修正后的各样地的植被覆盖度统计结果如图5所示。劣、低、中、高覆盖度样方分别占总样方的55.44%,13.49%,22.24%,8.83%,各样地内不同覆盖度样方均有分布。

图5 各样地植被覆盖度柱状堆积图

3.2 不同植被指数与植被覆盖度的相关性分析

将样方修正后的FVC与各植被指数进行Spearman相关性分析(图6)。由图6可以看出,在0.001水平(双侧)上,FVC与NDVI,DVI,RVI,SAVI,MSAVI,GEMI,ARVI植被指数间相关系数范围为0.66～0.72,为中度相关,且相关性均显著,说明使用无人机影像反演的FVC与各植被指数间建立估算模型是可行的。而各样方NDVI,DVI,RVI,SAVI,MSAVI,GEMI,ARVI均呈现正相关关系,且植被指数的相关系数超过0.9,自变量间表现出较强的多重共线性问题。MSR,RR,LR3种回归模型对变量间的多重共线性问题有较好的适用性,可以进行变量剔除。因此,本研究采用MSR,RR,LR回归模型估算研究区植被覆盖度。

图6 FVC与各植被指数相关分析图

3.3 FVC与各植被指数回归分析及模型验证

基于以上植被覆盖度与各植被指数的Spearman相关性分析结果,在R 4.1.1环境中,以FVC为因变量,NDVI,DVI,RVI,SAVI,MSAVI,GEMI,ARVI为自变量建立MSR,RR,LR回归模型。首先,利用caret包中的preProcess函数将原始数据进行标准化处理。使用step()功能建立逐步回归模型,通过AIC指标进行变量筛选,lm()重新拟合MSR回归模型。利用glmnet函数包建立RR回归和LR回归模型,通过lambda指标进行参数筛选,拟合出RR与LR回归模型。为了进一步验证已建立的MSR,RR与LR回归模型的可靠性和适应性,将1 000组检验数据代入模型中进行回判计算分析,通过因变量FVC的模型预测值与实际值的对比,以检验MSR,RR与LR回归模型的预测效果。同时利用验证数据进行结论验证,计算结果详见表4。

表4 植被覆盖度估算精度交叉验证

可以看出,由MSR所建立的回归模型的R2,TIC,RMSE,MAE在训练集与测试集上均优于由LR,RR所建立的回归模型。由分析可知,MSR,LR,RR由训练集与测试集所拟合的曲线变化较一致,拟合度相对较低,考虑到自变量共线问题,泛化能力相对较差。在验证数据中,MSR模型的训练集R2,TIC,RMSE,MAE分别为0.41,0.51,0.88%和0.63%,测试集上的R2,TIC,RMSE,MAE分别为0.40,0.56,0.92%和0.67%,均优于验证数据中的LR,RR模型。因此,已建立的MSR逐步回归模型具有一定的可靠性,较好地估算沙漠腹地达里雅布依天然绿洲的植被覆盖状况,故本研究采用MSR逐步回归模型来进行下一步的植被指数优选。

3.4 模型最优变量选择

在MSR逐步回归模型中,采用标准回归系数和确定系数两种评价指标进行最优植被指数筛选。由表5可知,在试验和验证数据中,SAVI,MSAVI两种植被指数的标准回归系数的绝对值远高于其他植被指数,确定系数R2均小于1。而试验数据中自变量间确定系数变化值范围为0.14～0.2,验证数据中自变量间确定系数变化值范围为0.03～0.37,确定系数结果在同一数量级内,差异并不显著,故以标准回归系数所得结果作为参考得出结论,在逐步多元回归模型中,自变量SAVI,MSAVI对因变量FVC变化影响程度最为显著,说明反演极端干旱区植被覆盖状况的最优植被指数为土壤调整型植被指数(SAVI,MSAVI)。

表5 最优变量评价表

4 讨论与结论

4.1 讨论

荒漠植被分布稀疏,地面信息提取难度相对较大。利用目视解译等的传统样方植被覆盖度获取方法耗时耗力,而基于无人机的航拍植被样方图像,经过Image J软件的批量提取,可以大幅提高植被覆盖度的获取效率,但仍可能存在植被与裸地间植被覆盖的提取混淆问题。参考本研究图像处理过程中出现的问题,太阳照射下植物在地面的投影冗长,加重在裸地上的投影像素,而造成样方植被覆盖度提取混淆。在无人机航拍样地时,应选择当地正午时间或者日出时进行拍摄,以减小样方植被覆盖度提取误差。在植被覆盖度提取过程中,还可以结合无人机航拍影像高精度的特点,利用稀疏植被纹理、形态、光谱等特征,提高沙漠地区植被信息的提取精度。

依据干旱、半干旱区地区常用植被指数及Sentinel-2B卫星影像数据特征,初步优选了NDVI,DVI,RVI,SAVI,MSAVI,GEMI,ARVI共7种较为典型的干旱区植被指数,是否有更合适反演干旱区天然绿洲植被覆盖状况的植被指数需要做进一步的研究。近年来,机器学习(machine learning,ML)十分流行,在数据处理上高效便捷,已经广泛应用到多个领域中。利用常见的随机森林、支持向量机等回归算法,建立起植被覆盖度与植被指数的关系,进一步提高回归模型的拟合精度。

本研究在Spearman相关性分析的基础上,运用逐步多元回归分析法得出结论,土壤调整型植被指数(SAVI,MSAVI)能够更好指示极端干旱区天然绿洲植被覆盖状况。从以往文献[40-44]中也可以看出,考虑了土壤背景影响的土壤调整型植被指数的估算结果与当地植被覆盖状况的空间解译结果高度契合。下一步将利用超低空遥感数据源,结合机器学习中的多种回归算法,建立起样方植被覆盖度与可见光、多光谱植被指数的回归模型,优选能够反演研究区地面植被覆盖状况的最佳模型。