CN101551455A - 干涉合成孔径雷达三维地形成像系统及其高程测绘方法 - Google Patents

Abstract

本发明公开了一种干涉合成孔径雷达三维地形成像系统及其高程测绘方法，主要解决现有干涉合成孔径雷达成像实效性差，并且不能对快速变化地形和突变地形进行三维高程测绘的问题。其系统包括：三个子孔径天线、雷达波发射机、雷达波接收机和成像数据处理器，该成像信号处理器包括SAR图像处理单元和InSAR图像处理单元。本发明通过三个子孔径接收雷达回波，然后对三个子孔径分别接收到的雷达回波进行SAR成像处理，再将得到的SAR复图像进行InSAR成像处理，其中包括：图像配准、相位滤波和基于聚类分析的相位展开。最后将处理好的InSAR展开相位图进行高程反演，恢复出三维数字高程图。本发明具有对测绘地形的适应性广，成像实效性强的优点，可用于对三维地形的测绘。

Description

干涉合成孔径雷达三维地形成像系统及其高程测绘方法

技术领域

本发明属于雷达技术领域，具体的说是一种干涉合成孔径雷达系统，可用于机载干涉合成孔径雷达对三维高程地形的测量。

背景技术

干涉合成孔径雷达InSAR是以合成孔径雷达SAR获得的复数据提取的相位信息为信息源来获取地表的三维信息和变化信息的一项技术。传统InSAR通过两副天线同时观测，或两次近平行的观测，获取同一地面场景的复图像对。由于地面目标与两天线位置的几何关系在复图像上产生了相位差，形成包含地面目标与两天线位置之差信息的干涉条纹，因此，利用传感器高度、雷达波长、波束视角及天线基线之间的几何关系，可以精确地测量出SAR图像上每一像素的三维位置和高度变化信息。

干涉合成孔径雷达测量发展30多年来，其应用领域不断地扩展。三维测高方面的应用包括：三维地形高度测量，DEM图的生成；地震学方面的应用，包括震前预测，抗震构造的机理研究；地面的下陷与抬升，海床的扩张等；冰川研究等地球动力学方面的应用，包括冰川地形测量，冰速测量；树高测量，林区地形，森林调查与制图等；海洋测绘方面的应用，包括海波与海流的测量等。

传统的两孔径机载干涉合成孔径雷达，由于只有两个接收孔径，其所能接收到的干涉相位图只有一幅，这也就导致了传统机载干涉合成孔径雷达的三维成像必须依赖干涉相位图内各相邻像素间的相位梯度，这就是造成传统机载干涉合成孔径雷达所能够观测的三维地形必须要求缓慢且连续变化的根本原因。事实上，这种对地形缓慢或连续变化的要求对于实际地形来说是十分苛刻的，这一要求导致传统机载干涉合成孔径雷达能够观测的地形十分有限。例如，在2008年汶川地震时由于汶川地区多陡峭的山峰，而且加上地震导致地貌剧烈形变，传统的机载干涉合成孔径雷达对上述这种地形是无能为力的。

传统的两孔径机载干涉合成孔径雷达，其三维高程成像是基于相邻像素的相位梯度信息，从而导致了相位展开算法复杂度高，运算时间长。由于突发事件例如地震、山体滑坡等对雷达测绘的实效性需求非常高，因而传统的两孔径机载干涉合成孔径雷达很难达到实时测绘成像的需求。

发明内容

本发明的目的在于避免上述传统两孔径机载干涉合成孔径雷达的不足，提供一

种三孔径机载干涉合成孔径雷达三维地形快视成像系统及其高程测绘方法，以实现对各类复杂三维地形的实时、快视观测。

发明实现手段

为实现上述目的，本发明提供的干涉合成孔径雷达三维地形成像系统，采用单航过干涉测量模型，包括：三个子孔径天线、雷达波发射机、雷达波接收机和成像数据处理器，该成像信号处理器包括SAR图像处理单元和InSAR图像处理单元，其中：三个子孔径天线垂直于飞机飞行方向放置，以实现三维高程地形的绘制；InSAR处理单元中包括：图像配准子单元、相位滤波子单元、高程反演子单元和多基线相位展开四个子单元，该多基线相位展开子单元固化有基于聚类分析的相位展开处理程序，实现对所有像素点的逐类相位展开。

