1.   引言

干涉合成孔径雷达(Interferometric Synthetic Aperture Radar, InSAR)是获取全球数字高程模型(Digital Elevation Model, DEM)的一种有效方法^[1]。DEM包含两种高程模型^[2]，一种是没有任何植被、建筑物或者其他物体的裸露地表的高程，称为数字地形模型(Digital Topography Model, DTM)。另一种是记录任何位置的表面高程，可以是植被区或者建筑区的表面高程，也可以是没有植被或其它物体的地面，称之为数字表面高程(Digital Surface Model, DSM)^[3]。极化干涉合成孔径雷达(Polarimetric Interferometry Synthetic Aperture Radar, PolInSAR)是获取高精度DSM的一种重要手段，对于地形复杂地区地理测绘具有重要意义。

文献[4,5]利用HH极化通道的干涉相位进行干涉处理得到植被区域DSM，在植被密集区域得到的结果与真实DSM较为吻合，但是对于植被较为稀疏的区域结果受植被形态影响较大，其反演误差较大；文献[6]利用机载激光雷达辅助获得植被下DTM，然后基于地上随机体散射模型(Random Volume Over Ground, RVOG)得到植被冠层高度，植被冠层高度与DTM结合得到植被区DSM，这种测量方法能得到较高精度的DSM，然而测量需要辅助激光雷达为前提，测量成本较大，且测量范围有限；文献[7]通过极化相干最优化中的相干性最优化和最大化相位中心分离方法获取极化干涉信息，然后进行干涉处理得到DSM，文中针对大范围山区对不同方法进行了对比分析，但对植被区DSM的进一步提升没有进行深入研究；文献[8]通过联合干涉相位和相干性幅度的最优相干性准则，有效区分了地表散射和森林冠层散射的相干性和散射中心，提高了植被参数反演的可靠性，但没有得到高精度DSM反演结果。以上针对植被区的DSM反演都是提高了电磁波在植被中的等效散射相位中心高度估计精度。而由于电磁波的穿透特性，等效散射相位中心与植被顶部之间仍存在一定距离，导致植被区DSM欠估计。

针对这一问题，本文提出一种植被区高精度DSM估计方法。首先估计出植被区不同极化散射机制下的干涉复相干系数，通过极化相干最优中最大化相位差最优化方法^[7]分离出电磁波在植被上、下两层散射区域的干涉相位。植被上层散射区域相位中心称之为高散射相位中心(High Scattering Phase Center, HSPC)；植被下层散射区域相位中心称之为低散射相位中心(Low Scattering Phase Center, LSPC)。然后提出一种归一化HSPC和归一化LSPC高度随消光系数变化的数学模型。基于该模型估计得到植被区HSPC与植被冠层顶部间距，即电磁波在植被中的最浅穿透深度。最后，在HSPC高度处补偿最浅穿透深度估计值，修正植被区由于散射相位中心欠估计导致的DSM估计值偏低的问题。实验表明该方法能有效提升植被覆盖区DSM反演精度。

2.   PolInSAR相干系数估计

植被区极化干涉复相干系数的相位是电磁波在植被中等效散射相位中心所处位置的重要体现^[9]。为了获取最优相干性，一种有效的方法是最大化复平面单位圆中相干系数相位差。该方法以相干区域内不同散射机制下相位中心为基础，寻找相位中心分离最大的极化散射机制作为最优极化状态，能够最大限度分离电磁波在不同植被层的相位中心。其本质是在极化状态空间中寻找一对散射机制${{{w}} _a}$和${{{w}} _b}$使极化干涉相干系数的相位差最大，即使得两个散射机制对应的HSPC和LSPC的相位差最大。首先建立特征方程，通过最大化特征值之差来获取最大相位差对应的散射机制^[10]：

其中，${{w}}$为任意散射机制；$\phi$是引入的一个旋转角；${{{Ω}} _{\rm{H}} } \!=\! ({{{Ω}} _{12}}{e^{{\rm i}\phi }} + {{Ω}} _{12}^{\rm{H} }{e^{ - {\rm i}\phi }})/2$; ${{T}} \!=\! ({{{T}}_{11}} + {{{T}}_{22}})/2$; ${{{Ω}} _{12}}$, ${{{T}}_{11}}$和${{{T}}_{22}}$为极化干涉相干矩阵^[11,12]。遍历$\phi$，可得到完整的极化干涉相干系数分布区域的边界。此时，即可得到相位差最大的两个边界点的相干系数值。

