一种InSAR非线性去平地效应方法

一种insar非线性去平地效应方法
技术领域
1.本发明涉及合成孔径雷达干涉处理领域，具体来说是一种insar非线性去平地效应方法，对干涉相位进行补偿，去除场景中呈非线性的平地相位。


背景技术：

2.合成孔径雷达干涉测量技术(interferometric synthetic aperture radar，insar)能够在全天时、全天候条件下大范围获取地面高程及形变信息，常用来取高精度的数字高程模型(digital elevation model，dem)数据，在地表形变探测、动目标检测、海洋测绘、森林制图、洪涝检测、交通监测和冰川研究等众多领域均具有重要的应用及研究价值。
3.根据空间几何关系，即便高度不变的平地在insar干涉相位图中也会产生距离向上呈周期性的密集条纹，致使干涉相位无法反映出实际的地面高程变化。因此有必要消除该平地相位，降低相位梯度，方便后续相位滤波和相位解缠处理。去除平地相位的过程即为去平地效应，或称平地相位补偿。
4.目前常用的去平地效应技术主要分为基于轨道参数、基于外部dem数据以及基于干涉条纹频率三种。在实际情景中，有时难以获取对应的精确轨道数据或外部dem数据，此时可采用第三种方法，通过估算出相位图中的最亮条纹频率并进行补偿，即频移以去除平地相位。此类方法无需提供干涉相位图之外的数据，且具有更快的计算速度，是最常用的去平地效应方法，但由于只补偿了固定为常数的线性相位，因而仅适用于估计均匀分布的平地条纹。
5.当下视角较小或坡度角较大时，干涉相位条纹将呈现明显的近端密集、远端稀疏现象，此时平地相位不再是均匀分布的线性相位，呈现近密远疏的情景。对此类场景，若依旧使用传统的去平地方法，由于估计的平地相位与真实相位存在较大差异，去平后的干涉相位图仍将在局部存在密集条纹，影响后续干涉处理。基于频移法进行改进的分块法虽然一定程度上能够改善上述问题，但针对不同场景分块的原则往往难以确定，且容易引入局部地形突变导致的局部估计错误，反而增加了人工工作量。


技术实现要素：

