CN101369019A - 基于极化数据融合的极化干涉合成孔径雷达三维成像方法 - Google Patents

Abstract

基于极化数据融合的极化干涉合成孔径雷达三维成像方法，属于极化雷达干涉测量数据处理技术领域。读入数据，对于每一个散射单元构造两个散射矢量k₁和k₂；对矩阵k₁ ^Hk₁－k₂ ^Hk₂进行特征值分解得到三个特征值λ₁＞0、λ₁＞0和λ₃＝0；分别计算见式(Ⅰ)的值以确定属于哪一种情况；属于第一种情况则按照式6)计算w＝w⁽¹⁾；第二种情况则按照式14)计算w＝w⁽²⁾；计算φ_s＝arg(w^Hk₁k₂ ^Hw)。根据式3)得到增强后的复图像；计算矩阵Ω₁₂＝<k₁k₂ ^H>，然后按照式22)计算φ_m＝arg(w^H [Ω₁₂]w)。对于所得到的整幅相位图(单视或多视)去平地效应；利用枝切法进行相位解缠。本发明适用范围广，计算时间比相干最优化方法的时间少；且融合后的相位质量要优于相干最优化方法，无论原信号强弱，均能有效地增强相位质量；为后续的精确的数字高程模型的产生奠定良好的基础。

Description

基于极化数据融合的极化干涉合成孔径雷达三维成像方法

技术领域

本发明涉及雷达***的极化干涉合成孔径雷达相位增强和三维成像方法。

背景技术

合成孔径雷达(Synthetic Aperture Radar，以下用SAR表示)是一种全天候、全天时、高分辨率雷达，通过距离向上对大时间带宽积的线性调频信号进行脉冲压缩和方位向的回波信号相干积累获得二维高分辨率的图像。

当飞行高度为H的雷达平台沿着直线飞行时，朝飞行方向的正侧方发射波束，如附图1所示，波束主瓣在地面覆盖一定面积，随着平台的运动，就会形成一条测绘带，定义雷达平台飞行的方向为方位向，与之垂直的方向为距离向，将雷达天线到测绘带内任意点目标C的距离称为斜距，其中，在雷达照射期间，天线到C点的最短距离称为该目标的最短斜距。雷达对点目标C的观测，如附图2所示，在某个时刻T₀，雷达波束在方位向上覆盖B～C的范围，此时点目标C刚刚进入波束；经过一段时间T_s后，波束在方位向上覆盖C～D段，C点刚刚脱离波束照射，雷达平台在T_s时间内走过的距离为L_s称为一个合成孔径长度，将T_s称为一个合成孔径时间，在此时间内，点目标C一直处于雷达波束照射下，点目标与雷达位置的几何关系，如附图3所示，其中，定义β为仰角，θ为斜视角。

雷达向地面发射脉冲，在一定时间延迟后，接收到地面场景内不同散射点反射的回波，经过距离向离散采样，积累多次回波得到数据阵列，如附图4所示，平行于距离向的数据代表一条距离线，平行于方位向的数据代表一条方位线，其中，距离向采样率对应的采样间隔构成的平行方位向的单元称为距离单元。在雷达平台直线运动过程中，同一目标到天线的斜距不断变化，如附图5所示，导致了回波的瞬时频率成线性调频特性，使SAR能够获得方位向的高分辨率。由于方位向和距离向的信号均为线性调频信号，对SAR信号成像的过程就是实现方位向和距离向二维脉冲压缩的过程。

SAR成像***，与其它成像***如光学成像***相比较，具有以下优点：

(1)SAR采用主动式微波成像，具有全天候、全天时成像的特点；

(2)选择合适波长，利用微波的穿透性，可以对被植被、沙漠、极冰或浅水覆盖的地域甚至地表下目标成像；

(3)SAR的分辨率与雷达工作波长、平台飞行高度、雷达作用距离无关，在太空或高空都能有效工作；

(4)合成孔径雷达作用距离远，测绘带宽；

(5)可以实现对地物进行多参数、多频段、多极化和多视角测绘。

目前，SAR***已广泛应用于资源勘探、战场侦察、环境保护、灾情检测、水文地质等国防和国民经济的重要领域，并在国民经济发展和军事领域中发挥着越来越重要的作用。

常规SAR采用单通道对场景大面积静止目标的成像，只能获取二维图像，不能得到目标的三维信息，为此，发展了干涉合成孔径雷达(Interferometric Synthetic ApertureRadar，以下用InSAR表示)。InSAR在常规SAR基础上添加了一个SAR数据获取渠道，使得对地面同一场景的两次观测存在着一个视角差异，通过提取复数据的相位信息为信息源，获取地表的三维信息和变化信息。InSAR通过两副天线同时观测或两次近似平行观测获取地面同一景观的复图像对。由于目标与天线位置的几何关系，在复图像上产生了相位差，形成了干涉相位纹图。干涉相位纹图中包含了斜距向上的点与两天线位置之差的精确信息。因此，利用传感器高度、雷达波长、视距视向以及天线基线之间的几何关系，可以精确的测量出图像上每一点的三维位置和变化信息，从而获得观测地区的三维图像。目前，InSAR已经广泛应用于地形测量、地球动力学应用、冰川研究、森林调查与制图、海洋测绘等方面。

