基于ALOS-2 PALSAR-2多基線極化干涉SAR數據反演森林高度

譚鵬源 朱建軍* 付海強 林 輝

①(中南大學地球科學與信息物理學院 長沙 410083)

②(中南林業科技大學林業遙感信息工程研究中心 長沙 410004)

1 引言

森林高度是估算森林蓄積量及生物量的重要基礎數據，對于研究森林資源狀況以及分析全球生態環境、氣候變化具有重要意義。極化合成孔徑雷達干涉測量技術(Polarimetric SAR Interferometry,PolInSAR)采用微波監測模式，其回波信號不僅記錄垂直結構及其屬性信息，且可以區分同一分辨單元內不同散射體高度的能力，已被視為大范圍、高分辨率、高精度反演森林高度的有效手段之一[1]。

為了實現利用PolInSAR觀測量準確地提取森林高度，Thrauhft等人[2]建立了隨機地體二層散射模型(Random Volume over Ground,RVoG)，該模型將森林散射場景抽象為兩層，即由隨機均勻分布的散射體組成的植被層，以及微波信號不可穿透的地表層。隨后，Papathanassiou等人[3,4]進一步分析PolInSAR復相干性與RVoG模型的關聯，建立了利用PolInSAR反演森林高度的框架。實質上，該框架是基于體散射去相干的模型表達來反演森林高度等參數的，并且利用不同PolInSAR數據都獲得了較高的反演精度[5–7]。

由于森林場景具有顯著的時變性，具有長時間間隔的星載重軌干涉SAR(如ALOS-1至少為46天，ALOS-2為14天)散射場景內介電常數變化(如降雨)和風動都會產生嚴重的時間去相干。因此，除了體去相干的影響，時間去相干也是星載重軌極化干涉SAR數據中不可忽略的去相干因素，決定了森林參數反演的精度，甚至是反演成敗的關鍵。為此，Yang等人[8]在隨機移動散射模型(Random Motion over Ground,RMoG)[9]和體時去相干散射模型 (Volume Temporal Decorrelation,VTD)模型[10]基礎上，提出一種時間去相干半經驗森林高度反演方法。該方法結合少量機載LiDAR森林高度數據輔助時間去相干半經驗模型解算，利用ALOS-1 PARSAR-1 HV極化相干幅度成功實現了大尺度森林高度反演。

然而，該方法需假設HV極化不包含地表散射回波能量貢獻，事實上，L波段SAR信號具有較強的穿透性，尤其當森林高度較低或密度較小時，HV極化方式會記錄顯著地表回波信號。此外，該方法只適用于單基線干涉數據，尚未考慮多基線條件下，如何充分利用觀測幾何的多樣性提升反演結果的可靠性。因此本文的目的是針對上述反演方法的限制，利用ALOS-2 PARSAR-2多基線PolInSAR數據更為準確地提取森林高度。主要思路如下：首先利用相干最大分離算法(Maximum Coherence Difference,MCD)在極化空間內尋求具有最少地面散射能量貢獻的極化方式，以獲得更為純凈的森林冠層散射貢獻。然后利用該極化方式的相干幅度，在少量森林高度地面調查數據輔助下基于時間去相干半經驗模型進行森林高度反演。在此基礎之上，結合多基線數據根據PolInSAR相干集在復數平面內的幾何表達，甄選最優觀測干涉數據的反演結果作為森林高度反演最終結果。

2 多基線PolInSAR森林高度反演策略

2.1 時間去相干半經驗模型

綜合顧及垂直方向上散射體分布產生的體去相干、散射場景內介電特性改變和植被風動引起的時間去相干同時占主導地位，星載重軌PolInSAR復相干系數一般形式表示為[8]

式中，φ0為地表相位；和分別表示植被體層和地面層介電特性改變引起的時間去相干復因子；μ(ω)為地體幅度比，與極化方式有關；γv/m為體散射去相干和時間去相干(植被風動引起)產生的耦合去相干

其中，h表示森林高度；f(z)為指數形式的垂直結構函數，描述垂直方向z上散射體的分布；σr(z)為散射體沿雷達視線方向的隨機運動標準差，假定與森林高度呈線性關系

式中，σr表示在參考高度hr(根據先驗信息一般設為15 m[8,9])處的運動標準差。

