CN103235301B - 基于复数域平差理论的POLInSAR植被高度反演方法 - Google Patents

Abstract

本发明涉及微波遥感在森林领域中的应用，公开了一种基于复数域平差理论的POLInSAR植被高度反演方法，该方法将复数域平差理论与POLInSAR植被高反演的物理模型进行有机结合，包括以下步骤：步骤1：对获取的N景全极化数据进行配准，获得M个基线干涉对；步骤2：对M个基线干涉对进行极化干涉，获得M*P个干涉条纹图，并估计其对应的复相干性；步骤3：基于复数域平差计算，求得植被高h_v。该基于复数域平差理论的POLInSAR植被高度反演方法原理直观，易于编程实现和扩展，是一种稳健可行的植被高反演算法，可以广泛应用于POLInSAR技术测量大范围区域乃至全球尺度植被高度。

Description

基于复数域平差理论的POLInSAR植被高度反演方法

技术领域

本发明涉及涉及微波遥感在森林领域中的应用，具体涉及一种基于复数域平差理论的POLInSAR植被高度反演方法。

背景技术

地形图是国家经济建设、社会发展和国家安全必不可少的基础图件，是基础设施规划、设计和施工及资源调查、开发所必须的基础资料。森林覆盖区域由于植被高度的影响，传统遥感手段测绘的地形图往往难以代表真实的地形图。因此，如何快速，准确地获取植被高度对于获取高精度的数字地形模型非常重要。另外，植被高度作为森林冠层结构的主要参数之一，影响着植物群落动态和构成以及边界层的气象和小气候，是研究、监测和管理森林资源的一个必不可少的信息来源，同时它还是森林蓄积量和森林地上生物量估计等生态系统模型的一个重要输入参数。

极化雷达干涉测量(Polarimetric SAR Interferometry，简称POLInSAR)技术集成了极化信息对植被散射体的形状、方位、介电特性敏感的优点和InSAR对植被散射体的空间分布和高度敏感的优势，具备全天侯、全天时、观测周期短、大范围乃至全球尺度监测植被高度的特点，在植被高测量领域体现出了巨大的潜力。

目前主要的植被高反演算法有：1)直接高度差法；2)基于相干散射模型的非线性迭代法和三阶段算法；3)类如旋转不变技术(Estimation of Signal Parametersvia Rotational Invariance Techniques，简称ESPRIT)等超分辨率方法。经过大量的实验验证，三阶段算法已经被证明是目前机载和室内极化干涉SAR数据(如EMSL数据)植被高反演精度最高的方法，是目前最常用的POLInSAR植被高度反演算法。它包括复平面直线拟合、植被偏差消除和植被高估计三个阶段，其反演精度主要受到地面相位估计误差和纯体相干性估计误差的两个关键问题的影响。

因此，有必要设计一种基于复数域平差理论的POLInSAR植被高度反演方法。

发明内容

本发明所要解决的技术问题是提供一种基于复数域平差理论的POLInSAR植被高度反演方法，该基于复数域平差理论的POLInSAR植被高度反演方法易于实施，应用广泛。

发明的技术解决方案如下：

一种基于复数域平差理论的POLInSAR植被高度反演方法，包括以下步骤：

步骤1：对获取的N景全极化数据进行配准，获得M个基线干涉对；

步骤2：对M个基线干涉对进行极化干涉，获得M*P(P为极化通道数)个干涉条纹图，并估计其对应的复相干性；

(植被高度反演模型的观测值是复相干性，所有前期处理(步骤1和2)是为了根据影像数据计算对应的相干性。然后利用步骤3中的反演模型和计算方法求得植被高度)

对于配准后的主影像S1和辅影像S2，复相干性计算公式如下：

$\mathrm{\γ}=\frac{|\⟨s_1{s}_2^{*}\⟩|}{\sqrt{\⟨s_1{s}_1^{*}\⟩\⟨s_2{s}_2^{*}\⟩}}$ 0≤γ≤1

式中s₁s₂*表示s₁与s₂的复共轭进行相乘，<·>表示求邻域平均。*是复数的共轭算子；

步骤3：基于复数域平差计算，求得植被高h_v；

所述的步骤1中：N为影像的数目，每景全极化数据即POLSAR数据，包含HH、HV、VH、VV四种极化数据；(由于N景全极化SAR数据是在不同条件下(天候、照度、摄像位置和角度等)获取的，故需要对其进行匹配将所有影像处理到同一坐标框架下，即配准。)所述的配准过程为：

选取其中一景POLSAR数据的HH极化SAR影像作为参考影像，将其余N-1景全极化数据中选取每景HH极化SAR影像采用GAMMA软件中灰度匹配方法配准到参考影像的坐标框架下；然后采用与HH极化组合方式相同的配准多项式对其它极化组合方式(HV、VH、VV)的SAR影像进行配准)；

任意两景全极化数据均进行干涉处理，从而得到个基线干涉对(每两景影像组合成一个基线)。

在所述的步骤2中，

对于每个基线干涉对，利用Pauli基合成技术或极化相干最优化技术生成Pauli极化或最优极化组合方式共P个，进而总的干涉对扩展为M*P个；(N景极化SAR影像数据可以生成M个基线干涉对，每个基线干涉对可以有P种极化方式组合，对应可以生成P个干涉对，故总干涉对个数为M*P个，生成M*P个干涉条纹图，每个干涉条纹图对可以获得一种复相干性观测量)

