CN111896955B - 一种船载sar交轨干涉处理方法 - Google Patents

Abstract

本发明公开一种船载SAR交轨干涉处理方法，包括以下步骤：前置滤波、块影像粗配准、块影像精配准、影像全局精配准、干涉图滤波、相位解缠、基于图像同名点的基线估计和基于控制点的基线估计；本发明包括主辅影像滤波、配准、干涉计算和干涉图的滤波、相位解缠等工作，针对船载双天线SAR的成像特点，图像配准采用了分块配准方法，利用块影像粗配准和块影像精配准，提高图像的相干系数，达到了较高的相干性；通过干涉图滤波确保了相干信息和滤波的结果在同一位置具有同样的大小，避免了滤波结果被重叠部分的相干值影响；通过基于图像同名点的基线估计方法和基于控制点的基线估计方法，实现毫米甚至亚毫米级的基线估计精度。

Description

一种船载SAR交轨干涉处理方法

技术领域

本发明涉及交轨干涉处理技术领域，尤其涉及一种船载SAR交轨干涉处理方法。

背景技术

根据InSAR平台和使用条件的不同，获取InSAR干涉数据的干涉模式主要有三种：交轨干涉测量、顺轨于涉测量和重复轨道干涉测量，顺轨干涉测量两天线所构成的直线方向与飞行方向平行，干涉相位差由地面目标的多普勒特性的差别引起，利用干涉测量技术就可以得到二维地面上的运动目标，实现动目标的监测，重复轨道干涉测量为同一幅天线不同时刻获得SAR数据进行干涉处理，利用该模式进行地形测绘或地表形变监测时，要求较精确的轨道信息和重复观测，这种模式最为成熟，应用也最为广泛，一般用于星载SAR平台，监测范围大，时间跨度较长，交轨干涉模式天线所构成的直线方向与飞行方向垂直，时间基线为零，能保证较强的相干性，干涉相位差由主辅天线与地面目标物之间的路径差造成，如果能够获取干涉系统的几何参数，就可以将相位信息转换成高程信息，从而得到地面点的三维坐标，建立相应地区的三维地形图，因此，该模式常用于地形制图和地形变化监测，既可运用于星载平台，也可用于机载平台以及位于更低成像高度的平台，如船载、船载以及无人机载等；

船载InSAR一般采用交轨干涉测量的形式，因为船载平台对应的成像视角变化一般较大，然而，现有技术中，难以保证主辅影像的强时间相干性，相干系数低，滤波结果容易重叠，基线估计精度低，影响后续计算精度，因此，本发明提出一种船载SAR交轨干涉处理方法以解决现有技术中存在的问题。

发明内容

针对上述问题，本发明的目的在于提出一种船载SAR交轨干涉处理方法，该船载SAR交轨干涉处理方法包括主辅影像滤波、配准、干涉计算和干涉图的滤波、相位解缠等工作，针对船载双天线SAR的成像特点，图像配准采用了分块配准方法，利用块影像粗配准和块影像精配准，提高图像的相干系数，达到了较高的相干性。

为实现本发明的目的，本发明通过以下技术方案实现：一种船载SAR交轨干涉处理方法，包括以下步骤：

步骤一：前置滤波

从频谱特性考虑干涉图的前置滤波，即滤波的操作针对影像的频谱进行，由于主辅影像方位向和距离向的频谱不完全重叠，在精配准之前对主辅影像进行距离向滤波和方位向滤波，降低频谱偏移去相干，由于成像时船载平台的位置和地面目标点之间存在着相对运动导致接收到的频率经历了多普勒频移，多普勒频率表达式为：

其中，Λ_s为传感器斜角，θ为雷达视角，v为雷达平台对地面速度，λ为波长，干涉处理中主辅图像对同一地面分辨率单元成像时，由于θ不同导致多普勒中心频率DCF_M与DCF_S不同，主辅影像的多普勒中心频率存在较大的差异，频谱偏移为：

Δf_Λ＝f_DCFm-f_DCFs (2)

频谱偏移导致主辅影像的去相干，利用B_w表示方位向带宽，利用f₁,f₂分别表示主辅图像考虑多普勒效应后的中心频率，利用B_new表示滤波后新的方位向带宽，则B_new满足:

B_new＝B_w-|f₁-f₂|＝B_w-|f_DC-f_D'_C| (3)

在式(3)中，当|f_DC-f_DC'|>B_W，则两幅图像方位向频谱没有重叠部分，两幅图像毫不相干，无法生成干涉图，由此得到两幅图像的多普勒质心对生成的干涉图的质量至关重要，采用矩形滤波器、Hamming滤波器进行方位向滤波；

矩形滤波器：

Hamming滤波器：

其中：α∈(0,1)，实际运算时常取α＝0.75；

步骤二：块影像粗配准