此外，常规SAR采用单一极化的发射与接收天线获得数据，只能获得散射波矢量的一个分量，不能完整描述目标的散射特性，在相当大程度上损失了回波中所包含的目标极化信息。极化合成孔径雷达(Polarimetric Synthetic Aperture Radar，以下用PolSAR表示)弥补了常规SAR这一不足。目标对不同的极化电磁波具有不同的散射特性。PolSAR通过使用不同极化的天线发射与接收电磁波，获得同一观测场景的VV、HH、VH、HV复图像，测量目标在不同极化通道下的散射特性，进而获得目标的复散射矩阵，能够完全描述目标散射回波的幅度与相位信息。由于不同极化通道的电磁回波在幅度和相位上均存在差异，目标复散射矩阵中所蕴含的信息量大大超过了目标雷达散射截面积中所蕴含的信息量，因此，PolSAR极大的提高了对目标信息的获取能力，从而为更加深入研究目标散射特性提供了重要依据。通过对目标的极化散射信息处理，例如散射矩阵元素统计特性分析，目标散射矩阵分解和目标散射成分分析，求解目标最优极化问题等，不仅可以分析散射单元内部的几何特性，而且可以调整天线极化状态使感兴趣的目标获得增强。目前，极化SAR主要应用于图像的滤波与分类研究，地形地表参数的获取，作物分析、水域观测、土地资源利用、环境监测、岩性分类等方面。

在SAR干涉测量中，由于电磁波在传播过程中存在衰减，且传播路径上的噪声与杂波干扰以及发射接收天线的热噪声干扰，使得信噪比很低，严重污染了干涉相位纹图，导致大量残差点的产生，使得相位解缠较为困难。极化干涉合成孔径雷达(PolarimetricInterferometric Synthetic Aperture Radar，以下用PolInSAR表示)综合了PolSAR和InSAR的优点，利用多个极化通道的复图像，通过复图像数据融合，增强相位信息，提高数据的信噪比，有效改善了干涉相位纹图质量，残差点数量大大减少，为最终得到高质量的三维图像起到重要的作用。PolInSAR三维成像处理过程的关键技术是相位增强。在PolInSAR成像中，每个像素的散射单元数据包含有两个极化散射矩阵或散射矢量，对应于两根空间分离的天线，因此，极化信息就可以用来增强相干性和改善两根天线接收信号之间的相位。在相位增强过程中，将PolInSAR每个散射单元的极化散射矩阵或者散射矢量，分为主、从散射矩阵或者散射矢量，进行优化处理，可以得到高质量的干涉相位纹图，然后根据InSAR三维成像方法得到优于常规InSAR的高质量三维图像。目前，PolInSAR的研究内容主要集中在树高反演，城市变化、矿区塌陷、山体滑坡地震和火山等引起的地表形变，海浪、海流运动等方面。

1998年S.R.Cloude和K.P.Papathanassiou在期刊IEEE Transactions on Geoscienceand Remote Sensing第36卷1551-1565页发表了论文《Polarimetric SAR interferometry》。文中提出基于干涉相关系数优化的双矢量的相干最优化方法，该算法能够最优化干涉相干性，能够有效改善相位质量，但是对于信号较弱的区域依然存在较多的残差点。

2003年POLinSAR 03 Proceedings会议上E.Colin，C.Titin-schnaider和W.Tabbara发表了论文《Investigation on different interferometric coherence optimizationmethods》。文中介绍了基于干涉相关系数优化的单矢量的相干最优化方法，并给出了理想情况下的解析解。但由于实际情况常常和理想情况有较大的偏差，因此并不能很好的改善相位质量。

2007年11月Progress in Electromagnetic Research Symposium会议上Tao XIONG，Jian YANG，Weijie ZHANG等人发表了论文《Coherence Enhancement for Polarimetric SARInterferometry》。文中讨论了基于幅度优化的PolInSAR的相位增强算法，主要讨论了在单视情况下干涉相位改善的情况。该算法得到的相位较原始各极化通道相位有一定的改善，但由于没有考虑多视信息，相位仍存在着大量的噪声，且文中没有给出解析解的完整表达形式。

通过对美国专利商标局USPTO，欧洲专利局EPO和日本专利局JPO的检索，没有找到类似的专利。

发明内容

本发明要解决的问题是提供一种基于信号幅度最大化干涉相位增强的方法的极化干涉合成孔径雷达三维成像方法。该方法利用复信号幅度和相位之间的关系，即较小幅度的信号具有不可靠的相位，相位的可信度可以通过最大化信号的幅度来改善。PolInSAR提供了两个散射矩阵或散射矢量，对于每一个像素，两个极化散射矢量同时向同一个最优方向投影，使得两个投影长度较小的信号幅度最大化，从而实现相位增强，改善干涉相位纹图，然后采用干涉SAR成像方法处理，得到高质量三维图像。

基于极化数据融合的极化干涉合成孔径雷达三维成像方法，该方法步骤如下：

1.对于极化合成孔径雷达干涉测量，对于每个散射单元，有两个极化散射矩阵[S1]和[S2]，或者两个散射矢量

$k_i=\frac{1}{\sqrt{2}}{(s_{\mathrm{HH}}+s_{\mathrm{VV}},s_{\mathrm{HH}}-s_{\mathrm{VV}},s_{\mathrm{HV}}+s_{\mathrm{VH}})}^T,i=\mathrm{1,2}---\left(1\right)$

其中^T表示矩阵转置操作，s_ij(i，j＝H or V)表示在HV极化基下j极化方式发射、i极化方式接收的复散射系数。这里仅考虑互易情况，即s_HV＝s_VH

定义6×6的相干矩阵[T]：

$\left[T\right]=\&lang;\left[\begin{array}{c}{k}_1\\ k_2\end{array}\right]\left[\begin{array}{cc}{k}_1^H& k_2^H\end{array}\right]\&rang;=\left[\begin{array}{cc}\left[T_{11}\right]& \left[{\mathrm{\&Omega;}}_{12}\right]\\ \left[{\mathrm{\&Omega;}}_{21}\right]& \left[T_{22}\right]\end{array}\right]---\left(2\right)$

其中^H表示复共轭转置。

2.为了将标量形式扩展到矢量形式，引入一个单位复矢量w，于是极化散射矢量k₁和k₂到w的投影η₁和η₂为如下两个复信号：

η₁＝w^Hk₁，η₂＝w^Hk₂                 (3)