所述的步骤3包括如下几个子步骤：

步骤a：采用粗差探测的方法，剔除观测误差大于2倍中误差的极化组合方式；(从而筛选出数据质量较高的极化组合方式)；

步骤b：建立复数域平差函数模型(测量平差中平差由函数模型和随机模型组成)和随机模型；(在此基础上，根据最小二乘原理求解待估参数的估值)

所述的函数模型为随机地体二层相干散射模型(RandomVolume Over Ground，RVOG)，其表达式为：

式中w为单位复数矢量，对应于某一种极化方式，γ(w)表示极化方式w对应的复相干性；【γ（w)即 $\mathrm{\&gamma;}=\frac{|\&lang;s_1{s}_2^{*}\&rang;|}{\sqrt{\&lang;s_1{s}_1^{*}\&rang;\&lang;s_2{s}_2^{*}\&rang;}}$ 0≤γ≤1中的r，r是相干性计算的通用公式，r(w)代表利用w极化数据计算的相干性，即极化方式w对应的复相干性。】代表植被下的地表相位，h_v代表植被高度，μ(w)为地体幅度比，γ_V表示只由植被层产生的纯体相干性，有：

${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)}$ 【步骤c中策略1是将γ_V的表达式(包含植被高和消光系数两个未知数)

${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)}$ 代入到γ(w)中求解植被高，策略2是将γ_V看为一个整体未知数，先将其整体值求出，之后再利用

${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)}$ 求出植被高度和消光系数)

式中σ为消光系数，为待求参数，用于描述植被层粒子对入射波的吸收和散射过程；h_v表示植被高度；在函数模型中，γ(w)为模型中的观测量，μ(w)，γ_V(σ，h_v)为待求参数，k_z表示有效垂直波数，为已知数，根据成像几何(垂直基线B_⊥，斜距R和入射角θ)和雷达波长λ求得，即：

B_⊥代表两次成像期间传感器上雷达天线间的距离；(垂直基线代表基线沿垂直于雷达视线方向(雷达天线到目标点)的分量)

任意极化通道观测值的权为： $P\left(w_i\right)=\frac{\underset{j=1}{\overset{N}{\min }}\left({\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}^2\right)}{{\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_i\right)\right|}^2}$

N代表观测值个数，根据极化SAR理论，不同极化的复相干性的模的中误差表达式如下：

${\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_i\right)\right|}=\frac{1-{\left|\mathrm{\&gamma;}\left(w_i\right)\right|}^2}{\sqrt{2N_{\mathrm{sample}}}};$

N_sample代表相干性估计时参与估计的独立像元样本总个数；

(函数模型和随机模型是测量数据处理领域中的两个基本概念。函数模型是指的观测方程(即观测值与未知数的函数关系)，随机模型是用未定义不同精度的观测值在求解参数时的权重，观测值好的权重给的大一点，对最后参数解算的贡献大一点。P(w_i)为w_i极化复相干性观测值对应的权，后面公式已经说明如何将P(w_i)引入解算方程中，在后面的最小二乘解算时，我们对观测值定权，即确定P矩阵。而P矩阵的确定就要用到上面两个公式。)

随机模型即权阵P，权阵P的表达式为：

P＝diag(P(w₁)，P(w₁)，P(w₂)，P(w₂)，...，P(w_N)，P(w_N))

其中，diag表示对角矩阵算子，N表示代表复相干性观测值个数，此时P为一2N*2N列的对角矩阵，即为随机模型；

步骤c：基于复数域最小二乘准则，即复观测值残差的实部和虚部同时最小，采用两种策略求解未知参数：

1)直接利用非线性最小二乘迭代算法，在解空间内逐步搜索寻找最优解，求解植被高度；

2)分两步求解植被高度：首先将观测模型线性化，构建误差方程，结合误差方程直接解算处理技术，采用修正的奇异值分解方法对误差方程进行求解，进而利用最小二乘原理直接求取纯体相干系数，(纯体去相干系数为RVOG模型中，在不考虑地表散射的条件下，植被层的在两次成像中的相似程度，其是植被高度和消光系数的函数，参见公式 ${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)};$ 然后利用查找表的思想迭代求取植被高度。

(这两种方法都是基于最小二乘原理的，都可以求得植被高度，不过采取的策略略有不同，在实际应用时任选一种都可以)

在步骤1中，对M个干涉对相应的主、辅影像采用方位向和距离向频谱滤波，分别消除多普勒中心不一致引起的多普勒去相关以及两次成像入射角不同引起的基线几何去相关；在步骤2中，利用不同基线条件下的卫星轨道信息分别滤除相应干涉条纹图中的平地相位分量，降低条纹频率，从而有利于准确估计相干性。

步骤c中，复数域最小二乘准则对应的表达式为：

$Q=V^T\mathrm{PV}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left(w_i\right){\left(\right|\mathrm{\&gamma;}\left(w_i\right)-\widehat{\mathrm{\&gamma;}}\left(w_i\right)\left|\right)}^2$

$=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left(w_i\right)[{(\mathrm{Re}\left(\mathrm{\&gamma;}\left(w_i\right)\right)-\mathrm{Re}\left(\widehat{\mathrm{\&gamma;}}\left(w_i\right)\right))}^2+{(\mathrm{Im}\left(\mathrm{\&gamma;}\left(w_i\right)\right)-\mathrm{Im}\left(\widehat{\mathrm{\&gamma;}}\left(w_i\right)\right))}^2]=\min$