对主辅影像A和B进行分块，得到块影像集合{A₁,A₂,…,A_m}和{B₁,B₂,…,B_m}，块的大小根据成像地区复杂性进行调整，基于块影像间的精准偏移量{Δ₁,Δ₂,…,Δ_m}得到全图整体偏移量，先基于块影像间的精准偏移量，得到全图整，偏移量，先，于轨道参数获得初始偏移量，配准初始偏移量为主图像块与辅图像块间得到的初始偏移量，将主辅影像块的中心点对应计算得到的偏移量作为配准初始偏移量即：

p_sla(m,n)＝p_mas(m,n)+offset(m,n) (6)

式中p_mas(m,n)表示主影像块中心点，p_sla(m,n)为辅影像块中心点，offset(m,n)表示两者偏移量，由于船载交轨干涉SAR得到的主辅影像及其相似，影像偏移量与基线相关，配准初始偏移量计算精度在5个像元以内，利用成像几何模型，首先计算主影像块中心点对应的地面坐标，然后求出地面点在辅图像上的对应点的坐标，从而解得主辅图像块的配准初始偏移量，根据成像模型计算主图像块的中心点对应的地面位置矢量p_sla(m,n)，即：

式中，f_D为多普勒频率，λ为波长，R为斜距，然后采用牛顿迭代法，求解出p(X_t,Y_t,Z_t)在辅图像上的对应点p_sla(m,n)所对应的方位向时间t，同样，初始值t₀用辅图像中心点所对应的时间表示，在得到点p_sla(m,n)所对应的方位向时间t后，得到辅图像上的行m和卫星的位置，再计算该点对应的斜距R，从而解求出距离向时间，最后得到辅图像上的列n，将同一点在两幅图像上的坐标相减即得到两幅图像的偏移；

然后基于匹配测度进行块影像的像元级配准，在影像的空间域、频率域进行，根据粗匹配结果在辅图像上选取搜索窗口，然后根据配准评价指标计算主影像块和窗口的相似程度，通过移动搜索窗得到最佳配准量，达到像元级的配准精度，每个主影像块计算得到辅影像对应的像元级偏移量，对于每一个待匹配的主影像块，在搜索范围内滑动窗口取得辅影像块，计算主辅影像块的配准评价参数，选择相关系数作为配准的指标，相关系数最大时得到像元匹配量：

其中：corr为相关系数；g_i，j；g_i+r，j+c为主辅影像对应点处的振幅强度；m,n为匹配窗的大小，r,c为粗匹配的偏移量；

步骤三：块影像精配准

块影像精配准是得到辅影像对应主块影像的子像素级偏移量的过程，先基于原始图像过采样数据进一步寻找精确的配准位置，称之像元过采样匹配法，将主辅影像本身进行过采样处理，插值的间隔决定了过采样的程度，对图像进行插值采样将图像放大，之后采用复影像的配准方法逐点计算匹配测度，寻找可靠的相对偏移量估算值，在这一步中选择小窗口，但窗口的数目增多，以进行相应的数据拟合，由于窗口大小变小，增大搜索窗的大小多次计算，来防止出现最大相干偏差，当相对偏移量保持稳定，表明这个估算值可靠，否则开始新一轮搜索；

步骤四：影像全局精配准

根据步骤三得到辅影像对应每块主影像子像元级配准偏移量，通过插值，到整个主辅，像的配准偏移量，差值通过一维三次样条函数插值条件实现：

其中x_i，Δ_i分别为块中心点像元坐标和配准偏移量，S(x)为得到的插值函数，偏移量是二维的，包含距离向和方位向偏移量，将二维配准偏移量的插值转化为二次一维三次样条插值即可；

步骤五：干涉图滤波

由干涉图的特点可知，由于受到各种因素的影响，干涉图中不可避免的含有大量的噪声，干涉图滤波就是去除和抑制干涉图噪声的过程，因此利用Goldstein和Werner针对干涉条纹特性提出的一种自适应滤波方法，表达式如下：

H(u,v)＝S{|Z(u,v)|}^αZ(u,v) (10)

其中，H(u,v)是滤波响应；S{}是平滑项；u和v是空间频率，α是滤波参数，取值范围为[0，1]，α值越大滤波效果越显著，

然后在Goldstein滤波基础上引入相干系数代替α作为滤波参数来控制滤波程度的强弱，单视情况下(L＝1)干涉相位标准差表达式：

其中

从式(11)看出干涉相位方差是相干系数γ的函数，利用γ作为滤波参数，滤波器函数式转化为：

H(u,v)＝S{|Z(u,v)|}^1-γZ(u,v) (12)；

步骤六：相位解缠

干涉测量的核心是获取与距离成比例的相位信息，然而干涉图包含的相位信息是经2π周期缠绕后的相位，即干涉相位值在(-π,π]之间，要获取正确的高程信息，必须进行相位解缠：

记φ_i,j为周期缠绕前的真实相位值，

为相应的缠绕相位，二维缠绕相位与解缠相位的关系为：

式中：k_i,j为第i行第j列个像元到卫星的整周模糊度，理想情况下干涉相位图中缠绕相位

在x方向和y方向差分与解缠相位φ_i,j在这两个方向的差分完全吻合，即

式中，