本方法的目的是找到一个最优矢量w来同时优化η₁和η₂的幅度，即最大化η₁和η₂中较小的幅度。如果两个幅度都变大，则干涉相位将被改善。

上述优化问题可以由如下数学模型描述：

$\underset{w}{\max }\left(\min (\left|w^H{k}_1\right|,\left|w^H{k}_2\right|)\right)---\left(4\right)$

s.t.‖w‖＝1

为了得到上述问题的解析解，将式(4)转化为如下等价问题。

max(a，b)

$s.t.a=\underset{w}{\max }{\left|w^H{k}_1\right|}^2,\mathrm{if}{\left|w^H{k}_1\right|}^2\&le;{\left|w^H{k}_2\right|}^2$

$b=\underset{w}{\max }{\left|w^H{k}_2\right|}^2,\mathrm{if}{\left|w^H{k}_1\right|}^2>{\left|w^H{k}_2\right|}^2$                   (5)

‖w‖＝1

根据二次规划理论，式(5)的解存在两种情况。

3.对于第一种情况，无约束的最大值|w^Hk₁|²或者|w^Hk₂|²位于约束区域内，则最优矢量w与k₁和k₂中具有较小长度的那个方向相同：

$w^{\left(1\right)}=\left\{\begin{array}{c}{k}_1/\left|\right|k_1\left|\right|,\mathrm{if}\left|\right|k_1\left|\right|\&le;\left|\right|k_2\left|\right|\\ k_2/\left|\right|k_2\left|\right|,\mathrm{if}\left|\right|k_1\left|\right|>\left|\right|k_2\left|\right|\end{array}\right.---\left(6\right)$

4.对于第二种情况，无约束的最大值|w^Hk₁|²和|w^Hk₂|²均不在约束区域内。此时a和b位于边界上。因此当w为最优投影方向时，两个投影具有相同的幅度：

|w^Hk₁|＝|w^Hk₂|             (7)

利用特征值分解来得到解析解。式(7)可以写成

$w^H(k_1{k}_1^H-k_2{k}_2^H)w=0---\left(8\right)$

如果k₁＝ck₂(c是任意非零复数)，则将归为第一种情况。

当k₁≠ck₂时，令[A]表示矩阵

于是[A]的秩为2。另外[A]是一个不定矩阵。因此它有三个特征值λ₁>0、λ₂<0以及λ₃＝0，对应的特征矢量分别是v₁、v₂和v₃。

如果w＝v₃，则完全符合式(8)，因为λ₃＝0。但是由k₁和k₂张成的空间，即Span{k₁，k₂}与Span{v₁，v₂}是等价的，且v₃与v₁和v₂正交，因此|w^Hk₁|＝|w^Hk₂|＝0。这不满足式(4)的目标。

因此w可以表示成为v₁和v₂的线性组合。

w＝β(αv₁+v₂)                    (9)

其中α是一个复系数，β是实的归一化系数。将式(9)代入式(8)，

$w^H(k_1{k}_1^H-k_2{k}_2^H)w={\mathrm{\&beta;}}^2{({\mathrm{\&alpha;v}}_1+v_2)}^H({\mathrm{\&lambda;}}_1{v}_1{v}_1^H+{\mathrm{\&lambda;}}_2{v}_2{v}_2^H)({\mathrm{\&alpha;v}}_1+v_2)$

                                                         (10)

$={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{v}_1^H{v}_1{v}_1^H{v}_1{\mathrm{\&lambda;}}_1+v_2^H{v}_2{v}_2^H{v}_2{\mathrm{\&lambda;}}_2)={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}_2^4{\mathrm{\&lambda;}}_1+{\left|\right|v_2\left|\right|}_2^4{\mathrm{\&lambda;}}_2)=0$

因此

$\left|\mathrm{\&alpha;}\right|=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\frac{{\left|\right|v_2\left|\right|}_2^2}{{\left|\right|v_1\left|\right|}_1^2}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}---\left(11\right)$

将式(9)代入 $w^H{k}_1{k}_1^{H}w,$

$w^H{k}_1{k}_1^{H}w={\mathrm{\&beta;}}^2{({\mathrm{\&alpha;v}}_1+v_2)}^H{k}_1{k}_1^H({\mathrm{\&alpha;v}}_1+v_2)$

$={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{\left|v_1^H{k}_1\right|}^2+{\mathrm{\&alpha;v}}_2^H{k}_1{k}_1{k}_1^H{v}_1+{\left({\mathrm{\&alpha;v}}_2^H{k}_1{k}_1^H{v}_1\right)}^{*}+{\left|v_2^H{k}_1\right|}^2)---\left(12\right)$

由柯西不等式，式(12)的最大值对应着 $\mathrm{\&alpha;}=d{v}_1^H{k}_1{k}_1^H{v}_2,$ 其中d是一个非零实数。因此α的幅角为

$\arg \left(\mathrm{\&alpha;}\right)=\arg \left(v_1^H{k}_1{k}_1^H{v}_2\right)---\left(13\right)$

在将式(11)和式(13)代入(9)并归一化之后，最后将得到：

$w^{\left(2\right)}=\frac{{\mathrm{\&alpha;v}}_1+v_2}{\left|\right|{\mathrm{\&alpha;v}}_1+v_2\left|\right|},$ $\mathrm{\&alpha;}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\exp \{-j\arg \left(v_1^H{k}_1{k}_1^H{v}_2\right)\}---\left(14\right)$

为了推导w⁽¹⁾和w⁽²⁾的判断条件，定义两个函数f₁(|α|)和f₂(|α|)，仅含有一个变量|α|