其中，表示平差后的复相干性观测值，V表示观测值残差，P表示观测值的权阵，P为2N*2N的对角阵，Re为取复数实部算子，Im为取复数虚部算子，N代表复相干性观测值个数(如说明书中我们用了5种复相干系数，这种条件下N＝5*2＝10.)

步骤c的策略2)中，采用下式求解：

X的表达式为：

$X=V\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& 0\\ 0& 0\end{array}\right]U^H{D}^{-1}L;$

X表示纯体去相干系数γ_V，地表相位有效地体散射幅度比μ(w)(这里面待估参数是测量平差中使用的术语，也就是方程中的未知数)的估计值；

式中，

∑^-1＝diag(1/σ₁，…，1/σ_s，...1)为∑的逆矩阵，其中1/σ_i(i＝1，2，...，S)为B₁矩阵奇异值的倒数，式中阂值S的确定原则为：当∑矩阵对角线上的前S个奇异值的和占总奇异值总和的90％以上时，则将∑的倒数矩阵∑^-1中前S个元素保留，剩下的N-S元素赋为1；

V和U为对B₁按下式进行SVD分解得到的：

$B_1=U\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]V^H;$ 其中B₁＝D^-1B，B为平差函数模型中的设计矩阵，其是通过对上述RVoG模型进行线性化得到的未知数改正数的系数矩阵；

D为对权阵P的逆阵进行Cholesky分解得到的三角阵；L为平差模型中的常数项矩阵；

求解得到去相干系数γ_V后，再利用查表方法反推求得植被高度h_v。

具体求解过程为：

根据先验信息给植被高度和消光系数一定的取值范围(如L波段数据消光系数一股不超过1dB/m，植被高度范围可根据NASA发布的全球植被高度图设定其范围)，将植被高度取值范围和消光系数取值范围各分成X等分(如1000等分)，并按照公式 ${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)}$ 计算出γ_V(σ，h_v)值，这样就建立了一个X*X大小的二维表，然后对比γ_V的估计值和查找表中的理论值，找到差异最小的那一组值，就得到了待求的植被高度h_V】

采用如下方法对方程进行求解：

其中，X表示纯体去相干系数γ_V，地表相位有效地体散射幅度比μ(w)(这里面待估参数是测量平差中使用的术语，也就是方程中的未知数)的估计值；

将γ_V写成γ_V＝a+bi，其中i代表单位虚数，a代表实部，b代表虚部；

(策略2是将γ_V看为一个整体未知数，先将其整体值求出，之后再利用γ_V求出植被高度和消光系数。所以这里首先将复数γ_V写成γ_V＝a+bi)

表示B₁的广义逆矩阵，L₁、B₁的具体形式为：

B₁＝D^-1B，L₁＝D^-1L；

其中D为对权阵P的逆阵进行Cholesky分解得到的三角阵；B为平差函数模型中的设计矩阵，其是通过对上述RVoG模型进行线性化得到的未知数改正数的系数矩阵(线性化操作为测量平差理论中数据处理的一种方法)；L为平差模型中的常数项矩阵，之后，对B₁按下式进行SVD分解得到U、V、 $\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]$ 三部分， $B_1=U\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]V^H;$ (这是两步运算，首先对设计矩阵B同D相乘得到B₁，之后对得到的乘积B₁再进行SVD分解，一个矩阵可以分解为多个矩阵的乘积，数学上有很多种分解方法，这里的SVD分解是矩阵分解的一种，是减弱解算过程病态问题的一种常用手段。B和D矩阵是已知的，由这两个可以求得B1，然后再将B1进行SVD分解成U，V， $\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]$ 三部分】；式中，U，V表示SVD矩阵分解结果中两个相互正交矩阵，H表示共轭转置算子；而 $\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]$ 代表SVD矩阵分解结果中一对角矩阵，∑表示非零对角阵，其表达式为：

∑＝diag(σ_１，σ_２，…，σ_n)，σ_１≥σ_２...≥σ_n

式中σ_i表示B₁矩阵的第i个奇异值，n为待估参数的个数；

则最终的参数估计X表达式为：

$X=V\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& 0\\ 0& 0\end{array}\right]U^H{D}^{-1}L;$

式中∑^-1表示对非零对角阵∑的逆阵；

采用SVD分解的修正方案，即将∑^-1修正为：

${\widetilde{\mathrm{\&Sigma;}}}^{-1}=\mathrm{diag}(1/{\mathrm{\&sigma;}}_1,\&CenterDot;\&CenterDot;\&CenterDot;,1/{\mathrm{\&sigma;}}_s,...1)$

其中1/σ_i(i＝1，2，...，S)为B₁矩阵奇异值的倒数，式中阂值S的确定原则为：当∑矩阵对角线上的前S个奇异值的和占总奇异值总和的90％以上时，则将∑的倒数矩阵∑^-1中前S个元素保留，剩下的N-S元素赋为1；】

计算植被高度是在权利3中的步骤c)，本专利在利用复数平差理论计算植被高时，采用了两种策略。

1)直接利用非线性最小二乘迭代算法求所有参数，在待求参数解向量空间内逐步搜索寻找使得权利5中满足复数域最小二乘准则(即Q最小)对应的解向量，解向量中的一个解对应要求的植被高度；