分别表示x,y方向的观测相位差，

为基于先验知识选取使Δ_i,j∈(-π,π](l＝1,2)成立的数组，由于积分路径不同，Δ_i,j不能和相邻点的差分保持一致，因而定义以下差分的残差：

和

用如下的最小化问题来估算：

其中，N表示图像的N个点，M表示剔除误匹配后的点数目，M-1表示最后一个点，M-2表示相邻的点，

和

是加权函数，同时满足下列约束条件：

以上所述问题为非线性最小化问题，转化为线性最小化问题提高计算效率：

那么上述问题就转化为求解

使得：

同时满足下列约束条件：

其中

均为大于等于0的实数，计算

和

后再计算相位梯度，最终通过下式计算出解缠后的相位得到最终的解缠结果。

步骤七：基于图像同名点的基线估计

基线是干涉测量中一个非常重要的参数，它决定了有效干涉相对的选择，在知道天线和地面点精确位置情况下容易计算基线和基线定向参数，当不知道天线位置就必须考虑采用特定的方法进行基线估计，先利用基于图像同名点的基线估计方法，通过得到足够数量的主辅影像同名点进行基线估计，图像同名点通过点匹配算法得到，并使用RANSAC算法以及相干系数图进一步消除误匹配点：

式中，k_i表示sift同名点对，i＝1,2…,N，Cor_i为该点对干涉处理得到的相干系数，P为RANSAC算法收敛后的满足模型的局内点集合，ε为相干性阈值，一般取0.9，剔除误匹配点后的同名点结果记为

i＝1,2…,M；最终得到的精确匹配点对可用于船载SAR基线的估计：

式中，t_i、t_i′分别为同名点对

在主辅影像的方位向时间，R_s,i、R′_s,i分别为同名点对

对应的主辅影天线相位中心的位置矢量,R₀、R′₀分别为主辅影像相位中心的初始位置矢量，B为基线长度估计结果，α为B取得最小值式对应的基线夹角；

步骤八：基于控制点的基线估计

然后利用基于控制点的基线估计方法，基线由主辅天线相位中心组成，与天线的物理中心不重合，此时根据R-D定位模型还原SAR成像几何得到精确的基线长度，R-D模型如下：

式中，R_S，V_S，R_G，V_G，f_D分别表示雷达成像中心的位置和速度矢量、地面点的位置和速度矢量和多普勒频率，

经运动补偿后，最终平台成像状态归化为匀速直线状态，同时平台的姿态保持不变，成像时的多普勒频率对于视场中的点为常量(f_D不变)，R_S与方位向成像时间t有关：

R_S＝V_S*t+R₀ (25)

R₀为起始雷达成像位置矢量，因此结合式(20)和(21)得到只与t有关的方程，通过给定初值迭代解出方位向时间t和控制点对应的天线相位中心的位置矢量，得到控制点对应的主辅天线相位中心的位置矢量R_S和R′_S后，基线估计结果为：

式中，B为基线长度估计结果，H_S、H_S分别为主辅天线相位中心的高程，α为基线倾角估计结果。

进一步改进在于：所述步骤一中，由于船载交轨干涉测量的基线较小，基线距导致的主辅影像地距频谱范围可忽略不计，只需考虑方位向滤波。

进一步改进在于：所述步骤二中，在搜索范围内滑动窗口取得辅影像块，搜索范围为10*10个像元。

进一步改进在于：所述步骤三中，过采样图像的实现有空间域和频率域法两种，频域法运用傅里叶变换将局部图像转换到频域，共轭相乘后再进行逆FFT变换得到相关函数，省去了循环搜索，提高了计算速度。

进一步改进在于：所述步骤五中，利用γ作为滤波参数在相干性好、噪声影响小的地方滤波强度弱，保留原始相位信息，对噪声影响大的地方滤波强度会变强。

进一步改进在于：所述步骤五中，式(12)确保了相干信息和滤波的结果在同一位置具有同样的大小，而且避免了滤波结果被重叠部分的相干值影响，防止goldstein滤波器在高相干区域进行过度滤波，同时保证在低相干地区实施较强滤波。

本发明的有益效果为：本发明包括主辅影像滤波、配准、干涉计算和干涉图的滤波、相位解缠等工作，针对船载双天线SAR的成像特点，图像配准采用了分块配准方法，利用块影像粗配准和块影像精配准，提高图像的相干系数，达到了较高的相干性；通过干涉图滤波确保了相干信息和滤波的结果在同一位置具有同样的大小，避免了滤波结果被重叠部分的相干值影响；通过相位解缠获取正确的高程信息；通过基于图像同名点的基线估计方法和基于控制点的基线估计方法，实现毫米甚至亚毫米级的基线估计精度，提高后续的DEM求解精度。

附图说明

图1为本发明的交轨干涉测量示意图；

图2为本发明的方位向滤波示意图；

图3为本发明的影像分块示意图；

图4为本发明的干涉相位标准差与相干系数的关系图；

图5为本发明的船载SAR成像结果图；

图6为本发明的两种配准方法相干系数分布图；

