基于連續小波變換的蒺藜植被覆蓋度預測

劉昕,楊光*,林群,張龍英,陳昊宇,王寧,劉峰,劉晨

(1.內蒙古農業大學沙漠治理學院，呼和浩特 010018；2.山西大學計算機與信息技術學院，太原 030006)

遙感技術的發展為植被覆蓋度提取提供了極大便利，可以快速、準確地獲取某一區域的植被覆蓋度，已被廣泛應用于生態學。目前，國內外學者在監測植被各項指標等方面做了較多研究。就植被覆蓋度而言，有學者提出了一種新的估算方法，利用紅邊斜率反演植被覆蓋度代替傳統的歸一化植被指數(normalized difference vegetation index，NDVI)反演植被覆蓋度模型[2]；以紅邊參數、藍邊參數以及吸收特征參數為輸入量來估算植被覆蓋度達到了良好的效果[3]；高永平等[4]采用無人機獲取可見光影像，用差異植被指數、歸一化綠化差異指數、歸一化綠藍差異指數等建立植被覆蓋度回歸模型，結果表明，利用無人機獲取數據提取出指數所建立的模型精度更高、更準確；上述研究都利用光譜特征參數作為自變量來反演植被覆蓋度的精度。近年來，由于小波變換具有多尺度、可以由粗到細逐步觀察信號的變化、豐富的小波基函數等優點[5]，在血壓、心電圖、氣象學、通用信號處理、語言識別、計算機圖形學等多領域都得到了廣泛應用。目前，小波變換大多數用于預測土壤有機質、植物葉綠素、植物冠層含量等方面，但未對植被覆蓋度進行深入研究。

本文以固沙植被蒺藜(TribulusterrestrisL.)為研究對象，探討了研究植被覆蓋度與兩波段原始光譜植被指數之間的關系，采用連續小波變換(continuous wavelet transform，CWT)對植被原始光譜反射率進行不同尺度的分解，分析各尺度小波系數與植被覆蓋度之間的相關性，找出每個尺度所對應的最優波段，以植被指數與最優小波系數為自變量，利用偏最小二乘法(partial least square，PLSR)、支持向量機(support vector machine，SVM)兩種算法建立植被覆蓋度的估算模型，對比分析選出最優模型，旨在為植被監測提供地面的理論支持。

1 材料與方法

1.1 植被選取

試驗區位于內蒙古呼和浩特市托克托縣(北緯40°5′—40°35′，東經111°2′—111°32′，圖1)。試驗于2019年7—8月進行，選取蒺藜為研究對象，土壤為沙壤土。將GIS影像與地形地貌、土壤質地和植被狀況相結合，進行實地調查并選擇植被類型(固沙)與土壤相適宜的試驗地，選取具有代表性的20個點作為地物光譜測試點，共測量70組植被光譜數據。

圖1 研究區地理位置圖及采樣點分布Fig.1 Geographical position map of study area and location of sampling points

1.2 光譜數據采集與處理

采用SVC HR-1024便攜式地物光譜儀(美國Spectra Vista公司)測定，光譜有效范圍為350～2 500 nm，波段個數為1 024，最小積分時間1 ms，信號采集方式為藍牙傳輸[6]。選擇晴朗無云或少云無風天氣在10:00—14:00進行采集，傳感器探頭距離蒺藜冠層的高度為30～60 cm，每個光譜數據的測定時間設置成5 s，每一組植物測定10條光譜曲線，取其平均值作為該樣點的光譜反射率值。

利用光譜儀自帶軟件SVC HR-1024 剔除光譜曲線變異較大的數據，進行融合。由于在測量時波段之間存在相應的能量差異和光譜儀器自身因素等影響，導致光譜曲線會產生一些噪聲，采用DB5與Savitzky-Golay方法進行去燥與平滑，并將光譜數據重采樣至10 nm。

1.3 植被覆蓋度測算

利用圖像提取植被覆蓋度，每測一組光譜，需在相同位置與高度進行垂直拍照，測量中使用的數碼相機型號為Canon EOS 90D，拍照時要將相機放在樣本的中心位置，每組拍攝2張照片，避免個別誤差對植被覆蓋度的提取造成影響，保證所有照片中植被覆蓋度都代表了實際的植被覆蓋度，在MATLAB R2016a軟件中利用照片中植物枝葉的像元數與照片中總像元數之比進行計算，再根據閾值法提取出植被覆蓋度(vegetation coverage，VC)。植被覆蓋度計算公式[7]如下。

(1)

