CN104007439A - 一种干涉圆迹sar高程估计处理方法 - Google Patents

Abstract

本发明公开了一种干涉圆迹SAR高程估计方法，包括：步骤1，粗估计出场景DEM；步骤2，在粗估计出的场景DEM的基础上，对InSAR平台主、副天线分别独立进行BP成像，其中采用分子孔径的方式对主副天线获取的圆迹SAR数据进行BP成像，每n度为一个子孔径，得到360/n个子孔径图像对，n是大于等于1的自然数；步骤3，对于每个子孔径图像对，进行干涉处理；步骤4，将步骤3获得的各子孔径图像对干涉处理估计的场景高程，根据其干涉系数的大小，进行加权平均，得到360度观测下的高程图；步骤5，将步骤4估计出的高程图作为步骤2中的参考DEM输入，重复步骤2到4若干次，直到干涉图像的干涉系数不再下降。

Description

一种干涉圆迹SAR高程估计处理方法

技术领域

本发明涉及高程估计技术领域，具体涉及一种干涉圆迹SAR高程估计方法。

背景技术

合成孔径雷达(SAR)具有全天时全天候，高分辨大区域成像的能力，已经成为军用和民用的重要工具。考虑到传统直线SAR模式成像几何下存在由于迭掩、透视缩短和阴影等SAR影像固有变形现象以及仅能获取观测对象在某小范围方向上的后向散射的问题，为了进一步提高SAR的观测能力，90年代初，美国学者Knaell等最早提出了圆迹合成孔径雷达(Circular SAR，CSAR)概念，此后CSAR因其特有的优势得到了广泛的研究。CSAR模式观测几何下，圆迹SAR通过对观测场景的360°全方位观测，能够获取目标在各方向的散射特征，同时其理论分辨率可以达到亚波长量级，以及具有3维成像能力，大地拓展了现有SAR成像系统的性能和应用领域。

如何利用圆迹SAR数据对场景进行精确重建，是圆迹SAR技术得以广泛应用的前提条件。有别于传统直线SAR成像，圆迹SAR成像时需要用到精确的场景DEM数据。这是由于，圆迹SAR的多角度成像几何通常采用正射成像的方法得到不同360°全方位观测的结果，若成像时所参考DEM数据与实际不符，回波在图像上的体现为朝观测方向发生偏移，尤其是对于各向同性的点目标而言，当参考高程不准确时，其在图像上将形成一个半径与参考高程的误差有关的圆环。因此，圆迹SAR成像时所采用的参考DEM数据的精确程度影响着最终的圆迹SAR图像质量。现有机载圆迹SAR实验中，成像时所参考的场景DEM主要通过激光雷达(Lidar)获取。由于激光雷达成本高、激光穿透性与SAR使用的微波波段不同、平台飞行高度有限等方面的局限性，仅仅依赖lidar获取场景DEM将会影响圆迹SAR技术应用的进一步推广(尤其是一些需要实时获取场景图像的军事应用)。若能够直接利用圆迹SAR数据本身进行场景高程信息提取，则可以有效实现成像处理的一体化，做到实时成像。在通过CSAR数据获取场景高程信息的研究方面，国内外仅对体现圆迹SAR三维分辨能力的立体测量方法进行了相关研究，2007年，法国宇航局利用SETHI机载SAR系统采集的X波段机载圆迹SAR数据获取了法国城区的全方位观测DEM。该试验中DEM的获取是通过立体测量方法实现，精度达2m。该试验主要是为了验证圆迹SAR的三维分辨能力，从而体现360度观测在城市测绘中的潜力，实际上获取的高程精度仍不能满足高分辨率CSAR精确成像需求。