所述的三个子孔径天线的任意两个子孔径组成一组，每一组对应一个基线长度，该两个基线长度不同，且两个长度之间的最大公约数为1，且两个基线长度同时满足极限基线、展宽模糊高度、高程精度、相位展开精度这四方面的要求。

为实现上述目的，本发明提供进行三维地形高程测绘的方法，包括如下步骤：

1)将雷达接收到的回波生成三幅二维复图像；

2)对三幅二维复图像进行距离向和方位向偏移的配准，并分别在不同基线长度下提取每两个孔径天线与地面同一目标之间的相位差，得到两幅干涉相位图；

3)对两幅干涉相位图像进行滤除噪声处理；

4)对滤除噪声后的干涉相位图进行基于聚类的多基线相位展开，得到真实相位图；

5)用真实相位图内每一像素点对应的真实相位值，乘以地形高度比例因子，得到三维地形高程图。

本发明具有如下优点

a)适用性广：由于本发明采用了三个子孔径天线，其中每一组对应一个基线长度，该两个基线长度不同，且两个长度之间的最大公约数为1。因而能够大大提高了模糊高度，实现对快速变化地形和突变地形的观测，如对悬崖和楼房的观测。

b)实效性强：由于本发明采用基于聚类分析的相位展开方法，其复杂度远远低于传统的相位展开方法，因此具有明显的实时和快视成像处理能力。

c)抗噪声能力强：由于本发明在InSAR处理单元的相位展开子单元中采用的了基于聚类分析的相位展开方法，该方法通过对像素点进行聚类后再逐类进行相位展开，因此有较强的噪声鲁棒性。

附图说明

图1是本发明的三维高程测绘原理图；

图2是本发明的干涉合成孔径雷达系统示意图；

图3是本发明的天线安装示意图；

图4是本发明的高程测绘流程图；

图5是本发明的仿真干涉条纹图；

图6是本发明的仿真干涉条纹图对应的成像标准结果；

图7是本发明与传统方法的成像效果比较图。

具体实施方式

以下结合附图详细说明本发明的系统结构及高程测量过程。

参照图2，本发明的干涉合成孔径雷达系统包括：三个子孔径天线、雷达波发射机、雷达波接收机和成像数据处理器。雷达在工作时首先由雷达波发射机通过整个孔径天线或一个子孔径发出雷达脉冲，雷达回波分别通过三个子孔径天线各自的雷达波接收机同时接收，并将接收到的信号传送给各自的数据处理器，即雷达在工作时首先由切换开关转换到雷达波发射模式，此时雷达波发射机通过整个孔径天线或者子孔径天线A、子孔径天线B、子孔径天线C中任意一个来发射雷达波，再发射工作完成后雷达通过切换开关切换到雷达波接收模式，此时雷达回波通过子孔径天线A，子孔径天线B，子孔径天线C由接收机A，接收机B，接收机C同时接收雷达。最终由数据处理器内固化的三维高程成像程序完成对三维地形的测绘。

参照图3，本发明将三孔径天线放在飞机腹部，摆放方向与飞机飞行方向垂直。其次，每两个子孔径天线对应一个基线长度，为了满足发明目的的需求，在设计基线长度时按以下要求进行：

1)展宽模糊高度的设计

对于SAR图像中的一个像素来说，真实高度与绝对干涉相位之间的关系为：

h＝a×ψ                                (1)

其中a是绝对干涉相位到高程的定标因子，ψ代表绝对干涉相位，同时

$a=\frac{\mathrm{\λ}\×r\×\sin \left(\mathrm{\θ}\right)}{B\×4\mathrm{\π}}---\left(2\right)$

