一种改进的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]对应的复相干系数作为观测值。可以组成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基彩色合成图。可以清楚看到在图2(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