胡楚鋒 周 洲 李南京 張麟兮 許家棟

①(西北工業(yè)大學(xué)無人機(jī)特種技術(shù)重點(diǎn)實(shí)驗(yàn)室 西安 710065)

②(西北工業(yè)大學(xué)航空學(xué)院 西安 710072 )

1 引言

地表植被覆蓋的估計(jì)與監(jiān)測(cè)是極化干涉雷達(dá)遙感的重要應(yīng)用領(lǐng)域[1,2]。自從 1997 年 Cloude和Papathanassiou發(fā)現(xiàn)干涉相干性強(qiáng)烈受極化的影響，國(guó)內(nèi)外學(xué)者廣泛開展了在極化干涉 SAR 測(cè)量條件下進(jìn)行植被參數(shù)提取的研究，建立了很多植被模型[3]并提出了大量的參數(shù)反演算法[4,5]。其中，由Treuhaft所建立的植被體及地表二層相干散射模型[6]，其在物理結(jié)構(gòu)及模型復(fù)雜度上進(jìn)行了很好的折衷，特別適用于地表植被參數(shù)反演。

根據(jù)入射電磁波極化狀態(tài)在體積層中的傳播性能，植被體及地表二層相干散射模型分為兩種：隨機(jī)體及地表二層模型(Random Volume over Ground model, RVoG)和有向體及地表二層模型(Oriented Volume over Ground model, OVoG)，兩者的區(qū)別在于：電磁波在前者中傳播時(shí)的衰減系數(shù)與極化狀態(tài)無關(guān)，而在后者傳播時(shí)衰減系數(shù)是極化狀態(tài)的函數(shù)。對(duì)于大多數(shù)植被來說，其物理特征滿足RVoG模型，在此基礎(chǔ)上采用不同的反演算法或改進(jìn)方法可獲得較好的參數(shù)反演精度[7,8]。2007年，研究人員發(fā)現(xiàn)對(duì)于一些具有極化特征的植被(如水稻等)，應(yīng)用RVoG模型具有局限性[9]，人們開始關(guān)注 OVoG模型。但到目前為止，研究主要集中于OVoG 模型的應(yīng)用范圍[10,11]，針對(duì)模型的反演算法研究較少[12]。

本文提出了一種基于OVoG模型的植被參數(shù)反演算法，根據(jù)模型中植被體衰減系數(shù)隨極化的不同而不同的特點(diǎn)，反演參數(shù)的有效解域由“線”變“面”，因而將同極化之差(如 HH-VV)所代表的二面角散射特征作為地形相位點(diǎn)的判別依據(jù)，并采取計(jì)算平均高度和標(biāo)準(zhǔn)差的方法進(jìn)行求解。對(duì)模型及算法的驗(yàn)證在環(huán)境可控的微波暗室中完成，通過構(gòu)建全極化干涉SAR半實(shí)物寬帶測(cè)量系統(tǒng)，獲取由多棵具有明顯極化特征南洋杉的極化干涉寬帶回波數(shù)據(jù)，反演結(jié)果表明：當(dāng)選取地體幅度比較大的頻點(diǎn)率時(shí)，植被高度的反演誤差僅為0.03 m，隨著頻率升高，地面的回波減少，地形分布的反演精度降低，植被高度的最大反演誤差仍能控制在0.2 m之內(nèi)，說明了反演算法的有效性。

2 OVoG模型

有向體(OV)是指植被體中的散射粒子具有極化相關(guān)性，例如樹干、農(nóng)作物等都有很明顯的方位取向，入射到有向體中的電磁波將變得各向異性。為了從極化干涉數(shù)據(jù)中提取出物理參數(shù)，必須要建立散射模型將測(cè)得的數(shù)據(jù)與需要反演的參數(shù)聯(lián)系起來，OVoG模型可表示為

式中ω為極化矢量，γ是由測(cè)量所獲得的復(fù)相干系數(shù)，γv是植被體相干函數(shù)，σ為植被衰減系數(shù)，μ為地體幅度比，φ0為地形相位。將測(cè)量所獲得的復(fù)相干系數(shù)γ畫在一個(gè)復(fù)單位圓內(nèi)，如圖1所示，兩條垂直的軸線分別代表γ的實(shí)部和虛部。與 RVoG模型不同，由于衰減系數(shù)與極化相關(guān)，所以式(1)在復(fù)平面內(nèi)不是一條直線，而是一個(gè)區(qū)域。