其中λ是雷达波长，r是天线相位中心到目标像素的斜距，θ是雷达波入射角，B是等效双程有效基线长度。直接测量的模糊干涉相位与绝对干涉相位之间的关系是

其中

k是整数。

由于本发明提出的机载干涉合成孔径雷达有三个孔径，因此三个孔径两两组合就可以得到不同长度的基线。对于同一地面单元的地形高度，则有

其中B₁，B₂代表不同的基线长度，并且k₁，k₂∈integer代表不同长度的基线所对应的高程模糊数。因此，当B₁，B₂的最大公约数为1的时候，模糊高度就会被展宽。

2)极限基线的的设计

在InSAR干涉测量的工作原理中，从测量角度讲基线与SAR成像的斜距构成了一组三角几何关系，而这组三角几何关系是InSAR能够进行高程测绘的基础，而组成这组三角几何关系的一个重要因素就是干涉基线。然而另一方面，基线又是导致图像对相干性损失的一个根源。

对于一个载频为f₀，带宽为B的雷达发射脉冲，其SAR回波从功率谱上看可以理解成一个带宽为

经过

平移后的信号从地面雷达后向散射系数频谱又称为波数谱上切下来的一小段，其中θ为雷达波束的地面入射角，c为光速。

本发明的机载InSAR系统三个孔径间的两两连线构成了基线，由于基线的存在，使得各天线与同一地面单元的入射角记为θ₁，θ₂，θ₃有了微小的差别。如上所述，不同的入射角对应于不同的波数谱移动量，所以三个孔径收集的地面信息的内容是从地面雷达后向散射系数频谱上切下来的不同频段。因此，InSAR回波谱的内容就可以分为各个孔径的公共部分和非公共部分。其中任意两个孔径的波数谱重叠部分的多少就取决于两个频谱移动量的差别Δf，而Δf取决于基线长度B_i，j i，j为两个孔径的编号，Δf与B_i，j的关系式如下：

$\mathrm{\Δf}=\frac{{2f}_0}{c}({\sin \mathrm{\θ}}_j-{\sin \mathrm{\θ}}_i)\≈\frac{{2f}_0\cos {\mathrm{\θ}}_c}{c\·r_0}{B}_{i,j}---\left(5\right)$

其中： ${\mathrm{\θ}}_c=\frac{{\mathrm{\θ}}_i+{\mathrm{\θ}}_j}{2},$ r₀为天线到地面目标的距离，i，j分别代表两个孔径的编号。

InSAR回波谱中，谱重叠的部分代表了SAR图像对之间的相干成分，重叠的部分越多，相干性就越高。

从(5)式中可以看出，基线的长度越长，重叠的部分就越少，因此当基线长度增大到使Δf等于整个带宽时，两孔径谱重叠的部分就消失，SAR图像对之间也就没有了相干信号成分，相干性彻底消失。这时InSAR测量就无法进行了，相应的基线值被称作极限基线B_i，j ^L，C，其中B_i，j ^L，C定义为

$B_{i,j}^{L,C}=r_0\frac{B}{f_0}\tan {\mathrm{\θ}}_c---\left(6\right)$

3)高程精度设计

在InSAR图像中像素高度和相位差的关系如下：

$\mathrm{\φ}=\frac{2\mathrm{\π}}{\mathrm{\λ}}[B_x\sin \mathrm{\θ}-B_y\cos \mathrm{\θ}]---\left(7\right)$

h＝H₀-r₀(cosθ)                                (8)

其中，λ代表雷达波长，H₀代表运动平台高度，B_x代表x方向基线长度分量，B_y代表y方向基线长度分量，θ代表雷达波入射角。

从干涉测量的角度出发，θ可根据式(7)计算出，地面高度可根据式(8)计算出。参数ρ，B_x，B_y，H₀和φ的偏差会导致计算高度的偏差。计算式(7)和式(8)对每个参量的导数，能够得到所述5个参数的偏差对地形三维高度造成的误差分别为：

${\mathrm{\σ}}_h^{\mathrm{\ρ}}={\mathrm{\σ}}_r\cos \mathrm{\θ}---\left(9\right)$