利用最大化相位差最优化方法进行极化干涉相干系数相位分离虽然得到了电磁波在植被中的HSPC, LSPC干涉相位，然而由于电磁波穿透特性，HSPC的干涉相位所对应的高程并没有达到植被冠层顶部高程。植被中HSPC位置与植被冠层顶部的仍存在一定差距，该差距即为某种散射机制下电磁波在植被中的最浅穿透深度。对于稠密植被区电磁波无法穿透植被到达地表并被反射接收导致植被下地表高程无法获取的情况下，通过HSPC相位干涉处理可以得到粗DSM。在得到的粗DSM中补偿估计的电磁波最浅穿透深度，则可以生成较高精度DSM，对于DSM反演具有重要的应用价值。

3.   高精度DSM反演方法

3.1   归一化散射相位中心变化模型

植被区SAR图像单个分辨单元内往往存在多种散射机制。由于植被分布为各向异性，体散射在SAR图像分辨单元中总是占据主导地位，体散射的存在引起严重去相干。文献[13]提出一种指数型垂直结构函数模型，即RVOG模型。植被区相干系数可表示为

(2)

式中，${\sigma _m} = 2\sigma /\cos \theta$; $\sigma$为消光系数；$\theta$为电磁波平均入射角；$M({{w}})$为由极化状态决定的地体幅度比；${h_v}$为植被体层的厚度(即植被的高度)；${\varphi _{\rm{0}}}$为地面相位；${\rm{i}}$为虚数单位；${\gamma _V\,}$为体散射复相干系数；${k_z}$为距离向谱滤波后的垂直有效波数。在没有地面散射分量、去相干因素为完全体散射去相干的情况下，植被区体散射复相干系数的相位为观测得到的相干系数相位与植被底部地形相位之差$\Delta \varphi$，该相位对应的高度称体散射相位中心(Volume Scattering Phase Center, VSPC)高度。在假设地面高程为零的情况下，VSPC高度对应最大化相位差方法分离得到的HSPC相位反演得到的高程。考虑体散射复相干系数相位与高度关系，将体散射复相干系数幅度和相位分离。对于分离出的体散射复相干系数相位$\Delta \varphi$，得到一个相对高度。

称为归一化VSPC高度。为了便于理解，将体散射复相干系数进行变换后分析两种不同消光系数下归一化VSPC高度的边界情况，一种是消光系数趋于0 dB/m，另一种是消光系数趋于无穷大。

当$\sigma \to 0$ dB/m时，由式(4)中间部分等式可以得到${\gamma _V} = {C_1}\exp \left(\displaystyle\frac{{{\rm{i}}{k_z}{h_v}}}{2}\right)$, ${C_1}$为幅度项。此时，由体散射复相干系数相位知，${P_c}$等于0.5，即VSPC位于植被体层的一半处。当$\sigma \to \infty$时，由式(4)右边部分等式可以得到${\gamma _V} = \exp ({\rm{i}}{k_z}{h_v})$，此时${P_c}$等于1, VSPC位于植被冠层顶部。

在给定垂直有效波束的情况下，根据式(3)可以得到不同植被高度归一化VSPC高度随消光系数变化的关系。可知，不同植被高度下归一化VSPC高度均大于0.5，而且随着消光系数的增大VSPC逐渐向植被冠层顶部移动，植被越高VSPC移动的越快。

由式(3)知，描述归一化VSPC高度与消光系数之间关系的数学表达式较为复杂。为了直观表达二者之间相对关系，利用关于消光系数的高阶多项式进行拟合近似得到归一化VSPC高度随消光系数变化的曲线模型。图1分别给出3阶、5阶和10阶多项式拟合VSPC模型曲线，并给出理论模型与拟合模型间的残差。由图1知，拟合多项式阶次越高，拟合模型越接近理论模型。考虑其极限情况，可得到拟合模型如式(5)：

图 1  不同阶次多项式拟合VSPC模型及残差曲线

下载: 全尺寸图片 幻灯片

其中，${a_0},{a_1}, ·\!·\!· ,{a_{{n}}}$为常数，$n$为阶次，$n \to \infty$。

实际中，式(5)所示无穷阶多项式模型在工程应用中难以实现。由理论归一化VSPC模型变化趋势结合指数函数泰勒级数展开形式，假设归一化VSPC高度与消光系数关系可近似为如式(6)所示指数模型：