為了解決上述模型過參數化問題，Yang等人[8]對散射場景做如下假設：(1) 散射場景內時間去相干與消光系數在空間上具有一致性；(2) HV極化方式具有較小地面散射能量貢獻，可假定其地體幅度比μmin=0，此時對于該極化方式可忽略地面層介電常數改變引起的時間去相干；(3) 假定干涉場景為零空間基線理想情況(忽略森林垂直結構引起的體散射去相干)，即垂直有效波束kz=0，此時對于式(2)適用積分第一中值定理(即對于在給定區間[a,b]有連續函數f(x)和同號可積函數g(x)，區間內存在一點ε滿足因此在上述假定條件下，式(1)可簡化為時間去相干半經驗模型[8,11]

其中，α為中值ε關于森林高度的比例因子，即ε=αh(0≤α ≤1);Sscene,Cscene分別與植被體層介電特性改變和風動引起的時間去相干有關

2.2 MCD相干優化算法

已有方法主要選用對森林冠層較為敏感的HV極化方式進行模型求解，但是ALOS-2 PALSAR-2發射具有較強穿透能力的L波段電磁波，HV極化方式回波信號中同樣會記錄顯著地表回波信號。鑒于此，本文利用ALOS-2 PALSAR-2全極化數據結合極化相干最優理論，盡可能抑制地表回波信號的干擾。具體方法如下：

在主輔極化SAR影像散射機制相同的情況下，極化干涉SAR的復相干系數表示為[12]

式中，自相關矩陣T11和T22都是標準Hermitian相干矩陣，分別描述主輔影像的極化特性，Ω12為互相關矩陣，不僅包含極化信息，還包含了主副影像不同極化通道間的干涉相位關系。ω為歸一化復投影矢量，通過轉換ω可以計算極化空間內任意極化基下對應散射機制的復相干系數，組成相干集。該復相干系數集合形成的區域邊界范圍可以看作將相干復平面旋轉任意角度，得到的實部最大和最小相干系數[13]

式中，φ為旋轉相位，在[0,π)范圍內等間隔采樣角度。式(7)求極值可以轉化為求解特征值問題即Aω=λT ω，進而得到最大和最小特征值分別對應的特征向量ω1和ω2，那么相干區域的一對邊界點可以表示為[14]

相比傳統InSAR技術只能獲取HH,HV或VV極化方式對應的復相干系數，相干集中包含了特定極化散射機理對應的復相干系數，為尋求極化空間內具有更為純凈森林冠層散射貢獻的極化方式提供了可能。相干區域范圍示意如圖1所示，其中在相干區域成對邊界點中距離最遠的一對相干系數點γA,γB(也就是相干區域長軸兩端點)，可以表征植被層和地表層有效相位中心的最大分離[15]。根據式(9)進一步確定體散射占優極化方式復相干γ(μmin)與地表散射占優極化方式的復相干γ(μmax)

圖1 不同干涉對相干區域在復平面上的示意圖Fig.1 Coherent regions on the complex plane for different interferometric pairs

其中，μmin表示具有最小地體幅度比的極化方式，對應體散射占優極化方式復相干；μmax表示具有最大地體幅度比的極化方式，對應表面散射占優極化方式復相干；kz為垂直有效波束，取決于成像相對幾何關系(垂直基線B⊥，斜距R，入射角θ和雷達波長λ)[16]

2.3 森林高度反演方法

即便簡化了模型參數和采用多基線PolIn-SAR數據增加了觀測量，利用傳統多維非線性迭代求解時間去相干半經驗模型仍存在秩虧問題。因此本文采用一種外部數據輔助反演法[8]，即先利用小范圍真實森林高度數據輔助解算出模型參數Sscene和Cscene，然后代入模型中即可得到整個散射場景范圍的森林高度結果。模型參數求解具體思路如下：對于給定模型參數初始值，利用訓練數據中的μmin極化方式相干幅度結合式(4)可以得到反演森林高度結果hinvert，它與對應真實森林高度數據hreal確定的散點圖如圖2所示。理想情況下，兩者數據散點應沿虛線y=x分布，但實際上在初始模型參數誤差存在情況下，兩者散點點陣橢圓主軸與y=x并非一致，而是存在一定的偏差。因此通過利用訓練數據對其調整來尋求散射場景最佳模型參數。

主成分分析思想[17]為實現上述思路提供了契機，即通過對訓練數據中hreal與hinvert這兩個二維數據的協方差矩陣進行特征值分解，可以確定該二維數據降維后的主軸(也就是散點點陣橢圓的長軸)斜率k

圖2 反演森林高度與真實森林高度散點點陣橢圓示意圖Fig.2 Scatters ellipse between invert height and real height

其中λ1和λ2為按降序排列的特征值，P為特征值對應的特征向量的元素。而點陣橢圓質心與虛線y=x的偏差b可以表示為

式中，M表示取平均運算。

散點點陣橢圓主軸確定后，顯然可以通過建立使逼近參數k,b分別趨近于1,0的目標函數

該目標函數可以利用高斯-牛頓迭代算法進行非線性最小二乘求解，如式(14)所示

式中，*表示最終迭代次數；通過給定模型參數初始值Sscene0,Cscene0，結合上述主成分思想可以得到初始點相應的k0,b0以及雅克比矩陣J0

然后將得到修正后的模型參數作為新的初始點進行下一次迭代，經過多次迭代后即可獲得最佳模型參數，迭代終止條件為(ε為經驗閾值，本文設為10-6)

在利用訓練數據求得時間去相干半經驗模型參數后，對每個像元求解一元非線性方程得到整個散射場景內的森林高度結果。

2.4 多基線融合策略

時間去相干、體去相干以及其他噪聲等因素會共同影響PolInSAR復相干性在復平面單位圓上的幾何表達[18]。在多基線配置下，不同干涉對在同一分辨單元內往往呈現出不同的相干區域結構(如圖1所示)。而相干特性P可以作為評價相干區域結構的指標[19]

式中，γ(μmin),γ(μmax)為2.2節所述相干區域長軸的兩端點，分別對應體散射極化通道與地表散射極化通道的復相干系數。|γ(μmin)-γ(μmax)|即為極化相干區域的長軸，反映了不同極化相干點在復數單位圓的分離程度；|γ(μmin)+γ(μmax)|為相干區域質心到坐標原點距離的2倍，反映了相干區域整體相干性的平均水平。因此，P值越大說明該干涉對具有更好的相干性質量與極化分離度，反演的結果更為可靠。

通過相干特性指標P甄選出不同干涉對在同一分辨單元內反演出的最優森林高度值作為多基線PolInSAR森林高度融合結果，多基線PolInSAR融合反演框架可表示為

式中，N為極化干涉SAR觀測基線數。

3 實驗與分析

3.1 實驗區及數據分析

研究區域黃豐橋國有林場(27°05′—27° 24′N,113° 35′—113° 55′E)呈帶狀分布，橫跨湖南省攸縣東西兩部(如圖3所示)。該林場屬亞熱帶季風濕潤性氣候區，年平均氣溫17.8 °C，年降水量1410.8 mm，大部分降雨發生于春、夏季。林場境內森林茂盛，擁有森林蓄積量90.12×104m3，森林覆蓋率達90%。林分類型以針葉林為主，包括杉木、油松、落葉松等。

地面實測數據由中南林業科技大學于2016年6～7月采集得到，通過在林區范圍內選取60個相互獨立的林分樣地以確保避免空間自相關，每個林分樣地規格為30×30 m。樹高則基于單木測高原理利用激光測高儀測得，林分高度范圍為4.60～20.20 m，平均高度為13.24 m。本文通過隨機采樣，將60個林分樣地數據隨機分為45個訓練數據(圖3黃點所示)和15個驗證數據(圖3紅點所示)兩組。

多基線星載重軌PolInSAR數據是利用日本宇航局(JAXA)提供的5景覆蓋研究區域的ALOS-2 PALSAR-2 L波段全極化數據。該SAR影像范圍如圖3藍色虛線所示，獲取時間為2016年6月至8月，獲取模式為StripMap2(SM2)，影像主要參數信息如表1所示。將5景SAR影像組成3個時間基線為14天的干涉影像對(BL1,BL2,BL3)，然后各自進行配準處理，并進行公共帶通濾波以確保去除幾何去相干。相干性以11×11窗口進行估計，并應用Boxcar濾波進行平滑處理以消除斑點效應。最后利用SRTM DEM對SAR影像進行地理編碼，并將其重采樣至與DEM空間分辨率一致(30×30 m)。圖4為3個干涉對的HV極化和μmin極化的相干性統計圖，由統計圖可見不同干涉對的相干性均較低，說明研究區域受時間去相干影響較為嚴重。

3.2 實驗結果與分析

以選取的15個驗證林分的實測森林高度(H-field)對反演結果(H-invert)進行分析評價，圖5為3個干涉對利用HV極化反演得到的散點圖結果，均方根誤差RMSE分別為：4.20 m,4.03 m和3.42 m。利用μmin極化方式反演的驗證結果如圖6所示，3個干涉對的反演精度分別提高了：20%,17%和12%，除此之外，相關系數R2也分別有所提高。分析認為采用全極化數據結合PolInSAR相干優化算法擴展了極化空間，相比已有方法中選用的HV極化，μmin極化含有更少地表散射貢獻，更貼近時間去相干半經驗模型推導過程中基于“零”地體幅度比的關鍵假設。

圖3 實驗區：綠線范圍為黃豐橋(HFQ)林場研究區域，藍色虛線為ALOS-2 PALSAR-2影像范圍，圓點為地面實測林分樣地Fig.3 The test site:the green line indicates the study area of HuangFengQiao (HFQ) forestry center,the blue dotted line indicates the ALOS-2 PALSAR-2 image range,and the dots are field measurement plots

表1 ALOS-2 PALSAR-2參數信息Tab.1 Parameter information of ALOS-2 PALSAR-2

圖4 相干性統計圖Fig.4 The histograms of coherence

圖5 單基線InSAR反演高度與驗證數據散點圖Fig.5 Scatterplot comparison between inversion height of single baseline InSAR and validation data

圖6 單基線PolInSAR反演高度與驗證數據散點圖Fig.6 Scatterplot comparison between inversion height of single baseline PolInSAR and validation data

從上述單基線森林高度反演結果看，不同干涉對反演整體精度較為接近，但是對于同一林分在利用不同干涉對反演的結果卻存在明顯差異。因此，當多基線數據可用時，我們進一步在單基線PolIn-SAR森林高度反演結果的基礎上挖掘PolInSAR數據本身特性并對其森林高度反演能力進行評判。與時間去相干相關的參數Sscene和Cscene共同反映了散射場景內的時間去相干影響水平，其中Sscene與植被體層介電特性變化相關，Cscene反映了植被體層隨機運動引起時間去相關水平。表2即為單基線PolInSAR模型參數解算結果，對于不同干涉對，Sscene越小，表明該基線在散射場景內植被體層介電變化(降水等引起)越顯著；Cscene越小，則表明植被體層隨機運動(風動等引起)越強烈。從圖4相干性統計圖也可以看出，干涉對BL1相干性相對更低，受時間去相干的影響更為嚴重。因此，在不同時間去相干以及其他噪聲影響下，每個干涉對在同一分辨單元內會具有不同的的相干特性，呈現出優劣不同的相干區域結構。

表2 單基線PolInSAR模型參數解算結果Tab.2 Model parameter results of single baseline PolInSAR inversion

3個干涉對在驗證林分的相干特性P值、反演森林高度值以及多基線融合森林高度值如表3所示，從整體看，根據相干特性P值大小從3個單基線PolInSAR反演結果中甄選出的森林高度結果更接近于實測真實森林高度。整個實驗區的多基線PolInSAR融合反演結果以及精度評定如圖7所示，均方根誤差RMSE為2.05 m，相比于已有的方法，本文提出的多基線PolInSAR融合反演策略精度至少提高了40%(與圖5中BL3基線結果對比)，同時，相關系數也提升至0.81。

表3 3個干涉對的相干特性P值以及森林高度值Tab.3 Coherence characteristicP-value and forest heights for three interferometric pairs

圖7 多基線PolInSAR融合反演Fig.7 Multi baseline PolInSAR fusion inversion

4 結束語

在多基線全極化數據可用條件下，彌補單基線InSAR觀測信息不足以及幾何結構單一的問題，對于反演結果整體精度提升具有重要作用。本文提出了一種星載重軌多基線PolInSAR反演森林高度的策略，對InSAR極化空間和觀測幾何空間進行擴展，主要結論如下：

(1) 該方法利用MCD相干優化算法獲得對體散射最為敏感的極化方式，并基于時間去相干半經驗模型進行森林高度反演，使每條單基線反演精度在一定程度上都有所提高。

(2) 利用由相干特性指標P確定的相干區域最優準則可以優選出同一分辨單元內最優的單基線森林高度反演結果。因此，相比僅利用單基線單一極化反演方法，多基線PolInSAR融合策略具有更好的穩定性，精度也更高。