譚鵬源 朱建軍* 付海強(qiáng) 林 輝

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

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

1 引言

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

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

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

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

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

2.1 時間去相干半經(jīng)驗(yàn)?zāi)Ｐ?/h3>

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

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

其中，h表示森林高度；f(z)為指數(shù)形式的垂直結(jié)構(gòu)函數(shù)，描述垂直方向z上散射體的分布；σr(z)為散射體沿雷達(dá)視線方向的隨機(jī)運(yùn)動標(biāo)準(zhǔn)差，假定與森林高度呈線性關(guān)系

式中，σr表示在參考高度hr(根據(jù)先驗(yàn)信息一般設(shè)為15 m[8,9])處的運(yùn)動標(biāo)準(zhǔn)差。

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

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

2.2 MCD相干優(yōu)化算法

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

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

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

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

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

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

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

2.3 森林高度反演方法

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

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

圖2 反演森林高度與真實(shí)森林高度散點(diǎn)點(diǎn)陣橢圓示意圖Fig.2 Scatters ellipse between invert height and real height

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

式中，M表示取平均運(yùn)算。

散點(diǎn)點(diǎn)陣橢圓主軸確定后，顯然可以通過建立使逼近參數(shù)k,b分別趨近于1,0的目標(biāo)函數(shù)

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

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

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

在利用訓(xùn)練數(shù)據(jù)求得時間去相干半經(jīng)驗(yàn)?zāi)Ｐ蛥?shù)后，對每個像元求解一元非線性方程得到整個散射場景內(nèi)的森林高度結(jié)果。

2.4 多基線融合策略

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

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

通過相干特性指標(biāo)P甄選出不同干涉對在同一分辨單元內(nèi)反演出的最優(yōu)森林高度值作為多基線PolInSAR森林高度融合結(jié)果，多基線PolInSAR融合反演框架可表示為

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

3 實(shí)驗(yàn)與分析

3.1 實(shí)驗(yàn)區(qū)及數(shù)據(jù)分析

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

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

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

3.2 實(shí)驗(yàn)結(jié)果與分析

以選取的15個驗(yàn)證林分的實(shí)測森林高度(H-field)對反演結(jié)果(H-invert)進(jìn)行分析評價，圖5為3個干涉對利用HV極化反演得到的散點(diǎn)圖結(jié)果，均方根誤差RMSE分別為：4.20 m,4.03 m和3.42 m。利用μmin極化方式反演的驗(yàn)證結(jié)果如圖6所示，3個干涉對的反演精度分別提高了：20%,17%和12%，除此之外，相關(guān)系數(shù)R2也分別有所提高。分析認(rèn)為采用全極化數(shù)據(jù)結(jié)合PolInSAR相干優(yōu)化算法擴(kuò)展了極化空間，相比已有方法中選用的HV極化，μmin極化含有更少地表散射貢獻(xiàn)，更貼近時間去相干半經(jīng)驗(yàn)?zāi)Ｐ屯茖?dǎo)過程中基于“零”地體幅度比的關(guān)鍵假設(shè)。

圖3 實(shí)驗(yàn)區(qū)：綠線范圍為黃豐橋(HFQ)林場研究區(qū)域，藍(lán)色虛線為ALOS-2 PALSAR-2影像范圍，圓點(diǎn)為地面實(shí)測林分樣地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參數(shù)信息Tab.1 Parameter information of ALOS-2 PALSAR-2

圖4 相干性統(tǒng)計圖Fig.4 The histograms of coherence

圖5 單基線InSAR反演高度與驗(yàn)證數(shù)據(jù)散點(diǎn)圖Fig.5 Scatterplot comparison between inversion height of single baseline InSAR and validation data

圖6 單基線PolInSAR反演高度與驗(yàn)證數(shù)據(jù)散點(diǎn)圖Fig.6 Scatterplot comparison between inversion height of single baseline PolInSAR and validation data

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

表2 單基線PolInSAR模型參數(shù)解算結(jié)果Tab.2 Model parameter results of single baseline PolInSAR inversion

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

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

(1) 該方法利用MCD相干優(yōu)化算法獲得對體散射最為敏感的極化方式，并基于時間去相干半經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行森林高度反演，使每條單基線反演精度在一定程度上都有所提高。

(2) 利用由相干特性指標(biāo)P確定的相干區(qū)域最優(yōu)準(zhǔn)則可以優(yōu)選出同一分辨單元內(nèi)最優(yōu)的單基線森林高度反演結(jié)果。因此，相比僅利用單基線單一極化反演方法，多基線PolInSAR融合策略具有更好的穩(wěn)定性，精度也更高。