蔡宗磊,苗正紅*,常 雪,劉艷慧,郝 剛,何龍濤
(1. 吉林省水利水電勘測設計研究院,吉林 長春 130021;2. 吉林省建筑科學研究設計院,吉林 長春 130021;3. 東北大學測繪遙感與數字礦山研究所,遼寧 沈陽 110819;4. 吉林省老龍口水庫管理局,吉林 琿春 133300)
植被覆蓋度是指植被在地面的垂直投影面積占研究區總面積的百分比[1-4],能直接反映陸地表面植被的生長情況,也是作為生態模型、水文模型等在內的一系列模型的輸入參數,在生態系統中起著重要的作用[5-7]。傳統的植被覆蓋度獲取方法一般選用地面測量法,如目測法、概率計算法等,該類方法受人力物力限制,只能獲取有限次數的觀測數據,難以呈現植被覆蓋度在時空上的變異性。遙感技術的特點為大范圍、實時觀測,可迅速簡便地獲得空間上植被覆蓋的變化,已成為估算區域植被覆蓋度的重要技術方法[8]。
像元二分模型、多端元混合像元分解法及回歸模型法等是遙感估算植被覆蓋度應用最廣泛的方法[9]。像元二分法以混合像元分解為基礎,對衛星信號中綠色植被與裸土信息做線性分解處理,分解后取得像元中植被信息,對植被覆蓋度進行反演,該方法可去除土壤等背景影響,但受到地物成分及類別的限制,此外要由純植被和裸土像元支撐[10-11]。多端元混合像元分解法可以彌補像元二分模型受兩種地物約束的缺點,在某種程度上能去掉土壤對植被覆蓋度反演的影響,削弱枯枝落葉層作用,在地物類型復雜的區域具有較高的反演精度,但在端元的選取上比較敏感,受端元類別、數量限制,易與相關性較強的地物產生“同譜異物”,從而影響植被覆蓋度精度[12-13]。回歸模型法的理論基礎是數理統計,創建遙感數據地表反射率與實測植被覆蓋度間的不同數學回歸模型,進行植被覆蓋度的反演。回歸模型法被普遍使用到提取植被覆蓋度中,且取得了良好效果,但需依靠較多實測數據進行植被覆蓋度的反演,普適性相對較差[14-19]。
由于草地下墊面植物群落的物種多樣性,以及生長環境復雜,應用回歸模型法中的支持向量機回歸估算草地植被覆蓋度,需實測數據進行建模,并且對模型進行精度評價。然而,在野外布設的1 m×1 m的小樣方,樣方大小不能與獲取的遙感影像的空間分辨率相匹配,從而影響建模的精度。無人機(Unmanned aerial vehicle,UAV)是一種由無線電遙控及本身程序裝置掌控,易于攜帶,為可以實施多項任務的無人駕駛飛行器[20]。無人機遙感具有靈活性強、可云下飛行、影像分辨率高、時效性強、成本低等眾多優勢,可以獲取野外大樣方的數據,可作為國產高分衛星數據估算草地植被覆蓋度的補充[21-22]。國產高分一號衛星,可以提供精度高、范圍廣的空間觀測服務,特別是搭載的PMS(Pantone matching system)相機,可以獲取8 m空間分辨率的遙感影像[23]。由于受到放牧、采礦及開墾等人類活動的影響,內蒙古自治區呼倫貝爾市伊敏露天煤礦區周邊草原具有較高的空間異質性,估算其植被覆蓋度,能為礦區復墾以及持續發展提供科學依據。
本文以伊敏露天煤礦附近的典型草甸草原為研究區,利用無人機獲取60 m×60 m大樣方數據與國產GF-1衛星PMS相機8 m空間分辨率數據構建多種植被指數作為自變量,應用支持向量機回歸,建立不同數據源(地面數據—無人機大樣方數據,無人機大樣方數據—GF-1數據)的植被覆蓋度估算模型,探討利用無人機大樣方數據結合國產GF-1數據反演中國北方草原礦區周邊草地覆蓋度的適用性與可行性,并且為國產衛星在草地生態系統中的應用提供基礎。
伊敏露天礦地處內蒙古自治區呼倫貝爾草原的中部,自1985年開采,并于1986年對礦區破壞土地進行了復墾及生態恢復,生態恢復后以針茅(Stipacapillata)、沙棘(Hippophaerhamnoides)、羊草(Leymuschinensis)、披堿草(Elymusdahuricus)等為主。其附近典型草甸草原區域以圍欄草場、退化草場、放牧草場以及灌木為主,植物群落主要分為:豬毛菜(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日進行地面采樣與觀測試驗,以伊敏露天煤礦區為中心向外10 km擴充為采樣區域,采樣涉及范圍達1600 km2。依照研究區植被特性,利用分區隨機取樣法布置樣點,為了保證樣方與像元的空間對應,在樣方周圍的一定范圍內,選擇在草地植被類型、分布以及生長狀態比較均勻的區域布設樣方,以此保證植被景觀的一致性。樣方大小為60 m×60 m,無人機(型號大疆Phantom 3 Professional)飛行在樣方上方100 m左右,獲取無人機航拍遙感數據。同時在無人機樣方四周布設9個1 m×1 m小樣方,樣方布設示意圖如圖1所示,從而控制大樣方的布設精度,使用索尼數碼相機(型號SELP 1650)在1 m×1 m小樣方上方1 m處垂直拍攝記錄樣方,通過GPS記錄每個樣方的地理坐標,人工記錄植被類型、地貌環境(表1)等,共計7個大樣方,63個小樣方(圖2)。

圖1 樣方布設示意圖

圖2 研究區及采樣點位置

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

表2 GF-1與無人機傳感器參數
1.3.1無人機數據植被指數的獲取 由于無人機數據只有紅色波段、綠色波段以及藍色波段,結合可見光植被指數與建模要求,無人機數據所用的植被指數為綠紅植被指數(Green Red Vegetation Index,GRVI)、可見光大氣阻抗植被指數(Visible Atmospherically Resistant Insex,VARI)、綠葉植被指數(Green Leaf Index,GLI)。以手持GPS記錄的樣方位置為中心,建立1 m×1 m的窗口,提取無人機數據3種植被指數的平均值作為建模的自變量。3種植被指數計算公式如下:
(1)
(2)
(3)
1.3.2國產衛星數據植被指數的獲取 由于國產GF-1衛星數據僅有紅、綠、藍及近紅外4個波段,將植被指數的優點與計算公式相結合,選擇與植被覆蓋度相關系數較高的歸一化植被指數(Normalized Difference Vegetation Index,NDVI)、增強型植被指數(Enhanced Vegetation Index,EVI)、土壤調整植被指數(Soil-Adjusted Vegetation Index,SAVI)與修正土壤調整植被指數(Modified Soil-Adjusted Vegetation Index,MSAVI)。

表3 植被指數
即便樣方的布設選擇在植被生長均勻的位置,但樣方在遙感影像上與鄰近像元的灰度值(Digital Number,DN)值及假彩色合成后的顏色仍有一定的不同,因而實地樣方對應的植被指數也不同。
為了保證手持GPS記錄的樣方位置提取的遙感影像植被指數在像元尺度上嚴格對應無人機實測植被覆蓋度,以小樣方所在的像元為中心,建立3×3的窗口,如圖3所示,同時記為位置C,其周邊鄰近像元依次為NW,N,NE,W,E,SW,S,SE。樣方植被指數(Vegetation Index,VI)的下標表示樣方所在像元位置,如VIC表示在位置C時,像元的植被指數。

圖3 3×3窗口
像元中樣方的位置可分成4種情況:
(1)C型(圖4):中心型(Central Type,C-Type),樣方的位置剛好基本在像元的中心,如圖所示,此時該像元的植被指數可以直接代表遙感植被信息與實測植被覆蓋度實施相應解析。

圖4 C-Type
(2)B型(圖5),邊界型(Border Type,B-Type),樣方的位置剛好基本在像元的邊緣,依東西南北4個方向再分為BE型、BS型、BW型、BN型。
根據像元中樣方的不同位置,分別對樣方所在像元及其鄰近像元的植被指數取平均值,作為樣方的植被指數結果。
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型的過度位置。根據樣方采樣中的五點法采樣法,對中心像元及其鄰近像元采用5次加權平均的方法。

圖5 B-Type

圖6 A-Type
a:計算公式為:
(8)
b:計算公式為:
(9)
T型的其他情況依次類推。

圖7 T型
A型的VI計算方式為:
(10)
(11)
(12)
(13)
照相法的原理為整個相片中綠度像元與像元總個數的比例即為植被覆蓋度。照相法估算植被覆蓋度由幾何校正、裁剪、計算相片的修正超綠特征、應用Otsu算法提取植被覆蓋度組成。校正點取相片中標準正方形樣方框架的4頂點,校正后相片邊框為正方形,而后依照邊框裁剪校正后的相片。由于健康的植被綠色波段(用G表示)的值大于紅色波段(用R表示)和藍色波段(用B表示),從而計算其修正超綠特征(用EXG表示)。當R 最大類間方差法,由日本學者大津(Nobuyuki Otsu)創建的一種自適應的閾值確定方法,又叫大津法,簡稱OTSU。OTSU算法根據計算的修正超綠特征圖像的灰度特性,將圖像分為植被與背景兩部分,并且植被與背景兩部分的類間方差很大,從而可以自動獲取植被與背景的分割閾值,從而計算植被覆蓋度[24]。 記t為植被與背景的分割閾值,植被像素占相片比例為W0,平均灰度為U0;背景像素占相片比例為W1,平均灰度為U1。 則相片的總平均灰度為: U=W0×U0+W1×U1 (14) 植被和背景相片的方差: g=W0×(U0-U)×(U0-U)+W1× (15) 當方差g最大時,可以認為此時植被和背景差異最大,此時的灰度t是最佳閾值: t=W0×W1×(U0-U1)×(U0-U1) (16) Vapnik等[25]創建了ε-不敏感損失函數: (17) 式中:ε控制擬合精度的不敏感系數,Lε(f(x),y)是損失函數,f(x)=ωx+b為用于擬合數據{xi,yi}i=1,2,…,m的回歸函數,其中ω,b是系數,xi∈Rd,yi∈R。 若訓練樣本采用線性SVM回歸時,設ε為全部訓練樣本數據的擬合誤差精度: (18) (19) (20) 式中s導出的是約束條件,C是平衡因子。 函數f(x)通過引入拉格朗日乘數,由拉格朗日對偶性把原始問題轉為對偶問題,原始問題的解由對偶問題的解求得。f(x)表示為: f(x)=∑nSV(αi-α*i)≤xi?x≥+b (21) 關于非線性回歸,是在線性回歸的基礎上引進核函數,輸入空間通過非線性映射到高維的特征空間,進而在高維空間上實施線性回歸。核函數種類較多,用的較廣泛的核函數是徑向基核函數(Radial Basis Function,RBF): (22) 式中xi表示核函數中心,p表示函數寬度參數,控制函數的徑向影響區域。 在引入核函數之后,原式將表達如下: f(x)=∑nSV(αi-α*i)×k(x,xi)+b (23) 相片進行幾何校正及裁剪等預處理后,通過計算相片的修正超綠特征應用OTSU算法提取植被覆蓋度,提取結果如圖8所示。由于相片的分辨率可以達到毫米級別,因此照相法估算植被覆蓋度的結果可以通過目視解譯來進行精度評價。照相法估算植被覆蓋度頻率分布直方圖如圖9所示,其中以植被覆蓋度為10%~30%為主,其次為30%~45%。 圖8 照相法裁剪圖(左)與估算結果圖(右) 無人機數據以GRVI,VARI,以及GLI作為自變量,共獲取50個數據,其中訓練樣本個數為35個,樣本在各個等級的植被覆蓋度都有選擇。結合照相法估算植被覆蓋度結果,應用SVM回歸模型建模,估算草地植被覆蓋度的精度如表4所示。模型的評價標準采用判定系數(Coefficient of Determination,記為R2)、均方根誤差(Root Mean Squared Error,RMSE)與相對分析誤差(Residual Predictive Deviation,RPD)來衡量。其中RPD(SD/RMSE,其中SD為標準層(Standard Deviation,SD))可以評價所建模型預測能力的優劣。其中SD為標準差,RMSE為均方根誤差。當RPD<1.75,認為所建立的模型不可用;當1.75≤RPD<2.25,認為所建立的模型基本可用;當2.25≤RPD<3,認為所建立的模型基本成功;當RPD≥3,認為所建立的模型非常成功。 圖9 照相法估算植被覆蓋度頻率分布直方圖 無人機數據以GRVI,VARI以及GLI為自變量與植被覆蓋度建模得到了較好的相關關系,SVM回歸建模估算植被覆蓋度的R2在0.83左右,RPD值均在2.50左右,RPD值均高于2.25,所建立模型比較成功(圖10),未知樣本驗證以GRVI為自變量建模總體精度略高,R2為0.75。 表4 數碼相片—無人機SVM建模精度 注:樣本總數50個,剔除異常值后,訓練樣本35個,驗證樣本15個 Note:After eliminating outliers,there are 50 samples,35 training samples and 15 validation samples 圖10 SVM回歸模型散點圖 GF-1衛星數據以NDVI,EVI,SAVI以及MSAVI作為自變量,共獲取54個數據,其中訓練樣本個數為38個,樣本在各個等級的植被覆蓋度都有選擇。結合無人機估算植被覆蓋度結果,應用SVM回歸模型建模,估算草地植被覆蓋度的精度如表5所示。 GF-1衛星以4種植被指數為自變量與植被覆蓋度建模均得到了較好的相關關系,SVM回歸建模估算植被覆蓋度的R2均在0.9以上,RPD值均在3以上,所建立模型非常成功(圖11),未知樣本驗證R2均在0.85以上,以SAVI建模總體精度以上要略優于其他3種植被指數。 圖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 注:樣本總數54個,剔除異常值后,訓練樣本38個,驗證樣本16個 Note:After eliminating outliers,54 samples,38 training samples and 16 validation samples were collected 圖12 基于SVM回歸的植被覆蓋度分級圖(SAVI) 草地生態系統在“大氣—水分—土壤”的碳轉化中有著重要的地位,并且在調節全球氣候及碳循環方面也發揮著重要作用。由于人類活動例如開墾、過度放牧、采礦等的影響,草地生態系統的平衡被打破,尤其是近年來,我國北方草原由于受到氣候條件與土地不合理利用等因素的影響,已經逐步出現土地干旱化、土壤沙化、草地荒漠化。因此,在我國北方開展草地植被覆蓋度的監測、預測與評估的研究迫在眉睫。然而,我國北方草原地區的時空異質性較強,傳統的測量方法監測大尺度范圍內的草地植被覆蓋度比較困難。衛星遙感可以獲取大區域的影像數據,能夠周期性地觀測,并且由于其空間分辨率、光譜分辨率及其時相特征的多樣性,遙感技術已逐步成為監測區域乃至全球植被覆蓋度的主要方法。地面數據以數碼相片布設的樣方因人力和物力等多種原因,布設的樣方大小無法與遙感影像的空間分辨率相匹配,而無人機可以獲取大尺度、寬范圍的數據,與遙感影像的空間分辨率相匹配,減少混合像元對植被覆蓋度估算的影響。本文通過應用SVM回歸模型,構建數碼相片—無人機大樣方數據植被覆蓋度估算模型與無人機大樣方數據—GF-1數據植被覆蓋度估算模型,探討國產GF-1衛星估算草原植被覆蓋度的方法。 本研究中,無人機數據以GRVI,VARI以及GLI 3種植被指數為自變量,以數碼相片獲取的植被覆蓋度為因變量,應用SVM回歸建立無人機植被覆蓋度反演模型,結果表明3種植被指數與植被覆蓋度擬合較好,其中以GRVI所建立的模型精度較高(R2=0.838,RMSE=7.795,RPD=2.534)。GF-1衛星PMS相機數據以NDVI,EVI,SAVI以及MSAVI4種植被指數為自變量,以無人機大樣方數據估算的植被覆蓋度為因變量,應用SVM回歸建立GF-1數據植被覆蓋度反演模型,結果表明4種植被指數均獲得較好的相關關系,其中以SAVI與草地植被覆蓋度相關系數最大(R2=0.956,RPD=4.857,RMSE=3.232)。伊敏露天煤礦區周邊草原,因2016年降水較少的原因,植被覆蓋度整體并不高,SAVI通過引入土壤調節系數,降低了土壤背景對于估算植被覆蓋度的精度。 通過應用SVM回歸模型分別構建地面數據—無人機大樣方數據植被覆蓋度與植被指數之間的關系、無人機大樣方數據—GF-1數據植被覆蓋度與植被指數之間的關系,建立不同尺度數據之間的植被覆蓋度反演模型,利用無人機大樣方數據以SAVI為自變量的SVM回歸模型反演國產高分衛星的草地植被覆蓋度,模型擬合較好,精度較高。 隨著我國高分衛星的研發,將無人機與高分衛星結合起來反演草地植被覆蓋度,可以推廣國產衛星遙感數據的應用,又能夠周期性地動態監測草原植被覆蓋度,為維護草地生態平衡與可持續發展奠定了基礎。
(U1-U)×(U1-U)=W0×
W1×(U0-U1)×(U0-U1)2.2 支持向量機回歸

3 結果分析
3.1 照相法估算植被覆蓋度結果

3.2 數碼相片—無人機數據建模估算植被覆蓋度結果及精度



3.3 無人機數據—國產衛星數據建模估算植被覆蓋度結果及精度



4 討論
5 結論