其中，$\sigma$为模型中的自变量，即消光系数；$a$, $b$为两个与植被高度有关的常数。简化的归一化VSPC模型取值范围满足$[0.5,1.0)$。为了区分后续所述LSPC，这里将等效后的归一化VSPC高度称归一化HSPC高度。

实际中，考虑归一化HSPC随着消光系数变化时是基于地体幅度比$M({ w})$为零的情况。对于LSPC，当消光系数为零时，地体幅度比$M({ w})$并不为零，即有部分电磁波能量到达地面。LSPC在消光系数为零时接近地面，即其归一化散射相位中心高度为0。而当消光系数无穷大时，电磁波无法穿透植被，此时植被中电磁波散射与极化状态无关，HSPC和LSPC相互重叠，归一化LSPC高度为1。也就是说，随着消光系数的增大，分离的两个散射相位中心逐渐靠近，并逐渐趋近于冠层顶部。因此，两种散射相位中心随消光系数增大有相似的变化趋势，但是由于取值范围不同，变化的快慢有区别。目前的研究中植被区的地体幅度比无法直接获取，分离的LSPC高度随消光系数变化关系也没有明确的模型。根据两种散射相位中心变化有相同趋势，同样假设归一化LSPC变化满足指数型变化模型。模型中固定植被高度下对应常数相同，只对指数项前的系数稍作修正使其满足归一化LSPC高度变化取值范围为$[0,1)$的情况。依据式(6)可得归一化LSPC高度变化模型如式(7)：

其中，$a$, $b$为给定植被高度下对应与归一化HSPC高度模型相同的常数。

对于构建的两个指数模型，关键在于如何准确求取模型中的常数$a$, $b$。为了使构建的归一化HSPC指数模型以较小的误差逼近归一化VSPC高度，建立如式(8)所示代价函数：

其中，angle(·)为取复数相位操作。式(8)代价函数的本质是一个多元非线性优化函数，消光系数范围为$\sigma \in [0, + \infty )$。由不同树高下归一化VSPC高度随消光系数变化关系知，对于常见的森林植被体中(植被高度20 m左右)当消光系数达到5 dB/m以上时VSPC几乎到达树冠顶部，此时相位中心高度达到饱和，继续增大消光系数，VSPC变化很小。因此，在求解代价函数时将消光系数范围设为$\sigma \in [0,5]$ dB/m，为了提高求解效率，可将消光系数以一定间隔离散化，重新构建代价函数如式(9)：

其中，$N$为消光系数的离散采样点数。通过优化^[14]式(9)最终可得到使指数模型近似逼近归一化VSPC模型的最优$a$, $b$值，从而得到归一化HSPC和归一化LSPC高度随消光系数变化模型。根据表1给出的雷达几何参数，由式(9)计算归一化HSPC在不同植被高度下的中误差(Root Mean Square Error, RMSE)，由此衡量拟合模型偏离理论模型的离散度。如图2所示，植被高度在1～30 m时，归一化HSPC模型的最大中误差为0.014，误差不超过2%，可见建立的指数模型引入的误差对DSM反演精度影响很小。

表 1  雷达几何参数

参数名	参数值
雷达平台高度(m)	3000
有效基线长度(m)	6.33
雷达下视角(°)	45
雷达工作频率(GHz)	1.3

下载: 导出CSV 

| 显示表格

图 2  归一化HSPC模型中误差

下载: 全尺寸图片 幻灯片

3.2   植被区电磁波最浅穿透深度估计

对于植被区，在植被平均高度已知而植被下地表高程无法获取情况下，如果能得到植被区的消光系数，则可以获取归一化HSPC高度，继而得到该消光系数下植被冠层顶部与HSPC高度之差，即植被中电磁波的最浅穿透深度：

对于由最大化相位差最优化方法分离的干涉相位，很容易得到分离的两种散射相位中心相对高度在不同相干系数幅度下的散点图，将两个散射相位中心散点相对高度差的平均值归一化后定义为相位中心高度差因子，如式(11)所示。

其中，${\varphi _h}$, ${\varphi _l}$分别为HSPC和LSPC干涉相位，mean(·)为求平均操作。

得到归一化HSPC和LSPC高度随消光系数变化的模型后，可以获得一系列对应消光系数下的归一化散射相位中心差。通过在归一化散射相位中心差中搜索得到的相位中心高度差因子，找出对应高度差因子下的消光系数${\sigma _{\rm ac}}$。至此，由搜索得到的消光系数可以求出植被区归一化HSPC高度${P_{\rm ac}}$，将${P_{\rm ac}}$代入式(10)中即可得到植被区电磁波的最浅穿透深度。

