沙漠腹地綠洲植被覆蓋度提取及植被指數優選

王 寧周明通魏 宣郭玉川

〔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]中也可以看出,考慮了土壤背景影響的土壤調整型植被指數的估算結果與當地植被覆蓋狀況的空間解譯結果高度契合。下一步將利用超低空遙感數據源,結合機器學習中的多種回歸算法,建立起樣方植被覆蓋度與可見光、多光譜植被指數的回歸模型,優選能夠反演研究區地面植被覆蓋狀況的最佳模型。

4.2 結論

(1)利用Image J軟件可以有效提高樣方植被覆蓋度解譯效率,平均總體精度OA達到88.67%,植被覆蓋度提取精度較高。提取大量樣方的植被覆蓋度,對于傳統的提取方式而言,Image J有較好的適用性。

(2)基于無人機航拍獲取植被覆蓋度,進而建立植被指數—覆蓋度相關關系來獲取植被覆蓋變化的研究方法是可行的。本研究旨在重要變量因素的篩選,即對最優植被指數的選擇,而不涉及地面植被覆蓋狀況的預測模型構建與模型優選,未來需引入多種機器學習回歸模型進行優選,提升反演精度。

(3)在極端干旱區綠洲生態系統中,土壤調整型植被指數(SAVI,MSAVI)能夠更好地反演天然綠洲植被覆蓋變化狀況。在不同環境下的荒漠地區是否能得到相同結果需要進一步的數據分析與驗證。