图7为本发明的基线长区间：0.16～0.3m统计直方图；

图8为本发明的基线长区间：0.166～0.18m统计直方图。

具体实施方式

为了加深对本发明的理解，下面将结合实施例对本发明做进一步详述，本实施例仅用于解释本发明，并不构成对本发明保护范围的限定。

根据图1、2、3、4所示，本实施例提供了一种船载SAR交轨干涉处理方法，包括以下步骤：

步骤一：前置滤波

从频谱特性考虑干涉图的前置滤波，即滤波的操作针对影像的频谱进行，如图1，由于主辅影像方位向和距离向的频谱不完全重叠，在精配准之前对主辅影像进行距离向滤波和方位向滤波，降低频谱偏移去相干，由于船载交轨干涉测量的基线较小，基线距导致的主辅影像地距频谱范围可忽略不计，只需考虑方位向滤波，由于成像时船载平台的位置和地面目标点之间存在着相对运动导致接收到的频率经历了多普勒频移，多普勒频率表达式为：

其中，Λ_s为传感器斜角，θ为雷达视角，v为雷达平台对地面速度，λ为波长，干涉处理中主辅图像对同一地面分辨率单元成像时，由于θ不同导致多普勒中心频率DCF_M与DCF_S不同，主辅影像的多普勒中心频率存在较大的差异，频谱偏移为：

Δf_Λ＝f_DCFm-f_DCFs (2)

频谱偏移导致主辅影像的去相干，船载SAR方位向滤波的过程如图2，利用B_w表示方位向带宽，利用f₁,f₂分别表示主辅图像考虑多普勒效应后的中心频率，利用B_new表示滤波后新的方位向带宽，则B_new满足:

B_new＝B_w-|f₁-f₂|＝B_w-|f_DC-f'_DC| (3)

在式(3)中，当|f_DC-f_DC'|>B_W，则两幅图像方位向频谱没有重叠部分，两幅图像毫不相干，无法生成干涉图，由此得到两幅图像的多普勒质心对生成的干涉图的质量至关重要，采用矩形滤波器、Hamming滤波器进行方位向滤波；

矩形滤波器：

Hamming滤波器：

其中：α∈(0,1)，实际运算时常取α＝0.75；

步骤二：块影像粗配准

对主辅影像A和B进行分块，得到块影像集合{A₁,A₂,…,A_m}和{B₁,B₂,…,B_m}，块的大小根据成像地区复杂性进行调整，基于块影像间的精准偏移量{Δ₁,Δ₂,…,Δ_m}得到全图整体偏移量，如图3，先基于块影像间的精准偏移量，得到全图整，偏移量，先，于轨道参数获得初始偏移量，配准初始偏移量为主图像块与辅图像块间得到的初始偏移量，将主辅影像块的中心点对应计算得到的偏移量作为配准初始偏移量即：

p_sla(m,n)＝p_mas(m,n)+offset(m,n) (6)

式中p_mas(m,n)表示主影像块中心点，p_sla(m,n)为辅影像块中心点，offset(m,n)表示两者偏移量，由于船载交轨干涉SAR得到的主辅影像及其相似，影像偏移量与基线相关，配准初始偏移量计算精度在5个像元以内，利用成像几何模型，首先计算主影像块中心点对应的地面坐标，然后求出地面点在辅图像上的对应点的坐标，从而解得主辅图像块的配准初始偏移量，根据成像模型计算主图像块的中心点对应的地面位置矢量p_sla(m,n)，即：

式中，f_D为多普勒频率，λ为波长，R为斜距，然后采用牛顿迭代法，求解出p(X_t,Y_t,Z_t)在辅图像上的对应点p_sla(m,n)所对应的方位向时间t，同样，初始值t₀用辅图像中心点所对应的时间表示，在得到点p_sla(m,n)所对应的方位向时间t后，得到辅图像上的行m和卫星的位置，再计算该点对应的斜距R，从而解求出距离向时间，最后得到辅图像上的列n，将同一点在两幅图像上的坐标相减即得到两幅图像的偏移；

然后基于匹配测度进行块影像的像元级配准，在影像的空间域、频率域进行，根据粗匹配结果在辅图像上选取搜索窗口，然后根据配准评价指标计算主影像块和窗口的相似程度，通过移动搜索窗得到最佳配准量，达到像元级的配准精度，每个主影像块计算得到辅影像对应的像元级偏移量，对于每一个待匹配的主影像块，在搜索范围(一般为10*10个像元)内滑动窗口取得辅影像块，计算主辅影像块的配准评价参数，选择相关系数作为配准的指标，相关系数最大时得到像元匹配量：

其中：corr为相关系数；g_i，j；g_i+r，j+c为主辅影像对应点处的振幅强度；m,n为匹配窗的大小，r,c为粗匹配的偏移量；

步骤三：块影像精配准

