基于無人機大樣方數(shù)據(jù)及國產(chǎn)衛(wèi)星反演草地植被覆蓋度方法研究

蔡宗磊，苗正紅*，常 雪，劉艷慧，郝 剛，何龍濤

(1. 吉林省水利水電勘測設計研究院，吉林 長春 130021；2. 吉林省建筑科學研究設計院，吉林 長春 130021；3. 東北大學測繪遙感與數(shù)字礦山研究所，遼寧 沈陽 110819；4. 吉林省老龍口水庫管理局，吉林 琿春 133300)

植被覆蓋度是指植被在地面的垂直投影面積占研究區(qū)總面積的百分比[1-4]，能直接反映陸地表面植被的生長情況，也是作為生態(tài)模型、水文模型等在內(nèi)的一系列模型的輸入?yún)?shù)，在生態(tài)系統(tǒng)中起著重要的作用[5-7]。傳統(tǒng)的植被覆蓋度獲取方法一般選用地面測量法，如目測法、概率計算法等，該類方法受人力物力限制，只能獲取有限次數(shù)的觀測數(shù)據(jù)，難以呈現(xiàn)植被覆蓋度在時空上的變異性。遙感技術(shù)的特點為大范圍、實時觀測，可迅速簡便地獲得空間上植被覆蓋的變化，已成為估算區(qū)域植被覆蓋度的重要技術(shù)方法[8]。

像元二分模型、多端元混合像元分解法及回歸模型法等是遙感估算植被覆蓋度應用最廣泛的方法[9]。像元二分法以混合像元分解為基礎，對衛(wèi)星信號中綠色植被與裸土信息做線性分解處理，分解后取得像元中植被信息，對植被覆蓋度進行反演，該方法可去除土壤等背景影響，但受到地物成分及類別的限制，此外要由純植被和裸土像元支撐[10-11]。多端元混合像元分解法可以彌補像元二分模型受兩種地物約束的缺點，在某種程度上能去掉土壤對植被覆蓋度反演的影響，削弱枯枝落葉層作用，在地物類型復雜的區(qū)域具有較高的反演精度，但在端元的選取上比較敏感，受端元類別、數(shù)量限制，易與相關(guān)性較強的地物產(chǎn)生“同譜異物”，從而影響植被覆蓋度精度[12-13]?；貧w模型法的理論基礎是數(shù)理統(tǒng)計，創(chuàng)建遙感數(shù)據(jù)地表反射率與實測植被覆蓋度間的不同數(shù)學回歸模型，進行植被覆蓋度的反演。回歸模型法被普遍使用到提取植被覆蓋度中，且取得了良好效果，但需依靠較多實測數(shù)據(jù)進行植被覆蓋度的反演，普適性相對較差[14-19]。

由于草地下墊面植物群落的物種多樣性，以及生長環(huán)境復雜，應用回歸模型法中的支持向量機回歸估算草地植被覆蓋度，需實測數(shù)據(jù)進行建模，并且對模型進行精度評價。然而，在野外布設的1 m×1 m的小樣方，樣方大小不能與獲取的遙感影像的空間分辨率相匹配，從而影響建模的精度。無人機(Unmanned aerial vehicle，UAV)是一種由無線電遙控及本身程序裝置掌控，易于攜帶，為可以實施多項任務的無人駕駛飛行器[20]。無人機遙感具有靈活性強、可云下飛行、影像分辨率高、時效性強、成本低等眾多優(yōu)勢，可以獲取野外大樣方的數(shù)據(jù)，可作為國產(chǎn)高分衛(wèi)星數(shù)據(jù)估算草地植被覆蓋度的補充[21-22]。國產(chǎn)高分一號衛(wèi)星，可以提供精度高、范圍廣的空間觀測服務，特別是搭載的PMS(Pantone matching system)相機，可以獲取8 m空間分辨率的遙感影像[23]。由于受到放牧、采礦及開墾等人類活動的影響，內(nèi)蒙古自治區(qū)呼倫貝爾市伊敏露天煤礦區(qū)周邊草原具有較高的空間異質(zhì)性，估算其植被覆蓋度，能為礦區(qū)復墾以及持續(xù)發(fā)展提供科學依據(jù)。