2)分两步求解植被高度，先γ_V看为一个整体未知数，此时求解γ_V是一个线性平差问题，利用复数域平差求得满足权利5中满足复数域最小二乘准则(即Q最小)时γ_V的估计值，事先利用γ_V(σ，h_v)理论计算公式建立二维查找表，然后对比γ_V的估计值和理论值，找到差异最小的哪一组值，其中一个即为所求植被高度

有益效果：

本发明的基于复数域平差理论的POLInSAR植被高度反演方法，该方法将复数域平差理论与POLInSAR植被高反演的物理模型进行有机结合：首先利用SAR灰度匹配技术将全极化数据影像配准；然后利用Pauli基合成技术和极化相干最优化技术扩展目标极化干涉观测量；接着利用不同极化方式对应的复相干性(包含相位信息及相干性信息)作为观测值，并对其采用粗差探测的方法剔除观测误差较大的观测量；之后利用POLInSAR植被高反演物理模型建立平差数学模型，采用复相干性及干涉相位的先验统计信息建立平差随机模型；最后基于复数域最小二乘理论，采用两种策略反演植被高度。该方法突破了传统三阶段植被高反演算法根据经验选取纯体相干性导致结果不准确的限制，且原理更为直观，易于编程实现和扩展，是一种稳健可行的植被高反演算法，可以广泛应用于POLInSAR技术测量大范围区域乃至全球尺度植被高度。

附图说明

图1为多基线POLInSAR观测几何模型。

图2为复数域平差法反演植被高流程；

图3为RVOG模型在复数平面内的表达。其中(a)图为不同极化方式对应复相干性在复平面上理论分布，(b)图为不同极化方式对应复相干性在复平面上实际分布。

图4为主、副影像利用Pauli基技术生成的三个Pauli极化RGB合成图。

图5为步骤4得到的不同极化方式的干涉条纹图。

图6为步骤5得到的去除平地相位后的不同极化方式的干涉条纹图。

图7为步骤6得到的不同极化方式的复相干系数的幅度图。

图8为基于复数域平差算法和实例数据反演的植被高度图。其中(a)图为步骤9中策略1反演的植被高度图，(b)图为步骤9中策略2反演的植被高度图。

图9为图2所述流程中极化干涉部分的具体实现流程。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，利用欧空局发布的POLSARpro软件模拟的2景全极化数据(模拟参数设置见附表1)，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明具体实施例以单基线线观测配置来阐述，即成像区域只有2景全极化SAR影像的情况。对应于每景全极化SAR影像，用极化散射矩阵(S矩阵)来表达这四个通道影像数据如公式(1)所示。该式中S_pq表示q极化发射、p极化接收时目标的后向散射，p、q为水平h或垂直v极化。一般而言，对互易媒质，根据单站互易原理，S矩阵将满足S_hv＝S_vh的条件，此时散射矩阵为对称矩阵：

$S=\left[\begin{array}{cc}{S}_{\mathrm{hh}}& S_{\mathrm{hv}}\\ S_{\mathrm{vh}}& S_{\mathrm{vv}}\end{array}\right]---\left(1\right)$

以S_I表示主影像，S₂表示辅影像，其对应的散射矩阵表示为：

$S_1=\left[\begin{array}{cc}{S}_{\mathrm{hh}}^1& S_{\mathrm{hv}}^1\\ S_{\mathrm{vh}}^1& S_{\mathrm{vv}}^1\end{array}\right]$ $S_2=\left[\begin{array}{cc}{S}_{\mathrm{hh}}^2& S_{\mathrm{hv}}^2\\ S_{\mathrm{vh}}^2& S_{\mathrm{vv}}^2\end{array}\right]---\left(2\right)$

本发明的具体实施步骤为：

步骤1：利用SAR影像匹配技术，将2景极化SAR影像配准。由于主影像4个极化通道的数据具有相同的成像几何结构，因此本发明首先利用GAMMA软件中基于相干性测度的配准方法对选取的两景极化影像进行配准，具体方法如下：利用主影像中的HH极化影像和辅影像中的HH极化影像以相干性为测度选取高相干性点作为同名点，然后利用多项式建立几何变换模型，最后将辅影像灰度重采样到主影像坐标框架下。接着，其他三个极化方式按照HH极化影像的几何变换模型，将其余3种极化数据影像配准到主影像的地理坐标框架下。相干性是指两幅影像的相似程度，对于如图1所示的InSAR成像几何中需要配准的主影像S1和辅影像S2，其相干性计算公式如下：

${\mathrm{\&gamma;}}_{\mathrm{Int}}=\frac{|\&lang;s_1{s}_2^{*}\&rang;|}{\sqrt{\&lang;s_1{s}_1^{*}\&rang;\&lang;s_2{s}_2^{*}\&rang;}}$ 0≤γ≤1      (3)

式中s₁s₂*表示s₁与s₂的复共轭进行相乘，<·>表示求邻域平均。【*是复数的共轭算子】

步骤2：对相应的主、辅影像采用GAMMA软件中方位向和距离向频谱滤波，分别消除多普勒中心不一致引起的多普勒去相关以及两次成像入射角不同引起的基线几何去相关。