${\mathrm{\σ}}_h^{B_x}=\frac{r_0{\left(\sin \mathrm{\θ}\right)}^2}{(B_x\cos \mathrm{\θ}+B_y\sin \mathrm{\θ})}{\mathrm{\σ}}_{B_x}---\left(10\right)$

${\mathrm{\σ}}_h^{B_y}=\frac{r_0\sin \mathrm{\θ}\cos \mathrm{\θ}}{(B_x\cos \mathrm{\θ}+B_y\sin \mathrm{\θ})}{\mathrm{\σ}}_{B_y}---\left(11\right)$

${\mathrm{\σ}}_h^{H_0}={\mathrm{\σ}}_{H0}---\left(12\right)$

${\mathrm{\σ}}_h^{\mathrm{\φ}}=\frac{\mathrm{\λ\ρ}\sin \mathrm{\θ}}{2\mathrm{\π}(B_x\cos \mathrm{\θ}+B_y\sin \mathrm{\θ})}{\mathrm{\σ}}_{\mathrm{\φ}}---\left(13\right)$

其中σ_h ^ρ代表着参数ρ对地形三维高度造成的误差；

代表着参数B_x对地形三维高度造成的误差；代表着参数B_y对地形三维高度造成的误差；

代表着参数H₀对地形三维高度造成的误差；σ_h ^φ代表着参数φ对地形三维高度造成的误差。显然当各类误差源一定时，基线越长，产生的测高误差就越小。因此在基线长度不超过极限基线的情况下，如果基线较长，InSAR系统的测高精度理论上也就越高。

4)相位展开精度设计

由于本发明采用的是基于聚类分析的相位展开方法，因此其测高精度主要是取决于不同基线的长度组合。不同长度的基线组合将决定距离最近两类像素的聚类中心之间的距离，显然这个距离越大，系统对噪声的鲁棒性也就越强；相反，这个距离越小，系统对噪声的鲁棒性也就越差。

多基线测高最基本的理论依据就是利用不同基线长度下得到的同一地形的高程信息，其数学表达式如公式(4)所示，即：

如果把k₁与k₂看做变量，显然(4)式是一个关于k₁与k₂的直线方程，而且B₁与B₂的基线长度不会变化，那么最近两类像素之间的距离显然就是当k₁与k₂取不同的整数组合时距离最近的两条平行线间的距离。根据平行线间的距离公式，平行线间的距离的数学表达式为：

$\frac{|B_1{k}_2-B_2{k}_1|}{\sqrt{B_1^2+B_2^2}}---\left(14\right)$

根据欧几里德定理的推论可以知道，如果a与b的最大公约数为d，则一定存在x和y使ax+by＝d。由于B₁与B₂的最大公约数为1，因此一定能够找到一组整数(k₁，k₂)使得|B₁k₂-B₂k₁|的值为1，因此最近距离的表达式应为：

$\frac{1}{\sqrt{B_1^2+B_2^2}}---\left(15\right)$

要使式(15)的值最大，假设机载平台能够承载总长度为a米的基线，也就是说：

B₁+B₂＝a                                (16)

且B₁与B₂的最大公约数为1

在(16)式的约束下对(15)式进行优化得：

$\min \frac{1}{\sqrt{B_1^2+B_2^2}}---\left(17\right)$

B₁+B₂＝a

s.t

且B₁与B₂的最大公约数为1

通过优化将满足约束条件的B₁与B₂优化出来，再通过比较代价函数值就可以确定高程精度最优的基线长度。

参照图4，本发明的高程测量方法包括如下步骤：

步骤1，将雷达接收到的回波生成三幅二维复图像。

雷达在接收模式下，雷达波通过子孔径天线A，子孔径天线B，子孔径天线C由接收机A，接收机B，接收机C接收到数据处理器，再由数据处理器里的SAR成像单元对雷达原始回波进行SAR成像处理，得到三幅二维复图像。

步骤2，对三幅二维复图像进行距离向和方位向偏移的配准。