实际上，通过SAR图像获取场景的DEM信息也是SAR的重要应用之一，其获取方法可以分为两种，即“合成孔径雷达干涉测量(InterferometricSAR,InSAR)”和“合成孔径雷达立体测量(Radargrammetry)”。SAR干涉测量利用雷达图像的相位信息而SAR立体测量利用其强度(灰度)信息获取场景DEM数据，因而，前者获取的DEM精度远高于后者。由于圆迹SAR技术的研究目前还处于实验初级阶段，国内外关于干涉圆迹SAR技术仅仅局限于小目标、小场景三维成像方面的研究，鲜有关于干涉圆迹SAR技术获取场景DEM的研究以及实验。圆迹SAR由于其全方位观测能力，可获得多角度乃至全方位的目标散射特性分布及其最大信噪比SAR回波信号，若利用圆迹SAR的不同子孔径观测目标区域，则得到不同视角方向的场景，其中的某个或某几个子孔径具有最大信噪比SAR回波信号。对于干涉SAR而言，信噪比越高，相干系数越高，从而干涉相位越精确。因此，通过干涉圆迹SAR观测，可以获取到目标的最高相位精度，相比于传统干涉SAR，圆迹干涉SAR方法获取的DEM将更为精确。因此，有必要针对CSAR数据，对理论测高精度可以达到波长量级的干涉SAR方法进行研究。现有关于圆迹干涉SAR的研究主要有小型目标三维成像以及地球表面地形形变估计，根据目前能够掌握的研究现状来看，基于圆迹干涉SAR的场景DEM估计还未开展深入的研究。

发明内容

为了克服现有技术的上述缺陷，本发明提出了一种干涉圆迹SAR高程估计方法。

根据本发明提出的干涉圆迹SAR高程估计方法，包括步骤：步骤1，粗估计出场景DEM；步骤2，在粗估计出的场景DEM的基础上，对InSAR平台主、副天线分别独立进行BP成像，其中采用分子孔径的方式对主副天线获取的圆迹SAR数据进行BP成像，每n度为一个子孔径，得到360/n个子孔径图像对，n是大于等于1的自然数；步骤3，对于每个子孔径图像对，进行干涉处理；步骤4，将步骤3获得的各子孔径图像对干涉处理估计的场景高程，根据其干涉系数的大小，进行加权平均，得到360度观测下的高程图；步骤5，将步骤4估计出的高程图作为步骤2中的参考DEM输入，重复步骤2到4若干次，直到干涉图像的干涉系数不再下降。

相比激光雷达测量成本高、平台飞行高度有限以及立体测量精度低的缺点，本发明的方法可以为圆迹SAR成像提供高精度的DEM，从而有效地实现圆迹SAR成像处理的一体化，做到实时的高精度高分辨率圆迹SAR成像。

附图说明

图1是干涉圆迹SAR几何模型。

图2是干涉圆迹SAR局部几何关系图。

图3是干涉圆迹SAR高程测量处理流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

本发明提出了一种干涉圆迹SAR高程估计方法。该方法包括的具体步骤如下：

步骤1，粗估计出场景DEM。在该步骤，可采用立体像对测量的方法粗估计出场景DEM。

圆迹SAR具有三维成像能力，在进行二维成像时，位于成像参考平面外的目标在不同方位角度观测时，将投影到成像平面的不同位置，因此可以通过立体像对的方法对场景的高程进行测量(具体参见Stephan Palm,Hélène M.Oriot,and Hubert M.Cantalloube.Radargrammetric DEMExtraction Over Urban Area Using Circular SAR Imagery.IEEETRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING,2012,VOL.50,NO.11,pp:4720-4725.)。立体测量方法的高程估计精度与高分辨率圆迹SAR精确成像的所需的高程精度相比略显粗糙(具体参见StephanPalm,Hélène M.Oriot,and Hubert M.Cantalloube.Radargrammetric DEMExtraction Over Urban Area Using Circular SAR Imagery.IEEETRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING,2012,VOL.50,NO.11:4720-4725.)，但该方法粗估计出的高程可以作为场景初始估计高度，从而可以避免由于成像参考高程与实际高程差距大导致的干涉图像对像素偏移大于分辨单元以及成像结果散焦现象。

利用圆迹SAR数据进行立体测量法进行DEM估计的过程如下：将360度干涉SAR(Interferometric SAR，INSAR)平台轨迹圆环划分为多段圆弧，对每段圆弧分割为多个子孔径，利用子孔径图像间的相似性，通过相关系数最大法(参见Tupin F,Nicolas J M.Matching criteria forradargrammetric Remote and Geoscience International IEEE[C]2002Sensing Symposium(IGARSS).2002,5:2608-2610.)，在其他子孔径图像上寻找与中心图像匹配的像素点，得到二者的位置偏移，进而提取地面场景的DEM。

