一種改進的PolInSAR PCT方法反演植被垂直結構

付海強，汪長城，朱建軍，解清華，趙 蓉

(中南大學 地球科學與信息物理工程學院，湖南 長沙 410083)

一種改進的PolInSAR PCT方法反演植被垂直結構

付海強，汪長城，朱建軍，解清華，趙 蓉

(中南大學 地球科學與信息物理工程學院，湖南 長沙 410083)

針對傳統PCT方法中“相干相位-幅度聯合反演算法”的缺點，采用RVOG模型，利用改進的非線性迭代算法反演植被高、地表相位。改進的非線性迭代算法不僅充分利用不同極化方式對應的復相干系數，同時兼顧復相干系數的先驗統計誤差，提高參數解算的可靠性，進而提高PCT結果的反演精度。最后，采用兩景德國E-SAR數據進行實驗，實驗結果表明：文中提出的方法能較好地反映植被的垂直結構信息；植被冠層對應的平均相對反射率函數近似服從高斯分布；反演的相對反射率值與植被的種類、密度存在一定關聯。

極化干涉SAR；極化干涉層析；植被垂直結構；非線性迭代；最小二乘

根據2007年ICPP的報告，全球98%的可利用土地變化體現在熱帶雨林的砍伐，這一人類活動將會導致以生物量方式存儲的碳以氣體形式排放到大氣當中[1]。而植被高度及植被垂直結構是計算生物量的重要輸入參數。已有研究表明，顧及植被垂直結構的生物量反演模型相比經驗模型[1]更為嚴謹，原因在于模型考慮了植被密度、種類等因素的影響[2]。故對植被垂直結構的研究具有重要意義。極化干涉SAR技術(PolInSAR)以其既具有對目標物高度分布敏感的特點又具有對目標物的形狀、方向敏感等特點，在研究植被垂直結構參數中具有巨大潛質[2]。

植被垂直結構可以體現在雷達相對反射率函數(Relative reflectivity function)f(z)中，對于f(z)的反演過程目前主要有兩種方法：① 2009年，Treuhaft等人假定f(z)服從高斯分布，建立了f(z)與極化相干系數之間的函數關系，利用多基線單極化數據建立參數解算模型，求解出高斯分布函數中的未知參數，實現了對植被垂直結構的重建[3]。2010年，Franck Garestier等人在Treuhaft的基礎上將該思想引入到全極化數據中，利用單基線全極化數據完成了對植被垂直結構的重建[4]。② 2006年，Cloude提出極化干涉層析技術(Polarization Coherence Tomography，PCT)，該方法利用傅立葉-勒讓德級數對f(z)進行展開，根據基線的數量，確定級數展開的次數，利用極化相干系數估計展開級數的系數，完成垂直結構的重建[5]。之后，2009年，Cloude將單基線PCT技術擴展到雙基線中，并對EMSL數據進行了成功實驗[6]。Cloude提出算法相比Treuhaft提出算法更為簡單，且不限定f(z)服從何種分布，得到了廣泛的應用研究。但是，該PCT技術反演的效果依賴于先驗信息(植被厚度、地表相位)的精度。Cloude提出利用“相干相位-幅度聯合反演算法”解算植被高度及地表相位。該算法的實質為對相干相位反演結果與相干幅度反演結果進行加權平均求和。反演結果一定程度上取決于權重因子的確定，Cloude提出將權重因子固定(如0.8)[7]。此外，該方法只采用兩種極化方式對應的復相干系數作為觀測量且沒有考慮觀測量的先驗統計誤差，這使得相干直線的擬合過程存在一定的不確定性，導致地表相位估計值存在偏差。針對這些問題，本文提出改進的非線性迭代算法解算植被厚度、地表相位，以此為輸入參數，利用PCT模型完成植被垂直結構的重建。最后，利用兩景E-SAR L波段數據對算法的可行性進行驗證。

1 一種改進的PCT算法原理

1.1 PCT原理

相對反射率函數f(z)與干涉相干系數之間的關系可以表達為[7]

(1)

式中:s1,s2分別為主、副影像的復數信號；φ0為地表相位；hv為植被高度；kz為有效垂直波數；f(w,z)表示極化狀態w下，隨高度z而變化的相對反射率函數，z的變化范圍為0≤z≤hv。對式(1)中f(w,z)進行傅立葉-勒讓德級數展開，可以得到f(w,z)歸一化后的估計值[5-7]

(2)

(3)

其中：

(4)

為了驗證2階勒讓德展開式的可行性，下面通過一組模擬數據實驗來說明該展開式能否表達2層RVOG模型。在假定消光系數為0時RVOG模型可以表達為[7]