OVoG模型的物理解釋如下：(1)由于衰減系數(shù)與極化狀態(tài)有關(guān)，因而對(duì)于每個(gè)極化通道來說，體相干系數(shù)也不一樣，它們分別對(duì)應(yīng)著不同的直線，根據(jù)地體幅度比從0到無窮大的變化范圍，可以看出OVoG模型的復(fù)相干系數(shù)分布在一個(gè)近似于三角形的區(qū)域內(nèi)。(2)有效區(qū)域的寬度除了與衰減系數(shù)的范圍有關(guān)，還與垂直波數(shù)和森林高度的乘積kz?hv有關(guān)，當(dāng)兩者乘積變大時(shí)，有效區(qū)域的寬度也隨之增加。(3)最大和最小的衰減系數(shù)分別對(duì)應(yīng)特征極化，因此，對(duì)應(yīng)特征極化的復(fù)相干系數(shù)是有效區(qū)域的邊界。

圖1 OVoG模型示意圖

3 參數(shù)反演算法

通過極化干涉雷達(dá)測(cè)量數(shù)據(jù)進(jìn)行植被參數(shù)反演的核心是森林頂部和真實(shí)地形高程的獲取。極化電磁波對(duì)散射體的電參數(shù)敏感，而干涉方式可以獲取地物散射特性的幾何參數(shù)，極化干涉SAR技術(shù)實(shí)現(xiàn)了極化和干涉技術(shù)的有效組合，既對(duì)地物散射體的空間分布和高度敏感，又對(duì)植被散射體的形狀和方向敏感。同時(shí)，極化干涉相干過程提供了相位散射中心的最優(yōu)分離。參數(shù)反演算法的具體步驟如下。

3.1 復(fù)相干系數(shù)的最小二乘直線擬合

該方法是在復(fù)平面單位圓內(nèi)尋找干涉相干系數(shù)的最優(yōu)擬合直線。根據(jù)單元圓上的兩個(gè)相位點(diǎn)變量φ1和φ2，每對(duì)相位點(diǎn)可畫出一條直線，然后選擇相干系數(shù)點(diǎn)與直線均方誤差最小的一對(duì)相位點(diǎn)。具體的實(shí)現(xiàn)方法是對(duì)相干系數(shù)點(diǎn)的實(shí)部和虛部使用最小二乘法進(jìn)行擬合，然后使用直線參數(shù)估計(jì)來獲得直線與單位圓的交點(diǎn)。首先將復(fù)相干系數(shù)表示成矩陣：

因此，擬合點(diǎn)與直線之間偏差的平方和為

確定a0,a1的值使得偏差的平方和L值最小，得到方程組的最小二乘解，從而完成復(fù)相干值的最小二乘直線擬合。

3.2 地形相位估計(jì)

與隨機(jī)體不同，有向體的散射主要來自于樹干和地面之間的相互作用，這種散射的貢獻(xiàn)類似于二面角的散射機(jī)理，即同極化之差所代表的散射特征，該極化通道的相位中心一般位于地表，因而選擇離HH-VV(或 VV-HH)較近的相干系數(shù)與單位圓的交點(diǎn)作為地形相位點(diǎn)。

在確定了地形相位點(diǎn)之后，將所有的相干點(diǎn)與地形相位點(diǎn)相連，構(gòu)成了復(fù)干涉相干值的有效區(qū)域。根據(jù)前面對(duì)相干點(diǎn)有效區(qū)域的解釋，右端的相干點(diǎn)具有最小的衰減系數(shù)σmin，而左端的點(diǎn)衰減系數(shù)最大。對(duì)于每一條直線來說，它們對(duì)應(yīng)于同一個(gè)森林高度hv。

3.3 森林高度和衰減系數(shù)估計(jì)

由于衰減系數(shù)隨極化的不同而不同，因而由先驗(yàn)知識(shí)確定一定的范圍，根據(jù)衰減系數(shù)的不同從單位圓上的地形相位點(diǎn)畫出一系列以hv為變量的體相干γv曲線。由確定的兩條邊線和這兩條曲線限定了反演植被高度hv有效的區(qū)域。由此得到許多可能的解hvi(i=1,2,…,n)，最后根據(jù)可能的解中，計(jì)算植被的平均高度和相應(yīng)的標(biāo)準(zhǔn)差ε，如式(5)，式(6)所示：