步骤2，在粗估计出的DEM的基础上，对InSAR平台主、副天线分别独立进行BP成像。

干涉圆迹SAR的几何模型如图1所示，两个雷达天线间的基线长度为b，与水平面的夹角为β。天线1为主天线，天线2为副天线，基线方向始终与雷达平台的速度方向保持垂直，即本发明采用交轨干涉模式。雷达平台在距离地面的高度平面作圆周运动，形成两条独立的圆形孔径。方位角θ∈[0,2π]，天线1的轨迹半径为R，高度为H。在本发明，假设两天线采用乒乓收发模式，即轮流自发自收，每个天线各自发射Chirp信号，中心频率为f_C，带宽为B_r。设场景中有一个任意点目标P，为了便于分析假设散射系数在方位角[θ_pc-Δθ/2,θ_pc+Δθ/2]范围内保持一致(其中θ_pc为目标零多普勒位置对应的方位角，Δθ为合成孔径的方位角度)，在其他方向，散射系数为0。对于机载单航过圆迹干涉SAR而言，由于基线长度有限，两幅天线视角十分接近，可以假设两幅天线视角下目标P的散射系数相同。

基于时域处理的后向投影(Back Projection,BP)成像算法具有逐点成像的特性，适用于进行圆迹SAR成像(Soumekh,M.Synthetic Aperture RadarSignal Processing with MATLAB Algorithms.Wiley.1999)。BP成像算法首先计算每个散射点到天线位置的距离，然后对距离压缩后的回波数据进行插值，得到各个脉冲对某个散射点所作出的贡献，最后通过相干累加对该散射点聚焦进而实现高分辨率成像。概括地说，BP成像算法是根据天线位置到成像区域中某一像素点的相对位置，取出对应回波信号中的相应时延处的能量，沿孔径方向叠加起来，使来自该像素点的回波信号为同相，从而信号得到加强，而来自其它点的回波信号则由于相位不同，叠加结果趋于零，因此可认为最终叠加结果即为该像素点的聚焦结果。对整个成像区域，按照上述过程逐像素点地进行相干叠加处理，即可获得最终图像。下面以对天线1获取的圆迹数据成像为例，简要叙述BP成像的主要实现步骤(Soumekh,M.Synthetic Aperture Radar Signal Processing with MATLABAlgorithms.Wiley.1999.)。

首先，对SAR原始回波信号进行距离向匹配滤波处理，如下所示，

$S_M(t,\mathrm{\θ})=s(t,\mathrm{\θ})\⊗p^{*}(-t)$

其中：表示时域卷积，*表示复数共轭，p(t)表示雷达反射信号，s(t,θ)＝p(t-t_θ)表示SAR原始回波信号，t_θ表示回波延时，θ表示天线相位中心方位角，θ∈(θ_pc-Δθ/2,θ_pc+Δθ/2)。

然后，将成像区域网格化，每个像素所在合成孔径的数据相干累加实现方位向聚焦。对于目标P而言，其成像结果可以表示为

$f(x_P,y_P,z_P)=\underset{\mathrm{\θ}\∈{\mathrm{\θ}}_{\mathrm{pc}}-\mathrm{\Δ\θ}/2,{\mathrm{\θ}}_{\mathrm{pc}}+\mathrm{\Δ\θ}/2}{\∫}{S}_M(t_{x_P,y_P,z_P}\left(\mathrm{\θ}\right),\mathrm{\θ})\mathrm{d\θ}$

其中，(x_P,y_P,z_P)表示目标P的距离、方位和高程坐标位置，表示天线相位中心到目标的回波时延， $R_{x_P,y_P,z_P}\left(\mathrm{\θ}\right)=\sqrt{{(R\cos \left(\mathrm{\θ}\right)-x_P)}^2+{(R\sin \left(\mathrm{\θ}\right)-y_P)}^2+{(H-z_P)}^2}$ 表示雷达到目标的距离。