3.3   植被区DSM估计

根据上述最大化散射相位中心差最优化方法分析和归一化HSPC, LSPC变化模型的构建，在植被平均高度已知而植被下地表高程无法获取情况下植被区高精度DSM估计步骤如下：

步骤1　最大化相位差最优化法分离植被区干涉相干系数：对植被区像素进行极化干涉相干矩阵估计，然后进行最大化相位差分离得到各个像素点处的HSPC干涉相位${\varphi _h}$和LSPC干涉相位${\varphi _l}$；

步骤2　归一化散射相位中心高度差因子计算：利用干涉相位(去平地效应后)和垂直有效波束之间的关系得到不同干涉相干系数幅度下分离的HSPC和LSPC相对高度的散点图。根据式(11)将两个散射相位中心散点的高度之差平均值归一化后得到相位中心高度差因子$f$；

步骤3　模型建立及模型参数估计：利用式(6)，式(7)建立HSPC和LSPC归一化高度随着消光系数变化的模型。对于某一给定植被高度，通过求解式(9)所示的代价函数得到模型中的常数$a$, $b$，从而得到归一化HSPC高度${P_h}$和归一化LSPC高度${P_l}$随消光系数变化模型；

步骤4　植被区消光系数确定：利用模型中归一化HSPC和LSPC相对高度差${P_h}-{P_l}$搜索归一化散射相位中心高度差因子$f$所在位置，该位置对应消光系数值即电磁波在植被区域的消光系数${\sigma _{\rm ac}}$；

步骤5　电磁波最浅穿透深度${h_p}$估计：将步骤4中得到的消光系数${\sigma _{\rm ac}}$代入式(6)中得到植被区归一化HSPC高度${P_{\rm ac}}$，再将${P_{\rm ac}}$代入式(10)中可得到植被区冠层顶部与HSPC高度差，即为电磁波在植被中的最浅穿透深度${h_p}$；

步骤6　DSM估计：利用HSPC干涉相位${\varphi _h}$进行干涉处理得到粗DSM^[15,16]，然后将得到的电磁波在植被中最浅穿透深度${h_p}$补偿到粗DSM中得到最终的高精度DSM估计结果。

4.   仿真实验结果与分析

利用欧空局研发的PolSARpr软件提供的L波段森林植被极化干涉SAR仿真数据对植被区高精度DSM反演方法进行验证。仿真数据植被类型选择阔叶林和针叶林代表不同树种类型。树木高度选择18 m和10 m代表典型的高植被和矮植被。雷达几何参数与表1相同，场景参数如表2所示。

表 2  场景参数

参数名	参数值
树种类型	阔叶林/针叶林
植被高度(m)	18/10
树种密度(株/Ha)	600
距离向地形坡度(°)	0

下载: 导出CSV 

| 显示表格

以平均高度为18 m的阔叶林和平均高度为10 m的针叶林为例，下面给出植被中电磁波最浅穿透深度获取过程。图3(a)，图3(b)分别为18 m植被高度下利用最大化相位差最优化方法得到的HSPC和LSPC干涉相位结果。由干涉相位与植被垂直有效波束关系得到不同相干系数幅度下两种散射相位中心相对高度散点图如图3(c)所示，图中散射相位中心散点的相干系数幅度大多集中在0.8～1.0之间，相干系数较低时，两相位中心差距较大，而构建的模型中理论上两散射相位中心差距不能大于植被高度的一半。因此，可舍弃相干系数小于0.8且两相位中心差距大于植被一半高度的散点。之后对两相位中心高度差求平均值，得到归一化相位中心高度差因子为0.16。图3(e)为在该植被高度下构建的归一化散射相位中心高度模型，优化得到参数$a$为6.85, $b$为0.73。此时，在高度差因子为0.16时两散射相位中心差位于图3(e)中黑色线段部分，其对应消光系数为0.6 dB/m。由图3(f)可以得到消光系数为0.6 dB/m，植被高度为18 m时，植被最浅穿透深度为2.48 m。由于仿真数据设置的地形高程为零，可以由HSPC干涉相位与植被垂直有效波束关系直接得到HSPC的实际高度。图3(d)为场景中心方位线上HSPC干涉相位得到的HSPC实际高度，由此可知HSPC平均高程为16.10 m。补偿由模型得到的电磁波最浅穿透深度得到最终植被表面高程为18.58 m。因此，反演得到的DSM估计误差为18–18.58=–0.58 (m)(负号表示补偿后冠层表面高程高于实际高程)。与实际高程相比，补偿后冠层高程误差为0.58/18=3.22%，误差小于4%。同样，经过上述最浅穿透深度估计过程可以得到植被高度为10 m时高度差因子为0.44，如图4(e)中黑色线段部分，对应消光系数为0.05 dB/m，由图4(f)得到消光系数为0.05 dB/m，植被高度为10 m时，电磁波最浅穿透深度为4.73 m。由图4(d)知HSPC平均高程为5.06 m。因此，由模型得到的电磁波最浅穿透深度补偿后的DSM估计误差为10–9.79=0.21 (m)。与实际比，DSM误差为0.21/10=2.1%，误差在3%内。仿真结果表明，两种植被高度下，通过建立的两个散射相位中心相对高度变化模型能够有效估计电磁波在植被中的最浅穿透深度，从而得到高精度DSM。