步骤3：基于不同基线下已有的3个线极化方式(HH、HV、VV，单站互易有HV＝VH)，利用Pauli基合成技术或极化相干最优化技术生成Pauli极化或最优极化方式的主辅影像对。对两景极化SAR数据S₁、S₂进行Pauli基转换得到两个Pauli基矢量k₁和k₂：

$k_1=V\left(\right[S_1\left]\right)=\frac{1}{\sqrt{2}}{[S_{\mathrm{hh}}^1+S_{\mathrm{vv}}^1,S_{\mathrm{hh}}^1-S_{\mathrm{vv}}^1,2S_{\mathrm{hv}}^1]}^T---\left(4\right)$

$k_2=V\left(\right[S_2\left]\right)=\frac{1}{\sqrt{2}}{[S_{\mathrm{hh}}^2+S_{\mathrm{vv}}^2,S_{\mathrm{hh}}^2-S_{\mathrm{vv}}^2,2S_{\mathrm{hv}}^2]}^T---\left(4\right)$

根据极化合成理论，用3个线极化观测量(单站互易有HV＝VH))可以进行极化合成，从而获取极化空间内任意一种极化观测量。具体表现为利用单位投影向量w_i，将Pauli基矢量元素线性组合得到任意极化观测量S_i，表达式如下：

S_i＝w_i ^*Tk      (6)

式中i表示任一极化方式(如HH，HV，VV等)，*T代表共轭转置算子，k代表Pauli基矢量

Pauli基合成极化方式：

对于Pauli基极化中的HH+VV极化方式有HH+VV和HH-VV分别代表一种极化方式；

对于Pauli基极化中的HH-VV极化

附图4为利用Pauli基合成技术分别对主、副影像进行极化合成，生成的RGB(R表示HH-VV极化方式的回波强度，G表示HV极化方式的回波强度，B表示HH+VV极化方式的回波强度)彩色合成图。图像中央较亮的“圆形”区域为植被覆盖区，其余区域为平地。

相干最优化极化方式：

所有极化相干最优算法都是寻求某些特殊的单位投影向量w_i，使得到的极化观测量能够满足一定的约束条件。本发明中极化相干最优化算法涵盖目前所有的极化相干最优化算法，如Cloude提出的相干性最大的相干最优算法，简称C&P算法，Tabb提出的相位最大分离的PD相干最优算法等等，并兼容未来提出其他极化相干最优化算法。本示例说明过程中，不采用相干最优的极化方式，故本发明在此略去相干最优的相关结果。

步骤4：将相同极化的主辅影像分别共轭相乘做干涉处理，生成不同极化方式的干涉条纹图(见附图5)。

步骤5：利用卫星轨道信息分别滤除相应干涉条纹图中的平地相位分量，降低条纹频率(见附图6)，从而有利于准确估计相干性。

步骤6：利用滤掉平地相位分量的干涉条纹图和主辅影像灰度图，利用PolSARpro软件中的boxcar算法估计不同极化方式的复相干性。图7为复相干系数的幅度。

步骤7：粗差探测。对于每一个像元，假设提取其5种极化组合方式(HH，HV，VV，HH-VV，HH+VV)，按照上面的步骤，获得复相干性：

γ_j＝Re(γ_j)+i*Im(γ_j)      (7)

式中j＝1、2、3、4、5，分别代表极化方式HH，HV，VV，HH-VV，HH+VV；Re(γ_j)表示复相干系数γ_j的实部，Im(γ_j)表示复相干系数γ_j的虚部；i表示复数虚部；*表示相乘。故可以组成坐标(Re(γ_j)，i*Im(γ_j))，将其表达为(x_j，γ_j)，选取模拟数据中的一个像元进行说明，将其在复数平面内进行表达(见附图3(b))。“相干范围”所在直线可以表达为y＝ax+b，按照整体最小二乘的方法，首先我们建立如下误差方程：

$y_j+v_{y_j}=a(x+v_{x_j})+b---\left(8\right)$

将上式整理得：

(A+E_A)x＝L+E_L      (9)

其中：

$A=\left[\begin{array}{cc}{x}_1& 1\\ x_2& 1\\ x_3& 1\\ x_4& 1\\ x_5& 1\end{array}\right]$ $E_A=\left[\begin{array}{cc}{v}_{x_1}& 0\\ v_{x_2}& 0\\ v_{x_3}& 0\\ {v_x}_4& 0\\ {v_x}_5& 0\end{array}\right]$ $E_A=\left[\begin{array}{cc}{v}_{x_1}& 0\\ v_{x_2}& 0\\ v_{x_3}& 0\\ {v_x}_4& 0\\ {v_x}_5& 0\end{array}\right]$ $L=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\\ y_5\end{array}\right]$ $E_L=\left[\begin{array}{c}{v}_{y_1}\\ v_{y_2}\\ v_{y_3}\\ v_{y_4}\\ v_{y_5}\end{array}\right]$

根据整体最小二乘准则利用SVD分解对待估参数a、b进行估计。接下来，计算上述5个复相干系数对应的坐标到拟合直线( 表示待估参数a、b的估计值)之间的距离d_j。

$d_j=\frac{|{\hat{a}x}_j-y_j+\hat{b}|}{\sqrt{1+a^2}}---\left(10\right)$

利用上式计算d_j的标准偏差：

$\mathrm{\&sigma;}=\sqrt{\frac{{(d_j-\stackrel{\&OverBar;}{d})}^T(d_j-\stackrel{\&OverBar;}{d})}{n-1}}$ $\stackrel{\&OverBar;}{d}=\frac{1}{n}\underset{j=1}{\overset{5}{\mathrm{\&Sigma;}}}{d}_j---\left(11\right)$