(5)

式中:φ0表示地表相位,γv表示純體去相干系數,μ(w)表示極化狀態w對應的有效地體散射幅度比,γ(w)表示極化狀態w對應的干涉相干系數。

圖1 RVOG模型的二次結構函數近似值隨地體散射幅度比變化關系

1.2 改進的非線性迭代法反演植被參數

由上面的分析可知，植被高度及地表相位是PCT模型必要的輸入參數。目前，用于植被高提取最為常用的模型為隨機地體二層散射模型(Random Volume over Ground，RVOG)[2,9-11]，“相干相位-幅度聯合反演算法”中的相干幅度反演算法的數學模型實際上是RVOG模型在不考慮消光系數影響時的一種特例，相比較而言，RVOG模型更為嚴謹。2001年，Papathanassiou以RVOG模型為基礎提出非線性迭代算法，并利用機載E-SAR數據驗證了算法的可靠性。然而，模型沒有考慮多余觀測信息及復相干系數的先驗統計誤差。針對這些問題，本文提出一種改進的非線性迭代算法。

RVOG模型的表達式為[9,12-13]

(6)

其中：γ(w)為極化狀態w對應的復相干系數，為已知量；φ0為地表相位，為未知實數；μ(w)為有效地體散射幅度比，與極化狀態有關，為未知實數；γv為“純”體去相干系數，其為消光系數σ、樹高hv、平均入射角θ及垂直向有效波數kz的函數。

當選擇大于3種極化方式的復相干系數時，便可以組成非線性方程組，求解未知參數。在測量平差領域，存在多余觀測時，即可采用最小二乘進行平差處理，其可有效抑制觀測噪聲，進而提高參數的估算精度。但是，注意到式(6)中的觀測量為復數，故該問題應該采用復數最小二乘進行解決。根據文獻[11，14]，以復數觀測值殘差的模的平方和最小作為平差準則：

(7)

其實質是對復數的實部、虛部進行聯合平差。按照三階段算法中相干直線擬合的方法，將復相干系數γ(w)拆分成實部、虛部[10]，之后，建立聯合平差模型。本文選用HH、HV、VV、HH+VV、HH-VV、PD (PDHigh、PDLow)、MCD (opt1、opt2、opt3)10種極化方式[7]對應的復相干系數作為觀測值?？梢越M成20個誤差方程(每個γ(w)可以組成兩個觀測方程)，13個未知數(φ0，σ，hv及10個地體幅度比μi)，多余觀測數為7。對于隨機模型，采用Cramer-Rao邊界計算復相干系數模的中誤差[15]，對不同極化方式的復相干系數進行定權，并假定實部為等權。最后，參數解算問題可歸結為非線性最小二乘問題[16]，參數解算的約束條件如下：

(8)

之后，將解算出來的地表相位φ0、樹高hv帶入PCT計算公式中。

2 實驗與結果分析

2.1 植被高度及地表相位的反演結果

為了驗證算法的可行性，本文采用德國DLR E-SAR數據對算法進行實驗驗證。圖2為植被高度及地表相位反演結果，圖2(a)為研究區域對應的Pauli基彩色合成圖?？梢郧宄吹皆趫D2(a)的左上角與右下角有兩塊林分。其中，林分1主要以針闊混交林為主[17]，林分2為人工種植管理的云杉[18]。剖面線A-A′與B-B′用于下文分析。

圖2(b)、圖2(c)為已有方法反演的植被高度及地表相位圖。圖2(d)、圖2(e)為新方法反演的植被高度及地表相位圖。對于林分1，通過查閱文獻[17]中的結果，發現兩種方法反演的植被高度變化范圍與文獻中的結果較為一致。但是已有方法反演的植被厚度在圖中右側出現一定程度的高估，見圖3(a)。對于地表相位，新方法反演的地表相位變化更加連續，主要在0 rad附近波動，見圖3(c),這與實際地形較為平坦這一事實較為相符。對于林分2，對照文獻[18]中的成果，由圖2(b)、圖2(d)可以明顯看出已有方法反演結果出現明顯偏差,其中植被高度出現高估，最大值達到55 m，見圖3(b),與實際情況不符。對于地表相位，已有方法相比新方法出現較多區域反演失敗,見圖3(d)。由于缺乏地面高分辨率實測數據，無法對結果進行定量分析，但是通過與已有結果對比分析可知，總體來講，新方法反演的植被參數更為可靠，原因在于[19-20]：

1)該算法采用平差的方法，充分利用了不同類型(10個)的觀測量，多余觀測量對參數解算起到可靠的約束作用，其相比“相干相位-幅度聯合反演算法”只利用兩個觀測量，更能準確“相干直線”，進而較為準確地反演地表相位及植被厚度。