上式中的S_M(t_x,y,z(θ),θ)表示带有载频的原始回波信号，而在实际的SAR系统中通常获取的是去载频的基带信号。因此，在BP成像处理中，需要先补偿载频再进行方位向相干累加，这样才能实现在目标点上同相相干累加，而在非目标点上非相干累加。针对实际应用情况，修改上述BP成像表达式如下

其中：表示补偿载频的相位，λ表示载频对应的波长。

由于场景的散射特性未知，多是各向异性的散射体，因此本发明采用分子孔径的方式对主副天线获取的圆迹SAR数据进行成像，每n度为一个子孔径，因此得到360/n个图像对。子孔径的合成孔径角度n由场景的目标特性来确定。

步骤3，对于每个子孔径图像对，进行干涉处理。该步骤进一步包括：步骤31，图像干涉；步骤32，干涉图像中相干系数估计；步骤33，干涉图像滤波去噪；步骤34，基线估计；步骤35，反演参考地形与原始地形的相对高度信息，然后进行地面位置校正。

在步骤31，将步骤2得到的子孔径复图像对进行共轭相乘，得到干涉图像。

下面对干涉图像的信号进行分析建模，即建立干涉相位与高程信息之间的关系。

干涉圆迹SAR的几何模型如图1所示，两个雷达天线间的基线长度为b，与水平面的夹角为β。天线1为主天线，天线2为副天线，基线方向始终与雷达平台的速度方向保持垂直，即这里采用交轨干涉模式。雷达平台在距离地面的高度平面作圆周运动，形成两条独立的圆形孔径。方位角θ∈[0,2π]，轨迹半径为R。在本发明，假设两天线采用乒乓收发模式，即轮流自发自收，每个天线各自发射Chirp信号，中心频率为f_C，带宽为B_r。设场景中有一个任意点目标P，为了便于分析假设散射系数在方位角[θ_pc-Δθ/2,θ_pc+Δθ/2]范围内保持一致(其中θ_pc为目标零多普勒位置对应的方位角，Δθ为合成孔径的方位角度)，在其他方向，散射系数为0。对于机载单航过圆迹干涉SAR而言，由于基线长度有限，两幅天线视角十分接近，可以假设两幅天线视角下目标P的散射系数相同。

为了使表述得更清楚，本发明建立了一个新的坐标系，如图2所示，其中，z′轴等同于z轴，y′表示零多普勒位置载机飞行方向，x′轴由右手坐标系确定，位置A、B对应天线1和天线2的零多普勒位置，ψ为天线关于目标的下视角(对于单航过干涉SAR而言，由于基线远远小于斜距，因此主副天线关于目标的下视角可以看作相等)。在三维空间中，对于天线1，目标P被重建为直线段PA′(实际上是弧线段，远场下可认为是直线段)，对于天线2，被重建为PB′。由于两干涉天线观测视角存在差异，两条直线段只有在目标P的真实高度处相交，而在其他高度平面，其真实高度与给定高度之间的高度差将引起切片对的图像偏移。

从图3中可以看出，在对于目标P进行BP成像时，在多普勒中心位置处，天线1、天线2精确的补偿相位分别为

干涉测量的目的是为了获取地形高程，因此在BP成像时地形高程未知，目标高度采参考高度。因此，补偿相位为

由于成像区域网格化，离散采样点中不一定存在位置，为了方便分析，假设离散采样点恰好位于A′，因此，实际的补偿相位为

基于前面的假设，对于天线1而言，A′处的相位为目标P散射系数的相位，而对于天线2而言，A′处的存在偏置相位

因此在A′处的干涉相位可以表示为

公式(5)体现了干涉相位与目标高程之间的关系。

步骤32，去除干涉图像中干涉系数低于0.5的区域。