$f_i\left(\left|\mathrm{\&alpha;}\right|\right)=\frac{{\left|w^H{k}_i\right|}^2}{{\left|\right|w\left|\right|}^2}=\frac{{\left|{({\mathrm{\&alpha;v}}_1+v_2)}^H{k}_i\right|}^2}{{\left|\right|{\mathrm{\&alpha;v}}_1+v_2\left|\right|}^2}=\frac{{\left|\mathrm{\&alpha;}\right|}^2{\left|v_1^H{k}_i\right|}^2+{\mathrm{\&alpha;}}^{*}{v}_1^H{k}_i{k}_i^H{v}_2+{\mathrm{\&alpha;v}}_2^H{k}_i{k}_i^H{v}_1+{\left|v_2^H{k}_i\right|}^2}{{\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}^2+{\mathrm{\&alpha;}}^{*}{v}_1^H{v}_2+{\mathrm{\&alpha;v}}_2^H{v}_1+{\left|\right|v_2\left|\right|}^2}$

                                                               (15)

$=\frac{{(\left|\mathrm{\&alpha;}\right|\left|v_1^H{k}_i\right|+\left|v_2^H{k}_i\right|)}^2}{{\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}^2+{\left|\right|v_2\left|\right|}^2},i=\mathrm{1,2}$

令

$\frac{{\mathrm{df}}_i\left(\left|\mathrm{\&alpha;}\right|\right)}{d\left|\mathrm{\&alpha;}\right|}=d\frac{{(\left|\mathrm{\&alpha;}\right|\left|v_1^H{k}_i\right|+\left|v_2^H{k}_i\right|)}^2}{{\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}^2+{\left|\right|v_2\left|\right|}^2}/d\left|\mathrm{\&alpha;}\right|=d\frac{z{\left|\mathrm{\&alpha;}\right|}^2+w\left|\mathrm{\&alpha;}\right|+u}{x{\left|\mathrm{\&alpha;}\right|}^2+y}/d\left|\mathrm{\&alpha;}\right|$

                                                      (16)

$=\frac{(2z\left|\mathrm{\&alpha;}\right|+w)(x{\left|\mathrm{\&alpha;}\right|}^2+y)-2\left|\mathrm{\&alpha;}\right|x(z{\left|\mathrm{\&alpha;}\right|}^2+w\left|\mathrm{\&alpha;}\right|+u)}{{(x{\left|\mathrm{\&alpha;}\right|}^2+y)}^2}=0$

于是

wx|α|²+2(xu-zy)|α|-wy＝0                         (17)

|α|对应的解为

$\left|\mathrm{\&alpha;}\right|=\frac{2\mathrm{zy}}{\mathrm{wx}}=\frac{\left|v_1^H{k}_i\right|}{\left|v_2^H{k}_i\right|}>0---\left(18\right)$

由于f₁和f₂的单调性，存在着如下两种情况(示意图如图6所示)：

如果 $\left|\mathrm{\&alpha;}\right|<\left|v_1^H{k}_1\right|/\left|v_2^H{k}_1\right|$ 以及 $\left|\mathrm{\&alpha;}\right|<\left|v_1^H{k}_2\right|/\left|v_2^H{k}_2\right|,$ 则式(5)中a大于b，且w＝w⁽¹⁾。

如果 $\left|\mathrm{\&alpha;}\right|>\left|v_1^H{k}_1\right|/\left|v_2^H{k}_1\right|$ 以及 $\left|\mathrm{\&alpha;}\right|>\left|v_1^H{k}_2\right|/\left|v_2^H{k}_2\right|,$ 则式(5)中b大于a，且w＝w⁽¹⁾。

否则w＝w⁽²⁾

因此，第二种情况下的解为

$w^{\left(2\right)}=\frac{{\mathrm{\&alpha;v}}_1+v_2}{\left|\right|{\mathrm{\&alpha;v}}_1+v_2\left|\right|},$ $\mathrm{\&alpha;}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\exp \{-j\arg \left(v_1^H{k}_1{k}_1^H{v}_2\right)\}---\left(19\right)$

其中λ₁>0和λ₂<0是矩阵

的两个非零特征值，v₁和v₂是对应的特征矢量。

5.综上所述，模型式(4)最终的解为