本文以伊敏露天煤礦附近的典型草甸草原為研究區(qū)，利用無人機獲取60 m×60 m大樣方數(shù)據(jù)與國產(chǎn)GF-1衛(wèi)星PMS相機8 m空間分辨率數(shù)據(jù)構(gòu)建多種植被指數(shù)作為自變量，應用支持向量機回歸，建立不同數(shù)據(jù)源(地面數(shù)據(jù)—無人機大樣方數(shù)據(jù)，無人機大樣方數(shù)據(jù)—GF-1數(shù)據(jù))的植被覆蓋度估算模型，探討利用無人機大樣方數(shù)據(jù)結(jié)合國產(chǎn)GF-1數(shù)據(jù)反演中國北方草原礦區(qū)周邊草地覆蓋度的適用性與可行性，并且為國產(chǎn)衛(wèi)星在草地生態(tài)系統(tǒng)中的應用提供基礎。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況及樣點選擇

伊敏露天礦地處內(nèi)蒙古自治區(qū)呼倫貝爾草原的中部，自1985年開采，并于1986年對礦區(qū)破壞土地進行了復墾及生態(tài)恢復，生態(tài)恢復后以針茅(Stipacapillata)、沙棘(Hippophaerhamnoides)、羊草(Leymuschinensis)、披堿草(Elymusdahuricus)等為主。其附近典型草甸草原區(qū)域以圍欄草場、退化草場、放牧草場以及灌木為主，植物群落主要分為：豬毛菜(Salsolacollina)、蒺藜(Tribulusterrestris)、灰綠藜(Chenopodiumglaucum)群落；沙生冰草(Agropyrondesertorum)、堿地風毛菊(Saussureajaponica)、蒲公英(TaraAacummongolicum)群落；大針茅(Stipagrandis)、差不嘎蒿(Artemisiahalodendron)、大籽蒿(Artemisiasieversiana)群落；羊草(Leymuschinensis)、披堿草(Elymusdahuricus)、差不嘎蒿(Artemisiahalodendron)群落。植被覆蓋度以中等蓋度為主，最低為2.7%，最高為88.7%。株高平均高度為11.2 cm，最低為5.1 cm，最高為21.3 cm。

2016年7月12—14日進行地面采樣與觀測試驗，以伊敏露天煤礦區(qū)為中心向外10 km擴充為采樣區(qū)域，采樣涉及范圍達1600 km2。依照研究區(qū)植被特性，利用分區(qū)隨機取樣法布置樣點，為了保證樣方與像元的空間對應，在樣方周圍的一定范圍內(nèi)，選擇在草地植被類型、分布以及生長狀態(tài)比較均勻的區(qū)域布設樣方，以此保證植被景觀的一致性。樣方大小為60 m×60 m，無人機(型號大疆Phantom 3 Professional)飛行在樣方上方100 m左右，獲取無人機航拍遙感數(shù)據(jù)。同時在無人機樣方四周布設9個1 m×1 m小樣方，樣方布設示意圖如圖1所示，從而控制大樣方的布設精度，使用索尼數(shù)碼相機(型號SELP 1650)在1 m×1 m小樣方上方1 m處垂直拍攝記錄樣方，通過GPS記錄每個樣方的地理坐標，人工記錄植被類型、地貌環(huán)境(表1)等，共計7個大樣方，63個小樣方(圖2)。

圖1 樣方布設示意圖

圖2 研究區(qū)及采樣點位置

表1 樣點地貌類型及植被類型

1.2 數(shù)據(jù)源及預處理

