本文内容来源于《测绘学报》2020年第10期，审图号GS（2020）5551号。解清华1 1.中国地质大学(武汉)地理与信息工程学院, 湖北 武汉 430074; 基金项目：国家自然科学基金（41804004；41820104005；41531068；41904004）；中国地质大学（武汉）中央高校基本科研业务费专项（CUG190633） 摘要：针对经典的PolInSAR森林高度三阶段几何反演算法在单基线条件容易受到地体幅度比假设以及地形坡度影响的问题，从测量平差角度提出了基于S-RVoG模型的PolInSAR非线性复数最小二乘森林高度反演算法。该算法不再需要假设某一个极化通道地体幅度比为零，且采用考虑地形坡度影响的S-RVoG模型作为平差模型。为了验证算法，本文采用欧空局BioSAR2008项目提供的3景P波段极化干涉SAR数据进行两组单基线森林高度反演试验。结果表明，在单基线条件下，基于RVoG模型的非线性复数最小二乘算法反演结果优于三阶段几何反演算法，而基于S-RVoG模型的非线性复数最小二乘算法进一步提高反演精度，对于坡度较大区域（坡度>10°），精度平均提高了18.48%。 关键词：极化干涉SAR 森林高度 地形坡度 S-RVoG模型 复数最小二乘   森林高度作为重要的森林特征参数之一，可以广泛应用于森林生物量估计、生态建模、森林管理、全球碳循环和气候变化研究等林业及环境相关领域[1-4]。极化干涉SAR(polarimetric interferometric synthetic aperture radar, PolInSAR)技术有效融合了雷达极化技术和雷达干涉测量技术的优点，在森林区域可以有效区分森林中植被冠层和地表层的散射贡献和高度，已经被证明是一项极具潜力的区域乃至全球尺度森林高度测绘技术手段[5-7]。 目前，森林参数反演中最常用的极化干涉散射模型是随机地体二层相干散射模型——RVoG(random volume over ground)模型[7-9]。该模型假设植被散射场景由地表层和地表上一定厚度的植被层构成。地表层不区分表面散射分量和地面-树干结构引起的二面角散射分量，并假定地表为平地，其散射中心固定在地表某一参考高度处。植被层假定由随机取向服从均匀分布的散射粒子组成，后向散射能量服从指数分布，并假设植被垂直结构为均匀同质，即微波信号在植被层传播过程中消光系数不变且与垂直位置无关。RVoG模型最终表达为极化复相干性与森林参数之间的函数关系。目前，三阶段几何反演算法[10]是基于RVoG模型最常用的森林高度反演算法，基于不同波段PolInSAR数据，已经成功地应用到多个试验区的森林高度反演[11-20]。然而该算法一方面需要假设某个极化通道对应的地体幅度比为零；另一方面，该算法测试的试验区通常集中在相对平坦的区域，主要是因为RVoG模型在建模时假设地表为平地。然而，从实际大范围应用和未来全球森林高度测图角度分析，森林覆盖区通常是在一些不平坦区域，甚至地形起伏较大区域。显然，利用RVoG模型在这些区域进行森林高度反演会存在地形坡度引起的误差。 针对三阶段几何反演算法中的地体幅度比假设限制，文献[21]首次从测量平差角度提出了一种基于RVoG模型的非线性复数最小二乘算法，根据观测值先验统计信息定权，建立平差函数模型、随机模型及平差准则，然后利用非线性最小二乘迭代反演森林高度。该方法反演过程简单直观，不需要分步计算和地体幅度比假设，然而该方法只采用模拟PolInSAR数据验证，且散射模型没有考虑地形坡度影响，从而限制了该算法在实际应用过程中的适用性。针对三阶段几何反演算法中的地形坡度的影响，文献[22]在RVoG模型基础上引入了距离向地形坡度参数，提出了改进的散射模型，即S-RVoG(sloped RVoG)模型[22]。文献[15]在此基础上，利用真实PolInSAR数据分析验证了S-RVoG模型的反演地形影响补偿的有效性，同时基于PD(phase diversity)相干最优化思想得到地体幅度比小的最优极化方式的相干性，从而削弱地体幅度比假设影响。文献[16]在此基础上提出了基于S-RVoG模型的双基线几何反演算法，不再需要地体幅度比假设，然而一方面要求有较大的基线比的双基线PolInSAR数据，另一方面由于基线观测误差的存在，一条基线误差较大将极大地影响反演精度。文献[23]针对基线问题进一步提出基于LiDAR数据辅助的双基线几何反演算法，基于支持向量机监督分类思想从多基线数据中进行双基线优化选择，从而得到更高精度的反演结果。上述分析表明，考虑坡度的S-RVoG模型目前在实际森林高度反演中仅应用于几何反演算法，且在单基线条件下受到地体幅度比假设限制。 针对这两方面研究的不足，本文提出基于S-RVoG模型的非线性复数最小二乘森林高度反演算法。该算法可以有效克服单基线几何反演算法的缺陷，一方面采用S-RVoG模型作为平差模型，补偿地形坡度影响；另一方面参数求解基本思想仍采用非线性复数最小二乘方法，从而不再需要假设某一个极化通道地体幅度比为零，且不需分步解算。为了验证算法的有效性，本文采用欧空局BioSAR2008项目提供的3景P波段极化干涉SAR数据进行森林高度反演试验，分别测试两组单基线配置下本文算法表现，同时测试三阶段几何反演算法和基于RVoG模型的非线性复数最小二乘算法用于方法对比验证和分析。最后，本文利用LiDAR数据获取的森林高度产品验证和对比不同算法获取的森林高度反演的精度。 1.1 三阶段几何反演算法 1.1.1 RVoG模型 文献[6]提出的随机地体两层植被相干散射模型——RVoG(random volume over ground)模型，将森林植被散射场景简化为地表层和植被层组成的两层散射问题。第1层考虑地表层贡献的地表散射分量，包含表面散射分量和地面-树干结构引起的二面角散射分量，并假定地表层为平底，其散射中心固定在z0处。第2层考虑位于地表层上面的一定厚度hv的植被层贡献的体散射分量，并假定植被体由随机取向的散射粒子组成，且整个植被冠层覆盖范围为从z0到z0+hv处。RVoG模型的最终表达式为[7-9]
式中，ϕ0代表地表相位，可以用于估计林下地形；μ(w)是地体幅度比，与极化有关；γV表示只有植被散射引起的相干性，又称为纯体相干性，与极化无关，依赖于森林高度hv和消光系数σ，对应表达式为[7-9]
式中，垂直波数
式中， 经过等价变换后，RVoG模型表达式实质为一个复线性模型，其几何特性表现为不同极化通道的复相干性在复平面上服从直线分布[10]。 1.1.2 反演算法流程 文献[10]根据RVoG模型的几何特性提出了三阶段算法，是目前最常用的基于RVoG模型的PolInSAR森林高度反演方法[10]。它的反演过程主要包括以下几个步骤。 (1) 直线拟合。基于RVoG模型的几何特性，选取一定数量极化通道对应的复相干性在复平面上直线拟合，利用总体最小二乘思想求取对应拟合直线参数。本文试验除了选取常用的线性极化和Pauli基极化方式外，即HH、VV、HV、HH+VV、HH-VV共5个极化通道的复相干性观测值，为提高直线拟合精度，另外增加目标函数为相位差异最大对应的PD(phase diversity)相干最优算法[24]获取的两种复相干性观测值γPDhigh和γPDlow。 (2) 地表相位ϕ0估计。拟合的直线与复平面上的单位圆有两个交点，其中地体幅度比最大的那个交点对应的相位即为地表相位。根据相干性与地面点的相对位置排序可以判断出正确的地表相位点。由于PD相干最优化算法满足相位中心最大分离的特点，一般γPDhigh的相位中心相比γPDlow离地表更远，因此本文采用的判别准则为
(3) 森林高度hv估计。假设某个极化通道完全没有地表散射贡献，对应的地体幅度比为零。由RVoG模型表达式可以看出，此时该极化通道复相干性对应于γVejϕ0，由于地表相位已知，剩余两个参数(植被高hv和消光系数σ)可以采用查找表方法求解。假设γPDhigh完全来自植被层散射贡献，相应的约束条件如下
1.2 基于RVoG模型的森林高度非线性复数最小二乘反演算法 文献[21]介绍了复数域中数据处理的最小二乘方法，并针对极化干涉SAR植被高度反演中观测值为复数的特点，从测量平差角度提出了基于RVoG模型的非线性复数最小二乘森林高度反演算法(RVoG-nonlinear complex least squares，RVoG-NCLS)，分别建立复数观测值条件下的平差函数模型、随机模型及平差准则，然后利用非线性最小二乘迭代反演森林高度。该方法不再需要假设某一个极化通道地体幅度比为零，克服了单基线三阶段几何反演算法的缺陷。 1.2.1 函数模型 关联不同极化通道复相干性观测值和森林参数的RVoG模型即为平差的函数模型，则森林高度求解可以看成是一个间接平差问题，对应函数模型可以表示为[21, 25]
式中，L表示不同极化通道复相干性观测向量；V表示改正数向量；F表示RVoG模型对应的函数关系，^表示附加改正数；N表示复相干性观测数量。由式(6)可知，只有部分参数与极化有关，随着极化通道增加，总实数观测量个数为2N，实数未知数个数为N+3，故为当观测量N>3时，存在多余观测。 1.2.2 随机模型 由极化SAR观测值统计特性可知，复相干性的模的中误差表达式为[26]
式中，A代表相干性估计窗口内像元总数。假定不同极化观测值之间互相独立，则任意极化通道观测值的权定义为[21, 25]
1.2.3 平差准则 利用建立的函数模型和随机模型，选取残差模的平方和最小作为平差准则构建复数最小二乘算法进行植被参数求解。由于观测模型为非线性模型，且模型的表达式十分复杂，线性化过程困难，因而采用非线性迭代法求解。具体地，对于N个复相干性极化观测量，则极化干涉SAR植被高反演的复数最小二乘准则表达式为[21, 25]
式中，[M]表示RVoG模型。为了降低参数反演的难度同时提高求解稳健性，一方面避免求解过多参数，另一方面尽量选择相关性小的极化通道观测值，本文在利用复数域最小二乘反演植被高时，选取了(phase diversity, PD)相干最优化[24]极化方式和C & P相干最优化[6]极化方式，即PDHigh、PDLow、Opt1、Opt2、Opt3，相应参数求解问题为基于10个实数观测值求解8个实数未知数。 1.3 基于S-RVoG模型的森林高度非线性复数最小二乘反演算法 1.3.1 S-RVoG模型 文献[22]针对距离向地形坡度影响，将RVoG模型进行了扩展，提出了考虑地形坡度的散射模型，即S-RVoG(Slope-RVoG)模型[22]。图 1为S-RVoG模型几何示意图，从中可以看出S-RVoG模型假设距离向上存在一个坡度角α，当地形坡度倾斜朝向雷达传感器时定义坡度为正，反之定义坡度为负。由于地面散射体沿着斜坡分布，沿着z轴描述其能量变化比沿着z′轴要更加复杂，故该模型是在沿斜坡方向和垂直斜坡方向的坐标系y′O′z′中建模，然后再参照RVoG模型的推导过程，可得S-RVoG模型的最终表达式为[22]
式中，垂直波数k′z=kzsinθ/sinθ′；入射角θ′=θ-α；γ′V和μ′(w)分别表示新坐标系统下的纯体相干性和地体幅度比。 1.3.2 反演算法流程 在1.2节算法基础上，为考虑地形坡度影响，本文提出基于S-RVoG模型的非线性复数最小二乘森林高度反演算法(S-RVoG-NCLS)。此时相应的非线性复数最小二乘森林高度反演算法需要更换测量平差的函数模型，而随机模型和平差准则保持不变。具体地，函数模型从式(1)表达的RVoG模型变为式(10)表达的S-RVoG模型。该模型可以仍然是以多极化复相干性为观测量，以森林高度，消光系数等森林参数作为未知数的观测方程。参数求解仍然采用PDHigh、PDLow、Opt1、Opt2、Opt3这5个最优极化方式的复相干性观测值，利用非线性迭代法求解。模型中的距离向坡度参数由外部DEM获取的地形坡度信息提供。 2.1 试验数据介绍与预处理 本文用于算法验证的试验数据采用欧空局BioSAR2008项目提供的E-SAR系统获取的瑞典北部Krycklan森林区域(64°14′N, 19°46′E)的3景机载P波段全极化SAR数据，其对应的主要参数表 1所示。该试验区是一片典型的北方森林区域，其中一部分长期被瑞典农业科学大学管理。该区域的主要树种类型为针叶林(苏格兰松和挪威云杉)及少部分的桦树。森林高度平均测量值约18 m, 其中最大测量高度值为30 m。该地区是一个中等地形坡度变化的丘陵区，高程变化在高于平均海平面150~380 m之间变化。覆盖试验区的1 m格网高精度DEM用于生成坡度信息。LiDAR观测获得森林高度产品用于验证森林高度反演结果。 表 1 试验区P波段机载全极化干涉SAR数据参数Tab. 1 Parameters of P-band airborne PolInSAR data over the study area
表选项 基于这3景极化干涉SAR数据，以影像1为主影像，影像2和影像3分别为辅影像，可以组成两个极化干涉对(影像1-2及影像1-3)。针对每个干涉对，数据预处理部分首先进行一系列极化干涉处理，包括影像配准、距离向频谱滤波、干涉、平地相位去除、多视处理。然后，采用11×11窗口获取不同极化方式(HH、HV、VV、HH+VV、HH-VV、PDHigh、PDLow、Opt1、Opt2和Opt3)对应的复相干性。 2.2 试验结果与讨论 基于两个极化干涉对获取的不同极化方式的复相干性观测值，针对每个干涉对，本文分别测试3种森林高度反演算法：三阶段几何反演算法，基于RVoG模型和基于S-RVoG模型的非线性复数最小二乘反演算法。在非线性优化计算过程中，为了兼顾计算效率和参数精度，本文试验中迭代的终止条件设定包括：迭代次数最大值为20，函数可容忍误差为0.01，参数可容忍误差为0.1。在计算量方面，非线性复数最小二乘平差算法，采用RVoG模型和S-RVoG模型都需要耗时大的非线性优化求解过程，两者的计算量相当，但要显著大于三阶段几何反演算法。图 2显示了最终的森林高度反演结果。具体的，图 2(a)和(b)显示了三阶段几何反演算法结果，图 2(c)和(d)显示了基于RVoG模型的非线性复数最小二乘反演算法结果，图 2(e)和(f)显示了基于S-RVoG模型的非线性复数最小二乘反演算法结果。从图 2可以看出，6个结果得到的森林高度空间分布趋势相似，但是基于RVoG模型和基于S-RVoG模型的非线性最小二乘算法相比三阶段几何反演得到的森林高度估计值更高(颜色更深)；基于S-RVoG模型的非线性最小二乘算法反演的结果与基于RVoG模型算法结果在局部有些差异，如图 2中红色矩形框标注的区域。 为了对比差异，图 3(a)和(b)展示了基于S-RVoG模型与基于RVoG模型的非线性复数最小二乘算法反演的森林高度差值，颜色比例限制在-5~5 m之间，在一些区域(图 3中红色矩形框内)可以明显看出差异。图 3(c)展示了根据LiDAR获取的DEM计算的试验区距离向地形坡度结果。通过对比，可以看出森林高度差值与地形坡度存在明显的相关性。当坡度为正时，改进算法得到森林高度结果相比原有算法结果更低。这是由于正的坡度会导致森林高度高估，而S-RVoG模型可以补偿坡度影响。反之，负的坡度会导致低估，S-RVoG模型同样可以补偿坡度影响。
由于LiDAR森林高度产品精度较高，常常作为参考用于验证极化干涉SAR森林高度反演结果精度。为了方便验证，首先将LiDAR森林高度产品经过地理编码处理转换到多视后的SAR成像几何坐标系下(方位向×距离向：6472× 1501像素)。为了进行样地交叉验证，本文采用一种基于窗口的方法从LiDAR森林高度产品中均匀选择森林样地。首先，一个均匀格网分布(方位向×距离向：30×15像素)被应用到LiDAR森林高度图用于选择森林样地的中心像素。接着，一个大小为51×51(即2601像素)的平均窗口应用到每个中心像素用于选择周围像素。然后，去除那些存在空值(即非森林像素)的样地。最后，共有362个森林样地被选择，相继计算每个样地对应的平均森林高度值用于算法反演结果精度验证。图 4给出了6种森林高度反演结果与LiDAR森林高度的样地交叉验证结果。均方根误差(RMSE)和相关系数(R2)被选择用于结果统计分析。图中每个点的颜色代表样地的平均距离向坡度，比例压缩在-15°~15°之间。图 4(a)和(b)可以看出，三阶段几何反演算法结果对应的RMSE都偏高，分别是7.88 m和8.63 m；对应的相关系数偏低，分别是0.26和0.33。精度偏低的原因是P波段数据所有极化通道都有显著的地面散射贡献，而假定的纯体相干性极化(PDHigh)也不能满足要求。图 4(c)和(d)可以看出单基线条件下非线性复数最小二乘算法森林高度反演结果相比几何反演算法明显改善，RMSE分别是5.68 m和6.05 m，精度平均提高了28.95%，对应的相关系数也提高到了0.45和0.54，平均提高了67.8%。这表明非线性复数最小二乘算法相比几何反演算法可以有效补偿地体幅度比假设引起的偏差，获取更高的森林高度反演精度。从图 4(e)和(f)可以看出采用S-RVoG模型后，森林高度反演精度进一步提高，其对应的RMSE平均降低了4.43%，相关系数平均增加8.08%。精度改善较少的原因是所有森林样地的平均坡度不大，大约3.54°。为了更好地展示地形坡度校正的效果，以坡度大于10°为标准从所有样地中筛选出23个森林样地结果去比较森林高度精度。图 5显示了23个坡度大于10°的样地的交叉验证结果，可以看出考虑坡度校正的复数非线性最小二乘算法对应的RMSE平均降低了18.48%，相关系数平均增长88.68%。
针对经典的三阶段几何反演算法在单基线条件容易受到地体幅度比假设和地形坡度影响的问题，本文从测量平差角度提出了基于S-RVoG模型的PolInSAR非线性复数最小二乘森林高度反演算法。该算法一方面采用非线性复数最小二乘解算思想，不再需要假设某一个极化通道地体幅度比为零；另一方面采用S-RVoG模型作为平差模型，能够补偿地形坡度影响，从而有效弥补了单基线几何反演算法在这两方面的缺陷。本文采用欧空局BioSAR2008项目提供的3景P波段瑞典森林地区极化干涉SAR数据进行两组森林高度反演试验。试验结果表明，在单基线条件下，基于RVoG模型的非线性复数最小二乘森林高度反演算法优于几何反演算法结果，森林高度反演精度平均提高了28.95%；基于S-RVoG模型的非线性复数最小二乘算法可以进一步提高森林高度反演精度，对于坡度较大区域(坡度>10°)，精度提高了18.48%。这些试验结果有效验证了算法的有效性，并体现了算法在地形起伏区域森林高度反演的潜力。 第一作者简介：解清华(1989-), 男, 博士, 副教授, 研究方向为PolSAR/PolInSAR植被生物物理参数反演及林下地形测绘。 |
|
|
,朱建军2
(1)
(2)
，依赖于成像几何(垂直基线B⊥，斜距R和入射角θ)和雷达波长λ。经过数学变换，式(2)表达的RVoG模型可以等价为[10]
(3)
。
