一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法与流程

1.本发明属于测绘科学技术领域，尤其涉及一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法。


背景技术：

2.sar(合成孔径雷达)作为一种全天时、全天候、高分辨率的成像雷达，被广泛应用在遥感领域。由于sar图像具有相位信息，因此可以通过对sar图像序列的干涉处理，实现对成像区域形变信息的非接触式的高精度提取，鉴于sar系统的上述优势，它成为了地面形变监测的重要手段之一。传统的sar多搭载在卫星和飞机等平台上，这使得雷达的重访时间很长，不能对特定目标持续监测。地基sar(ground based sar，gbsar)是一种小型化、低成本的sar系统，可实现对特定区域的持续观测，它的重访周期短，监测精度高，结合干涉和差分干涉技术，可以对一个特定场景长时间连续形变监测，被广泛应用于塌陷、滑坡等灾害预警。然而，通常gbsar都安装在直线轨道上，由于轨道长度有限，gbsar在方位向的观测范围受到了限制，另一种在上世纪90年代初提出的旋转式合成孔径雷达(rotor sar，rosar)，通过将雷达置于旋转机械臂上，从而获得了更广的观测范围，rosar通常安装在直升机旋翼上，并使用脉冲形式的信号，这使得rosar的系统十分复杂。基于rosar的运动形式，一种使用调频连续波(frequency modulated continuous wave，fmcw)体制的改进的旋转式sar在2012年被提出，并命名为圆弧sar(arc-sar)。arc-insar是一种将干涉测量与arc-sar相结合的高精度的形变测量技术，它通过对同一位置、不同时刻获取的两幅sar图像进行差分干涉处理，基于相位信息来实现形变测量，能够实现在有限空间内对地面大范围、高精度的持续监测。arc-insar一般工作在x或者ku波段，形变测量精度可以达到毫米或者亚毫米量级。
3.大气干扰是由于雷达电磁波在空气传播过程中，受时空变化的大气压力、温度、湿度等因素导致不同相位延迟所引起的相位误差。大气干扰严重影响着insar干涉测量的形变反演精度和可靠性，是阻碍insar干涉测量应用发展的重要限制因素之一，而对于具有大范围形变测量能力的arc-insar而言，大气干扰更是需要亟待解决的实际问题。目前，insar大气干扰相位校正方法主要分为两大类：(1)基于外部辅助数据(2)基于干涉相位本身。基于外部辅助数据的方法主要有基于地表气象数据的建模、基于gps数据反演大气干扰相位、基于空间辐射计测量法等，这些基于外部辅助数据的相位校正方法主要受辅助数据低空间分辨率的影响，校正效果往往不是十分理想。基于干涉相位本身的校正方法有干涉相位滤波法、基于时间序列sar图像的校正方法和基于稳定控制点的大气干扰相位校正算法，其中干涉相位滤波法、基于时间序列sar图像的校正方法又面临着分离大气干扰相位与非线性形变的难题，故目前应用最多的insar大气干扰相位校正方法是基于稳定控制点的大气干扰相位校正，但该方法虽然有良好的校正效果，但它是在建立假定大气干扰相位延迟在方位向一致的条件下提出的，未考虑复杂情况下大气干扰相位在方位向的空间变化特性，会使得空变大气干扰相位校正过程中出现方位向欠校正的问题。
4.因此，如何提出一种不引入外部辅助数据而仅需依靠干涉相位本身就能进行大气
干扰相位校正的方法是本技术领域人员亟待解决的问题。


技术实现要素：