遙感數(shù)據(jù)源包括2016年7月21日獲取的國產(chǎn)GF-1衛(wèi)星PMS相機數(shù)據(jù)，無人機樣方大小為60 m×60 m，像元個數(shù)在77萬左右，因此計算無人機數(shù)據(jù)空間分辨率在7 cm左右。參數(shù)如表2所示。根據(jù)數(shù)據(jù)源以及建模的要求，因使用的飛行平臺沒搭載POS系統(tǒng)，為降低由平臺本身的不穩(wěn)固性(如俯仰、側(cè)滾、偏轉(zhuǎn)等)引起的圖像畸變，由目視選擇地面變形不大的影像。無人機數(shù)據(jù)的預處理主要有幾何校正，GF-1衛(wèi)星影像數(shù)據(jù)分輻射定標、大氣校正、幾何校正等3步處理。本文基于野外實測的GPS點進行遙感影像及無人機的幾何糾正，幾何校正的模型使用一階多項式(仿射)，校正精度控制在半個像元以內(nèi)。

表2 GF-1與無人機傳感器參數(shù)

1.3 植被指數(shù)

1.3.1無人機數(shù)據(jù)植被指數(shù)的獲取 由于無人機數(shù)據(jù)只有紅色波段、綠色波段以及藍色波段，結(jié)合可見光植被指數(shù)與建模要求，無人機數(shù)據(jù)所用的植被指數(shù)為綠紅植被指數(shù)(Green Red Vegetation Index，GRVI)、可見光大氣阻抗植被指數(shù)(Visible Atmospherically Resistant Insex，VARI)、綠葉植被指數(shù)(Green Leaf Index，GLI)。以手持GPS記錄的樣方位置為中心，建立1 m×1 m的窗口，提取無人機數(shù)據(jù)3種植被指數(shù)的平均值作為建模的自變量。3種植被指數(shù)計算公式如下：

(1)

(2)

(3)

1.3.2國產(chǎn)衛(wèi)星數(shù)據(jù)植被指數(shù)的獲取 由于國產(chǎn)GF-1衛(wèi)星數(shù)據(jù)僅有紅、綠、藍及近紅外4個波段，將植被指數(shù)的優(yōu)點與計算公式相結(jié)合，選擇與植被覆蓋度相關(guān)系數(shù)較高的歸一化植被指數(shù)(Normalized Difference Vegetation Index，NDVI)、增強型植被指數(shù)(Enhanced Vegetation Index，EVI)、土壤調(diào)整植被指數(shù)(Soil-Adjusted Vegetation Index，SAVI)與修正土壤調(diào)整植被指數(shù)(Modified Soil-Adjusted Vegetation Index，MSAVI)。

表3 植被指數(shù)

即便樣方的布設選擇在植被生長均勻的位置，但樣方在遙感影像上與鄰近像元的灰度值(Digital Number，DN)值及假彩色合成后的顏色仍有一定的不同，因而實地樣方對應的植被指數(shù)也不同。

為了保證手持GPS記錄的樣方位置提取的遙感影像植被指數(shù)在像元尺度上嚴格對應無人機實測植被覆蓋度，以小樣方所在的像元為中心，建立3×3的窗口，如圖3所示，同時記為位置C，其周邊鄰近像元依次為NW，N，NE，W，E，SW，S，SE。樣方植被指數(shù)(Vegetation Index，VI)的下標表示樣方所在像元位置，如VIC表示在位置C時，像元的植被指數(shù)。

圖3 3×3窗口

像元中樣方的位置可分成4種情況：

(1)C型(圖4)：中心型(Central Type，C-Type)，樣方的位置剛好基本在像元的中心，如圖所示，此時該像元的植被指數(shù)可以直接代表遙感植被信息與實測植被覆蓋度實施相應解析。

圖4 C-Type

(2)B型(圖5)，邊界型(Border Type，B-Type)，樣方的位置剛好基本在像元的邊緣，依東西南北4個方向再分為BE型、BS型、BW型、BN型。

根據(jù)像元中樣方的不同位置，分別對樣方所在像元及其鄰近像元的植被指數(shù)取平均值，作為樣方的植被指數(shù)結(jié)果。

B型的VI計算方式為：

(4)

(5)

(6)

(7)

(3)A型(圖6)，頂點型(Angle Type，A-Type)，樣方的位置剛好基本在像元的頂點附近，按西北、東北、東南、西南四個方向又分為ANW型、ANE型、ASE型、ASW型。