块影像精配准是得到辅影像对应主块影像的子像素级偏移量的过程，先基于原始图像过采样数据进一步寻找精确的配准位置，称之像元过采样匹配法，将主辅影像本身进行过采样处理，插值的间隔决定了过采样的程度，对图像进行插值采样将图像放大，之后采用复影像的配准方法逐点计算匹配测度，寻找可靠的相对偏移量估算值，在这一步中选择小窗口，但窗口的数目增多，以进行相应的数据拟合，由于窗口大小变小，增大搜索窗的大小多次计算，来防止出现最大相干偏差，当相对偏移量保持稳定，表明这个估算值可靠，否则开始新一轮搜索；过采样图像的实现有空间域和频率域法两种，上述为空间域法，而频域法运用傅里叶变换将局部图像转换到频域，共轭相乘后再进行逆FFT变换得到相关函数，省去了循环搜索，提高了计算速度；

步骤四：影像全局精配准

根据步骤三得到辅影像对应每块主影像子像元级配准偏移量，通过插值，到整个主辅，像的配准偏移量，差值通过一维三次样条函数插值条件实现：

其中x_i，Δ_i分别为块中心点像元坐标和配准偏移量，S(x)为得到的插值函数，偏移量是二维的，包含距离向和方位向偏移量，将二维配准偏移量的插值转化为二次一维三次样条插值即可；

步骤五：干涉图滤波

由干涉图的特点可知，由于受到各种因素的影响，干涉图中不可避免的含有大量的噪声，干涉图滤波就是去除和抑制干涉图噪声的过程，因此利用Goldstein和Werner针对干涉条纹特性提出的一种自适应滤波方法，表达式如下：

H(u,v)＝S{|Z(u,v)|}^αZ(u,v) (10)

其中，H(u,v)是滤波响应；S{}是平滑项；u和v是空间频率，α是滤波参数，取值范围为[0，1]，α值越大滤波效果越显著，

然后在Goldstein滤波基础上引入相干系数代替α作为滤波参数来控制滤波程度的强弱，单视情况下(L＝1)干涉相位标准差表达式：

其中

从式(11)看出干涉相位方差是相干系数γ的函数，干涉相位标准差与相干系数之间的关系如图4，利用γ作为滤波参数在相干性好、噪声影响小的地方滤波强度弱，保留原始相位信息，对噪声影响大的地方滤波强度会变强，滤波器函数式转化为：

H(u,v)＝S{|Z(u,v)|}^1-γZ(u,v) (12)

该做法不但确保了相干信息和滤波的结果在同一位置具有同样的大小，而且它还避免了滤波结果被重叠部分的相干值影响。这种改进可以防止goldstein滤波器在高相干区域进行过度滤波，同时保证在低相干地区实施较强滤波

步骤六：相位解缠

干涉测量的核心是获取与距离成比例的相位信息，然而干涉图包含的相位信息是经2π周期缠绕后的相位，即干涉相位值在(-π,π]之间，要获取正确的高程信息，必须进行相位解缠：

记φ_i,j为周期缠绕前的真实相位值，

为相应的缠绕相位，二维缠绕相位与解缠相位的关系为：

式中：k_i,j为第i行第j列个像元到卫星的整周模糊度，理想情况下干涉相位图中缠绕相位

在x方向和y方向差分与解缠相位φ_i,j在这两个方向的差分完全吻合，即

式中，

分别表示x,y方向的观测相位差，Δk_i,j为基于先验知识选取使Δ_i,j∈(-π,π](l＝1,2)成立的数组，由于积分路径不同，Δ_i,j不能和相邻点的差分保持一致，因而定义以下差分的残差：

和

用如下的最小化问题来估算：

其中，N表示图像的N个点，M表示剔除误匹配后的点数目，M-1表示最后一个点，M-2表示相邻的点，

和

是加权函数，同时满足下列约束条件：

以上所述问题为非线性最小化问题，转化为线性最小化问题提高计算效率：

那么上述问题就转化为求解

使得：

同时满足下列约束条件：

其中

均为大于等于0的实数，计算

和

后再计算相位梯度，最终通过下式计算出解缠后的相位得到最终的解缠结果。

步骤七：基于图像同名点的基线估计

基线是干涉测量中一个非常重要的参数，它决定了有效干涉相对的选择，在知道天线和地面点精确位置情况下容易计算基线和基线定向参数，当不知道天线位置就必须考虑采用特定的方法进行基线估计，先利用基于图像同名点的基线估计方法，通过得到足够数量的主辅影像同名点进行基线估计，图像同名点通过点匹配算法得到，并使用RANSAC算法以及相干系数图进一步消除误匹配点：

式中，k_i表示sift同名点对，i＝1,2…,N，Cor_i为该点对干涉处理得到的相干系数，P为RANSAC算法收敛后的满足模型的局内点集合，ε为相干性阈值，一般取0.9，剔除误匹配点后的同名点结果记为

i＝1,2…,M；最终得到的精确匹配点对可用于船载SAR基线的估计：

式中，t_i、t_i′分别为同名点对

在主辅影像的方位向时间，R_s,i、R′_s,i分别为同名点对

