基于機載P波段全極化SAR數據的森林地上生物量估測

姬永杰，楊叢瑞，張王菲，曾 鵬，張甫香，屈亞妮

（1. 西南林業大學 地理與生態旅游學院，云南 昆明 650224；2. 西南林業大學 國家林業和草原局西南生態文明研究中心，云南 昆明 650224；3. 云南省紅河州測繪地理信息服務中心，云南 紅河 661199；4. 西南林業大學林學院，云南 昆明 650224）

森林在陸地生態系統中發揮著重要的碳匯作用。森林生物量的空間精準量化對了解陸地碳儲量、碳收支、碳平衡，以及它們與全球氣候變化之間的關系具有重要意義[1-3]。森林地上生物量(AGB)為林木干、枝、葉的有機體干質量，不含林木根部、灌木、草本及林下枯枝落葉，是森林生物量的主體部分。近年來，采用遙感技術反演森林AGB展露出巨大的潛力，光學、激光雷達(LiDAR)、被動微波雷達和合成孔徑雷達(SAR)等多種遙感數據用于森林AGB的估測，使得森林AGB估測的區域尺度性、經濟性及精度等不斷提升[3-7]。其中，SAR不受天時、天氣的影響，且微波波長較長，在森林中具有較強的穿透能力，因此在獲取森林密度、樹高、AGB等森林垂直結構因子方面極具潛力[2,8-10]。SAR估測森林AGB的能力有賴于其工作頻率的高低，大量研究證實了低頻波段后向散射對森林AGB變化的敏感性更強。在常用于森林監測的微波波段中，頻率較低的P波段SAR后向散射系數對森林AGB變化最為敏感[11-13]。目前P波段的研究多基于機載數據展開，且研究區多位于國外，而國內相關研究則開展較少。CARTUS等[1]基于AfriSAR、BioSAR和TropiSAR機載飛行試驗，使用P、C、L波段聯合進行了森林AGB反演，結果表明P波段與森林結構因子的樹干、樹枝的相關性較高，P波段的加入有力地提高了反演的準確度和精度。LIAO等[14]使用TropiSAR機載P波段數據，采用層析方法(TomoSAR)將相干幅度、干涉相位和后向散射特征建模對法屬圭亞那熱帶雨林AGB進行了反演，結果表明樹高特征的引入可有效提高森林AGB反演精度，決定系數(R2)最高可達0.7。馮琦等[10]使用國產機載P波段數據結合坡度因子，在考慮當地入射角和坡度的情況下建立對數統計模型，對內蒙古根河市生態站的寒溫帶針葉林進行森林AGB的反演，最高反演精度R2為0.634、均方根誤差(RMSE)為12.07 t·hm-2。

現有采用P波段SAR數據進行的森林AGB估測，多集中在采用P波段的后向散射信息、相位和相干性信息，而對于極化信息的利用則較少。由于SAR信息的差異，使得其用于森林AGB估測的方法也差異明顯，如采用后向散射信息估測森林地上生物量，多采用線性回歸參數模型；而采用相位和相干性信息利用層析技術進行森林AGB估測則多基于植被微波散射模型或電磁波信號模型進行反演。線性回歸模型簡單靈活，但通常無法表征P波段特征與森林AGB變化之間的復雜關系；植被散射模型或者電磁波信號模型能夠體現森林與P波段電磁波之間的部分物理作用機制，但模型較復雜，應用推廣困難。

近年來，隨著極化SAR數據的豐富，可提取的極化SAR特征涌現，非參數模型被廣泛應用于SAR極化特征農作物生長參數的定量反演，并表現出較強的反演能力和較好的應用推廣性。鑒于P在森林監測中的潛力，其極化特征在森林AGB估測中并未深入探索，參數模型過于簡單、植被微波散射模型理解過于困難，本研究以中國北方典型寒溫帶森林作為研究對象，使用機載P波段SAR數據，提取多種極化特征，在分析其P波段極化散射特征的基礎上，探索采用參數和非參數模型進行森林AGB估測的可行性，旨在明確P波段森林的極化散射特征，探索采用P波段極化SAR數據進行森林AGB估測具體應用的有效方法。

1 研究區概況

研究區為位于內蒙古呼倫貝爾盟的根河市大興安嶺森林生態系統國家野外科學觀測研究站(大興安嶺生態站)，50°49′～50°51′N，121°30′～121°31′E。該研究站是中國緯度最高的森林生態研究站，面積為102 km2。研究區氣候類型為典型的寒溫帶大陸性季風氣候。研究區地勢相對平緩，區域內80%的坡度小于15°，海拔高度分布在800～1 200 m，森林覆蓋率大于75%，主要樹種為興安落葉松Larix gmelinii、白樺Betula platyphylla、樟子松Pinus sylvestrisvar.mongolica等。

2 數據獲取與預處理