(4)T型(圖7)，即過度型，除了C型、B型、A型，樣方位置位于像元的過度位置，如圖7a為C型與B型的過度位置，圖7b為C型與A型的過度位置。根據(jù)樣方采樣中的五點法采樣法，對中心像元及其鄰近像元采用5次加權(quán)平均的方法。

圖5 B-Type

圖6 A-Type

a：計算公式為：

(8)

b：計算公式為：

(9)

T型的其他情況依次類推。

圖7 T型

A型的VI計算方式為：

(10)

(11)

(12)

(13)

2 方法

2.1 照相法

照相法的原理為整個相片中綠度像元與像元總個數(shù)的比例即為植被覆蓋度。照相法估算植被覆蓋度由幾何校正、裁剪、計算相片的修正超綠特征、應用Otsu算法提取植被覆蓋度組成。校正點取相片中標準正方形樣方框架的4頂點，校正后相片邊框為正方形，而后依照邊框裁剪校正后的相片。由于健康的植被綠色波段(用G表示)的值大于紅色波段(用R表示)和藍色波段(用B表示)，從而計算其修正超綠特征(用EXG表示)。當R

最大類間方差法，由日本學者大津(Nobuyuki Otsu)創(chuàng)建的一種自適應的閾值確定方法，又叫大津法,簡稱OTSU。OTSU算法根據(jù)計算的修正超綠特征圖像的灰度特性，將圖像分為植被與背景兩部分，并且植被與背景兩部分的類間方差很大，從而可以自動獲取植被與背景的分割閾值，從而計算植被覆蓋度[24]。

記t為植被與背景的分割閾值，植被像素占相片比例為W0，平均灰度為U0；背景像素占相片比例為W1，平均灰度為U1。

則相片的總平均灰度為：

U=W0×U0+W1×U1

(14)

植被和背景相片的方差：

g=W0×(U0-U)×(U0-U)+W1×
(U1-U)×(U1-U)=W0×
W1×(U0-U1)×(U0-U1)

(15)

當方差g最大時，可以認為此時植被和背景差異最大，此時的灰度t是最佳閾值：

t=W0×W1×(U0-U1)×(U0-U1)

(16)

2.2 支持向量機回歸

Vapnik等[25]創(chuàng)建了ε-不敏感損失函數(shù)：

(17)

式中：ε控制擬合精度的不敏感系數(shù)，Lε(f(x),y)是損失函數(shù)，f(x)=ωx+b為用于擬合數(shù)據(jù){xi,yi}i=1,2,…,m的回歸函數(shù)，其中ω，b是系數(shù)，xi∈Rd，yi∈R。

若訓練樣本采用線性SVM回歸時，設ε為全部訓練樣本數(shù)據(jù)的擬合誤差精度：

(18)

(19)

(20)

式中s導出的是約束條件，C是平衡因子。

函數(shù)f(x)通過引入拉格朗日乘數(shù)，由拉格朗日對偶性把原始問題轉(zhuǎn)為對偶問題，原始問題的解由對偶問題的解求得。f(x)表示為：

f(x)=∑nSV(αi-α*i)≤xi?x≥+b

(21)

關(guān)于非線性回歸，是在線性回歸的基礎上引進核函數(shù)，輸入空間通過非線性映射到高維的特征空間，進而在高維空間上實施線性回歸。核函數(shù)種類較多，用的較廣泛的核函數(shù)是徑向基核函數(shù)(Radial Basis Function，RBF)：

(22)

式中xi表示核函數(shù)中心，p表示函數(shù)寬度參數(shù)，控制函數(shù)的徑向影響區(qū)域。

在引入核函數(shù)之后，原式將表達如下：

f(x)=∑nSV(αi-α*i)×k(x,xi)+b

(23)

3 結(jié)果分析

3.1 照相法估算植被覆蓋度結(jié)果

相片進行幾何校正及裁剪等預處理后，通過計算相片的修正超綠特征應用OTSU算法提取植被覆蓋度，提取結(jié)果如圖8所示。由于相片的分辨率可以達到毫米級別，因此照相法估算植被覆蓋度的結(jié)果可以通過目視解譯來進行精度評價。照相法估算植被覆蓋度頻率分布直方圖如圖9所示，其中以植被覆蓋度為10%～30%為主，其次為30%～45%。