利用现有的互相关图像配准方法或者控制点图像配准方法，但不拘泥于这些方法，对两幅二维图像之间距离向和方位向偏移进行纠正，使得被纠正后的两幅二维复图像相同位置上的像素点是由地面同一位置的回波得到。

步骤3，分别在不同基线长度下提取每两个孔径天线与地面同一目标之间的相位差，得到两幅干涉相位图。

首先，在第一复图像与第二复图像对应的子孔径天线(A)与子孔径天线(B)之间的基线长度l_A，B下，提取子孔径天线(A)与子孔径天线(B)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(B)对应的干涉相位图；

然后，在第一复图像与第三复图像对应的子孔径(天线A)与子孔径天线(C)之间的另一基线长度l_A，C下，提取子孔径天线(A)与子孔径天线(C)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(C)对应的干涉相位图。

或者先在第一复图像与第三复图像对应的子孔径天线(A)与子孔径天线(C)之间的另一基线长度l_A，C下提取子孔径天线(A)与子孔径天线(C)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(C)对应的干涉相位图；再在第一复图像与第二复图像对应的子孔径天线(A)与子孔径天线(B)之间的基线长度l_A，B下提取子孔径天线(A)与子孔径天线(B)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(B)对应的干涉相位图。

步骤4，对两幅干涉相位图像进行滤除噪声处理。

由步骤3直接得到的干涉相位图信噪比很低，因此在这里利用现有的均值滤波方法或中值滤波方法，但并不拘泥于这些方法，来对干涉相位图进行滤除噪声处理，从而提高三维高程图的测量精度。

步骤5，对滤除噪声后的干涉相位图进行基于聚类的多基线相位展开，得到真实相位图。

首先，对滤波后的干涉相位图中的像素点进行聚类，并找出每一聚类的中心；

然后，利用中国剩余定理求出聚类中心对应的干涉相位的高程模糊数；

最后，用上述高程模糊数代表该聚类内其他像素的高程模糊数，并通过式(18)得到每一聚类像素点的展开相位：

其中k_i，j代表坐标为(i，j)位置上的像素高程模糊数，

代表坐标为(i，j)位置上的像素的干涉相位。

步骤6，用真实相位图内每一像素点对应的真实相位值，乘以地形高度比例因子，得到三维地形高程图。

首先，通过干涉合成孔径雷达的几何关系，利用基线长度，雷达脉冲入射角，雷达波长，远景中心距来按公式(19)确定地形高度比例因子：

$a=\frac{\mathrm{\λ}\×r\×\sin \left(\mathrm{\θ}\right)}{B\×4\mathrm{\π}}---\left(19\right)$

其中λ是雷达波长，r是天线相位中心到目标像素的斜距，θ是雷达波入射角，B是基线长度。

然后，再通过地形高度比例因子与相对应的展开相位的乘积来恢复出高程图。每个像素的高程比例因子是随着像素位置的不同而变化的。

本发明的效果可以通过以下仿真进一步说明：

本仿真假定的雷达运动平台为飞机，仿真的基线长度l_A，B为-0.638米；基线长度l_A，C为2.815米，仿真测量地形为一个山丘和两座楼房。

本发明测高效果的仿真干涉条纹如图5，其中，图5(a)为在基线长度为l_A，B下对仿真地形测量后的干涉条纹图，图5(b)为在基线长度为l_A，C下对仿真地形测量后的干涉条纹图。

图5中所给出的干涉条纹图对应的标准高程模糊数如图6所示。其中，图6(a)为图5(a)对应的标准高程模糊数，图6(b)为图5(b)对应的标准高程模糊数。

通过本发明方法得到的仿真地形的高程模糊数与传统方法所得到的仿真地形高程模糊数的比较结果如图7所示。其中图7(a)和图7(b)为本发明方法得到的高程模糊数，图7(c)和图7(d)为传统CRT方法得到的高程模糊数，图7(e)和图7(f)为Projection方法得到的高程模糊数。通过图7和图6比较可以看出，本发明方法比起传统方法得到的结果更加接近图6中所给出标准高程模糊数。