$w=\left\{\begin{array}{c}{w}^{\left(2\right)},\mathrm{if}(\sqrt{\frac{-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}>\max (\frac{\left|v_1^H{k}_1\right|}{\left|v_2^H{k}_1\right|},\frac{\left|v_1^H{k}_2\right|}{\left|v_2^H{k}_2\right|}))\mathrm{or}(\sqrt{\frac{{-\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}<\min (\frac{\left|v_1^H{k}_1\right|}{\left|v_2^H{k}_1\right|},\frac{\left|v_1^H{k}_2\right|}{\left|v_2^H{k}_2\right|}))\\ w^{\left(1\right)},\mathrm{else}\end{array}\right.---\left(20\right)$

6.确定w后，即可计算干涉相位。在单视情况下，融合的相位为

${\mathrm{\&phi;}}_s=\arg \left({\mathrm{\&eta;}}_1{\mathrm{\&eta;}}_2^{*}\right)=\arg \left(w^H{k}_1{k}_2^{H}w\right)---\left(21\right)$

多视情况下，融合的相位为

${\mathrm{\&phi;}}_m=\arg ({\&lang;\mathrm{\&eta;}}_1{\mathrm{\&eta;}}_2^{*}\&rang;)=\arg \left(w^H\left[{\mathrm{\&Omega;}}_{12}\right]w\right)---\left(22\right)$

与现有技术相比，本发明的优点如下：

(1)其关键在于最大化信号的幅度，这是基于复信号的幅度与相位之间的关系的，即在统计意义下，信号的幅度越大，则对应的相位越可信。图7给出了信号幅度和反映相位质量的相干的二维统计直方图。在单视情况下，增强的相位与两个散射矢量的相似性参数的相位是等价的。它可以认为是全极化信息的加权平均，权重与各个通道的信号幅度成正比。

(2)利用极化干涉数据对改善相位质量的性能进行了验证。本方法的计算时间比相干最优化方法的时间少，且融合后的相位质量要优于相干最优化方法，尤其是在弱信号区域。

(3)对于多视情况，利用本方法，干涉相位中超过99％的残差点被去除，且地形变化的细节信息能够更加清晰的观察到。这大大降低了相位解缠的难度，并在高精度数字高程模型反演中将起到重要的作用。

附图说明

图1 为合成孔径雷达工作示意图；

图2 为合成孔径雷达对点目标的观测示意图；

图3 为合成孔径雷达与点目标的位置几何关系示意图；

图4 为合成孔径雷达采样得到的数据阵列示意图；

图5 为合成孔径雷达到目标的斜距变化示意图；

图6 为模型解的两种情况的示意图；

图7 为相干和幅度的二维统计直方示意图；

图8(a) 为植被和裸地区域HH通道的幅度示意图；

   (b) 为植被和裸地区域的光学示意图；

图9(a) 为山地测试区域的HH通道幅度示意图；

   (b) 为山地测试区域的光学示意图；

   (c) 为山地测试区域的HH通道的相位示意图；

   (d) 为由幅度最优化方法得到的相位示意图；

图10(a) 为图8(a)中方框A放大区域的HH通道的幅度示意图；

    (b) 为HH通道的相位示意图；

    (c) 为由幅度最优化方法得到的幅度示意图；

    (d) 为由幅度最优化方法得到的相位示意图；

    (e) 为由相干最优化方法得到的相干示意图；

    (f) 为由相干最优化方法得到的相位示意图；

图11(a) 为图8(a)中方框B放大区域的HH通道的幅度示意图；

    (b) 为HH通道的相位示意图；

    (c) 为HH通道的残差点示意图；

    (d) 为由幅度最优化方法得到的相位示意图；

    (e)为由幅度最优化方法得到的残差点示意图；

    (f)为由相干最优化方法得到的相位示意图；

图12(a)为由幅度最优化方法得到的对应图9(a)中方框区域的相位示意图；

    (b)为由相干最优化方法得到的相位示意图；

    (c)为由幅度最优化方法得到的残差点示意图；

    (d)为由相干最优化方法得到的残差点示意图；

图13(a)为图8对应的植被地区的去平地效应后的相位示意图；

    (b)为反映地形的解缠相位示意图；

图14(a)为图9对应的山地地区的去平地效应后的相位示意图；

    (b)为反映地形的解缠相位示意图；

图15是发明的流程示意图。

附表说明

表1图10(a)中区域的比较。

表2图9(a)中区域的比较。

具体实施方式

下面结合附图对本发明作进一步说明。

为验证本方法的有效性，本发明采用SIR-C/X-SAR雷达***的L波段全极化单视复图像数据对算法进行验证。

该数据是由SIR-C/X-SAR雷达***于1994年10月8日到9日采集的位于天山测试区域的L波段全极化单视复图像对。测试区域位于俄罗斯的贝加尔湖的东南岸，包含不同的地物，如树林、庄稼地、裸地和山地。未经任何预处理之前，干涉相位被强噪声严重污染，产生大量的残差点。

图15是发明的流程示意图。在实例中，对数据中的植被和裸地、山地分别进行处理，通过与相干最优方法的对比说明本发明的方法普适性和有效性。根据上述说明，该方法步骤如下：

(1)读入极化干涉数据。

图8(a)所示的是植被和裸地测试区域(1000×1000像素)的HH通道的幅度|sHH|，包含了多种不同的地物，如树林(F)、道路(R)、裸地(BG)以及庄稼地(C)。由Google地图得到的对应的相同分辨率的光学图像由图8(b)所示。图9(a)所示的是一块包含山地的测试区域，其对应的光学图像由图9(b)所示，来源同样是Google Map。

(2)对于每一个散射单元，利用极化干涉数据按式(1)构造散射矢量k₁和k₂。

(3)对矩阵

进行特征值分解，得到三个特征值λ₁>0、λ₁>0和λ₃＝0以及对应的特征向量v₁、v₂和v₃。

(4)分别计算以及的值。根据式(20)比较

和

以及

之间的关系，判断解属于哪一种情况。若属于第一种情况，则按照式(6)计算w＝w⁽¹⁾；若属于第二种情况，则按照式(14)计算w＝w⁽²⁾。

(5)为得到增强后的单视相位，按照式(21)计算 ${\mathrm{\&phi;}}_s=\arg \left(w^H{k}_1{k}_2^{H}w\right).$ 并根据式(3)得到增强后的复图像。

对于植被和裸地，图8(a)中方框A的放大图如图10(a)所示。地面大部分被树林覆盖，其中有三条道路通过。黑色区域为裸地。由于噪声和去相干的影响，HH通道的噪声非常严重，导致地形细节被淹没(图10(b))。HH、HV和VV通道的平均较小幅度分别为0.3604、0.1401和0.2578。利用本方法，融合后的图像对的较小平均幅度达到0.4768，且式(3)中η1的幅度如图10(c)所示。它明显比HH通道的幅度图要“亮”一些。

图11(a)所示的是图8(a)中方框B的放大区域。大部分地面被矮植被如庄稼和草地所覆盖。其中平行的直线看起来像田埂。田埂和地面之间形成的二面角导致强烈的回波信号。HH通道的相位(如图11(b)所示)噪声很大以至于有6550个残差点存在，如图11(c)。

对于山地HH通道的干涉相位如图9(c)所示，在1000×1000个像素中，对应着95329个残差点，因此其密度接近于10％。将全极化通道的信息融合，图像对的平均较小幅度可以从HH、HV和VV通道的0.3250、0.1310和0.2831增大到0.4574。

由于增强后的单视相位仍存在着较大的噪声，因此这里不作考虑。

(6)为得到增强后的多视相位，首先计算矩阵 ${\mathrm{\&Omega;}}_{12}=\&lang;k_1{k}_2^H\&rang;,$ 然后按照式(22)计算φ_m＝arg(w^H[Ω₁₂]w)。

对于植被和裸地的区域A，融合后相位(图10(d))噪声被明显去除，相位的跳变能够清晰的分辨出来。99.76％的残差点被去除。值得注意的是，增强后，树林和裸地区域的相位没有明显的边界。这意味着树林地区的相位可以被认为是地面相位，因为裸地的相位肯定对应着地面相位。

作为比较，由相干最优化方法得到的相位如图10(f)所示。虽然对应的相干值被优化接近于1(如图10(e)所示)，且树林区域大部分噪声被去除，但裸地的相位依然受到噪声的影响，这与由光学图像观察到的以及雷达图像种地幅度(一次散射)所反映的裸地的平坦性并不相符。

表1列出了更多关于原始数据、由幅度最优化方法和相干最优化方法增强的数据之间的比较，包括平均较小幅度，平均相干值以及单视和多视情况下的残差点数。

对于植被和裸地的区域B，利用幅度最优化和相干最优化方法增强的相位分别如图11(d)和11(f)所示。明显地，由幅度最优化方法得到的相位具有最好的质量以及最小的噪声。99.63％的残差点被成功的去除，如图11(e)，证明了本方法的优势。虽然相干最优化方法得到的相干值最高，并接近于1，但对应的相位仍然受到噪声和去相干的影响。

对于山地，融合的相位如图9(d)所示。99.96％的残差点被去除，只剩下34个。虽然对于具有中等和强信号的山地，相干最优化方法能够有效的改善相位质量，但对于具有微弱信号的平地，融合的相位仍然存在着大量的残差点。注意图9(a)中的白框中的区域，右半部分幅度较小。由幅度最优化和相干最优化得到的增强相位分别示于图12(a)和12(b)。对应于低幅度的区域，图12(d)的右半部分存在大量的残差点，图12(b)中对应的相位受到噪声的干扰。

基于最大化信号的幅度，由本方法得到的相位无论在强信号区域还是弱信号区域都存在很少的残差点，如图12(c)所示。这表明本方法的稳健性。

更多的比较列于表2中。

(7)对于所得到的整幅相位图(单视或多视)去平地效应。高度相同的平地在干涉相位图中形成的干涉条纹随距离向和方位向的变化而呈周期性变化的现象称为“平地效应”现象。“平地效应”使干涉相位图不能直观反映地形变化，同时给降噪和相位解缠带来了困难，因此必须进行“平地效应”消除处理。“平地效应”可以通过对干涉纹乘以复相位函数来去除。

设基线距为B，波长为λ，下视角为θ₀，基线与水平方向的夹角为α，斜距为R，地距方向两点到天线的距离差为ΔR，则这两点之间的相位差为

${\mathrm{\&Delta;\&phi;}}_R=-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{B\cos ({\mathrm{\&theta;}}_0-\mathrm{\&alpha;})\mathrm{\&Delta;R}}{R\tan {\mathrm{\&theta;}}_0}$

即由平地引起的相位。以距离向某个点为基准相位，依上式计算整个区域的平地相位φ_flat，则去平地效应后的相位为

φ_m，f＝W(φ_m-φ_flat)或φ_s，f＝W(φ_s-φ_flat)

图13(a)和图14(a)分别给出了植被和裸地区域、山地区域去平地效应之后的相位图。

(8)对去平地效应后的相位进行相位解缠，由于增强后的相位噪声被大大减小，因此对应的参差点数非常少，因此利用Goldstein提出的枝切法进行相位解缠。大致分为四个步骤：1.计算相位的参差点；2.根据参差点的位置和电荷连接枝切线；3.进行相位梯度积分；4.进行后处理，将枝切线上的像素进行解缠得到最终的结果。

图13(b)和图14(b)分别给出了植被和裸地区域、山地区域最终的解缠相位图。从图中可以看到，地形变化能够更清晰的被反映出来，且细节信息保持得很好。

表1图10(a)中区域的比较

 

	平均较小幅度	平均相干值	残差点数
				HH通道	0.3604	0.7843	8331
HV通道	0.1401	0.7077	12033
				VV通道	0.2578	0.7450	10279
幅度最优化，单视	0.4768	0.8523	1699
				相干最优化，单视	0.2795	0.8858	1889
幅度最优化，多视	--------	0.8228	20
				相干最优化，多视	--------	0.9634	173

表2图9(a)中区域的比较

 

	平均较小幅度	平均相干值	残差点数
				HH通道	0.3250	0.7632	95329
HV通道	0.1310	0.6974	125954
				VV通道	0.2831	0.7536	100533
幅度最优化，单视	0.4574	0.8309	21851
				相干最优化，单视	0.2535	0.8623	24788
幅度最优化，多视	--------	0.7998	34
				相干最优化，多视	--------	0.9551	2078

Claims (1)

1.基于极化数据融合的极化干涉合成孔径雷达三维成像方法，其特征在于，该方法步骤如下：

(1)对于极化合成孔径雷达干涉测量，每个散射单元有两个极化散射矩阵[S₁]和[S₂]，或者两个散射矢量

$k_i=\frac{1}{\sqrt{2}}{[s_{\mathrm{HH}}+s_{\mathrm{VV}},s_{\mathrm{HH}}-s_{\mathrm{VV}},s_{\mathrm{HV}}+s_{\mathrm{VH}}]}^T,i=\mathrm{1,2}---1)$

其中，^T表示矩阵转置操作，s_ij(i，j＝H or V)表示在HV极化基下j极化方式发射、i极化方式接收的复散射系数；这里我们仅考虑互易情况，即s_HV＝s_VH，

定义6×6的相干矩阵[T]：

$\left[T\right]=\&lang;\left[\begin{array}{c}{k}_1\\ k_2\end{array}\right]\left[\begin{array}{cc}{k}_1^H& k_2^H\end{array}\right]\&rang;=\left[\begin{array}{cc}\left[T_{11}\right]& \left[{\mathrm{\&Omega;}}_{12}\right]\\ \left[{\mathrm{\&Omega;}}_{21}\right]& \left[T_{22}\right]\end{array}\right]---2)$

其中^H表示复共轭转置；

(2)引入单位复矢量w，极化散射矢量k₁和k₂到w的投影η₁和η₂为如下两个复信号：

η₁＝w^Hk₁，η₂w^Hk₂                            3)

建立数学模型：

$\underset{w}{\max }\left(\min (\left|w^H{k}_1\right|,\left|w^H{k}_2\right|)\right)$

$\begin{array}{cc}s.t.& \left|\right|w\left|\right|=1\end{array}$                                   4)