圖8 照相法裁剪圖(左)與估算結(jié)果圖(右)

3.2 數(shù)碼相片—無人機數(shù)據(jù)建模估算植被覆蓋度結(jié)果及精度

無人機數(shù)據(jù)以GRVI,VARI,以及GLI作為自變量，共獲取50個數(shù)據(jù)，其中訓練樣本個數(shù)為35個，樣本在各個等級的植被覆蓋度都有選擇。結(jié)合照相法估算植被覆蓋度結(jié)果，應用SVM回歸模型建模，估算草地植被覆蓋度的精度如表4所示。模型的評價標準采用判定系數(shù)(Coefficient of Determination，記為R2)、均方根誤差(Root Mean Squared Error，RMSE)與相對分析誤差(Residual Predictive Deviation，RPD)來衡量。其中RPD(SD/RMSE，其中SD為標準層(Standard Deviation，SD))可以評價所建模型預測能力的優(yōu)劣。其中SD為標準差，RMSE為均方根誤差。當RPD<1.75，認為所建立的模型不可用；當1.75≤RPD<2.25，認為所建立的模型基本可用；當2.25≤RPD<3，認為所建立的模型基本成功；當RPD≥3，認為所建立的模型非常成功。

圖9 照相法估算植被覆蓋度頻率分布直方圖

無人機數(shù)據(jù)以GRVI,VARI以及GLI為自變量與植被覆蓋度建模得到了較好的相關(guān)關(guān)系，SVM回歸建模估算植被覆蓋度的R2在0.83左右，RPD值均在2.50左右，RPD值均高于2.25，所建立模型比較成功(圖10)，未知樣本驗證以GRVI為自變量建?？傮w精度略高，R2為0.75。

表4 數(shù)碼相片—無人機SVM建模精度

注：樣本總數(shù)50個，剔除異常值后，訓練樣本35個，驗證樣本15個

Note:After eliminating outliers,there are 50 samples,35 training samples and 15 validation samples

圖10 SVM回歸模型散點圖

3.3 無人機數(shù)據(jù)—國產(chǎn)衛(wèi)星數(shù)據(jù)建模估算植被覆蓋度結(jié)果及精度

GF-1衛(wèi)星數(shù)據(jù)以NDVI,EVI,SAVI以及MSAVI作為自變量，共獲取54個數(shù)據(jù)，其中訓練樣本個數(shù)為38個，樣本在各個等級的植被覆蓋度都有選擇。結(jié)合無人機估算植被覆蓋度結(jié)果，應用SVM回歸模型建模，估算草地植被覆蓋度的精度如表5所示。

GF-1衛(wèi)星以4種植被指數(shù)為自變量與植被覆蓋度建模均得到了較好的相關(guān)關(guān)系，SVM回歸建模估算植被覆蓋度的R2均在0.9以上，RPD值均在3以上，所建立模型非常成功(圖11)，未知樣本驗證R2均在0.85以上，以SAVI建?？傮w精度以上要略優(yōu)于其他3種植被指數(shù)。

圖11 SVM回歸模型散點圖(SAVI)

方法Method傳感器Sensor訓練樣本(n=38)Knownsamplevalidation(n=38)驗證樣本(n=16)Unknownsamplevalidation(n=16)NDVIEVISAVIMSAVINDVIEVISAVIMSAVISVMPMS相機PMSCameraRMSE3.354.503.233.405.007.244.735.38R20.950.920.960.950.930.860.950.93RPD4.633.444.864.554.032.764.263.74

注：樣本總數(shù)54個，剔除異常值后，訓練樣本38個，驗證樣本16個

Note:After eliminating outliers,54 samples,38 training samples and 16 validation samples were collected

圖12 基于SVM回歸的植被覆蓋度分級圖(SAVI)

4 討論