2.1 P波段SAR數據獲取與預處理

采用的P波段SAR數據由機載CAMSAR系統獲得[15]。通過協議方式獲取根河實驗區機載P波段全極化(HH/HV/VV/VH) SLC數據。該數據是2013年9月以“獎狀Ⅱ”飛機平臺，以CASMSAR系統右視觀測獲取的全極化SAR數據。獲取過程中飛行方向為自西向東，飛行高度為5 807 m，獲取時間為2013年9月13—16日，極化方式為HH/HV/VH/VV，產品模式為SLC，幅寬為6 km × 7 km，中心入射角為55.058°，距離向分辨率為0.666 m，方位向分辨率為0.625 m。機載P波段SAR數據的預處理關鍵步驟包括多視處理、正射校正、入射角校正和極化方位角校正。多視視數在距離向和方位向均為3，數據的正射校正、入射校正參考文獻[10]；極化方位角校正過程和算法見文獻[16]。P波段SAR數據預處理結果、森林AGB抽樣點及林分概況見圖1。

圖 1 P波段SAR數據(A)及其覆蓋區樣點分布(B)、林分概況(C)Figure 1 P band SAR data (A), the distribution of samples (B) and stand examples (C)

2.2 激光雷達(LiDAR)數據獲取與預處理

LiDAR獲取的數字表面模型(DSM)、數字高程模型(DEM)用于P波段SAR數據的地理編碼，由冠層高度模型(CHM)獲取的LiDAR森林AGB數據用于反演模型的訓練及驗證。本研究獲取的機載LiDAR數據是將Leica機載雷達系統荷載于“運-5”飛機平臺上，于2012年8—9月在根河實驗區開展飛行任務。該原始數據密度為5.6個·m-2的點云數據，激光中心波譜值為1 550 nm，在初始點云數據的基礎上，提取了研究區域高精度的DEM (圖2A)、CHM (圖2B)和森林AGB(圖2C)等衍生產品。DSM、DEM、CHM數據的詳細生成方法參考文獻[17]，高精度LiDAR森林AGB的詳細提取過程與方法參考文獻[10,18]。

圖 2 LiDAR衍生數據Figure 2 Lidar derived data

為保證建模樣本和驗證樣本能夠代表整個研究區的森林AGB水平，以高精度LiDAR森林AGB圖為基礎，按照750 m的空間采樣間隔，在ArcGIS中采用交互人工干預(去除道路及裸露地物)的方法選取113個樣點(圖1B)作為森林AGB反演模型的訓練與驗證。機載P波段SAR數據覆蓋區森林林相不均勻，平均AGB較低，約46.7 t·hm-2，且大于100 t·hm-2的采樣點僅有5個。113個樣點AGB以10 t·hm-2間 隔 得 分 布 情 況(圖1B)：0～10 t·hm-2，3個；＞10～20 t·hm-2，14個；＞20～30 t·hm-2，17個；＞30～40 t·hm-2，17個；＞40～50 t·hm-2，18個；＞50～60 t·hm-2，14個；＞60～70 t·hm-2，10個；＞70～80 t·hm-2，8個；＞80～90 t·hm-2，7個；＞90 t·hm-2，5個。

3 方法

3.1 森林P波段極化散射特征分析

極化目標分解方法是從全極化數據中提取地物極化信息的有效方法，目前多種極化分解參數在森林類型識別中表現出巨大的潛力。本研究采用目前常用的極化后向散射系數、常用的3種極化分解方法提取極化SAR特征，并基于此分析森林的極化散射特征。提取P波段HH、HV和VV等3個極化的后向散射系數，基于3個后向散射系數的雷達植被指數(RVI)、極化辨別率參數(PDR)；基于Freeman-Durden三分量分解的體散射分量(FVOL)、單次散射分量(FODD)、二次散射分量(FDBL)、體地散射比分量(FD1/FD2，1表示地散射分量是ODD和DBL的和，2表示地散射分量僅為ODD)；基于Yamaguchi的體散射分量(YVOL)、單次散射分量(YODD)、二次散射分量(YDBL)、螺旋體散射分量(YHLX)；基于H-A-ALPHA極化分解的極化散射熵(entropy)、反熵(anisotropy)、散射角(alpha)、目標方位向角(beta)、相位差角1(gamma)、相位差角2(delta)。3種極化分解方法參考文獻[19-20]。由于本次飛行試驗時并未布設角反射器用于定標，因此本研究使用的后向散射系數值僅有相對含義。