对应的主辅影天线相位中心的位置矢量,R₀、R′₀分别为主辅影像相位中心的初始位置矢量，B为基线长度估计结果，α为B取得最小值式对应的基线夹角；

步骤八：基于控制点的基线估计

然后利用基于控制点的基线估计方法，基线由主辅天线相位中心组成，与天线的物理中心不重合，此时根据R-D定位模型还原SAR成像几何得到精确的基线长度，R-D模型如下：

式中，R_S，V_S，R_G，V_G，f_D分别表示雷达成像中心的位置和速度矢量、地面点的位置和速度矢量和多普勒频率，

经运动补偿后，最终平台成像状态归化为匀速直线状态，同时平台的姿态保持不变，成像时的多普勒频率对于视场中的点为常量(f_D不变)，R_S与方位向成像时间t有关：

R_S＝V_S*t+R₀ (25)

R₀为起始雷达成像位置矢量，因此结合式(20)和(21)得到只与t有关的方程，通过给定初值迭代解出方位向时间t和控制点对应的天线相位中心的位置矢量，得到控制点对应的主辅天线相位中心的位置矢量R_S和R′_S后，基线估计结果为：

式中，B为基线长度估计结果，H_S、H_S分别为主辅天线相位中心的高程，α为基线倾角估计结果。

试验例：