對(duì)于每種極化下的衰減系數(shù)，若復(fù)相干點(diǎn)的有效區(qū)域比較窄小，則算法對(duì)該參數(shù)不敏感，也就是說，當(dāng)衰減系數(shù)的差異較大時(shí)(如1 dB/m)，體相干函數(shù)在單位圓中的偏移可以忽略，因而不同極化下的衰減系數(shù)也不明顯。

4 室內(nèi)寬帶實(shí)驗(yàn)驗(yàn)證

系統(tǒng)在西北工業(yè)大學(xué)無人機(jī)特種技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室的綜合性能微波屏蔽暗室內(nèi)搭建完成[12]，如圖2所示。矢量網(wǎng)絡(luò)分析儀和4個(gè)雙脊喇叭天線構(gòu)成整個(gè)全極化寬帶雷達(dá)系統(tǒng)，分時(shí)發(fā)射和接收兩路寬帶信號(hào)，進(jìn)行準(zhǔn)單站測(cè)試。天線向下傾斜，保持45?的入射角，距離場(chǎng)景中心的距離為7 m。為了獲取干涉數(shù)據(jù)，在取樣架的不同位置進(jìn)行了測(cè)量，為保證在高頻段不出現(xiàn)相位模糊，兩次移動(dòng)之間的距離(即基線距)很小。測(cè)試樣品由9棵高1.5 m的小南洋杉構(gòu)成，它們種在一個(gè)邊長(zhǎng)1.5 m的方形容器內(nèi)，土壤層的厚度為0.2 m，從土壤表面到植被頂端的距離為1.3 m。南洋杉樹干直徑約為0.02 m，樹干上每隔0.2 m有5-6個(gè)長(zhǎng)約0.5 m的樹枝構(gòu)成，每個(gè)樹枝的直徑大約為0.01 m，頂端的樹枝與樹干有約60?的夾角，越接近底端時(shí)，樹枝越接近水平分布。

實(shí)驗(yàn)中，測(cè)試樣品繞著垂直軸在方位向上旋轉(zhuǎn)360?，獲得了120個(gè)角度間隔3?的全極化后向散射頻域數(shù)據(jù)。系統(tǒng)采用矢量網(wǎng)絡(luò)分析儀的步進(jìn)頻率體制進(jìn)行測(cè)量，測(cè)試的頻率范圍是2-11 GHz，頻率間隔為10 MHz，原始頻域數(shù)據(jù)是通過與空暗室的頻域數(shù)據(jù)相減所獲得的，使用時(shí)域門來隔離目標(biāo)響應(yīng)與天線耦合及暗室雜波。然后，采用一個(gè)金屬球和兩種不同旋轉(zhuǎn)角度下的二面角作為參考目標(biāo)進(jìn)行全極化校準(zhǔn)。該測(cè)量系統(tǒng)的校準(zhǔn)精度在0 dBsm時(shí)幅度誤差優(yōu)于0.1 dB，相位誤差優(yōu)于1?。

原始數(shù)據(jù)獲取了兩條基線下場(chǎng)景的極化散射矩陣SHH,SHV,SVH,SVV，根據(jù)復(fù)相干系數(shù)的定義(如式(7))，可以計(jì)算出γHH,γHV(假設(shè)γHV=γVH),γVV,γHH+VV,γHH-VV的大小。

圖2 室內(nèi)極化干涉寬帶測(cè)量系統(tǒng)及植被照片

式中S1和S2分別表示上基線和下基線的后向散射測(cè)量值。對(duì)于整個(gè)場(chǎng)景來說，假定它由許多個(gè)散射體構(gòu)成，并且每一個(gè)散射體都具有獨(dú)立的相干特性，則式(7)中的相干系數(shù)估計(jì)可以用空間平均來代替集合平均。因此，通過轉(zhuǎn)臺(tái)旋轉(zhuǎn)一圈后對(duì)樣品的多次采樣進(jìn)行空間平均獲得了不同極化通道的復(fù)相干系數(shù)。首先以發(fā)射頻率為4.5 GHz的原始數(shù)據(jù)進(jìn)行分析，對(duì)獲得的不同極化通道的數(shù)據(jù)計(jì)算其在不同角度平均情況下的相干值，計(jì)算結(jié)果如表1所示。

表1 不同角度平均下的相干值(L代表角度平均次數(shù))