式中，VC為植被覆蓋度；Nveg為提取出植被枝葉的像元數；Nsum為整幅圖像的像元數。

在整個試驗期間各處理對黃瓜未見藥害現象，黃瓜生長正常，對黃瓜安全，黃瓜葉色青綠，黃瓜果實鮮靚，對產量無不良影響。

1.4 植被指數選取

目前，利用植被指數與光譜數據進行相關性分析以及反演植被各種生理化參數是常用方法。本研究選取350～2 500 nm波段范圍內原始光譜反射率中任意2個波段組合的差值指數(difference index，DI)、歸一化指數(normalized index，NI)、比值指數(ratio index,RI)、重新歸一化差異指數(re-normalize difference index，RDI)和修正簡單比率(modified simple ratio，MSR)。計算公式[8]如下。

DI(Rλ1,Rλ2)=Rλ1-Rλ2

(2)

NI(Rλ1,Rλ2)=|Rλ1-Rλ2|/Rλ1+Rλ2

(3)

RI=(Rλ1,Rλ2)=Rλ1/Rλ2

(4)

(5)

(6)

式中，Rλ1、Rλ2為任意兩波段光譜反射率。

1.5 連續小波變換

小波變換包括離散型小波變換(discrete wavelet transform，DWT)和連續型小波變換(CWT)，由于離散型小波變換會隨著分解尺度的增加而減少樣本的光譜維數，不能對植被光譜進行有效的特征信息提取。因此，本研究選用CWT對植被光譜數據數據進行不同尺度的分解(Ψa,b)，生成在不同分解尺度與不同波長上的小波系數[Wf(a,b)]，由波長和分解尺度構成的二維矩陣，將一維高光譜反射數據經過連續小波變換轉變成二維小波系數。小波系數可以提取出植被反射光譜吸收特征的形狀，并提取出特征參數與敏感波段[9]；連續小波變換有許多小波基函數，經過對比與分析后，本研究選擇Daubechies小波系中的db4母小波基函數，為了避免過多的分解尺度，選取1、2、3、4、5、6、7、8、9、10為小波系數的分解尺度，通過MATLAB R2016a對植被光譜數據進行連續小波變換。

(7)

式中，Ψa,b為連續小波分解，a為伸縮因子，b為平移因子，λ為植物高光譜數據的波段數。

(8)

式中，Wf(a,b)為不同分解尺度的小波系數；f(λ)為植被光譜反射率。

1.6 模型構建與精度驗證

PLSR最開始來源于分析化學，現在已經在許多領域應用，屬于多元統計回歸方法之一，是相關分析、主成分分析、多元線性回歸分析的集合體[10]。PLSR主要適用于含有多個自變量或多個因變量的回歸分析，可以將高光譜數據降維，利用有效數據進行建模；PLSR還會在建模的同時考慮自變量和因變量中主成分的提取，建模精度更高。

SVM方法是基于統計學基礎的一種新型機器學習方法，被用來解決分類和非線性回歸等問題，通過多次訓練使得模型最優化，選用RBF(radial basis function)為核函數，采用網絡搜索優化思想對c和g值進行優化，從2-10到210進行優化，找出最優解c和g，并得出最小交叉驗證均方誤差(mean square erro，MSE)。此方法結構簡單，可以較好地預測植被覆蓋度。

本文將植被指數與小波系數作為自變量建立PLSR回歸模型和SVM模型，將70組樣本中的52組作為訓練樣本，18組作為驗證樣本；高光譜預測模型的準確性評估主要基于決定系數(R2)和總均方根差(root mean square error，RMSE)，即R2值越大，則建立模型的穩定性越好；RMSE值越小，建立模型的準確度越高，預測能力越強。

2 結果與分析

2.1 不同植被覆蓋度光譜特征分析

單株蒺藜與蒺藜群落所表現出來的光譜曲線會隨著植被覆蓋度發生變化，圖2為不同植被覆蓋度(10%～50%)原始光譜反射率曲線。從整體來看，隨著植被覆蓋度的增加，光譜曲線的“五谷四峰”特征越來越明顯。在可見光400～760 nm波段范圍內，紅邊位置附近的“兩谷一峰”特征逐漸突出，主要由于葉綠素含量也隨植被覆蓋度發生變化；在近紅外760～1 300 nm波段范圍內，主要受植被內部結構的影響，導致光譜反射率出現變化；在短波紅外1 350～2 500 nm波段范圍內，光譜反射率受植物中水分與CO2的影響，光譜曲線呈下降趨勢。

圖2 不同植被覆蓋度原始光譜反射率曲線Fig.2 Original spectral reflectance curves of different vegetation coverage