式中，n为选取的极化方式个数，对于本示例，n＝5。当d_j>2σ时，认为其对应的复相干系数误差较大，被认为含有粗差，将其其对应的复相干性剔除，不参与接下来的平差计算。

步骤8：对于每一个像元而言，根据平差模型的构成条件，需利用相应的复相干性观测值和随机地体二层相干散射模型(Random Volume Over Ground，简称RVOG)，建立复数域平差数学模型和随机模型，在此基础上，根据最小二乘原理求解待估参数的估值。

POLInSAR植被高反演的函数模型是基于目前POLInSAR在森林植被高度估计应用中最常用且有效的散射模型，即随机地体二层相干散射模型(RandomVolume Over Ground，RVOG)。该模型从复相干性出发，假定植被层中散射能量随着高度增加呈指数变化，在进行距离向频谱滤波消除基线几何去相干并且不考虑时间去相干的复相干性可以表示为：

0≤L(w)≤1

式中w为单位复数矢量，代表某一种物理散射机制，代表植被下的地表相位，h_v代表植被高度，μ(w)为地体幅度比，γ_V表示只由植被层产生的纯体相干性，其一般表达式为：

${\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+i{k}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+i{k}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)}---\left(13\right)$

式中σ为消光系数，为待求参数，用于描述植被层粒子对入射波的吸收和散射过程；h_v表示植被高度，为待求参数，k_z表示有效垂直波数，为已知数，根据成像几何(垂直基线B_⊥，斜距R和入射角θ)和雷达波长λ求得，其一般表达式为：

$k_z=\frac{4\mathrm{\&pi;\&Delta;\&theta;}}{\mathrm{\&lambda;}\sin \mathrm{\&theta;}}\&ap;\frac{4\mathrm{\&pi;}{B}_{\&perp;}}{\mathrm{\&lambda;}R\sin \mathrm{\&theta;}}---\left(14\right)$

式(12)即为RVOG模型的最终表达式。观察RVOG模型的表达式可以看出，植被高反演是基于复数观测值，即复相干性γ(w)。当具备一定观测数量时，可以采用复数域最小二乘求解。

POLInSAR植被高反演的随机模型，其特征在于，根据极化SAR理论，不同极化的复相干性的模的中误差表达式如下：

${\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}=\frac{1-{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}^2}{\sqrt{2L}}---\left(15\right)$

其中，L代表相干性估计时参与估计的独立像元样本总个数，对于本发明中的实例选取7×7大小的窗口进行估计。本发明的植被高反演平差策略中，假定不同极化观测值之间互相独立，选取相干值中误差最小对应的极化方式的中误差为单位权中误差，则根据权的定义可以得到任意极化通道观测值的权为：

$P\left(w_j\right)=\frac{\underset{j=1}{\overset{N}{\min }}\left({\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}^2\right)}{{\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}^2}---\left(16\right)$

步骤9：由于该模型自身的特点，基于复数域最小二乘准则，即复观测值残差的实部和虚部同时最小，我们可以采用两种策略求解未知参数：

(1)策略1：直接利用非线性迭代算法求解植被高度

假定选取的极化通道数为N，则POLInSAR植被高反演的复数最小二乘准则表达式为：

其中，[M]表示观测模型(12)，然后采用矩阵运算软件Matlab软件自带的最小二乘迭代算法寻找最优解，从而求得植被高度h_v(见附图8(a))。

(2)策略2：分两步求解植被高度。首先将观测模型线性化，构建误差方程，采用修正的奇异值分解法直接对误差方程进行解算，进而直接得到纯体相干系数γ_v，然后利用查找表的思想迭代求取植被高度h_v。

具体地，首先令γ_v＝a+bi，同时将变换为：带入(12)式，将实部、虚部拆分并分别进行泰勒级数展开，得到平差模型中的数学模型为：

(18)

其中：

式中下标1代表实部；2代表虚部；i代表极化方式；a₀、b₀、μ_i0、分别为a、b、u_i、的近似值；da、db、du_i、为对应待估参数的改正值。本实施例中采用γ_HH、γ_HV、γ_VV、γ_HH-VV、γ_HH+VV5种线性极化方式，可以得到10个观测量(每个极化方式组合对应的相干系数有实部和虚部)，8个待估参数：γ_v(含有2个未知参数：a、b、μ₁、μ₂、μ₃、μ₄、μ₅，多余观测数为2。对每种极化方式分别按照(18)、(19)式列误差方程，用矩阵形式表示为如下：

$\underset{10\&times;1}{V}=\underset{10\&times;8}{B}\underset{8\&times;1}{X}-\underset{10\&times;1}{L}---\left(20\right)$

其中， $\underset{10\&times;1}{V}={\left[\begin{array}{ccccc}{v}_{1\mathrm{HH}}& v_{2\mathrm{HH}}& \&CenterDot;\&CenterDot;\&CenterDot;& v_{1\mathrm{HH}+\mathrm{VV}}& v_{2\mathrm{HH}+\mathrm{VV}}\end{array}\right]}^T---\left(21\right)$

$\underset{10\&times;1}{L}={\left[\begin{array}{ccccc}{l}_{1\mathrm{HH}}& l_{2\mathrm{HH}}& \&CenterDot;\&CenterDot;\&CenterDot;& l_{1\mathrm{HH}+\mathrm{VV}}& l_{2\mathrm{HH}+\mathrm{VV}}\end{array}\right]}^T---\left(24\right)$

由于对平差模型的数学模型过度参数化，导致误差方程设计矩阵病态，无法按照间接平差解算方法，解算得到的最小二乘解。故本发明采用一种修正的奇异值分解方法直接对误差方程(20)进行求解。具体地，对于观测量不等精度的情况下，采用如下方法对法方程进行求解得参数的估计值为：

$X=B_1^{+}{L}_1---\left(25\right)$

其中，X表示待估参数的估计值；表示B₁的广义逆矩阵，L₁、B₁的具体形式为：

B₁＝D^-1B L₁＝D^-1L      (26)

其中D为对权阵P的逆阵进行Cholesky分解得到的三角阵；B为平差模型中的设计矩阵；L为平差模型中的常数项矩阵。之后，对B₁进行SVD分解可以得到；

$B_1=U\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]V^H---\left(27\right)$