草地生態(tài)系統(tǒng)在“大氣—水分—土壤”的碳轉(zhuǎn)化中有著重要的地位，并且在調(diào)節(jié)全球氣候及碳循環(huán)方面也發(fā)揮著重要作用。由于人類活動例如開墾、過度放牧、采礦等的影響，草地生態(tài)系統(tǒng)的平衡被打破，尤其是近年來，我國北方草原由于受到氣候條件與土地不合理利用等因素的影響，已經(jīng)逐步出現(xiàn)土地干旱化、土壤沙化、草地荒漠化。因此，在我國北方開展草地植被覆蓋度的監(jiān)測、預測與評估的研究迫在眉睫。然而，我國北方草原地區(qū)的時空異質(zhì)性較強，傳統(tǒng)的測量方法監(jiān)測大尺度范圍內(nèi)的草地植被覆蓋度比較困難。衛(wèi)星遙感可以獲取大區(qū)域的影像數(shù)據(jù)，能夠周期性地觀測，并且由于其空間分辨率、光譜分辨率及其時相特征的多樣性，遙感技術(shù)已逐步成為監(jiān)測區(qū)域乃至全球植被覆蓋度的主要方法。地面數(shù)據(jù)以數(shù)碼相片布設的樣方因人力和物力等多種原因，布設的樣方大小無法與遙感影像的空間分辨率相匹配，而無人機可以獲取大尺度、寬范圍的數(shù)據(jù)，與遙感影像的空間分辨率相匹配，減少混合像元對植被覆蓋度估算的影響。本文通過應用SVM回歸模型，構(gòu)建數(shù)碼相片—無人機大樣方數(shù)據(jù)植被覆蓋度估算模型與無人機大樣方數(shù)據(jù)—GF-1數(shù)據(jù)植被覆蓋度估算模型，探討國產(chǎn)GF-1衛(wèi)星估算草原植被覆蓋度的方法。

本研究中，無人機數(shù)據(jù)以GRVI，VARI以及GLI 3種植被指數(shù)為自變量，以數(shù)碼相片獲取的植被覆蓋度為因變量，應用SVM回歸建立無人機植被覆蓋度反演模型，結(jié)果表明3種植被指數(shù)與植被覆蓋度擬合較好，其中以GRVI所建立的模型精度較高(R2=0.838，RMSE=7.795，RPD=2.534)。GF-1衛(wèi)星PMS相機數(shù)據(jù)以NDVI，EVI，SAVI以及MSAVI4種植被指數(shù)為自變量，以無人機大樣方數(shù)據(jù)估算的植被覆蓋度為因變量，應用SVM回歸建立GF-1數(shù)據(jù)植被覆蓋度反演模型，結(jié)果表明4種植被指數(shù)均獲得較好的相關(guān)關(guān)系，其中以SAVI與草地植被覆蓋度相關(guān)系數(shù)最大(R2=0.956，RPD=4.857，RMSE=3.232)。伊敏露天煤礦區(qū)周邊草原，因2016年降水較少的原因，植被覆蓋度整體并不高，SAVI通過引入土壤調(diào)節(jié)系數(shù)，降低了土壤背景對于估算植被覆蓋度的精度。

5 結(jié)論

通過應用SVM回歸模型分別構(gòu)建地面數(shù)據(jù)—無人機大樣方數(shù)據(jù)植被覆蓋度與植被指數(shù)之間的關(guān)系、無人機大樣方數(shù)據(jù)—GF-1數(shù)據(jù)植被覆蓋度與植被指數(shù)之間的關(guān)系，建立不同尺度數(shù)據(jù)之間的植被覆蓋度反演模型，利用無人機大樣方數(shù)據(jù)以SAVI為自變量的SVM回歸模型反演國產(chǎn)高分衛(wèi)星的草地植被覆蓋度，模型擬合較好，精度較高。

隨著我國高分衛(wèi)星的研發(fā)，將無人機與高分衛(wèi)星結(jié)合起來反演草地植被覆蓋度，可以推廣國產(chǎn)衛(wèi)星遙感數(shù)據(jù)的應用，又能夠周期性地動態(tài)監(jiān)測草原植被覆蓋度，為維護草地生態(tài)平衡與可持續(xù)發(fā)展奠定了基礎。