采用干涉进行高程测量的前提是主从影像是相干的，而在InSAR处理过程中存在许多去相干因素，包括时间去相干、空间去相干、热噪声去相干、数据处理去相干和大气效应引起的去相干等。这些去相干因素的综合作用，造成InSAR相干退化，甚至完全不相干，其也最终决定了干涉SAR测高或测形变成像的质量。对于机载单航过圆迹干涉SAR而言，热噪声去相干是去相干因素的一个重要方面。热噪声去相干主要取决于目标的后向散射强度。当目标的后向散射强度很弱时，热噪声去相干因素十分明显。由于圆迹SAR具有全方位观测的优势，可获得多角度乃至全方位的目标散射特性分布。因此，当某个方位观测角度下出现由于目标散射强度低而导致低相干区域时，我们在该方位观测角度下的高程恢复时不采用该目标区域的干涉相位信息，采用其它高相干区域的干涉相位信息，从而可以避免低相干相位信息对高程估计精度的影响。

干涉系数的计算如下：

对于复干涉图像对s₁(由天线1获取)和s₂(由天线2获取)，其相干系数定义为：

$\mathrm{\γ}=\frac{\left|E\right[s_1\·{s_2}^{*}\left]\right|}{\sqrt{E\left[\right|{s_1|}^2]}E\left[{\left|s_2\right|}^2\right]}$

其中γ＝E[]表示统计平均操作。在实际估算干涉相干的过程中，我们采用最大似然(MLS)相干估算器[R Touzi,A Lopes,J Brumquet,and P WVachon.1999Coherence estimation for SAR imaginery.IEEE Trans onGeosci and Remote Sens,37(1):135-149]

假设在估算窗口内，散射体是各态历经的，因此采用空间平均代替总体平均，同时考虑到图像间的相位差异，则MLS估算表达式为：

其中是每个像元的系统相位分量。

步骤33，干涉图像滤波去噪。

原始干涉图中存在着大量的相位噪声,这会导致高程估计精度的下降,为了降低相位噪声影响,通常需对干涉图进行噪声抑制，常用干涉图滤波方法有均值滤波[[I.Lim,T.S.Yeo,C.S.Ng,et al,Phase NoiseFilter for Interferometric SAR,1997InternationalGeoscience and Remote Sensing Symposium Proceedings(IGRASS.97),pp:445-447,1997.]、中值滤波[J.S.Lee,Intensity andphase statistics of multi-look polarimetric and Interferometricimagery,IEEE Transaction on Geoscience and Remote Sensing,Vol.32,No.5,pp:1017-1028,1994.]、Goldstein滤波[R.M.Goldstein and C.L.Werner,Radar interferogram filtering for geophysicalapplication Geophysical Research Letters,Vol.25,No.21,pp:4035-4038,1998]等。

由于本发明发明采用在粗估计的DEM上进行成像，参考成像高程较实际场景高程而言，误差不大，干涉图像的相位不会发生缠绕现象。在进行滤波去噪时，无需考虑滤波方法对干涉图条纹边缘的破坏。本发明采用简单快速均值滤波方法，在滤波窗口(2M+1)×(2N+1)的均值滤波器为：

其中：表示干涉图在像素坐标(i,j)的相位值,表示滤波后像素坐标(m,n)的相位值。滤波窗口内尽量选择独立同分布的样本，即尽量选择干涉图中幅度接近或是相干系数接近的区域样本。

步骤34，在步骤2获取的复图像对的基础上进行基线估计。

最优的垂直轨迹基线长度是获得高精度地形高度的基础,而准确的基线参数则是非模糊相位图向高度图转换的必要条件。现有的基线估计方法主要有两种:一种是利用星历参数,即卫星的位置信息,根据干涉SAR的空间几何关系计算基线；另一种则是从SAR干涉图像出发,利用SAR干涉机理来估计基线。第一种方法实现简单,但是它的精度受平台位置测量精度和测量采样频率的限制。第二种方法在实现中需要地面参考点信息,它利用回波信号本身来估计基线,因此精度较高。

考虑到实际因素，本发明采用第二种估计基线的方法，具体方法采用Mohr和Madsen[Mohr J J,Madsen S N.Parametric Estimation of TimeVarying Baselines in Airborne Interferometric SAR[C].IGARSSp1996.Piscataway,NJ,USA:IEEE,1996:677-679.]给出了一种依据图像对间偏移参数的基线估计方法(在步骤2获取的复图像对的基础上进行基线估计)，这里不再详细叙述。

步骤35，反演参考地形与原始地形的相对高度信息，然后进行地面位置校正。

对于，反演参考场景DEM与真实场景DEM的相对高度信息，根据InSAR高程测量的机理,可以求出步骤33获得的干涉相位图像中每个像素点所对应的地面样点的高度值，这一过程就是干涉相位值到地面高度值的转换。需要首先计算成像参考地形到雷达的斜距信息以及参考DEM到雷达的下视角，然后根据步骤1中分析得到的干涉相位与目标高程之间的关系式，反演出参考地形与原始地形的相对高度信息。

对于地面位置校正，由于干涉图像上的像素点对应的高度不同,使得间隔分布的样点,在地面上呈不等间隔分布，因此,地面位置校正的目的,就是恢复地形在地面上的等间隔采样。地面位置校正的处理步骤如下：

1)结合各像素点的高度值,求出它们在地面上的具体位置(求解过程参考文献[付涛.InSAR高保相成像及关键技术研究.[D].成都:电子科技大学,2013.来源cnki全文数据库]),生成准确的地貌曲面。