式中，U，V表示SVD矩阵分解结果中两个相互正交矩阵；而 $\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]$ 代表SVD矩阵分解结果中一对角矩阵，∑表示非零对角阵；H表示共轭转置算子。

求解(22)式M-P逆，并带入(20)式有：

$X=V\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& 0\\ 0& 0\end{array}\right]U^H{D}^{-1}L---\left(28\right)$

式中，∑^-1表示对非零对角阵∑的逆阵。当系数矩阵严重病态时，采用一种对SVD分解的修正方案。即将(28)式中∑^-1修正为：

${\widetilde{\mathrm{\&Sigma;}}}^{-1}=\mathrm{diag}(1/{\mathrm{\&sigma;}}_1,\&CenterDot;\&CenterDot;\&CenterDot;,1/{\mathrm{\&sigma;}}_s,...1)---\left(29\right)$

其中1/σ_i为B₁矩阵的奇异值的倒数。利用修正的奇异矩阵分解法对无差方程进行求解，根据(28)式求得a、b，从而得到利用体散射相干系数γ_v＝a+bi。

在求得纯体相干性的估计值后，接着利用查找表思想，在消光系数和植被高一定的解空间内，搜索到满足复数模差异最小的哪一组解作为参数的树高和消光系数的最佳估计值，从而求得植被高hv(见附图8(b))。实例中的约束条件表达式为：

$\underset{\mathrm{\&sigma;},h_v}{\min }\left|\right|{\widehat{\mathrm{\&gamma;}}}_V-{\mathrm{\&gamma;}}_V(\mathrm{\&sigma;},h_V)\left|\right|---\left(30\right)$

表1模拟数据参数设置

Claims (5)

1.一种基于复数域平差理论的POLInSAR植被高度反演方法，其特征在于，包括以下步骤：

步骤1：对获取的N景全极化数据进行配准，获得M个基线干涉对；

步骤2：对M个基线干涉对进行极化干涉，获得M*P个干涉条纹图，并估计其对应的复相干性；P为极化通道数；

对于配准后的主影像S1和辅影像S2,复相干性计算公式如下：

$\begin{array}{cc}\mathrm{\&gamma;}=\frac{|\&lang;s_1{s}_2^{*}\&rang;|}{\sqrt{\&lang;s_1{s}_1^{*}\&rang;\&lang;s_2{s}_2^{*}\&rang;}}& 0\&le;\mathrm{\&gamma;}\&le;1\end{array}$

式中s₁s₂*表示s₁与s₂的复共轭进行相乘，〈·〉表示求邻域平均，*是复数的共轭算子；

步骤3：基于复数域平差计算，求得植被高h_v；

所述的步骤1中：N为影像的数目，每景全极化数据即POLSAR数据，包含HH、HV、VH、VV四种极化数据；所述的配准过程为：

选取其中一景POLSAR数据的HH极化SAR影像作为参考影像，将其余N-1景全极化数据中选取每景HH极化SAR影像采用GAMMA软件中灰度匹配方法配准到参考影像的坐标框架下；然后采用与HH极化组合方式相同的配准多项式对HV、VH和VV极化组合方式的SAR影像进行配准；

任意两景全极化数据均进行干涉处理，从而得到个基线干涉对；

在所述的步骤2中，

对于每个基线干涉对，利用Pauli基合成技术或极化相干最优化技术生成Pauli极化或最优极化组合方式共P个，进而总的干涉对扩展为M*P个。

2.根据权利要求1所述的基于复数域平差理论的POLInSAR植被高度反演方法，其特征在于，所述的步骤3包括如下几个子步骤：

步骤a：采用粗差探测的方法，剔除观测误差大于2倍中误差的极化组合方式；

步骤b：建立复数域平差函数模型和随机模型；

所述的函数模型为随机地体二层相干散射模型(Random Volume Over Ground，RVOG)，其表达式为：

式中w为单位复数矢量，对应于某一种极化方式，γ(w)表示极化方式w对应的复相干性；代表植被下的地表相位，h_v代表植被高度，μ(w)为地体幅度比，γ_V表示只由植被层产生的纯体相干性，有:

${\mathrm{\&gamma;}}_V=(\mathrm{\&sigma;},h_v)=\frac{2\mathrm{\&sigma;}(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;}+{\mathrm{ik}}_z{h}_v)-1)}{(2\mathrm{\&sigma;}+{\mathrm{ik}}_z\cos \mathrm{\&theta;})(\exp (2\mathrm{\&sigma;}{h}_v/\cos \mathrm{\&theta;})-1)};$

式中σ为消光系数，为待求参数，用于描述植被层粒子对入射波的吸收和散射过程；h_v表示植被高度；在函数模型中，γ(w)为模型中的观测量，μ(w)，γ_V(σ,h_v)为待求参数，k_z表示有效垂直波数，为已知数，根据成像几何和雷达波长λ求得，即：

B_⊥代表两次成像期间传感器上雷达天线间的距离；(垂直基线代表基线沿垂直于雷达视线方向(雷达天线到目标点)的分量)

任意极化通道观测值的权为： $P\left(w_i\right)=\frac{\underset{j=1}{\overset{N}{\min }}\left({\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_j\right)\right|}^2\right)}{{\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_i\right)\right|}^2}$

N代表观测值个数，根据极化SAR理论，不同极化的复相干性的模的中误差表达式如下：

${\mathrm{\&sigma;}}_{\left|\mathrm{\&gamma;}\left(w_i\right)\right|}=\frac{{1-\left|\mathrm{\&gamma;}\left(w_i\right)\right|}^2}{\sqrt{{2N}_{\mathrm{sample}}}};$

N_sample代表相干性估计时参与估计的独立像元样本总个数；

随机模型即权阵P，权阵P的表达式为：

P＝diag(P(w₁),P(w₁),P(w₂),P(w₂),...,P(w_N),P(w_N))

其中，diag表示对角矩阵算子，N表示代表复相干性观测值个数，此时P为一2N*2N列的对角矩阵，即为随机模型；

步骤c：基于复数域最小二乘准则，即复观测值残差的实部和虚部同时最小，采用两种策略求解未知参数：

1)直接利用非线性最小二乘迭代算法，在解空间内逐步搜索寻找最优解，求解植被高度；

2)分两步求解植被高度：首先将观测模型线性化，构建误差方程，结合误差方程直接解算处理技术，采用修正的奇异值分解方法对误差方程进行求解，进而利用最小二乘原理直接求取纯体相干系数，然后利用查找表的思想迭代求取植被高度。

3.根据权利要求1所述的基于复数域平差理论的POLInSAR植被高度反演方法，其特征在于，在步骤1中，对M个干涉对相应的主、辅影像采用方位向和距离向频谱滤波，分别消除多普勒中心不一致引起的多普勒去相关以及两次成像入射角不同引起的基线几何去相关；在步骤2中，利用不同基线条件下的卫星轨道信息分别滤除相应干涉条纹图中的平地相位分量，降低条纹频率，从而有利于准确估计相干性。

4.根据权利要求2所述的基于复数域平差理论的POLInSAR植被高度反演方法，其特征在于，步骤c中，复数域最小二乘准则对应的表达式为：

$\begin{array}{c}Q=V^T\mathrm{PV}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left(w_i\right){\left(\right|\mathrm{\&gamma;}\left(w_i\right)-\widehat{\mathrm{\&gamma;}}\left(w_i\right)\left|\right)}^2\\ =\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}P\left(w_i\right)[{(\mathrm{Re}\left(\mathrm{\&gamma;}\left(w_i\right)\right)-\mathrm{Re}\left(\widehat{\mathrm{\&gamma;}}\left(w_i\right)\right))}^2+{(\mathrm{Im}\left(\mathrm{\&gamma;}\left(w_i\right)\right)-\mathrm{Im}\left(\widehat{\mathrm{\&gamma;}}\left(w_i\right)\right))}^2]=\min \end{array}$

其中，表示平差后的复相干性观测值，V表示观测值残差，P表示观测值的权阵，P为2N*2N的对角阵，Re为取复数实部算子，Im为取复数虚部算子，N代表复相干性观测值个数。

5.根据权利要求2所述的基于复数域平差理论的POLInSAR植被高度反演方法，其特征在于，步骤c的策略2)中，采用下式求解：

$X=V\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& 0\\ 0& 0\end{array}\right]U^H{D}^{-1}L;$

X表示纯体去相干系数γ_V，地表相位有效地体散射幅度比μ(w)(这里面待估参数是测量平差中使用的术语，也就是方程中的未知数)的估计值；

式中，

                         ∑^-1＝diag(1/σ₁,…,1/σ_s,...1)为∑的逆矩阵，其中1/σ_i(i＝1,2,…,S)为B₁矩阵奇异值的倒数，式中阈值S的确定原则为：当∑矩阵对角线上的前S个奇异值的和占总奇异值总和的90％以上时，则将∑的倒数矩阵∑^-1中前S个元素保留，剩下的N-S元素赋为1；

V和U为对B₁按下式进行SVD分解得到的：

$B_1=U\left[\begin{array}{cc}\mathrm{\&Sigma;}& 0\\ 0& 0\end{array}\right]V^H;$ 其中B₁＝D^-1B，B为平差函数模型中的设计矩阵，其是通过对上述RVoG模型进行线性化得到的未知数改正数的系数矩阵；

D为对权阵P的逆阵进行Cholesky分解得到的三角阵；L为平差模型中的常数项矩阵；

求解得到去相干系数γ_V后，再利用查表方法反推求得植被高度h_v。