2.2 兩波段原始光譜植被指數與植被覆蓋度的相關性

圖3顯示,5種指數中，比值指數(RI)表現形式最好，其次是修正簡單比率(MSR)、歸一化指數(NI)、差值指數(DI)、重新歸一化差異指數(RDI),其中DI形式中組合較好的波段是560～620 nm與560～1 200 nm的組合和630～1 250 nm與1 300～2 350 nm的組合；RI中組合較好的波段是570～600 nm與570～1 770 nm的組合和1 275～1 550 nm與1 370～1 770 nm的組合；NI中組合較好的波段是570～600 nm與570～1 780 nm的組合和620～1 750 nm與620～2 350 nm的組合；RDI中組合較好的波段是350～750 nm與700～1 300 nm的組合和710～1 800 nm與1 300～2 150 nm的組合；MSR中組合較好的波段是700～750 nm與700～1 900 nm的組合和1 500～1 850 nm與1 500～1 850 nm的組合。5種指數形式中，與植被覆蓋度相關性最好的波段分別是DI(2 260 nm，2 210 nm)、RI(1 410 nm，660 nm)、NI(1 470 nm，670 nm)、RDI(1 770 nm，670 nm)、MSR(1 410 nm，660 nm)。圖4為5種參數與植被覆蓋度的相關性趨勢圖，5種參數中的最優波段與植被覆蓋度相關性最好的植被參數為比值植被指數，相關系數為0.636 3，其次為歸一化植被指數和修正簡單比率指數，相關系數均達到0.6以上。

圖3 兩波段原始光譜指數與植被覆蓋度的相關性Fig.3 Correlation between two-band original spectral index and vegetation coverage

圖4 最優波段與植被覆蓋度的相關性Fig.4 Correlation between optimal band and vegetation coverage

2.3 植被覆蓋度與不同分解尺度小波系數的相關性分析

2.3.1連續小波變換光譜曲線 蒺藜高光譜數據經過連續小波變換后光譜曲線變化情況見圖5。不同分解尺度的光譜曲線表現的形態基本一致,具有明顯的規律，但變化起伏較大；分解后不同尺度的小波系數將原始光譜中吸收峰與吸收谷的位置放大；與原始光譜曲線對比后可知，光譜曲線特征隨著尺度的增加而增強，小波系數與光譜信號也逐漸增大；與原始光譜曲線特征位置進行比較后，可以看出：隨著分解尺度的增加，在500、600、670、900、1 000、1 200、1 450、1 950 nm的波谷位置均發生“藍移”的現象，向藍波方向持續移動；在550、750、950、1 100、1 500、1 870、2 010 nm的波峰位置出現向右“紅移”的現象，逐漸向紅光波段移動。從結果來看，利用連續小波變換不但可以加強信號反射的強度，也可以將光譜特征參數放大，有利于對蒺藜光譜特征進行深入分析。

圖5 1～10尺度分解光譜曲線對比Fig.5 Comparison of 1～10 scale decomposition spectral curves

2.3.2不同分解尺度小波系數與植被覆蓋度的相關性 圖6為不同分解尺度(原始光譜、一階微分光譜、二階微分光譜)小波系數與植被覆蓋度的相關系數矩陣。圖中紅色的區域均代表相關性強的地方，顏色越深，代表相關性越大。原始光譜的植被覆蓋度敏感波段主要集中在可見光波段(500～800 nm)、近紅外波段(1 200～1 350 nm)和短波紅外波段(1 800～2 050 nm、2 200～2 300 nm)。每個分解尺度對應的相關系數分別為0.720 1、0.735 7、0.741 6、0.724 8、0.763 9、0.788 9、0.769 1、0.780 7、0.761 4、0.760 1，其中尺度6相關系數最高，對應的波段為630 nm。一階微分的植被覆蓋度敏感波段主要集中在可見光波段(450～950 nm)、近紅外波段(1 200～1 300 nm)和短波紅外波段(2 200～2 350 nm)。每個分解尺度對應的相關系數分別為0.728 5、0.738 0、0.709 2、0.707 0、0.768 5、0.806 9、0.762 3、0.770 3、0.735 2、0.774 5，尺度6的相關系數最高，對應的波段是590 nm。二階微分的植被覆蓋度敏感波段主要集中在可見光波段(350～859 nm)和短波紅外波段(1 800～2 100 nm)。每個分解尺度對應的相關系數分別為0.691 1、0.632 6、0.615 1、0.642 6、0.763 8、0.781 8、0.744 0、0.765 9、0.730 3、0.756 3，尺度6的相關系數最高，對應的波段是510 nm。