图 3  18 m植被高度下电磁波最浅穿透深度估计过程

下载: 全尺寸图片 幻灯片

图 4  10 m植被高度下电磁波最浅穿透深度估计过程

下载: 全尺寸图片 幻灯片

为了验证本文方法的有效性，将文献[7]中利用数值半径法求解极化相干最优化相干系数从而生成DSM方法和文献[8]中利用联合干涉相位和相干性幅度的极化SAR最优相干性估计从而生成DSM方法以及本文方法得到的DSM进行对比。图5分别给出了针叶林和阔叶林中18 m和10 m两种植被高度情况下3种方法反演DSM结果。可以看出，文献[7]和文献[8]两种DSM反演方法在不同树种不同树高下得到的结果与真实值相比均偏低；相同树高下这两种方法在针叶林中反演的DSM结果与在阔叶林中反演的DSM结果相比较差，阔叶林反演结果更接近真实值，这是由于阔叶林与针叶林相比其结构比例更接近波长，电磁波在阔叶林中穿透性更差，在植被中的有效散射相位中心更加接近植被表面。对比图5(a)和图5(c)以及图5(b)和图5(d)，可以看出无论针叶林和阔叶林本文方法得到的DSM结果均更加接近真实DSM。

图 5  针叶林和阔叶林中两种植被高度下不同DSM反演结果对比

下载: 全尺寸图片 幻灯片

为了进一步验证上述方法的有效性，本文利用X波段下电磁波在稀疏森林区域的传播过程等效于L波段下电磁波在稠密森林区域的传播过程进行DSM反演验证，即电磁波较难穿透植被层到达地表时用机载实测数据进行DSM反演验证。试验数据为中国渭南东部黄河滩地区机载N-SAR系统录取的全极化干涉数据，该地区地势较为平坦，随机采样测得地表平均海拔高度为371 m，该地区植被为大片人工种植树林，较为稀疏，随机采样测得树林平均高度为19.29 m，即植被区DSM平均高度为390.29 m。图6(a)为该地区Google Earth光学图，图6(b)对应全极化SAR数据Pauli基极化分解图。数据大小为8192×4096像素。图7为图6(b)中红色线段部分数据不同方法下机载实测数据生成DSM结果的截面图。可以看出，不同方法反演得到的DSM结果较实际均偏低，文献[7]与文献[8]方法反演生成的DSM误差较大，而本文方法得到的DSM结果更加接近真实值。

图 6  机载全极化试验数据

下载: 全尺寸图片 幻灯片

图 7  不同方法机载全极化实测数据DSM反演结果对比

下载: 全尺寸图片 幻灯片

5.   结论

在植被平均高度已知而植被下地表高程无法获取情况下，为了解决植被区传统DSM反演方法精度低的问题，本文提出一种高精度DSM反演方法。本文以最大化相位差最优化方法估计相干系数为基础，通过归一化体散射相位中心随着消光系数变化趋势的分析，构建了一种电磁波在植被中高散射相位中心和低散射相位中心随消光系数变化的指数模型，由该模型得到植被中电磁波的消光系数从而获得电磁波在植被中的最浅穿透深度。将得到的电磁波最浅穿透深度补偿到高散射相位中心干涉相位处理得到的高程中可获得最终DSM结果。试验结果表明，本文方法较传统方法能有效提高植被区DSM反演精度。在植被高度变化不大的情况下，由此也可以反演出植被下高精度DTM结果^[17]，为植被区高精度地形测绘提供了一种可能。