由表1的結(jié)果可以看出，隨著角度平均數(shù)L的減少，各極化通道的相干值大體上也隨之降低。從表中還可以看出，HH-VV信號(hào)的相干值較小，并且在平均角度減少時(shí)，相干值降低的幅度也不大，說明該信號(hào)的散射具有較強(qiáng)的隨機(jī)性，更接近于地表相位點(diǎn)。除了以上獲得的5個(gè)復(fù)相干系數(shù)，還可以采用相干最優(yōu)化算法來尋找極化空間中具有最佳相干性的復(fù)相干系數(shù)，通過增加復(fù)相干系數(shù)的個(gè)數(shù)可以提高直線擬合的精度。

由垂直波數(shù)的表達(dá)式kz=4πB/λRs inθ可知，干涉基線距B直接關(guān)系到垂直波數(shù)的大小，因此在做極化干涉測(cè)量時(shí)，需要估算出基線距，本實(shí)驗(yàn)綜合考慮不同極化通道的影響，采用高度差為0.5 m的兩個(gè)三面角進(jìn)行基線距估計(jì)。雙天線與目標(biāo)之間的關(guān)系如圖3所示，從而，兩條基線下天線的雙程相位差為

當(dāng)a1=a2時(shí)，得到三面角1的高度為

圖3 干涉基線距估計(jì)示意圖 (單位: m)

同樣地，可以得到三面角2的高度z2，由于兩個(gè)二面角之間的高度差已知，根據(jù)測(cè)量得到的相位差，便可由高度差推算出上、下兩條基線之間的距離d1+d2。在兩種不同基線下對(duì)兩個(gè)三面角進(jìn)行 2維成像，對(duì)兩幅圖像進(jìn)行像素元匹配，在三面角所對(duì)應(yīng)像素點(diǎn)進(jìn)行高度估計(jì)，結(jié)果如圖4所示。從4個(gè)不同通道獲得的干涉基線距分別為：0.073 m,0.061 m, 0.055 m, 0.072 m，對(duì)其取算數(shù)平均，得到基線距為0.065 m。

圖4 兩個(gè)三面角的高度 (單位：m)

采用有向體及地表模型進(jìn)行參數(shù)反演。測(cè)試頻率為4.5 GHz，入射角為45?，天線距離樣品中心距離為7 m，由此計(jì)算得到垂直波數(shù)kz等于2.47。對(duì)地形相位和植被高度進(jìn)行估計(jì)，反演過程如圖5所示。圖5(a)給出了由8個(gè)不同的干涉相干值以及根據(jù)最小二乘法進(jìn)行擬合的直線，該直線與單位圓相交于兩點(diǎn)。如圖5(b)所示，由γHH-VV的位置可知，與單位圓的交點(diǎn)(0.994, -0.09909)是真實(shí)的地形相位點(diǎn)，這一點(diǎn)對(duì)應(yīng)的地形為-0.1186 m。從該點(diǎn)出發(fā)，與最大和最小相干點(diǎn)相連，由此確定了相干系數(shù)的取值范圍。根據(jù)體相干函數(shù)的表達(dá)式，衰減系數(shù)的范圍從0.3 dB/m到2 dB/m，如圖5(c)所示，畫出

圖5 參數(shù)反演過程

以這兩個(gè)衰減系數(shù)為邊界、植被高度hv為變量的兩條曲線，與上一階段中的兩條直線相交，形成了一個(gè)類似于四邊形的有效解域。在解域中取9個(gè)有效值，計(jì)算其平均值和標(biāo)準(zhǔn)差，得到的結(jié)果分別為1.33 m和0.2 m，估算的結(jié)果與實(shí)際高度誤差為0.03 m。

接下來討論反演地形和植被高度隨頻率的變化情況。分別采用基于RVoG模型和OVoG模型的參數(shù)反演方法對(duì)原始數(shù)據(jù)進(jìn)行處理，圖 6給出了從 3 GHz到10 GHz反演出的地形和植被高度。從圖中可以看出，采用OVoG模型的反演精度較高。當(dāng)發(fā)射頻率為3 GHz時(shí)，由于低頻信號(hào)對(duì)植被的穿透力較強(qiáng)，來自于地面的回波較多，因而地形的反演精度較高，但卻降低了植被高度的反演精度。隨著頻率的增加，信號(hào)對(duì)植被的穿透越來越弱，因而地面的回波也越來越少，對(duì)應(yīng)的地體幅度比越來越小，地形反演的誤差也逐漸增大。

5 結(jié)論

圖6 地形及植被高度反演結(jié)果