为了得到上述问题的解析解，将式4)转化如下：

$\max (a,b)$

$s.t.a=\underset{w}{\max }{\left|w^H{k}_1\right|}^2,\mathrm{if}{\left|w^H{k}_1\right|}^2\&le;{\left|w^H{k}_2\right|}^2$

                                             5)

$b=\underset{w}{\max }{\left|w^H{k}_2\right|}^2,\mathrm{if}{\left|w^H{k}_1\right|}^2>{\left|w^H{k}_2\right|}^2$

$\left|\right|w\left|\right|=1$

根据二次规划理论，式5)的解存在两种情况：

(2.1)对于第一种情况，无约束的最大值|w^Hk₁|²或者|w^Hk₂|²位于约束区域内，则最优矢量w与k₁和k₂中具有较小长度的那个方向相同：

$w^{\left(1\right)}=\left\{\begin{array}{c}{k}_1/\left|\right|k_1\left|\right|,\mathrm{if}\left|\right|k_1\left|\right|\&le;\left|\right|k_2\left|\right|\\ k_2/\left|\right|k_2\left|\right|,\mathrm{if}\left|\right|k_1\left|\right|>\left|\right|k_2\left|\right|\end{array}\right.---6)$

(2.2)对于第二种情况，无约束的最大值|w^Hk₁|²和|w^Hk₂|²均不在约束区域内；此时a和b位于边界上；因此当w为最优投影方向时，两个投影具有相同的幅度：

