光学遥感或摄影测量手段绘制大范围地形图时无法准确去除植被高，直接影响了DTM的提取精度。准确获得植被高对获得高精度DTM具有重要意义[1, 2, 3, 4, 5]。极化干涉合成孔径雷达(polarimetric SAR interferometry,PolInSAR)集PolSAR和InSAR技术于一体，既具有InSAR对地表植被散射体的空间分布和高度敏感的特性，又具有PolSAR对植被散射体的形状和方向敏感的特性[2]。PolInSAR的这一特点为大范围获取植被高度带来了契机[3, 4]。
目前，利用PolInSAR技术提取植被高最为常用的模型为随机地体二层散射模型(random volume over ground,RVoG)。该模型将植被覆盖区抽象为植被层和地表层两层[5, 6]。2001年,文献[7]将该模型扩展到极化干涉技术中，并提出六维非线性模型来反演植被高。2003年，文献[9]在文献[7]的基础之上提出三阶段法反演植被高度，该算法简化了参数解算过程的复杂程度。大量试验证明了这两种算法的可行性[7, 8, 9, 10, 11]。然而，这两种算法并没有考虑复数观测量的先验统计误差，这可能导致参数解算存在较大偏差。针对这一问题，文献[12]将该模型的解算概括为复数非线性最小二乘问题，对复数最小二乘的平差准则及随机模型的建立进行了讨论。但是，文献[12]只是将植被高反演问题作为验证复数最小二乘的实例，并没有从植被高度反演的角度深入剖析算法的有效性。当将该方法推广至真实SAR数据反演植被高度时，存在以下问题：①该方法由于采用的函数模型过于复杂，不易进行线性化处理，需要通过迭代运算进行解算，这就使得参数的解算精度受限于迭代初值选取的可靠性；②该算法的运算效率较低，对于大范围植被参数反演存在难度；③该方法没有考虑时间去相干这一因素，限制了该方法在真实全极化数据中的应用。
本文以文献[12]为基础，提出PolInSAR植被参数反演复数最小二乘平差法。该方法首先将RVoG模型扩展为考虑时间去相干影响的RVoG+VTD(volume temporal decorrelation,VTD)模型[13]，以削弱时间去相干的影响；之后，提出了RVoG+VTD模型线性化的方法；随后，采用DEM差分法的结果作为模型解算的初值，利用截断奇异值修正方法对误差方程进行解算[21]；最后，利用德国Oberpfaffenhofen地区的E-SAR数据进行试验验证。
2 RVoG+VTD散射模型在顾及时间去相干影响时，极化相干系数可以通过RVoG+VTD模型进行表达[13, 14, 15]
式中，w表示极化状态，表征某种特定目标的散射特性；γ(w)表示极化状态w对应的复相干系数，为已知复数； φ0表示地表相位，为未知实数；μ(w)表示有效地体散射幅度比，与极化状态有关，为未知实数；γt表示时间去相干系数，与极化状态无关，为未知实数；γv表示“纯”体去相干系数[13, 14, 15] 式中，θ为平均入射角；σ为植被的平均消光系数；hv为植被高度；kz为垂直向有效波数。
在考虑植被层时间去相干时，并没有改变RVoG模型在复数平面内成线性分布的特性，尽管相干直线的斜率发生了改变，但该相干直线与复数单位圆用于确定地表相位的交点并没有发生改变[13]。在实际应用中，通常获取的复相干系数由于观测噪声的影响，其在复平面上的分布并不位于同一相干直线上(见图 1(a))。利用等权最小二乘直线拟合可以确定相干直线[9, 10]。当采用非线性迭代算法[7]时，其对应的模糊空间为一条线段(见图 1(a))，该模糊空间主要由时间去相干所引起。而当采用三阶段算法[9]时，由于该算法要假定某种极化方式对应的地体散射幅度比为0，加之时间去相干的影响，其模糊空间扩散为“扇形”(见图 1(a)阴影区域)。本文算法就如何利用RVoG+VTD模型削减模糊范围将在本文的第3部分进行详细说明。
3 POLInSAR植被高反演的复数域最小二乘平差法根据文献[16]，不同的极化散射机理，可以得到多种复相干系数(线性极化方式，无限制条件的极化相干最优，有限制条件的极化相干最优)，进而可以得到多个不同的式(1)。按照测量平差中的做法，存在冗余观测时，可以利用平差方法解算待估参数。
3.1 平差准则与定权方案测量平差领域中，处理对象多为实数，还没有完整的理论体系来阐述复数域平差。文献[17, 18, 19]阐述了两种复数最小二乘平差准则：一种是将复数的实部、虚部拆分，分别进行实数域最小二乘处理；另一种是以复数观测值残差的模的平方和最小作为平差准则，同时考虑复数观测值的实部和虚部误差，整体求解得到复数最小二乘解。文献[12]已经证明，第2种准则比第1种准则更优。故本文采用第2种平差准则如下
式中，Re()、Im()表示取复数实部、虚部算子；
(w)表示γ(w)的估计值。对该平差准则进一步分析
可知，该平差准则的实质是复数实部、虚部联合平差。对于随机模型，文献[12]中提出利用相干幅度标准差(Cramer-Rao边界)对复数的模进行定权。同时假定实部、虚部为独立等精度观测量，这与文献[9]中直线拟合过程提出的策略是一致的。本文采用相同的定权方案对观测量进行定权。
3.2 函数模型及线性化方法
由式(1)知平差数学模型为非线性模型，但是如果将γtγv简化成a+bi，可以看到式(1)的形式大大简化，易于进行线性化处理。将eiφ0变换为eiφ0=cos φ0+isin φ0，代入式(1)，将实部、虚部拆分并分别进行泰勒级数展开
式中，下标1、2代表实部、虚部；
分别为a、b、uw、φ0的近似值；
为对应待估参数近似值的改正数。本文采用γHH、γHV、γVV、γHH－VV、γHH+VV、MCD (相干最大分离算法：opt1、opt2、opt3)与PD (相位分离算法：PDhigh、PDlow)[16]10种极化方式对应的复相干系数作为观测值，可以得到20个观测量(每个极化方式组合对应的相干系数有实部和虚部)、13个待估参数(φ0与γv含有2个未知参数a、b)、10个地体散射幅度比，多余观测数为7。对每种极化方式分别按照式(5)列出误差方程，用矩阵形式表示如下
式中，B为设计矩阵；X为待求参数；L为常数矩阵；V为残差矩阵。
3.3 误差方程的解算
式(6)中，设计矩阵B可以归结为稀疏矩阵。为了克服法方程系数阵(BTPB)极易出现的病态问题，本文采用截断奇异值分解法直接对误差方程进行求解，具体公式参见文献[20]。该方法优点在于对于良态、病态、秩亏的方程均适用[20, 21]。此外，该方法可在一定程度上有效抑制由观测值较差的几何结构性所引起的病态问题，如对于稀疏及低矮植被区，极化干涉表现出相似的相干性，在复数平面内表现出较差的几何结构性。在进行平差计算之前，采用DEM差分法进行确定系数矩阵B的初值[12, 16]。由于待估参数的初值较为粗略，故采用迭代平差法计算。
3.4 顾及时间去相关的植被高度反演通过解算得到的a、b，即可得到“纯”体去相干系数γv与时间去相干系数γt的乘积。由于γtγv的解算过程没有假定有效地体散射幅度比为0[9]，故经过有效地体散射幅度比的补偿，使得模糊区域由扇形(图 1(a))退化为一条线段(图 1(b))。对于图 1(b)中的模糊区域主要由时间去相干的不确定所引起[13]。本文采用的RVoG+VTD模型，可使模糊区域由线段退化为一个点(图 1(b))。在此过程中，γtγv=a+bi，该等式中含有实部、虚部2个观测量，但是含有3个未知参数γt、h、σ，无法对参数进行直接解算。为了解决该问题，可以根据先验信息将消光系数σ设为定值[8]。已有研究表明，该方法可以成功地用于植被参数的反演[4, 5]。为了说明该问题，假设kz=0.3、σ=0.15 dB、h=10 m、γt=0.3∶0.1∶0.9做模拟试验，其中设定σ=0.1∶0.01∶0.2为σ的估值。研究h的估值误差与σ、γt的变化关系(见图 2)。由模拟试验可以看出，σ30%的相对误差引起h最大相对误差为4%(见图 2中实线)，这对树高反演是可以接受的。而时间去相干则引起较大的h变化，相对误差最大达到34%(见图 2中短虚线)。由此可见在能获得研究区相对可靠的σ均值估值的条件下，可以将其固定，而考虑时间去相干是必要的。根据[7, 21]研究成果，令σ=0.2 dB。最后，建立二维查找表，以理论计算值与真实解算值差异最小为约束条件，进而得到植被高度h，相应的约束条件如下
|
| 图 1 RVoG模型的几何解释 Fig. 1 Geometric explanation of RVoG model |
|
| 图 2 植被高度估计值随消光系数、时间失相关变化图 Fig. 2 The figure of tree height estimations change with extinction and temporal decorrelation |
综上所述，本文算法整体流程如图 3所示。
|
| 图 3 算法流程 Fig. 3 Flow chart of the algorithm |
为了验证CLSA法的反演性能，本文使用覆盖德国Oberpfaffenhofen 地区E-SAR L波段机载全极化数据进行试验(图 4)。图 4(a)①为试验区数据光学影像(下载自谷歌地球)，红色矩形为研究区域；图 4(a)②为极化数据的Pauli基彩色合成图；图 4(a)③为德国行政区划图，红点为研究区域地理位置。干涉基线为10 m，平均入射角为40°，空间分辨率为3 m×3 m。图 4(b)为不同极化方式的相干性统计图。由统计图可见不同极化方式相干性均较低，可以推断研究区受时间去相干影响较为严重，按照3.4节分析，有必要考虑时间去相干的影响。对照光学遥感影像，该区域主要地物为森林、居民地、农田、机场等。在数据处理之前，对非植被区进行掩膜处理。虚线为用于下文分析的剖面线，其在图中的位置与文献[21]中的树高剖面线位置较为相近，有一定参照价值。由于缺乏研究区域地面实测数据，难以对本文反演的结果进行定量精度评定。为了验证本方法的优越性，本文同时给出非线性迭代及三阶段算法的反演结果，用于对比分析。算法特性对比见表 1。
|
| 图 4 试验区范围及其相干性统计图 Fig. 4 The range of the test area and its histograms of coherence |
| 算法名称 | 初值确定方法 | 输入参数 | 输出参数 |
| 非线性迭代算法 | DEM差分法 | opt | 3个 | uw | φ0 | h | σ |
| 三阶段算法 | 不需要 | 线性极化线性极化+相干最优 | φ0 | h | σ | ||
| 新方法 | DEM差分法 | 线性极化+相干最优+消光系数 | 11个 | uw | φ0 | h | σ |
图 5(a)、(b)、(c)分别为非线性迭代法、三阶段法和本文方法反演的地表相位图。图 5(d)为沿图 4中虚线段所做的剖面对比图。结果表明：相比之下，非线性迭代法与三阶段法反演的地表相位波动较大，其中三阶段法反演的地表相位跳变最大。本文方法(图 5(c))算得的地表相位整体变化趋势相比已有方法变化更为平缓，主要集中在-0.5~0.5 rad之间，这与实际地表相对平坦这一事实更为相符[7]。理论上，根据RVoG+VTD模型的定义，3种方法反演的地表相位应该相等[13]。本文方法优于其他两种方法的原因在于：①本文方法考虑了观测量的先验统计误差，通过平差方法降低了观测噪声对参数估计的影响；②采用修正奇异值分解法可以一定程度上克服由较差的观测量几何结构性所引起的病态问题。
|
| 图 5 不同方法地表相位反演结果 Fig. 5 The ground phase inversion results of different methods |
图 6为非线性迭代算法、三阶段算法和本文方法反演的植被高度图。整体来看，3种方法反演的结果差别较大。非线性迭代算法在低矮植被覆盖区呈现高估，在高植被区，出现低估，植被高主要集中在5~10 m之间。分析原因在于初值的难以确定[9, 10]、高度与消光系数的模糊问题[7]及时间去相干的影响[13, 14, 15]，这些因素都会不同程度地限制参数的解算精度。通过试验发现，初值的不确定性是造成结果出现偏差的主要原因。三阶段算法出现明显高估，主要集中在25~30 m之间。分析原因在于体去相干系数误差较大[9]及时间去相干[13, 14, 15]使得结果出现明显偏差。本文方法反演的植被高度主要集中15~25 m之间，与文献[7, 21]中反演的结果及实测数据范围较为一致，其结果明显优于其他两种算法。进一步对比，对照文献[21]中树高反演结果的剖面，植被的变化范围为：低矮植被覆盖区0~5 m；高植被区15~25 m。对照图 6(d)，无论是低矮植被区还是高植被区，本文方法反演的剖面线与文献[21]中的剖面线变化范围及趋势更为相近。进一步说明了本文方法的优越性。原因在于本文方法获得更为可靠的体去相干系数，有效改善了已有算法中树高解算时出现的明显的高估与低估问题。此外，采用复数平差的方法有效抑制了观测误差，提高了参数的估算精度。在运算效率方面，使用Pentium(R)双核处理器(主频2 GHz)及2 GB内存的计算机对该研究区域的数据进行解算，3种方法的运算时间分别为：7440 s、1050 s、1210 s。本文方法运算效率明显优于非线性迭代算法，略低于三阶段算法。
|
| 图 6 不同方法的植被高度反演结果 Fig. 6 The vegetation height inversion results by different methods |
(1) 本文提出的RVoG+VTD线性化的方法及解算方法，尽管没有对非线性模型的非线性强度做定量分析后进行线性化处理，但通过试验结果的分析表明方法是行之有效的。
(2) 真实数据试验表明本文算法相比已有算法在反演精度上具有明显优越性。在时间效率方面，略低于三阶段算法，明显高于非线性迭代算法。由反演的地表相位结果来看，植被层时间去相干对RVoG模型确定地表相位影响较小，这为今后利用星载全极化数据提取森林覆盖区域的DTM提供了新思路。在反演树高方面，本文提出的方法反演的结果更贴合实际。这种顾及时间去相干的模型对于即将发射的ALOS-2 PALSAR数据提取植被参数具有重要意义。
(3) 本文提出的方法还需要进一步完善，未来工作中将研究复数的实部、虚部如何进行合理定权；本文提出的线性化平差模型为参数解算的精度评定提供了可能，这对于定量评价参数反演的效果具有重要研究价值，对于建立植被覆盖区高精度DTM及林业普查具有重要意义；将该方法扩展到多基线，大大增加了观测量，可以避免对消光系数的固定。总之，本文提出的植被高度反演的复数最小二乘法，将平差理论引入到复数数据处理中，这为PolInSAR数据处理及平差理论的完善提供了一种新的思路。
致谢：感谢欧洲空间局(ESA)免费发布的POLSARpro 4.2软件及德国宇航局发布的E-SAR数据。
| [1] | ZHANG Qiyong, CEN Minyi, ZHOU Guoqing, et al. Extracting Trees from LiDAR Data in Urban Region[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(4): 330-336. (张齐勇,岑敏仪, 周国清, 等. 城区LiDAR点云数据的树木提取[J]. 测绘学报, 2009, 38(4): 330-336.) |
| [2] | NEUMANN M. Estimation of Forest Structure, Ground, and Canopy Layer Characteristics from Multibaseline Polarimetric Interferometric SAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(3): 1086-1103. |
| [3] | TAN Lulu, YANG Libo, YANG Ruliang. Investigation of Tree Height Retrieval with Polarimetric SAR Interferometry Based on ESPRIT Algorithm[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(3): 296-300. (谈璐璐, 杨立波, 杨汝良. 基于ESPRIT算法的极化干涉SAR植被高度反演研究[J]. 测绘学报, 2011, 40(3): 296-300.) |
| [4] | GARESTIER F, FERNANDEZ D, PAPATHANASSIOU K P. Pine Forest Height Inversion Using Single-pass X-band PolInSAR Data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(1): 59-68. |
| [5] | TREUHAFT R N, MADSEN S N, MOGHADDAM M, et al. Vegetation Characteristics and Underlying Topography from Interferometric Radar[J]. Radio Science, 1996, 31(6): 1449-1485. |
| [6] | TREUHAFT R N, SIQUEIRA P R. Vertical Structure of Vegetated Land Surfaces from Interferometric and Polarimetric Radar [J]. Radio Science, 2000, 35(1): 141-177. |
| [7] | PAPATHANASSIOU K P, CLOUDE S R. Single-baseline Polarimetric SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(11): 2352-2362. |
| [8] | LI X W, GUO H D, LIAO J J, et al. Inversion of Vegetation Parameters Using Spaceborne Polarimetric SAR Interferometry[J]. Journal of Remote Sensing, 2002, 6(6): 424-429. |
| [9] | CLOUDE S R, PAPATHANASSIOU K P. Three-stage Inversion Process for Polarimetric SAR Interferometry[J]. IEEE Proceedings on Radar, Sonar and Navigation, 2003, 1503(3): 125-134. |
| [10] | CHEN E X, LI Z Y, PANG Y, et al. Polarimetric Synthetic Aperture Radar Interferometry Based Mean Tree Height Extraction Technique[J]. Scientia Silvae Sinicae, 2007, 43(4): 66-70. |
| [11] | LUO H M, CHEN E X, CHENG J, et al. Forest Height Estimation Methods Using Polarimetric SAR Interferometry[J]. Journal of Remote Sensing, 2010, 14(4): 806-813. |
| [12] | ZHU Jianjun, XIE Qinghua, ZUO Tingying, et al. Criterion of Complex Least Squares Adjustment and Its Application in Tree Height Inversion with PolInSAR Data[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(1): 41-51. (朱建军,解清华,左廷英,等. 复数域最小二乘平差及其在POLInSAR植被高反演中的应用[J]. 测绘学报, 2014, 43(1):41-51.) |
| [13] | PAPATHANASSIOU K P, CLOUDE S R. The Effect of Temporal Decorrelation on the Inversion of Forest Parameters from PolInSAR Data[C]// Proceedings of IEEE International Geoscience and Remote Sensing Symposium: 3. New York: IEEE, 2003: 1429-1431. |
| [14] | SEUNG-KUK L, KUGLER F, PAPATHANASSIOU K P, et al. Quantifying Temporal Decorrelation over Boreal Forest at L-and P-band[C]// Proceedings of the 7th European Conference on Synthetic Aperture Radar(EUSAR). Friedrichshafen: VDE, 2008: 1-4. |
| [15] | LAVALLE M, SIMARD M, HENSLEY S. A Temporal Decorrelation Model for Polarimetric Radar Interferometers[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(7): 2880-2888. |
| [16] | CLOUDE S R. Polarisation: Applications in Remote Sensing[M]. London: Oxford University Press, 2009. |
| [17] | GU Xiangqian, KANG Hongwen, CAO Hongxing. The Least-square Method in Complex Number Domain[J]. Progress in Natural Science, 2006, 16(1): 49-54. (谷湘潜, 康红文, 曹洪兴. 复数域内的最小二乘法[J]. 自然科学进展, 2006, 16(1): 49-54.) |
| [18] | LI Mengxia, CHEN Zhong. The Modification of Least Square Method(LSM) in the Complex Field[J]. Journal of Yangtze University: Natural Science Edition, 2008, 5(3): 7-8. (李梦霞, 陈忠. 复数域内最小二乘法的一种改进[J]. 长江大学学报: 自然科学版, 2008, 5(3): 7-8.) |
| [19] | MILLER K S. Complex Linear Least Squares[J]. SIAM Review, 1973, 15(4): 706-726. |
| [20] | WANG Zhenjie. Research on the Regularization Solutions of Ill-posed Problems in Geodesy[D]. Wuhan: Institute of Geodesy and Geophysics of Chinese Academy of Sciences, 2003. (王振杰. 大地测量中不适定问题的正则化解法研究[D]. 武汉: 中国科学院测量与地球物理研究所, 2003. |
| [21] | PAPATHANASSIOU K P, CLOUDE S R, REIGBER A, et al. Multi-baseline Polarimetric SAR Interferometry for Vegetation Parameters Estimation[C]// Proceedings of IEEE 2000 International Geoscience and Remote Sensing Symposium. Honolulu: IEEE, 2000, 6: 2762-2764. |