2)根据前一步生成的曲面,计算等间距差值,进而得到等间距采样的参考地形与原始地形的相对高度信息。

步骤4，将步骤3获得的各子孔径干涉处理估计的场景高程，根据其干涉系数的大小，进行加权平均，得到360度观测下的高程图。

参考下式：

$\stackrel{\‾}{\mathrm{\Δh}(m,n)}=\frac{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\γ}}_1(m,n)\mathrm{\Δ}{h}_l(m,n)}{\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\γ}}_l(m,n)}$

其中：表示全方位观测下的像素坐标(m,n)处参考地形与原始地形的相对高度的加权平均，Δh_l(m,n)表示第l个子孔径下像素坐标(m,n)处参考地形与原始地形的相对高度，γ_l(m,n)表示第l个子孔径下像素坐标(m,n)处的相干系数，L表示子孔径的个数。

最后将全方位观测得到的参考地形与原始地形的相对高度加上参考地形高程，最终得到场景地形高程。

步骤5，将步骤4估计出的场景地形高程作为步骤2中的参考DEM输入，重复步骤二到四若干次，直到干涉图像的干涉系数(或图像对比度)不再下降。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种干涉圆迹SAR高程估计方法，该方法包括：

步骤1，粗估计出场景DEM，作为首次干涉处理时成像采用的参考场景DEM；

步骤2，在参考场景DEM的基础上，对InSAR平台主、副天线分别独立进行BP成像，其中采用分子孔径的方式对主副天线获取的圆迹SAR数据进行BP成像，每n度为一个子孔径，得到360/n个子孔径复图像对，n是大于等于1的自然数；

步骤3，对于每个子孔径复图像对，进行干涉处理；

步骤4，将步骤3获得的各子孔径复图像对干涉处理估计的场景DEM，根据其干涉系数的大小，进行加权平均，得到360度观测下的场景DEM；

步骤5，将步骤4估计出的场景DEM，作为步骤2中的参考场景DEM，重复步骤2到4若干次，直到干涉图像的干涉系数不再下降。

2.根据权利要求1所述的方法，其特征在于，步骤1采用立体像对测量法粗估计出场景DEM。

3.根据权利要求2所述的方法，其特征在于，步骤1进一步包括：将360度雷达平台轨迹圆环划分为多段圆弧，对每段圆弧分割为多个子孔径，利用子孔径图像间的相似性，通过相关系数最大法，在其他子孔径图像上寻找与中心图像匹配的像素点，得到二者的位置偏移，进而提取地面场景的DEM。

4.根据权利要求3所述的方法，其特征在于，步骤3进一步包括：

步骤31，将步骤2得到的子孔径复图像对进行共轭相乘，得到干涉图像；

步骤32，去除干涉图像中干涉系数低的区域；

步骤33，对经步骤32处理后的干涉图像进行滤波去噪；

步骤34，在步骤2获取的复图像对的基础上进行基线估计；

步骤35，基于所述基线估计，反演参考地形与真实地形的相对高度信息，然后进行地面位置校正。