本文提出了一種基于OVoG模型的植被參數(shù)反演算法，通過對(duì)測(cè)量的極化干涉相干值反演出植被的物理參數(shù)，逆算法中將HH-VV(或VV-HH)極化通道所代表的二面角散射機(jī)理作為地形相位點(diǎn)的判別依據(jù)，更加真實(shí)地反映了模型中的主要散射源，結(jié)合計(jì)算平均值和標(biāo)準(zhǔn)差的方法準(zhǔn)確地反演了植被高度信息。并通過室內(nèi)寬帶極化干涉測(cè)量對(duì)模型及算法進(jìn)行了驗(yàn)證。通過對(duì)具有極化特征的南洋杉場(chǎng)景的回波數(shù)據(jù)進(jìn)行處理，反演出整個(gè)寬帶頻率范圍內(nèi)植被的高度和地形分布，反演結(jié)果具有較高的精度，說明了反演算法的有效性。文中的算法及實(shí)驗(yàn)系統(tǒng)具有很強(qiáng)的工程實(shí)用價(jià)值，可應(yīng)用于極化干涉SAR領(lǐng)域。

[1]吳一戎, 洪文, 王彥平. 極化干涉 SAR的研究現(xiàn)狀與啟示[J].電子與信息學(xué)報(bào), 2007, 29(5): 1258-1262.Wu Yi-rong, Hong Wen, and Wang Yan-ping. The current status and implications of polarimetric SAR interferometry[J].Journal of Electronics&Information Technology, 2007, 29(5):1258-1262.

[2]Lavalle M, Solimini D, Pottier E,et al.. Compact polarimetric SAR interferometry[J].IET Radar,Sonar&Navigation, 2010,4(3): 449-456.

[3]Jin Y Q. Polarimetric scattering modeling and information retrieval of SAR remote sensinga review of FDU work[J].Progress in Electromagnetics Research, 2010, 104(5):333-384.

[4]Pottier E. Recent advances in the development of the open source toolbox for polarimetric and interferometric polarimetric SAR data processing: the POLSARPRO V4.1.5 software[C]. IEEE International Geoscience and Remote Sensing Symposium, Honolulu, HI, June 25-30, 2010:2527-2530.

[5]羅環(huán)敏, 陳爾學(xué), 程建, 等. 極化干涉SAR森林高度反演方法研究[J]. 遙感學(xué)報(bào), 2010, 14(4): 814-821.Luo Huan-min, Chen Er-xue, Cheng Jian,et al.. Forest height estimation methods using polarimetric SAR interferometry[J].Journal of Remote Sensing, 2010, 14(4): 814-821.

[6]Treuhaft R N and Siqueira P R. Vertical structure of vegetated land surfaces from interferometric and polarimetric radar[J].Radio Science, 2000, 35(1): 141-177.

[7]李廷偉, 黃海風(fēng), 梁甸農(nóng), 等. 基于Freeman分解的植被參數(shù)反演新方法[J]. 電子與信息學(xué)報(bào), 2011, 33(4): 781-786.Li Ting-wei, Huang Hai-feng, Liang Dian-nong,et al.. A novel vegetation parameters inversion method based on the freeman decomposition[J].Journal of Electronics&Information Technology, 2011, 33(4): 781-786.

[8]談璐璐, 陳兵, 楊汝良. 利用POLInSAR數(shù)據(jù)反演植被高度的改進(jìn)三階段算法[J]. 系統(tǒng)仿真學(xué)報(bào), 2010, 22(4): 996-999.Tan Lu-lu, Chen Bing, and Yang Ru-liang. Improved three-stage algorithm of tree height retrieval with PolInSAR[J].Journal of System Simulation, 2010, 22(4):996-999.

[9]Lopez J M, Ballester J D, and Marquez Y. Model limitations and parameter-estimation methods for agricultural applications of polarimetric SAR interferometry[J].IEEETransactions on Geoscience Remote Sensing, 2007, 45(11):3481-3493.

[10]Garestier F, Dubois-Fernandez P C, and Champion I. Forest height inversion using high-resolution P-band Pol-InSAR data[J].IEEE Transactions on Geoscience Remote Sensing,2008, 46(11): 3544-3559.

[11]Sanjuan M J, Lopez J M, and Ballester J D. Microwave scattering profiles of a rice sample by means of polarization coherence tomography[C]. IEEE International Geoscience and Remote Sensing Symposium, Boston, MA, USA, July 7-11, 2008: 554-557.

[12]Neumann M, Ferro L, and Reigber A. Estimation of forest structure, ground, and canopy layer characteristics from multibaseline polarimetric interferometric SAR data[J].IEEE Transactions on Geoscience Remote Sensing, 2010, 48(3):1086-1104.