2)該方法考慮了觀測量的先驗統計誤差，利用平差的思想可以盡可能地削弱誤差對待估參數的影響。而“相干相位-幅度聯合反演算法”在觀測噪聲較大的條件下，通過誤差傳遞，引起待估參數估值出現偏差。

2.2 PCT反演結果

通過對植被高度、地表相位反演結果的對比分析，可見，本文提出的算法更為可靠。這為獲得準確的PCT結果打下良好的基礎。由于已有方法反演的地表相位及植被高度出現較大的偏差，再將其作為PCT的輸入參數已經失去參考價值，故下文不給出以其作為輸入參數的PCT結果。利用本文方法反演的植被厚度、地表相位作為PCT的輸入參數進行分析。圖3為新方法和已有方法反演植被高度、地表相位的剖面圖。

圖2 植被高度及地表相位反演結果

圖3 新方法和已有方法反演植被高度、地表相位的剖面圖

圖4(a)、圖4(b)為圖3中沿剖面線A-A′和B-B′做HV極化方式PCT結果剖面線。HV極化方式對于體散射較為敏感[10]，可以看到對于植被冠層有明顯的回波信號。圖4(a)通過對照光學遙感影像，發現在像素100附近，植被密度相對較小，微波信號衰減較弱，植被層相對反射率較小，到達地面的信號強度較強。新方法反演的相對散射率與先驗知識較為相符，間接證明了新方法的可靠性。此外，注意到在林分1植被底層，仍有明顯的回波信號，說明對于針闊混交林，在植被層底部容易發育低矮的植被。對于林分2，主要以人工種植管理的云杉為主，由于該林分中樹種較為單一，呈現出較為均勻的結構特性。由圖4(b)看到利用本文提出的方法反演的相對反射率在植被層變化均一，這與林分2人工管理、樹種單一、植被冠層結構特性均一的信息較為相符。按照文獻[8]中提到森林密度在PCT結果中反映的特征，可以判斷林分2相比林分1具有較大的密度。植被底層沒有出現明顯的體散射，分析原因在于：林分2覆蓋的地表主要以裸地為主，體散射回波信號較弱，這與云杉林區底層不適于其他植被生長這一林業先驗知識較為相符；該區域植被密度較大，信號衰減嚴重，到達地面的回波信號較弱。

圖5(a)、圖5(b)分別為林分1與林分2植被冠層對應的平均相對反射率，二者分布有明顯差別，整體來看林分1回波信號更靠近植被冠層頂部。這與兩個林分樹種的不同有關，對此進行如下分析可以看到：針闊混交林，回波信號主要來自植被冠層頂部，說明該區域植被冠層厚度占總樹高的比重較小，靠近樹高頂部冠層密度較大。對圖5(a)中靠近頂部的反射率曲線進行高斯函數擬合，其相關系數達到0.91，說明對于林分1，其冠層相對反射率函數服從高斯分布。對于植被底部，反射率明顯增強，反射率曲線近似服從線性分布。對于林分2，云杉植被冠層厚度占總高度的比重較大，其最大反射率均值對應高度相比林分1要偏低，平均相對反射率曲線同樣符合高斯分布，其相關系數為0.95。林分2底部對應的相對反射率較低。以上說明對于針葉林，植被冠層厚度占樹高比例較大，植被冠層密度較大的區域靠近樹干底部?？梢姼鶕撉€，可以大致判斷植被冠層的形態特征，這對于生物量模型的建立，以及利用極化數據判別樹的種類具有重要指導意義。本文方法反演的相對反射率函數變化特點能較為準確地反映植被的真實結構特性，進一步說明了方法的可行性。

3 結束語

本文提出一種改進的PolInSAR PCT技術反演植被垂直結構。采用改進的非線性迭代算法對植被高及地表相位進行解算。最后將計算得到的樹高、地表相位帶入PCT模型中，對相對反射率函數進行解算。實驗結果表明，本文提出的方法能較準確描述植被的垂直結構信息。通過對反演結果的分析發現，PCT技術反演的相對反射率能區分出植被密度、形態特征的區別。這一點對于深入研究生物量具有重要指導意義及研究價值。

[1]LETOAN T,et al. The BIOMASS mission: mapping global forest biomass to better understand the terrestrial carbon cycle[J]. Remote Sensing of Environment, 2011, 115(11): 2850-2860.

[2]羅環敏. 基于極化干涉SAR的森林結構信息提取模型與方法[D]. 成都：電子科技大學，2011.