5.针对现有技术的不足，本发明的目的是提供一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法，以解决现有技术中地基圆弧干涉合成孔径雷达大气干扰相位校正方法需要引入外部数据同时精确度不高的问题。
6.为了解决上述技术问题，本发明采用了如下的技术方案：
7.本发明提供了一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法，包括以下步骤：
8.s10、构建地基圆弧干涉合成孔径雷达差分干涉相位模型，所述相位模型具体可表示为：
[0009][0010]
其中，为主图像u与从图像v的干涉相位，为形变相位，为大气干扰相位，为噪声相位；
[0011]
s20、对所述干涉相位进行滤波处理，消除相位噪声；
[0012]
s30、根据统计特征量从所述干涉相位中筛选出永久散射点，其中所述永久散射点为所述干涉相位近似为零的控制点；
[0013]
s40、对所述步骤s20获得的所述干涉相位在方位向上进行子空间划分，再根据设定窗大小对子空间划分网格，根据相关系数和幅度离差对网格内的所述永久散射点进行相位加权平均处理后，作为网格中心点的大气干扰相位样本点；
[0014]
s50、对所述步骤s40获得的各子空间中的所述大气干扰相位样本点进行筛选，得到可用样本点；
[0015]
s60、利用所述可用样本点估算对应子空间的所述大气干扰相位误差模型参数，进而获取对应探测子空间的所述大气干扰相位误差；
[0016]
s70：获得补偿后的全部探测空间的所述形变相位。
[0017]
进一步的，所述步骤s20中的滤波处理为加权圆周期中值滤波，加权系数为：
[0018][0019]
其中，(kn，km)为干涉图像中所述永久散射点的坐标序号，n为待滤波处理的所述干涉相位的横坐标，m为待滤波处理的所述干涉相位的纵坐标，其中，m为滤波窗在横轴上的长度，n为滤波窗在纵轴上的长度，
[0020]
进一步的，所述步骤s30中筛选出永久散射点的步骤包括：
[0021]
s31、根据主从复散射图像计算相关系数γ，公式如下：
[0022][0023]
其中，m(i，j)和所述s(i，j)分别为构成所述干涉相位的主复散射图像和从复散射图像，*表示复共轭，w1为第一滑动窗大小，w2为第二滑动窗大小；
[0024]
s32、利用n
p
幅时序复散射图像，可计算得到平均相关系数公式如下：
[0025][0026]
其中，q表示第q帧所述干涉相位，γq表示第q帧的相关系数；
[0027]
s33、设置第一阈值门限t1，筛选出所述的像素点作为所述永久散射点的候选点；
[0028]
s34、计算时序复散射图像的幅度离差da，公式如下：
[0029][0030]
其中，δa表示像素点幅度时间序列的标准差，ma像素点幅度时间序列的均值；
[0031]
s35、设置第二阈值门限t2，从所述步骤s31所得所述永久散射点的候选点集中筛选出所述da≤t2的所述永久散射点。
[0032]
进一步的，所述步骤s32的所述第一阈值门限t1设置为0.9，所述第二阈值门限t2设置为0.1。
[0033]
进一步的，所述步骤s40中获取所述大气干扰相位样本点的步骤包括：
[0034]
s41、对所述干涉相位在方位向按角度θ
azm
对图像探测范围进行子空间划分；
[0035]
s42、对子空间按照w
x
*wy的窗大小进行空间网格划分，其中w
x
与wy分别代表子空间划分时窗的横纵长度，利用平均相关系数对各网格内所有所述永久散射点的所述干涉相位进行加权平均，得到第一均值相位其中加权值为r；
[0036]
s43、利用幅度离差对子网格内所有所述永久散射点的所述干涉相位再次进行加权平均，得到第二均值相位其中加权值为1-da；
[0037]
s44、将作为网格中心点的所述大气干扰相位近似值，同时计算该网格中心点的与雷达的径向距离rw，即为所述大气干扰相位样本点，其中，所述w为空间网格编号；
[0038]
s45、按所述步骤s42至所述步骤s44遍历所有已划分的子空间。
[0039]
进一步的，所述步骤s41的所述角度θ
azm
设置为30
°
。
[0040]
进一步的，所述步骤s50中得到所述可用样本点的步骤包括：
[0041]
s51、统计各子空间已划分网格内所述永久散射点数量占当前网格总像素点的百分比pw(％)；
[0042]
s52、设置第三阈值门限t3，当所述pw≥t3时，即为可用的所述大气干扰
相位样本点。
[0043]
进一步的，所述步骤s52中所述第三阈值门限t3设置为10。
[0044]
进一步的，所述步骤s60包括：
[0045]
s61、将所述大气干扰相位建模成随斜距变化的线性模型：
[0046][0047]
其中，β0为常数分量，β1为与斜距相关的线性系数，λ为系统发射信号波长；
[0048]
s62、建立线性方程组：
[0049][0050]
其中，β＝[β
0 β1]
t
，ε＝[ε1，ε2，
…
，εw]
t
为随机误差向量，t表示矩阵转置；
[0051]
s63、利用最小二乘法对β进行估算，可得：
[0052][0053]
则所述大气干扰相位估值为：
[0054][0055]
其中，与之间的差值即为补偿后的所述干涉相位；
[0056]
s64、重复所述步骤s61至所述步骤s63，即可完成对全部探测空间的所述大气干扰相位的估算。
[0057]
本发明提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法与现有技术相比，至少具有如下有益效果：
[0058]
本发明引入了永久散射点的概念，通过相关系数均值和幅度离差对永久散射点进行筛选，提取出形变相位近似为零的稳定永久散射点，并以相关系数和幅度离差对稳定永久散射点干涉相位进行加权平均作为大气干扰相位样本点，同时对探测区域按方位向进行子空间划分，并对各子空间再次进行空间网格划分，利用网格内的大气干扰相位样本点的确定等效大气干扰相位样本点，进而完成大气干扰相位模型参数的最小二乘估算，在不引入外部辅助数据的情况下，依靠干涉相位本身进行了大气干扰相位的有效校正；充分考虑了大气干扰相位二维空变特性，通过子空间和网格划分解决了大气干扰相位方位向欠校正问题；方法流程简单易实现，且具有较强的适应性与鲁棒性，对大气干扰相位的校正效果良好，符合工程实际需求。
附图说明
[0059]
为了更清楚地说明本发明的方案，下面将对实施例描述中所需要使用的图作一个简单的介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0060]
图1为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法的流程图；
[0061]
图2为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法中永久散射点筛选方法流程图；
[0062]
图3为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法中大气干扰相位样本点筛选流程图。
具体实施方式
[0063]
为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的较佳实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。
[0064]
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。
[0065]
如图1所示为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法流程图，该方法考虑了大气干扰相位的二维空间变化特性，通过引入了永久散射点(permanent scatters，ps)的概念，利用相关系数均值和幅度离差对永久散射点进行筛选，提取出形变相位近似为零的稳定永久散射点，并以相关系数和幅度离差对稳定永久散射点干涉相位进行加权平均作为大气干扰相位样本点，同时对探测区域按方位向进行子空间划分，并对各子空间再次进行空间网格划分，利用网格内的大气干扰相位样本点的确定等效大气干扰相位样本点，进而完成大气干扰相位模型参数的最小二乘估计，在不引入外部辅助数据的情况下，依靠干涉相位本身进行了大气干扰相位的有效校正。
[0066]
具体步骤如下：
[0067]
s10、地基圆弧干涉合成孔径雷达系统实际工作时，受大气干扰相位误差和噪声相位误差影响，干涉相位满足如下模型：
[0068][0069]
其中，为主图像u与从图像v的干涉相位，为形变相位，为大气干扰相位，为噪声相位；
[0070]
s20、对干涉相位进行滤波处理，消除相位噪声；
[0071]
鉴于地基圆弧干涉合成孔径雷达形变测量对噪声抑制的较高要求，噪声相位滤波算法需要兼具效率与性能，加权圆周期中值滤波算法是目前被广泛采用的地基合成孔径雷达相位滤波算法，表达式(1)干涉相位经滤波处理后可的到如下输出：
[0072][0073]
其中，(kn，km)为干涉图像中ps点的坐标序号，n为待滤波处理的所述干涉相位的横坐标，m为待滤波处理的所述干涉相位的纵坐标，加权系数为
其中，m为滤波窗在横轴上的长度，n为滤波窗在纵轴上的长度，(median为中值滤波)，
[0074]
s30、根据统计特征量从干涉相位中筛选出永久散射点(permanent scatters，ps)，其中永久散射点为干涉相位近似为零的控制点，即因此经过步骤s20滤波后的干涉相位可表示为：
[0075][0076]
该步骤中永久散射点的筛选特征量主要包括相关系数和幅度离差，通过二级筛选完成稳定永久散射点的筛选；
[0077]
s40、对步骤s20获得的干涉相位在方位向上进行子空间划分，再根据设定窗大小对子空间划分网格，根据相关系数和幅度离差对网格内的永久散射点进行相位加权平均处理后，作为网格中心点的大气干扰相位样本点；
[0078]
s50、对步骤s40获得的各子空间中的大气干扰相位样本点进行筛选，得到可用样本点；
[0079]
s60、利用可用样本点估计对应子空间的大气干扰相位误差模型参数，进而获取对应探测子空间的大气干扰相位误差；
[0080]
s70、根据步骤s10至步骤s60即可求得补偿后的全部探测空间的形变相位。
[0081]
如图2所示为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法中步骤s30永久散射点筛选方法流程图，具体包括如下步骤：
[0082]
s31：由于相关系数能够反映干涉相位图的信噪比，因此可通过设置合理的相关系数阈值门限(本实施例中设定为[0.9，1])，选取均值相关系数较大的散射点作为永久散射点的候选点集，相关系数可根据待计算像素点与周围邻域点的幅相信息计算得到，如表达式(4)所示：
[0083][0084]
其中，m(i，j)和s(i，j)分别为构成干涉相位的主从复散射图像，*表示复共轭，w1为第一滑动窗大小，w2为第二滑动窗大小；
[0085]
s32、利用n
p
幅时序复散射图像，可计算得到n
p-1个相关系数矩阵，则平均相关系数可由表达式(5)计算：
[0086][0087]
其中，q表示第q帧所述干涉相位，γq表示第q帧的相关系数；
[0088]
s33：对于高质量的像素点，时序幅度信息可以表征其辐射稳定性，而幅度离差则是时序幅度信息的直接表征量，故可设置合理的阈值门限(本实施例中设定为[0，0.1])，利
用幅度离差来实现永久散射点的二次筛选；
[0089]
s34、计算时序复散射图像的幅度离差da：
[0090][0091]
其中，δa表示像素点幅度时间序列的标准差，ma像素点幅度时间序列的均值。由表达式(6)可知，像素点的幅度离差越小，其辐射稳定性越高，对应的相位变化就越小，即干涉相位约接近于零。
[0092]
如图3所示为本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法中步骤s40大气干扰相位样本点筛选流程图，具体包括如下步骤：
[0093]
s41：由于地基圆弧干涉合成孔径雷达通常工作于大范围扫描模式(大角度、远距离)，使得大气干扰相位误差具有二维空变特性，但在局部小区域内，大气的影响依然具有一致性，因此可以考虑对探测空间在方位向以角度间隔θ
azm
(本实施例中设定为30
°
)进行子空间划分，每个子空间单独进行距离向的一维大气误差估计与补偿；
[0094]
s42：对子空间按照w
x
*wy(本实施例中设定为30*30)的窗大小进行空间网格划分，其中w
x
与wy分别代表子空间划分时窗的横纵长度，利用平均相关系数对各网格内所有永久散射点的干涉相位进行加权平均，得到均值相位其中加权值为r；
[0095]
s43：利用幅度离差对子网格内所有永久散射点的干涉相位再次进行加权平均，得到均值相位其中加权值为1-da；
[0096]
s44：将作为网格中心点的大气干扰相位近似值，同时计算该网格中心点的与雷达的径向距离rw，即为大气干扰相位样本点，其中，w为空间网格编号；
[0097]
s45：按步骤s42至步骤s44遍历所有已划分的子空间，获得各子空间对应的大气干扰相位点集。
[0098]
为提高大气干扰相位参数的精度，需要对等效大气样本点的质量进行筛选，本发明实施例提供的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法中步骤s50给出了可用的等效大气干扰相位样本点筛选过程，具体包括如下步骤：
[0099]
s51：统计各子空间已划分网格内永久散射点数量占当前网格总像素点的百分比pw(％)，算式如表达式(7)所示：
[0100]
pw＝n
ps
/(w
x
*wy)*100
ꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0101]
s52、设置阈值门限t3(本实施例中设定为10)，当pw≥t3时，即为可用的等效大气干扰相位样本点；否则，直接丢弃。
[0102]
遍历当前子空间所有的空间网格，可以得到可用的等效大气干扰相位样本点集。
[0103]
由于在方位向进行了子空间划分，因此在各子空间内，大气在方位向上不存在空变性，大气干扰相位误差只与斜距有关。那么通过采取本实施例中步骤s60即可完成大气干扰相位误差的估计与补偿，具体步骤如下：
[0104]
s61：构建大气干扰相位随斜距变化的线性模型，如表达式(8)所示：
[0105][0106]
其中，β0为常数分量，β1为与斜距相关的线性系数；
[0107]
s62、建立线性方程组：
[0108][0109]
其中，β＝[β
0 β1]
t
，ε＝[ε1，ε2，
…
，εw]
t
为随机误差向量，t表示矩阵转置；
[0110]
s63、利用最小二乘法对β进行估算，可得：
[0111][0112]
则大气干扰相位估值为：
[0113][0114]
其中，与之间的差值即为补偿后的干涉相位。
[0115]
以上实施例所述的一种地基圆弧干涉合成孔径雷达大气干扰相位校正方法，引入了永久散射点的概念，通过相关系数均值和幅度离差对永久散射点进行筛选，提取出形变相位近似为零的稳定永久散射点，并以相关系数和幅度离差对稳定永久散射点干涉相位进行加权平均作为大气干扰相位样本点，同时对探测区域按方位向进行子空间划分，并对各子空间再次进行空间网格划分，利用网格内的大气干扰相位样本点的确定等效大气干扰相位样本点，进而完成大气干扰相位模型参数的最小二乘估算，在不引入外部辅助数据的情况下，依靠干涉相位本身进行了大气干扰相位的有效校正；充分考虑了大气干扰相位二维空变特性，通过子空间和网格划分解决了大气干扰相位方位向欠校正问题；方法流程简单易实现，且具有较强的适应性与鲁棒性，对大气干扰相位的校正效果良好，符合工程实际需求。
[0116]
显然，以上所描述的实施例仅仅是本发明较佳实施例，而不是全部的实施例，附图中给出了本发明的较佳实施例，但并不限制本发明的专利范围。本发明可以以许多不同的形式来实现，相反地，提供这些实施例的目的是使对本发明的公开内容的理解更加透彻全面。尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来而言，其依然可以对前述各具体实施方式所记载的技术方案进行修改，或者对其中部分技术特征进行等效替换。凡是利用本发明说明书及附图内容所做的等效结构，直接或间接运用在其他相关的技术领域，均同理在本发明专利保护范围之内。