付海強(qiáng)，汪長城，朱建軍，解清華，趙 蓉

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

一種改進(jìn)的PolInSAR PCT方法反演植被垂直結(jié)構(gòu)

付海強(qiáng)，汪長城，朱建軍，解清華，趙 蓉

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

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

極化干涉SAR；極化干涉層析；植被垂直結(jié)構(gòu)；非線性迭代；最小二乘

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

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

1 一種改進(jìn)的PCT算法原理

1.1 PCT原理

相對反射率函數(shù)f(z)與干涉相干系數(shù)之間的關(guān)系可以表達(dá)為[7]

(1)

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

(2)

(3)

其中：

(4)

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

(5)

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

圖1 RVOG模型的二次結(jié)構(gòu)函數(shù)近似值隨地體散射幅度比變化關(guān)系

1.2 改進(jìn)的非線性迭代法反演植被參數(shù)

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

RVOG模型的表達(dá)式為[9,12-13]

(6)

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

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

(7)

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

(8)

之后，將解算出來的地表相位φ0、樹高h(yuǎn)v帶入PCT計算公式中。

2 實驗與結(jié)果分析

2.1 植被高度及地表相位的反演結(jié)果

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

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

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

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

2.2 PCT反演結(jié)果

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

圖2 植被高度及地表相位反演結(jié)果

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

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

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

3 結(jié)束語

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

[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]羅環(huán)敏. 基于極化干涉SAR的森林結(jié)構(gòu)信息提取模型與方法[D]. 成都：電子科技大學(xué)，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]朱建軍，解清華，左延英，等. 復(fù)數(shù)域最小二乘平差及其在PolInSAR植被高反演中的應(yīng)用[J].測繪學(xué)報，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]谷湘潛,康紅文，曹洪興. 復(fù)數(shù)域內(nèi)的最小二乘法[J]. 自然科學(xué)進(jìn)展, 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]王新洲. 非線性模型參數(shù)估計理論與應(yīng)用[M]. 武漢：武漢大學(xué)出版社, 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相干最優(yōu)理論及其驗證分析[J].測繪與空間地理信息，2014，37(1)：58-62.

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

[責(zé)任編輯：劉文霞]

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；補(bǔ)充更新日期：2014-06-08

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

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

TN951

：A

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