[3]TREUHAFT R N,et al. Vegetation profiles in tropical forests from multibaseline interferometric synthetic aperture radar, field, and lidar measurements [J]. Journal of Geophysical Research, 2009, 114(D23).

[4]GARESTIER F, LETOAN T. Estimation of the backscatter vertical profile of a pine forest using single baseline P-Band PolInSAR data [J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(9): 3340-3348.

[5]CLOUDE S R. Polarization coherence tomography [J]. Radio Sci., 2006, 41(4).

[6]CLOUDE S R. dual-baseline coherence tomography [J].IEEE Geoscience and Remote Sensing Letters, 2007, 4(1): 127-131.

[7]CLOUDE S R. Polarisation: Applications in remote sensing [M]. New York: Oxford University Press, 2009.

[8]PRAKS J,et al. SAR coherence tomography for boreal forest with aid of laser measurements[C]. Geoscience and Remote Sensing Symposium 2008.IGARSS 2008. IEEE International. IEEE, 2008, 2: II-469- II-472.

[9]PAPATHANASSIOU K P, CLOUDE S R. The effect of temporal decorrelation on the inversion of forest parameters from Pol-InSAR data[C]. International Ernational Geoscience and Remote Sensing Symposium. 3. 2003.

[10]CLOUDE S R, PAPATHANASSIOU K P. Three-stage inversion process for polarimetric SAR interferometry[J]. IEE Proc.-Radar Sonar Navig (S1350-2395), 2003, 1503(3): 125-134.

[11]朱建軍，解清華，左延英，等. 復數域最小二乘平差及其在PolInSAR植被高反演中的應用[J].測繪學報，2014，43(1)：45-51.

[12]SEUNG-KUK L, KUGLER F, PAPATHANASSIOU K P, et al. Quantifying temporal decorrelation over boreal forest at L-and P-band[C]. In Synthetic Aperture Radar (EUSAR), 2008 7th European Conference on, 2008, 1-4.

[13]LAVALL M, SIMARD M, HENSLEY S. A temporal decorrelation model for polarimetric radar interferometers. Geoscience and Remote Sensing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(7): 2880-2888.

[14]谷湘潛,康紅文，曹洪興. 復數域內的最小二乘法[J]. 自然科學進展, 2006,16(1): 49-54.

[15]SEYMOUR M S, CUMMING I G. Maximum likelihood estimation for SAR interferometry[C]. IEEE Geoscience and Remote Sensing Symposium, 1994, 4: 2272-2275.

[16]王新洲. 非線性模型參數估計理論與應用[M]. 武漢：武漢大學出版社, 2002.

[17]PAPATHANASSIOU K P, CLOUDE S R. Single-baseline polarimetric SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001,39(11):2352- 2362.

[18]REIGBER A, NEUMANN M, GUILLASO S, et al. Evaluating PolInSAR parameter estimation using tomographic imaging results[C]. In Radar Conference, 2005, 189-192.

[19]尚雙林，劉國林，陶秋香.極化干涉SAR相干最優理論及其驗證分析[J].測繪與空間地理信息，2014，37(1)：58-62.

[20]程雪姣，徐佳，劉慶群，等.面向城市地物分類的SAR圖像紋理特征提取與分析[J].測繪與空間地理信息，2014，37(4)：47-50.

[責任編輯：劉文霞]

A modified PolInSAR PCT method to invert vegetation vertical structure

FU Hai-qiang, WANG Chang-cheng, ZHU Jian-jun, XIE Qing-hua, ZHAO Rong

(School of Geosciences and Info-Physics, Central South University, Changsha 410083,China)

There are shortcomings for “phase-amplitude joint inversion algorithm” of PCT method. A modified nonlinear iteration method is applied to inverting the vegetation height, ground phase and temporal decorrelation, based on RVOG model. The modified method can not only take full advantage of different complex coherence values of polarizations, but also can consider the prior statistics errors of different complex coherence values. This makes parameter estimations more reliable and products more accuracy result. The modified method can prove two necessary parameter, vegetation height and ground phase. Finally, the new approach is validated on E-SAR data of Germany. It demonstrates that this method can reflect vegetation vertical structure well. The mean relative reflection ratio function follows the Gaussian distribution. And the mean relative reflection ratio values have relationship with vegetation species and density.

polarimetric interferometric SAR (PolInSAR); polarization coherence tomography (PCT); vegetation vertical structure; nonlinear iteration；least squares

2013-08-12；補充更新日期：2014-06-08

國家自然科學基金資助項目(41274010, 41371335)；湖南省自然科學基金資助項目(12JJ4035)；國家863計劃資助項目(2012AA121301)

付海強(1987-)，男，碩士研究生.

TN951

：A

：1006-7949(2014)11-0056-06