森林在P波段極化特征響應分析的目的是為了確定P波段對森林AGB動態變化敏感的極化特征參數，從而確定有效的極化特征參數進行森林AGB的估測。將研究區的森林AGB劃分為A (表示生物量在0～30 t·hm-2變化時對應后向散射系數的箱線變化，均值約20 t·hm-2)，B (31～50 t·hm-2，40 t·hm-2)，C (51～70 t·hm-2，60 t·hm-2)，D (71～90 t·hm-2，80 t·hm-2)，E (＞91 t·hm-2，100 t·hm-2)等5個變化等級，分別制作各等級相應P波段SAR提取參數值的箱線圖，分析研究區各極化特征參數對森林AGB動態變化的響應，進而分析其極化散射特征。為了定量的分析各極化特征與森林AGB變化的關系，計算了它們之間的皮爾遜相關系數及顯著性水平。

圖 3 P波段森林散射機制分析Figure 3 Forest scattering mechanism analysis at P band

3.2 森林AGB估測方法

3.2.1 多元線性逐步回歸模型 多元線性逐步回歸模型(MLR)是森林AGB估測中最常用、最經典的方法之一。與常規的線性回歸模型相比，MLR可同時完成模型輸入參數的優選，進而提高森林AGB估測的效率和精度。MLR是將自變量逐個引入，每次判斷自變量對因變量影響的顯著性，并對模型中的自變量進行檢驗，逐個從模型中剔除不顯著的變量，從而得到最優模型。即在保證顯著性值在0.05以下的情況下，篩選相關性高的特征變量，進而得到因變量的最優估計。MLR的實現算法詳見文獻[21]。

3.2.2 KNN、SVR、RF非參數模型 KNN[22-25]、SVR[22,24]和RF[24,26]是森林AGB估測中常用的非參數模型，與參數模型相比，無固定的模型結構，通常通過數據驅動的方法來確定模型結構并用于森林AGB的估測，因此也稱為機器學習方法。在森林AGB估測中，這3種方法各有優勢，因此選取這3種方法來探究非參數方法在P波段森林AGB估測中的潛力。

3.3 P波段森林AGB估測結果精度驗證

反演結果精度的定量評價通過反演結果與真值之間的決定系數(R2)、均方根誤差(RMSE)、平均絕對誤差(MAE)、估測精度(Acc)來表征。這4個參數的計算公式參考文獻[27]。

4 結果與分析

4.1 森林P波段極化散射特征分析

本研究獲取的P波段中心波長為64 cm，該波段在森林中的穿透性較強。由圖3可知：P波段20個SAR特征均表現出對研究區森林AGB變化的敏感性。P波段3種極化方式后向散射系數對森林AGB在0～110 t·hm-2內動態變化的響應與已有研究[28]相似，即后向散射系數的值均隨著森林AGB的升高而緩慢增長。在各森林AGB水平中，HH、HV和VV散射能量強度相差并不明顯，并且從箱線圖(圖3)的中值可以看出：HH的單調增長趨勢較其他2個極化明顯，且異常值較少。這可能是由于P波段波長較長，比較粗樹干等圓柱體形狀的散射體為森林中的散射體，使得去極化特征明顯降低，因此同極化的HH和VV后向散射能量對AGB變化的敏感性明顯高于HV[29]。此外，由后向散射系數計算的RVI也隨著森林AGB增加而單調增加且未出現飽和現象，并且通過箱線圖的寬度變化可以看出：各水平RVI箱寬變化不明顯，說明該值對森林結構的變化不敏感，而對森林AGB的變化比較敏感。相比RVI，PDR非單調增長，在不同森林AGB水平，占主導散射的散射體會發生明顯變化[30-31]，該現象也可由Freeman-Durden中各個參數隨森林AGB變化的響應加以說明。P波段Freeman-Durden分解相關特征中，除FODD和FDBL外，均隨著森林AGB的增加而增加，并且2種體地散射比的值均遠遠大于1。Yamaguchi分解特征對森林AGB的響應趨勢與Freeman-Durden分解相關特征的趨勢基本一致。P波段H-A-ALPHA分解提取的參數中，散射熵(entropy)隨著森林AGB增加而單調增加，而alpha則敏感性不明顯。但beta、gamma和delta對森林AGB變化的敏感性則高于散射角alpha。P波段森林的極化散射特征隨森林AGB的變化與已有研究中短波長的研究差異明顯，說明現有的極化分解建模方法與該森林AGB水平P波段的散射特征有明顯差異。

由表1可知：極化分解中表征體散射分量的極化特征與森林AGB變化的顯著性較低或不存在顯著性，而二次散射特征分量則與森林AGB變化具有強相關性且顯著性水平最高。這說明在本研究區，樹干與地表形成的二次散射對森林AGB的變化最為敏感，且目前極化分解中的體散模型并不適合P波段研究區森林AGB水平下森林的散射機制。表1的定量分析結果與圖4的定性分析結果對比可知：由于各極化參數的值域范圍差異明顯，圖3對森林AGB敏感性較明顯的參數在定量分析時相關性并不一定最高，因此圖4的分析僅可作為初步參考，對于森林AGB敏感的極化特征參數的分析，仍需要采用皮爾遜系數進行相關性定量分析。