Claims (8)

1、一种干涉合成孔径雷达三维地形成像系统，包括：三个子孔径天线、雷达波发射机、雷达波接收机和成像数据处理器，该成像信号处理器包括SAR图像处理单元和InSAR图像处理单元，其特征在于：

三个子孔径天线垂直于飞机飞行方向放置，以实现三维高程地形的绘制；

InSAR处理单元中包括：图像配准子单元、相位滤波子单元、高程反演子单元和多基线相位展开四个子单元，该多基线相位展开子单元固化有基于聚类分析的相位展开处理程序，实现对所有像素点的逐类相位展开。

2、根据权利要求1所述的三维地形成像系统，其特征在于：三个子孔径天线的任意两个子孔径组成一组，每一组对应一个基线长度，该两个基线长度同时满足极限基线、展宽模糊高度、高程精度、相位展开精度这四方面的要求。

3、根据权利要求2所述的三维地形成像系统，其特征在于：两个基线长度不同，且两个长度之间的最大公约数为1。

4、一种干涉合成孔径雷达进行三维地形高程测绘的方法，包括如下步骤：

1)将雷达接收到的回波生成三幅二维复图像；

2)对三幅二维复图像进行距离向和方位向偏移的配准，并分别在不同基线长度下提取每两个孔径天线与地面同一目标之间的相位差，得到两幅干涉相位图；

3)对两幅干涉相位图像进行滤除噪声处理；

4)对滤除噪声后的干涉相位图进行基于聚类的多基线相位展开，得到真实相位图；

5)用真实相位图内每一像素点对应的真实相位值，乘以地形高度比例因子，得到三维地形高程图。

5、根据权利要求4所述的三维地形高程测绘的方法，其中步骤2)所述的对三幅二维复图像进行距离向和方位向偏移的配准，是先将第一复图像与第二复图像进行距离向和方位向偏移的纠正；再将第一复图像与第三复图像进行距离向和方位向偏移的纠正。

6、根据权利要求4所述的三维地形高程测绘的方法，其中步骤2)所述的分别在不同基线长度下提取每两个孔径天线与地面同一目标之间的相位差，得到两幅干涉相位图，是先在第一复图像与第二复图像对应的子孔径天线(A)与子孔径天线(B)之间的基线长度l_A，B下，提取子孔径天线(A)与子孔径天线(B)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(B)对应的干涉相位图；再在第一复图像与第三复图像对应的子孔径天线(A)与子孔径天线(C)之间的另一基线长度l_A，C下，提取子孔径天线(A)与子孔径天线(C)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(C)对应的干涉相位图。

7、根据权利要求4所述的三维地形高程测绘的方法，其中步骤2)所述的分别在不同基线长度下提取每两个孔径天线与地面同一目标之间的相位差，得到两幅干涉相位图，是先在第一复图像与第三复图像对应的子孔径天线(A)与子孔径天线(C)之间的另一基线长度l_A，C下提取子孔径天线(A)与子孔径天线(C)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(C)对应的干涉相位图；再在第一复图像与第二复图像对应的子孔径天线(A)与子孔径天线(B)之间的基线长度l_A，B下提取子孔径天线(A)与子孔径天线(B)这两个天线与地面同一目标之间的相位差，得到复图像(A)与复图像(B)对应的干涉相位图。

8、根据权利要求4所述的三维地形高程测绘的方法，其中步骤4)所述的对滤除噪声后的干涉相位图进行基于聚类的多基线相位展开，按如下步骤进行：

8a)对滤波后的干涉相位图中的像素点进行聚类，并找出每一聚类的中心；

8b)利用中国剩余定理求出聚类中心对应的干涉相位的高程模糊数；

8c)用上述高程模糊数代表该聚类内其他像素的高程模糊数，并通过如下公式得到每一聚类像素点的展开相位：

其中k_i，j代表坐标为(i，j)位置上的像素高程模糊数，

代表坐标为(i，j)位置上的像素的干涉相位。

Images