|w^Hk₁|＝|w^Hk₂|                         7)

利用特征值分解来得到解析解，式7)写成

$w^H(k_1{k}_1^H-k_2{k}_2^H)w=0---8)$

如果k₁＝ck₂，c是任意非零复数，则将归为第一种情况；

当k₁≠ck₂时，令[A]表示矩阵于是[A]的秩为2；另外[A]是一个不定矩阵；因此它有三个特征值λ₁>0、λ₂<0以及λ₃＝0，对应的特征矢量分别是v₁、v₂和v₃；

如果w＝v₃，则完全符合式8)，因为λ₃＝0；但是由k₁和k₂张成的空间，即Span{k₁，k₂}与Span{v₁，v₂}是等价的，且v₃与v₁和v₂正交，因此|w^Hk₁|＝|w^Hk₂|＝0；这不满足式4)的目标；

因此w表示成为v₁和v₂的线性组合为

w＝β(αv₁+v₂)                            9)

其中，a是一个复系数，β是实的归一化系数；将式9)代入式8)，

$w^H(k_1{k}_1^H-k_2{k}_2^H)w={\mathrm{\&beta;}}^2{(\mathrm{\&alpha;}{v}_1+v_2)}^H({\mathrm{\&lambda;}}_1{v}_1{v}_1^H+{\mathrm{\&lambda;}}_2{v}_2{v}_2^H)(\mathrm{\&alpha;}{v}_1+v_2)$

                                        10)

$={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{v}_1^H{v}_1{v}_1^H{v}_1{\mathrm{\&lambda;}}_1+v_2^H{v}_2{v}_2^H{v}_2{\mathrm{\&lambda;}}_2)={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}_2^4{\mathrm{\&lambda;}}_1+{\left|\right|v_2\left|\right|}_2^4{\mathrm{\&lambda;}}_2)=0$

因此

$\left|\mathrm{\&alpha;}\right|=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\frac{{\left|\right|v_2\left|\right|}_2^2}{{\left|\right|v_1\left|\right|}_1^2}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}---11)$

将式9)代入

$w^H{k}_1{k}_1^{H}w={\mathrm{\&beta;}}^2{(\mathrm{\&alpha;}{v}_1+v_2)}^H{k}_1{k}_1^H(\mathrm{\&alpha;}{v}_1+v_2)$

$={\mathrm{\&beta;}}^2({\left|\mathrm{\&alpha;}\right|}^2{\left|v_1^H{k}_1\right|}^2+\mathrm{\&alpha;}{v}_2^H{k}_1{k}_1^H{v}_1+{\left(\mathrm{\&alpha;}{v}_2^H{k}_1{k}_1^H{v}_1\right)}^{*}+{\left|v_2^H{k}_1\right|}^2)---12)$

由柯西不等式，式12)的最大值对应着 $\mathrm{\&alpha;}=d{v}_1^H{k}_1{k}_1^H{v}_2,$ 其中d是一个非零实数；因此a的幅角为

$\arg \left(\mathrm{\&alpha;}\right)=\arg \left(v_1^H{k}_1{k}_1^H{v}_2\right)---13)$

在将式11)和式13)代入9)并归一化之后，最后将得到：

$w^{\left(2\right)}=\frac{\mathrm{\&alpha;}{v}_1+v_2}{||\mathrm{\&alpha;}{v}_1+v_2||},\mathrm{\&alpha;}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\exp \{-j\arg \left({v_1}^H{k}_1{k_1}^H{v}_2\right)\}---14)$