6.本发明为了在无监督条件下，尽可能精确去除小下视角或大坡度角场景的平地相位，降低后续干涉处理难度，提供了一种insar非线性去平地效应方法，以在缺少额外轨道和dem数据的情况下，能够无监督、高效并且准确去除场景的平地相位。
7.所述insar非线性去平地效应方法，包括以下几个步骤：
8.步骤一：合成孔径雷达对同一待测目标进行两次成像，得到两幅复图像并进行配准，对两幅配准后的sar复图像进行干涉，得到对应的干涉相位图；
9.在干涉相位图中，将含噪相位分解为平地相位与反应细节变化的残余相位与反应细节变化的残余相位
10.步骤二：对干涉复图像进行预滤波，并对预滤波后的干涉复图像距离向数据进行离散傅里叶变换并求和，从而得到对应于距离向的频谱幅值序列及相位序列。
11.根据离散傅里叶变换，沿距离向分布的预滤波后的干涉复图像复干涉相位分解公式：
[0012][0013]
其中，是含噪相位的复数形式，为虚数单位，m为坐标位置，g代表组成复干涉相位的矢量个数，ci为矢量的权值，代表对应矢量的距离向频率。
[0014]
在每组距离向上对复干涉相位序列进行离散傅里叶变换并求和，得到频谱幅值序列si的求和结果，获得反映整体距离向分布的频谱sr以及对应的频谱幅值序列p和相位序列φ：
[0015][0016]
其中代表第i个距离向上的复干涉相位序列，u为频谱序列的坐标位置，
⊙
表示shur乘积算子。
[0017]
步骤三：根据频谱幅值序列构造归一化直方图，根据直方图分布自适应选取阈值。
[0018]
根据直方图自适应选择阈值的实现过程为：
[0019]
(1)依据频谱幅值序列构造归一化直方图，对归一化直方图作二阶差分运算，并寻找其中所有的极大值点；
[0020]
对直方图序列h作二阶差分：
[0021]
d(k)＝[h(k+1)-h(k)]-[h(k)-h(k-1)],k＝2,3,...,n-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0022]
其中d为二阶差分序列，n为直方图组数，通过判断d(k)≥0即可找到直方图中的所有极大值点。
[0023]
(2)从左端起依次判断极大值点处对应的频率与直方图均匀分布下的频率之间的关系，若则代表噪声的频率水平；若则代表平地相位的频率水平，在满足平地相位频率条件的极大值点中选择对应直方图中幅值最小的极大值点作为阈值t。
[0024]
步骤四：针对于频谱幅值序列中低于阈值的部分，取其均值作为背景噪声功率估计值，并将该部分置零；对于高于阈值的部分，则减去背景噪声功率估计值以进一步抑制噪声的影响，经过以上操作后得到处理后的频谱幅值序列。
[0025]
背景噪声功率估计值为：noise＝mean(p(u)＜t)，mean(
·
)表示取均值操作。
[0026]
处理后的频谱幅值序列为：
[0027][0028]
步骤五：将处理后的频谱幅值序列与频谱相位序列重新结合，并作离散傅里叶逆变换后，取相角即可得到非线性平地相位序列，并将其扩展为原始尺寸。
[0029]
非线性平地相位序列为：
[0030]
φ
p
＝arg(ifft(p
′⊙
φ))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0031]
非线性平地相位序列扩展为原始尺寸的公式为：
[0032][0033]
步骤六：从干涉相位图中减去扩展为原始尺寸后的非线性平地相位序列，得到最终的去平地结果。
[0034]
本发明与现有技术相比的优点在于：
[0035]
(1)相较于传统频移法，本发明能够估计非线性平地相位，在小视角系统或地形复杂情况下去平地结果更准确，去平地后的干涉条纹更稀疏，有利于后续的相位滤波和解缠处理；
[0036]
(2)相较于分块法，本发明能够自适应选取阈值，适用于不同场景；且本发明由于考虑到整个距离向上的频谱，能够有效排除局部地形突变导致的频谱低峰，不易受局部区域影响，整个过程无需人工介入。
附图说明
[0037]
图1是本发明insar非线性去平地效应方法的流程图；
[0038]
图2是本发明中实施例生成的含噪干涉相位图；
[0039]
图3是本发明中实施例生成的真实平地相位图；
[0040]
图4是本发明中实施例生成的距离向频谱图；
[0041]
图5是本发明中实施例生成的归一化直方图；
[0042]
图6是本发明中实施例生成的非线性平地相位图；
[0043]
图7是本发明中实施例生成的去平地结果图；
[0044]
图8是本发明中实施例生成的去平地结果距离向频谱图；
[0045]
图9是传统频移法实施例生成的平地相位图；
[0046]
图10是传统频移法实施例生成的残余相位图；
[0047]
图11是传统频移法实施例生成的去平地结果距离向频谱图；
[0048]
图12是分块法实施例生成的平地相位图；
[0049]
图13是分块法实施例生成的残余相位图；
[0050]
图14是分块法实施例生成的去平地结果距离向频谱图。
具体实施方式
[0051]
下面将结合附图和实施例对本发明作进一步的详细说明。
[0052]
一种insar非线性去平地效应方法，流程图如图1所示，包括以下几个步骤：
[0053]
步骤一，合成孔径雷达对同一待测目标进行两次成像，得到两幅复图像并进行配准，对两幅配准后的sar复图像进行干涉，得到对应的干涉相位图；
[0054]
按表1所示雷达参数仿真得到小视角的海洋场景含噪干涉相位图如图2所示，其真实平地相位如图3所示。可以看出，平地相位左端条纹相对密集，右端条纹相对稀疏。
[0055]
表1雷达参数
[0056][0057]
在干涉相位图中，条纹的密集程度由占主要作用的平地相位以及次要作用的局部地形细节相位和噪声决定，因此可基于加性噪声模型，将含噪相位分解为平地相位与反应细节变化的残余相位
[0058][0059]
步骤二，对干涉复图像进行预滤波，对预滤波后的干涉复图像在距离向进行离散傅里叶变换并求和，从而得到对应于距离向的频谱幅值序列及频谱相位序列。
[0060]
根据离散傅里叶变换，复干涉相位可分解为多个线性相位的组合。考虑到平地相位一般为沿距离向的平行条纹，不失一般性给出沿距离向分布的复干涉相位分解公式：
[0061][0062]
其中，是含噪相位的复数形式，为虚数单位，m为坐标位置，g代表组成复干涉相位的矢量个数，ci为矢量的权值，代表对应矢量的距离向频率。
[0063]
在估计过程中，为防止图像中的局部区域地形的突变对估计平地相位造成干扰，in不应简单选择某一距离向数据作为样本，而应综合考虑整幅图像所有方位向上的距离向频率特征。因此需要对每组距离向上的频谱序列si进行求和，从而获得可以反映整体距离向分布的频谱sr以及对应的频谱幅值序列p和相位序列φ：
[0064][0065]
其中代表第i个距离向上的复干涉相位序列，u为频谱序列的坐标位置，
⊙
表示shur乘积算子。
[0066]
仿真图像的频谱幅值序列p如图4所示，可以看出，其谱峰相较于零频存在偏移，且具有一定宽度。
[0067]
而平地相位作为相位中的主要成分，可根据矢量权值的大小，选取前g个权值较大的矢量组合用来估计：
[0068][0069]
反映在sr中，矢量权值的大小即为p的大小，距离向频率特征则可由φ确定，后续的关键便在于如何确定g以估计平地相位。
[0070]
步骤三，根据频谱幅值序列构造归一化直方图，根据直方图分布自适应选取阈值。
[0071]
对于一般场景的干涉相位图，其频谱往往呈单一谱峰状，而当地形复杂时可能会出现多个次峰。传统频移法仅依赖于峰值点的位置，当谱峰较宽或存在次峰时容易丢失部分平地相位信息，导致估计的平地相位与真实平地相位无法拟合。由前面的分析可知，平地相位作为干涉相位图的主要影响因子，应在频谱中占主导地位，分布在低频区域，频谱幅值较高；而残余相位起次要作用，分布在高频区域，频谱幅值一般较低。因此，可以根据频谱幅值选取阈值进行区分。
[0072]
根据频谱幅值选取阈值的具体过程为：
[0073]
首先，通过构建归一化直方图的方法，将整个距离向上的频谱幅值进行归并，高频但幅值较低的噪声分量将会聚集在直方图低端形成谱峰，低频但幅值较高的平地相位则会比较零散地分布在高端。需要注意：为了确保平地相位分量与残余相位分量被有效区分，直方图的组数不应过小。在选择阈值时，若t过小，则筛选出的幅值序列中会包含地形细节、噪声等造成干扰；t过大，则仅有频谱峰值被保留，整个算法退化为传统频移法，导致平地相位损失。考虑到平地相位多分布于谱峰附近，相邻频点幅值差异较大，且同幅值水平的频点较少，映射在直方图中将呈一个个频率较低的凸点，此凸点在直方图中的位置即噪声频谱与平地相位频谱的分割点，因此本发明选择将直方图中对应幅值最小的极大值点作为t的选择标准。
[0074]
为计算直方图极大值点的位置，需要对直方图序列h作二阶差分：
[0075]
d(k)＝[h(k+1)-h(k)]-[h(k)-h(k-1)],k＝2,3,...,n-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0076]
其中d为二阶差分序列，n为直方图组数，通过判断d(k)≥0即可找到直方图中的所有极大值点。
[0077]
需要注意的是，由于噪声的随机性，当直方图组数较多，且噪声较严重时，幅值最小的极大值点可能出现在直方图左端以噪声为主的区域干扰阈值的选择，若选此点作为阈值，估计出的平地相位中将包含大量噪声成分。为避免此现象，可根据该点的频率水平进一步判断是否适合作为阈值：
[0078]
对于噪声的频率水平，由于其代表了整个直方图的主要部分，通常即在直方图均匀分布下对应的频率以上，相反平地相位的频率水平则在该频率以下。综上所述，只有当幅值最小的极大值点其频率小于或等于时才将其作为阈值，否则应继续判断剩余极大值点是否满足阈值条件并将其中幅值最小的点作为阈值。至此本发明可根据直方图的组数自适应选择出合适的阈值。
[0079]
仿真图像的距离向频谱幅值序列所构建的归一化直方图如图5所示，可以看出频谱幅值主要集中在直方图左侧，且呈下落式分布，右端则分布着多个零散的凸点，即为所需寻找的极大值点。
[0080]
步骤四，针对于频谱幅值序列中低于阈值的部分，取其均值作为背景噪声功率估计值，并将该部分置零；对于高于阈值的部分，则减去背景噪声功率估计值以进一步抑制噪声的影响，得到处理后的频谱幅值序列：
[0081]
处理后的频谱幅值序列为：
[0082][0083]
其中，noise＝mean(p(u)＜t)为背景噪声估计值，mean()表示取均值操作。需要注意的是，由于噪声均匀分布在整个频谱上，因此对大于阈值的幅值部分，也应当将该部分噪声减去，减小噪声产生的影响。
[0084]
步骤五，将处理后的频谱幅值序列与频谱相位序列重新结合，并作离散傅里叶逆变换后，取相角即可得到非线性平地相位序列，将其扩展为原始尺寸，并从原干涉相位图中减去即为最终的去平地结果。
[0085]
针对处理后的频谱幅值序列，还需要进一步将其还原为平地相位形式，将其重新与相位序列结合后，即可得到以平地相位为主所对应的频谱，对其进行离散傅里叶逆变换后，得到对应的距离向平地复相位形式，取相位主值即得到平地相位序列：
[0086]
φ
p
＝arg(ifft(p
′⊙
φ))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)
[0087]
为获得与原干涉相位图对应的平地相位图，还需进一步将其该序列扩展，由于平地相位表现为在距离向上平行的干涉条纹，因此可将其与全1列矩阵相乘：
[0088][0089]
原干涉相位图与平地相位图相减并缠绕后，即得到最终的去平地结果。
[0090]
对图2采用本发明提出的方法进行去平地效应处理，得到的平地相位图如图6所示，去平地结果如图7所示，以及去平地后距离向频谱如图8所示。
[0091]
为了验证本发明提出的非线性去平地效应方法的有效性，分别采用传统频移法与分块法对同样的干涉相位图进行去平地处理，其结果分别如图9-11和图12-14所示。
[0092]
将三种方法得到的平地相位与真实平地相位进行对比，并引入相位标准差(phase standard deviation,psd)指标进行比较。psd以真实相位作为参考，其值越小，说明估计结果与真值的偏差越小，可用下式计算：
[0093][0094]
其中，p,q为相位图的尺寸，为估计的平地相位，为真实平地相位，w[
·
]用于对相位进行缠绕，取[-π,π]的相位主值。
[0095]
统计结果如表2所示：
[0096]
表2仿真数据评估结果
[0097]
[0098]
从表2以及图6-14所示三种方法的去平地结果可以看出：
[0099]
利用传统频移法估计的平地条纹为分布均匀的周期性条纹，去平地后残余相位的左侧仍存在明显的条纹变化，距离向频谱的峰值被搬移到零频附近；
[0100]
分块法估计的平地相位条纹有了明显的粗细变化，使得残余相位图相较传统频移法颜色变化更加缓慢，谱峰相比传统频移法更窄；
[0101]
相较于现有两种方法，本发明提出的方法估计的平地相位条纹从左至右宽度缓慢增大，残余相位中颜色变化更加平缓，仅在局部区域表现出地形细节导致的相位变化，谱峰相较于前两种方法，更贴合零频。
[0102]
另一方面，从psd值看，本发明提出的非线性去平地效应方法psd最小，证明与真实平地相位拟合的最好。
[0103]
综合以上分析可知，本发明提出的insar非线性去平地效应方法针对不同场景的适应性更强，去平地结果相较现有方法有显著提升。
[0104]
以上所述，仅为本发明的部分实施示例，并非用于限定本发明的保护范围。任何不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