圖6 不同分解尺度小波系數與植被覆蓋度相關系數矩陣Fig.6 Correlation coefficient matrix of wavelet coefficients and vegetation coverage at different decomposition scales

2.4 植被覆蓋度高光譜反演模型的優選

以原始光譜反射率提取的植被指數以及經過CWT提取的10組最優小波系數為自變量，建立與植被覆蓋度PLSR回歸模型、SVM回歸模型。為了獲得最佳的支持向量機模型，通過支持向量機函數的訓練，得到最優的c與RBF核函數中的g參數，圖7為支持向量機中最優c、g選擇圖。植被指數最優c=4，g=0.125；原始光譜小波系數最優c=102 4，g=0.003 9；一階微分光譜小波系數最優c=32，g=8；二階微分光譜小波系數最優c=102 4，g=32。由表1可以看出,植被指數、小波系數植被覆蓋度都有一定的相關性，經過CWT后所建立的模型比原始光譜反射率建立的模型精度高，更具有可行性。以植被指數為自變量建立的模型中，SVM回歸模型的相關系數最高，模型精度達到0.859 9；以不同導數變換后提取的小波系數為輸入量建立的模型中，SVM模型的相關系數與檢驗精度沒有PLSR模型效果顯著，模型不穩定。經過以上對比與分析，以二階微分小波系數為自變量的PLSR模型精度最高，可以很好地實現植被覆蓋度的估算，而且方法簡單，便于操作。經過CWT處理后的光譜數據不僅可以放大植被的特征信息，也可以提高建模時的精度，使模型更加穩定，有利于植被覆蓋度以及其他生理化參數的估算研究。

表1 植被覆蓋度估測模型的建模集與預測集結果Table 1 Vegetation coverage estimation model modeling and prediction set results

圖7 SVM最優c、g選擇Fig.7 SVM optimal c,g selection

3 討論

不同植被類型光譜反射率曲線由于生長環境、內部因素等差異表現出的光譜曲線并不完全相同，但是表現形式都為“五谷四峰”形式。因此，本研究不只是針對蒺藜，其他植被均適用。對植被原始光譜反射率進行植被指數提取、導數變換并提取小波系數，有助于凸顯出植被光譜特征，利于不同植被覆蓋度含量的光譜曲線特征變化分析。這與李棟[11]的研究結果一致，構建連續小波變換反演模型代替原始反射率構建的反演模型，反演結果更為突出，顯著地提高了模型的精度。經過對比與分析可知：以不同導數變換小波系數為自變量比植被指數為自變量的建模精度最高；2種建模方法中，二階微分變換的小波系數與PLSR結合后的現象最好，可以有效地解決噪聲以及外界因素等影響，將植被覆蓋度很好地估算出來。國內外學者利用CWT對固沙植被的研究相對較少，研究主要集中在土壤有機質、農作物等方面，因此本研究采用CWT對蒺藜進行光譜處理具有新意。

本文所優選的植被指數都是較常見的，但是該方法有一定的局限性，不具有普適性，隨著操作方式和測量條件的變化，優選植被指數的波長也會發生變化，但是經過CWT可以使整個光譜波段有序地分解，提取出的小波系數建立的模型適用性較強。劉小軍等[12]與Maire等[13]均選取了幾種典型植被指數，通過任意波段組合建立植被指數空間，在監測水稻等植物均有較高的監測精度。利用CWT處理光譜信息不但可以突出植被光譜特征，也可以將光譜隱藏信息挖掘出來，更有利于植被覆蓋度的預測，具有十分重要的意義，這與Cheng等[14]的研究結果一致；由于具有較多的連續小波基函數，在今后研究中將進一步深入分析不同小波基函數對植被光譜的變化情況。

本研究結果表明，原始光譜植被指數與植被覆蓋度呈顯著相關；利用CWT變換不僅可以增強光譜的特征信號，提取出的小波系數也與植被覆蓋度之間有良好的相關性，3種光譜變換后的小波系數與植被覆蓋度均具有良好的相關性，經過一階微分變換后小波系數與植被覆蓋度相關性最大，具有良好的效果；以原始光譜植被指數和(原始光譜、一階微分光譜、二階微分光譜)CWT提取出的小波系數為自變量建立的4種模型中，經過CWT處理后的光譜數據建立的模型精度最高；利用CWT可以提高植被光譜中的有效特征信息，PLSR可以作為反演植被中其他生理化參數的方法，為植被監測提供更有效的參考價值。