表 1 P波段SAR極化特征與森林AGB相關性分析Table 1 Correlation coefficients of polarization feature and forest AGB

4.2 森林AGB估測結果分析

從表2可知：4種模型中，MLR模型估測結果精度最低，RF精度最高；而KNN和SVR模型的估測精度則相近，R2相同，RMSE、MAE和Acc也基本相同，比RF估測結果的最高值僅低約2%。

表 2 基于4種模型的P波段SAR 森林AGB估測結果Table 2 P band SAR forest AGB inversion using four models

從圖4可看出：在整個森林AGB分布中，MLR估測結果分布較其他3種非參數方法分散，在森林AGB大于80 t·hm-2時出現了明顯的飽和低估現象。其他3種方法盡管也有低估的現象，但是飽和現象并不明顯。此外，本研究中4種方法估測結果的相對誤差約30%，而以往區域性森林AGB的估測誤差為37%～67%[32]，采用P波段HV后向散射系數的估測結果中，同質性森林地區的相對誤差約13%，而異質性地區相對誤差則約60%[28]。在采用L-波段極化分解參數進行森林AGB反演研究中，在森林AGB水平低于120 t·hm-2時，后向散射的特征優于3種極化分解的特征，這與本研究中P波段的研究結果基本一致[33]。4種方法的估測結果中均出現了低值高估和高值低估的現象，這與以往采用不同遙感數據源進行森林AGB估測的結果一致[32]。圖5中，RF的殘差分布最接近高斯分布，盡管峰值出現在10 t·hm-2左右，由于其值分布較窄，因此具有較高的估測精度。KNN和SVR的殘差分布圖均較為連續，并且分布區間明顯高于MLR，但由于正值占比較大，因此總體上呈現高估現象。MLR估測結果的殘差分布圖在10和20 t·hm-2出現了2個明顯的峰值，且所有殘差值的分布范圍較寬，解釋了其估測精度低于其他3種方法的原因。在ENGHART等[34]基于C和L-波段參數和非參數方法森林AGB估測的研究中，也發現MLR的估測結果要略低于非參數模型。然而參數和非參數模型的適用性還受到訓練樣本大小的影響，因此在訓練樣本較大時可選擇非參數模型進行估測，當樣本較小時則可優先選擇MLR方法進行森林AGB估測[32]。

圖 4 4種模型森林AGB估測與LiDAR抽樣點散點圖Figure 4 Scatter plots of estimated and LiDAR forest AGB of four models

圖 5 森林AGB估測結果與LiDAR抽樣點差值直Figure 5 Histograms of difference between forest AGB estimation results and LiDAR sampling points

為了進一步分析P波段極化特征對森林AGB估測的潛力，本研究將研究區森林AGB劃分為不同的等級，然后采用4種估測模型中估測結果最優的RF模型進行估測，估測結果見表3。在本研究中將森林AGB劃分為3組，其中第1組分別以30、60 t·hm-2為邊界劃分為3個子組；而第2組和第3組分別以40和50 t·hm-2為界劃分為2個子組。由表3可知：分組界限不同，估測精度有明顯的差異，分組劃分越詳細，估測精度(Acc)越高；此外，3種分組情況的估測結果均表明：在森林AGB平均值約45 t·hm-2，最高值不超過120 t·hm-2時，P波段在森林AGB水平較高的分組估測精度較高，如在第1組中，森林AGB大于30 t·hm-2分組的估測精度比小于30 t·hm-2分組的估測精約高5%；而在以50 t·hm-2為分組界限的2組中，森林高AGB組的估測精度比低AGB組的估測精度約高出6%。表3的結果分析表明：待估森林的AGB水平對P波段極化特征進行森林AGB估測的估測精度有明顯影響，且P波段極化信息更適合森林AGB較高區域森林AGB的估測。

表 3 基于4種模型的P波段SAR 森林AGB反演情況Table 3 P band SAR forest AGB inversion based on four models

5 結論

針對P波段極化SAR數據在森林AGB估測中的潛力，提取了20個P波段極化SAR參數，探索了森林AGB估測中常用的MLR、KNN、SVR和RF方法在使用P波段極化SAR數據進行森林AGB估測的潛力。結果表明：P波段極化SAR信息在森林AGB估測中具潛力，但估測精度受到待估區域森林AGB水平高低的影響；4種估測方法中，非參數方法的估測結果明顯優于MLR估測。P波段森林AGB估測結果中，同樣存在AGB低值高估和高值低估的現象，其原因仍需要進一步探索。此外由于本研究區的森林AGB均值為45 t·hm-2，最高值低于120 t·hm-2，所以P波段極化信息在AGB高于120 t·hm-2的森林覆蓋區中對森林AGB的估測能力仍有待進一步研究。