为了推导w⁽¹⁾和w⁽²⁾的判断条件，定义两个函数f₁(|α|)和f₂(|α|)，仅含有一个变量|α|

$f_i\left(\left|\mathrm{\&alpha;}\right|\right)=\frac{{\left|w^H{k}_i\right|}^2}{{\left|\right|w\left|\right|}^2}=\frac{{\left|{({\mathrm{\&alpha;v}}_1+v_2)}^H{k}_i\right|}^2}{{\left|\right|{\mathrm{\&alpha;v}}_1+v_2\left|\right|}^2}=\frac{{\left|\mathrm{\&alpha;}\right|}^2{\left|v_1^H{k}_i\right|}^2+{\mathrm{\&alpha;}}^{*}{v}_1^H{k}_i{k}_i^H{v}_2+{\mathrm{\&alpha;v}}_2^H{k}_i{k}_i^H{v}_1+{\left|v_2^H{k}_i\right|}^2}{{\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}^2+{\mathrm{\&alpha;}}^{*}{v}_1^H{v}_2+{\mathrm{\&alpha;v}}_2^H{v}_1+{\left|\right|v_2\left|\right|}^2}---15)$

$=\frac{{(\left|\mathrm{\&alpha;}\right|\left|v_1^H{k}_i\right|+\left|v_2^H{k}_i\right|)}^2}{{\left|\mathrm{\&alpha;}\right|}^2{\left|\right|v_1\left|\right|}^2+{\left|\right|v_2\left|\right|}^2},i=\mathrm{1,2}$

令

$\frac{{\mathrm{df}}_i\left(\left|\mathrm{\&alpha;}\right|\right)}{d\left|\mathrm{\&alpha;}\right|}=d\frac{{(\left|\mathrm{\&alpha;}\right|\left|{v_1}^H{k}_i\right|+\left|{v_2}^H{k}_i\right|)}^2}{{\left|\mathrm{\&alpha;}\right|}^2{||v_1||}^2+{||v_2||}^2}/d\left|\mathrm{\&alpha;}\right|=d\frac{z{\left|\mathrm{\&alpha;}\right|}^2+w\left|\mathrm{\&alpha;}\right|+u}{x{\left|\mathrm{\&alpha;}\right|}^2+y}/d\left|\mathrm{\&alpha;}\right|$

$=\frac{(2z\left|\mathrm{\&alpha;}\right|+w)(x{\left|\mathrm{\&alpha;}\right|}^2+y)-2\left|\mathrm{\&alpha;}\right|x(z{\left|\mathrm{\&alpha;}\right|}^2+w\left|\mathrm{\&alpha;}\right|+u)}{{(x{\left|\mathrm{\&alpha;}\right|}^2+y)}^2}=0---16)$

于是

wx|α|²+2(xu-zy)|α|-wy＝0            17)

|α|对应的解为

$\left|\mathrm{\&alpha;}\right|=\frac{2\mathrm{zy}}{\mathrm{wx}}=\frac{\left|v_1^H{k}_i\right|}{\left|v_2^H{k}_i\right|}>0---18)$

由于f₁和f₂的单调性，存在着如下两种情况：

如果 $\left|\mathrm{\&alpha;}\right|<\left|v_1^H{k}_1\right|/\left|v_2^H{k}_1\right|$ 以及 $\left|\mathrm{\&alpha;}\right|<\left|v_1^H{k}_2\right|/\left|v_2^H{k}_2\right|,$ 则式(5)中a大于b，且w＝w⁽¹⁾，

如果 $\left|\mathrm{\&alpha;}\right|>\left|v_1^H{k}_1\right|/\left|v_2^H{k}_1\right|$ 以及 $\left|\mathrm{\&alpha;}\right|>\left|v_1^H{k}_2\right|/\left|v_2^H{k}_2\right|,$ 则式(5)中b大于a，且w＝w⁽¹⁾；

否则w＝w⁽²⁾

因此，第二种情况下的解为

$w^{\left(2\right)}=\frac{\mathrm{\&alpha;}{v}_1+v_2}{\left|\right|\mathrm{\&alpha;}{v}_1+v_2\left|\right|},\mathrm{\&alpha;}=\sqrt{-\frac{{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}\exp \{-j\arg \left(v_1^H{k}_1{k}_1^H{v}_2\right)\}---19)$

其中λ₁>0和λ₂<0是矩阵

的两个非零特征值，v₁和v₂是对应的特征矢量；

(3)因此，模型式4)最终的解为

$w=\left\{\begin{array}{c}{w}^{\left(2\right)},\mathrm{if}(\sqrt{\frac{-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}>\max (\frac{\left|v_1^H{k}_1\right|}{\left|v_2^H{k}_1\right|},\frac{\left|v_1^H{k}_2\right|}{\left|v_2^H{k}_2\right|}))\mathrm{or}(\sqrt{\frac{-{\mathrm{\&lambda;}}_2}{{\mathrm{\&lambda;}}_1}}<\min (\frac{\left|v_1^H{k}_1\right|}{\left|v_2^H{k}_1\right|},\frac{\left|v_1^H{k}_2\right|}{\left|v_2^H{k}_2\right|}))\\ w^{\left(1\right)},\mathrm{else}\end{array}\right.---20)$

(4)确定w后，即可计算干涉相位；

在单视情况下，融合的相位为

${\mathrm{\&phi;}}_s=\arg \left({\mathrm{\&eta;}}_1{\mathrm{\&eta;}}_2^{*}\right)=\arg \left(w^H{k}_1{k}_2^{H}w\right)---21)$

多视情况下，融合的相位为

${\mathrm{\&phi;}}_m=\arg (\&lang;{\mathrm{\&eta;}}_1{\mathrm{\&eta;}}_2^{*}\&rang;)=\arg \left(w^H\left[{\mathrm{\&Omega;}}_{12}\right]w\right)---22)$

Images