试验使用中国科学院电子学研究所研制的Geo-MiniSAR设备，搭建了船载双天线InSAR系统，结合高精度POS系统和聚焦算法，可同时得到双通道高分辨率SAR影像。采集时间2018年6月，成像区域为湖北省武汉市某地区,船载运行轨道为东风大道高架桥上一段(，船载SAR影像主要参数为表1，成像效果为图5。

表1船载SAR影像的主要参数

为船载影像分块配准方法的精度进行分析，将其与星载SAR配准经常使用的由粗到精的多项式配准方法进行对比实验。对船载SAR主辅影像分别进行由粗到精的多项式配准和分块配准得到相干系数，如图6，低相干区域(0～0.5)相干性低，无比较意义，中高相干区(0.5～1)分块配准在高相干区占占比很高，而由粗到精的多项式配准得到相干系数分布较为均一，可知分块配准得到的配准结果更为可靠。

经过船载SAR影像配准后，进行干涉计算，对干涉图使用Goldstein进行滤波得到，滤波后高相干性区域的干涉条纹较清晰，低相干性区域，如植被和和海拔突变的地方，干涉条纹模糊，没有进行滤波。综合分析船载SAR主辅影像的相干系数及解缠相位，可知阴影区和植被区相干性较低，相位变化复杂；道路、裸地和草地等区域相干性高，相位变化平滑，相位解缠结果符合实际情况。

然而进行船载SAR基线估计的试验，根据剩余匹配同名点计算得到基线，大概统计结果如图7，进一步减小统计区间为0.166～0.180m，得到图8，可知出现频率最高的基线长为基线的最小值，所以将估计值取所有基线结果的最小值，即式(23)是合理的，表2为前述两种基线估计方法的结果，两种方法的结果接近，结果误差在亚毫米级别，在没有控制点的情况下，基于图像同名点的船载SAR基线估计方法是可靠的。

表2基线估计结果

该船载SAR交轨干涉处理方法包括主辅影像滤波、配准、干涉计算和干涉图的滤波、相位解缠等工作，针对船载双天线SAR的成像特点，图像配准采用了分块配准方法，利用块影像粗配准和块影像精配准，提高图像的相干系数，达到了较高的相干性；通过干涉图滤波确保了相干信息和滤波的结果在同一位置具有同样的大小，避免了滤波结果被重叠部分的相干值影响；通过相位解缠获取正确的高程信息；通过基于图像同名点的基线估计方法和基于控制点的基线估计方法，实现毫米甚至亚毫米级的基线估计精度，提高后续的DEM求解精度。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (6)

1.一种船载SAR交轨干涉处理方法，其特征在于，包括以下步骤：

步骤一：前置滤波

从频谱特性考虑干涉图的前置滤波，即滤波的操作针对影像的频谱进行，由于主辅影像方位向和距离向的频谱不完全重叠，在精配准之前对主辅影像进行距离向滤波和方位向滤波，降低频谱偏移去相干，由于成像时船载平台的位置和地面目标点之间存在着相对运动导致接收到的频率经历了多普勒频移，多普勒频率表达式为：

其中，Λ_s为传感器斜角，θ为雷达视角，v为雷达平台对地面速度，λ为波长，干涉处理中主辅图像对同一地面分辨率单元成像时，由于θ不同导致多普勒中心频率DCF_M与DCF_S不同，主辅影像的多普勒中心频率存在较大的差异，频谱偏移为：

频谱偏移导致主辅影像的去相干，利用B_w表示方位向带宽，利用f₁,f₂分别表示主辅图像考虑多普勒效应后的中心频率，利用B_new表示滤波后新的方位向带宽，则B_new满足:

B_new＝B_w-|f₁-f₂|＝B_w-|f_DC-f'_DC| (3)

在式(3)中，当|f_DC-f_DC'|>B_W，则两幅图像方位向频谱没有重叠部分，两幅图像毫不相干，无法生成干涉图，由此得到两幅图像的多普勒质心对生成的干涉图的质量至关重要，采用矩形滤波器、Hamming滤波器进行方位向滤波；

矩形滤波器：

Hamming滤波器：

其中：α∈(0,1)，实际运算时常取α＝0.75；

步骤二：块影像粗配准

对主辅影像A和B进行分块，得到块影像集合{A₁,A₂,…,A_m}和{B₁,B₂,…,B_m}，块的大小根据成像地区复杂性进行调整，基于块影像间的精准偏移量{Δ₁,Δ₂,…,Δ_m}得到全图整体偏移量，先基于轨道参数获得初始偏移量，配准初始偏移量为主图像块与辅图像块间得到的初始偏移量，将主辅影像块的中心点对应计算得到的偏移量作为配准初始偏移量即：

p_sla(m,n)＝p_mas(m,n)+offset(m,n) (6)

式中p_mas(m,n)表示主影像块中心点，p_sla(m,n)为辅影像块中心点，offset(m,n)表示两者偏移量，由于船载交轨干涉SAR得到的主辅影像及其相似，影像偏移量与基线相关，配准初始偏移量计算精度在5个像元以内，利用成像几何模型，首先计算主影像块中心点对应的地面坐标，然后求出地面点在辅图像上的对应点的坐标，从而解得主辅图像块的配准初始偏移量，根据成像模型计算主图像块的中心点对应的地面位置矢量p_sla(m,n)，即：

式中，f_D为多普勒频率，λ为波长，R为斜距，然后采用牛顿迭代法，求解出p(X_t,Y_t,Z_t)在辅图像上的对应点p_sla(m,n)所对应的方位向时间t，同样，初始值t₀用辅图像中心点所对应的时间表示，在得到点p_sla(m,n)所对应的方位向时间t后，得到辅图像上的行m和卫星的位置，再计算该点对应的斜距R，从而解求出距离向时间，最后得到辅图像上的列n，将同一点在两幅图像上的坐标相减即得到两幅图像的偏移；

然后基于匹配测度进行块影像的像元级配准，在影像的空间域、频率域进行，根据粗匹配结果在辅图像上选取搜索窗口，然后根据配准评价指标计算主影像块和窗口的相似程度，通过移动搜索窗得到最佳配准量，达到像元级的配准精度，每个主影像块计算得到辅影像对应的像元级偏移量，对于每一个待匹配的主影像块，在搜索范围内滑动窗口取得辅影像块，计算主辅影像块的配准评价参数，选择相关系数作为配准的指标，相关系数最大时得到像元匹配量：

其中：corr为相关系数；g_i，j；g_i+r，j+c为主辅影像对应点处的振幅强度；m,n为匹配窗的大小，r,c为粗匹配的偏移量；

步骤三：块影像精配准

块影像精配准是得到辅影像对应主块影像的子像素级偏移量的过程，先基于原始图像过采样数据进一步寻找精确的配准位置，称之像元过采样匹配法，将主辅影像本身进行过采样处理，插值的间隔决定了过采样的程度，对图像进行插值采样将图像放大，之后采用复影像的配准方法逐点计算匹配测度，寻找可靠的相对偏移量估算值，在这一步中选择小窗口，但窗口的数目增多，以进行相应的数据拟合，由于窗口大小变小，增大搜索窗的大小多次计算，来防止出现最大相干偏差，当相对偏移量保持稳定，表明这个估算值可靠，否则开始新一轮搜索；

步骤四：影像全局精配准

根据步骤三得到辅影像对应每块主影像子像元级配准偏移量，通过插值，到整个主辅，像的配准偏移量，差值通过一维三次样条函数插值条件实现：

其中x_i，Δ_i分别为块中心点像元坐标和配准偏移量，S(x)为得到的插值函数，偏移量是二维的，包含距离向和方位向偏移量，将二维配准偏移量的插值转化为二次一维三次样条插值即可；

步骤五：干涉图滤波

由干涉图的特点可知，由于受到各种因素的影响，干涉图中不可避免的含有大量的噪声，干涉图滤波就是去除和抑制干涉图噪声的过程，因此利用Goldstein和Werner针对干涉条纹特性提出的一种自适应滤波方法，表达式如下：

H(u,v)＝S{|Z(u,v)|}^αZ(u,v) (10)

其中，H(u,v)是滤波响应；S{}是平滑项；u和v是空间频率，α是滤波参数，取值范围为[0，1]，α值越大滤波效果越显著，

然后在Goldstein滤波基础上引入相干系数代替α作为滤波参数来控制滤波程度的强弱，单视情况下(L＝1)干涉相位标准差表达式：

其中

从式(11)看出干涉相位方差是相干系数γ的函数，利用γ作为滤波参数，滤波器函数式转化为：

H(u,v)＝S{|Z(u,v)|}^1-γZ(u,v) (12)；

步骤六：相位解缠

干涉测量的核心是获取与距离成比例的相位信息，然而干涉图包含的相位信息是经2π周期缠绕后的相位，即干涉相位值在(-π,π]之间，要获取正确的高程信息，必须进行相位解缠：

记φ_i,j为周期缠绕前的真实相位值，

为相应的缠绕相位，二维缠绕相位与解缠相位的关系为：

式中：k_i,j为第i行第j列个像元到卫星的整周模糊度，理想情况下干涉相位图中缠绕相位

在x方向和y方向差分与解缠相位φ_i,j在这两个方向的差分完全吻合，即

式中，

分别表示x,y方向的观测相位差，Δk_i,j为基于先验知识选取使Δ_i,j∈(-π,π](l＝1,2)成立的数组，由于积分路径不同，Δ_i,j不能和相邻点的差分保持一致，因而定义以下差分的残差：

和

用如下的最小化问题来估算：

其中，N表示图像的N个点，M表示剔除误匹配后的点数目，M-1表示最后一个点，M-2表示相邻的点，

和

是加权函数，同时满足下列约束条件：

以上所述问题为非线性最小化问题，转化为线性最小化问题提高计算效率：

那么上述问题就转化为求解

使得：

同时满足下列约束条件：

其中

均为大于等于0的实数，计算

和

后再计算相位梯度，最终通过下式计算出解缠后的相位得到最终的解缠结果：

步骤七：基于图像同名点的基线估计

基线是干涉测量中一个非常重要的参数，它决定了有效干涉相对的选择，在知道天线和地面点精确位置情况下容易计算基线和基线定向参数，当不知道天线位置就必须考虑采用特定的方法进行基线估计，先利用基于图像同名点的基线估计方法，通过得到足够数量的主辅影像同名点进行基线估计，图像同名点通过点匹配算法得到，并使用RANSAC算法以及相干系数图进一步消除误匹配点：

式中，k_i表示sift同名点对，i=1,2…,N，Cor_i为该点对干涉处理得到的相干系数，P为RANSAC算法收敛后的满足模型的局内点集合，ε为相干性阈值，一般取0.9，剔除误匹配点后的同名点结果记为

最终得到的精确匹配点对可用于船载SAR基线的估计：

式中，t_i、t_i′分别为同名点对

在主辅影像的方位向时间，R_s,i、R′_s,i分别为同名点对

对应的主辅影天线相位中心的位置矢量,R₀、R′₀分别为主辅影像相位中心的初始位置矢量，B为基线长度估计结果，α为B取得最小值式对应的基线夹角；

步骤八：基于控制点的基线估计

然后利用基于控制点的基线估计方法，基线由主辅天线相位中心组成，与天线的物理中心不重合，此时根据R-D定位模型还原SAR成像几何得到精确的基线长度，R-D模型如下：

式中，R_S，V_S，R_G，V_G，f_D分别表示雷达成像中心的位置和速度矢量、地面点的位置和速度矢量和多普勒频率，

经运动补偿后，最终平台成像状态归化为匀速直线状态，同时平台的姿态保持不变，成像时的多普勒频率对于视场中的点为常量(f_D不变)，R_S与方位向成像时间t有关：

R_S＝V_S*t+R₀ (25)

R₀为起始雷达成像位置矢量，因此结合式(20)和(21)得到只与t有关的方程，通过给定初值迭代解出方位向时间t和控制点对应的天线相位中心的位置矢量，得到控制点对应的主辅天线相位中心的位置矢量R_S和R′_S后，基线估计结果为：

式中，B为基线长度估计结果，H_S、H_S分别为主辅天线相位中心的高程，α为基线倾角估计结果。

2.根据权利要求1所述的一种船载SAR交轨干涉处理方法，其特征在于：所述步骤一中，由于船载交轨干涉测量的基线较小，基线距导致的主辅影像地距频谱范围可忽略不计，只需考虑方位向滤波。

3.根据权利要求1所述的一种船载SAR交轨干涉处理方法，其特征在于：所述步骤二中，在搜索范围内滑动窗口取得辅影像块，搜索范围为10*10个像元。

4.根据权利要求1所述的一种船载SAR交轨干涉处理方法，其特征在于：所述步骤三中，过采样图像的实现有空间域和频率域法两种，频域法运用傅里叶变换将局部图像转换到频域，共轭相乘后再进行逆FFT变换得到相关函数，省去了循环搜索，提高了计算速度。

5.根据权利要求1所述的一种船载SAR交轨干涉处理方法，其特征在于：所述步骤五中，利用γ作为滤波参数在相干性好、噪声影响小的地方滤波强度弱，保留原始相位信息，对噪声影响大的地方滤波强度会变强。

6.根据权利要求1所述的一种船载SAR交轨干涉处理方法，其特征在于：所述步骤五中，式(12)确保了相干信息和滤波的结果在同一位置具有同样的大小，而且避免了滤波结果被重叠部分的相干值影响，防止goldstein滤波器在高相干区域进行过度滤波，同时保证在低相干地区实施较强滤波。

Images