技术特征：
1.一种insar非线性去平地效应方法，其特征在于，包括以下几个步骤：步骤一，合成孔径雷达对同一待测目标进行两次成像，得到两幅复图像并进行配准，对两幅配准后的sar复图像进行干涉，得到对应的干涉相位图；步骤二，对干涉复图像进行预滤波，并对预滤波后的干涉复图像距离向数据进行离散傅里叶变换并求和，从而得到对应于距离向的频谱幅值序列及相位序列；预滤波后的干涉复图像的复干涉相位分解公式：其中，是含噪相位的复数形式，为虚数单位，m为坐标位置，g代表组成复干涉相位的矢量个数，c
i
为第i个矢量的权值，代表对应第i个矢量的距离向频率；然后，在每组距离向上对复干涉相位序列进行离散傅里叶变换并求和，得到频谱幅值序列s
i
的求和结果，获得频谱幅值序列p和相位序列φ：其中代表第i个距离向上的复干涉相位序列，u为频谱序列的坐标位置，
⊙
表示shur乘积算子；步骤三，根据频谱幅值序列构造归一化直方图，根据直方图分布自适应选取阈值；根据直方图自适应选择阈值的实现过程为：步骤301，依据频谱幅值序列构造归一化直方图，对归一化直方图序列h作二阶差分运算，并寻找其中所有的极大值点；步骤302，从左端起依次判断极大值点处对应的频率与直方图均匀分布下的频率之间的关系，若则代表噪声的频率水平；若则代表平地相位的频率水平，n为直方图组数；步骤303，在满足平地相位频率条件的极大值点中选择对应直方图中幅值最小的极大值点作为阈值t；步骤四，针对于频谱幅值序列中低于阈值的部分，取其均值作为背景噪声功率估计值，并将该部分置零；对于高于阈值的部分，则减去背景噪声功率估计值以进一步抑制噪声的影响，经过以上操作后得到处理后的频谱幅值序列；步骤五，将处理后的频谱幅值序列与频谱相位序列重新结合，并作离散傅里叶逆变换后，取相角即可得到非线性平地相位序列，并将其扩展为原始尺寸；步骤六，从干涉相位图中减去扩展为原始尺寸后的非线性平地相位序列，得到最终的去平地结果。2.根据权利要求1所述的一种insar非线性去平地效应方法，其特征在于，步骤一中，所述的干涉相位图中，将含噪相位分解为平地相位与反应细节变化的残余相位与反应细节变化的残余相位3.根据权利要求1所述的一种insar非线性去平地效应方法，其特征在于，步骤三中，所
述的对归一化直方图序列h作二阶差分运算：d(k)＝[h(k+1)-h(k)]-[h(k)-h(k-1)],k＝2,3,...,n-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)其中d为二阶差分序列，通过判断d(k)≥0即可找到直方图中的所有极大值点。4.根据权利要求1所述的一种insar非线性去平地效应方法，其特征在于，步骤四中，所述的背景噪声功率估计值为：noise＝mean(p(u)＜t)，mean()表示取均值操作；处理后的频谱幅值序列为：5.根据权利要求1所述的一种insar非线性去平地效应方法，其特征在于，步骤五中，所述的非线性平地相位序列为：φ
p
＝arg(ifft(p
′⊙
φ))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)非线性平地相位序列扩展为原始尺寸的公式为：

技术总结
本发明公开了一种InSAR非线性去平地效应方法，涉及合成孔径雷达干涉处理领域。合成孔径雷达对待测目标进行检测，得到干涉相位图和干涉复图像，对干涉复图像进行预滤波，并进行离散傅里叶变换并求和，得到对应于距离向的频谱幅值序列及相位序列。然后频谱幅值序列构造归一化直方图，根据直方图分布自适应选取阈值，在频谱幅值序列高于阈值的部分减去低于阈值的部分的均值，得到处理后的频谱幅值序列，并与频谱相位序列重新结合，作离散傅里叶逆变换，取相角，得到非线性平地相位序列，并将其扩展为原始尺寸。从干涉相位图中减去扩展为原始尺寸后的非线性平地相位序列，得到最终的去平地结果。本发明去平地结果更准确，适用于不同场景。场景。场景。


技术研发人员：徐华平 梁若斌 李硕
受保护的技术使用者：北京航空航天大学
技术研发日：2023.07.04
技术公布日：2023/9/25