一种针对矿区大梯度形变区域的相位解缠方法及系统

1.本发明属于合成孔径雷达干涉数据处理领域，具体涉及一种针对矿区大梯度形变区域的相位解缠方法及系统。


背景技术：

2.合成孔径雷达干涉测量(interferometric syntheticaperture radar,insar)是一种用于获取高精度数字高程模型(digital elevation model,dem)、监测地表形变与预警滑坡的关键技术。相位解缠是insar数据处理中的关键步骤。大梯度变化区域相位解缠结果精度将直接影响最终地表形变监测精度。然而在实际情况中，大梯度变化区域相邻像素的绝对相位跳变很可能会大于π，常用的单基线相位解缠方法受到相位连续性假设的条件限制，相位解缠结果精度不甚理想，从而影响到insar获取最终产品的精度。相位解缠一直是insar数据处理流程中的研究热点，常用的相位解缠方法主要可以分为两类，基于路径的相位解缠与基于最小范数的相位解缠方法。前者以枝切法、质量图法和最小不连续法为代表，该类方法通过计算残差点、优化积分路径等手段来减小误差的传播。但，在干涉图的低质量区域，该类方法会产生误差传递现象，导致累积误差，甚至出现孤岛现象，严重影响到相位解缠结果精度。后者以最小费用流(minimum-cost flow，mcf)方法、统计费用流(statistical costnetwork flow,snaphu)及l
p
范数法为代表，该类方法通过使真实相位梯度与相位梯度估计值之间差异最小化来获取最终解缠结果。然而，该类方法易产生全局误差传递且鲁棒性较差。矿区的大梯度区域由于受到相位连续性假设的限制，传统的相位解缠方法往往难以得到较好的解缠结果，极大影响了矿区形变监测的精度。针对上述问题，多基线相位解缠方法被提出，多基线相位解缠绕方法可以利用干涉相位图的基线多样性来实现相位解缠，该类方法可以摆脱相位连续性假设的束缚，在矿区大梯度变形区域也可以获得较好的解缠结果。然而，多基线相位解缠方法所需要的数据更多且其存在噪声鲁棒性较差的问题。因此，如何通过单基线相位解缠方法更为精确的获取大梯度变化区域的相位解缠结果显得尤为重要。


技术实现要素：

3.针对现有技术的不足，本发明提出了一种针对矿区大梯度形变区域的相位解缠方法及系统，能够解决大梯度变化区域解缠精度低的问题，有效的提高相位解缠精度。
4.为实现上述目的，本发明提供了如下方案：
5.一种针对矿区大梯度形变区域的相位解缠方法，包括以下步骤：
6.s1：对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，对所述相位梯度进行相位差分求得二阶差分结果，对所述二阶差分结果进行fft2变换；
7.s2：通过对原干涉相位图二阶差分的fft2变换结果进行泊松方程的求解，对求解结果进行二维傅里叶反变换ifft2得到原干涉相位图的预解缠结果；
8.s3：通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条
纹，得到噪声相位图；
9.s4：通过对所述噪声相位图进行l1范数相位解缠，得到噪声相位中残留的有效相位信息；
10.s5：将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值。
11.优选的，所述s1中，
12.对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的方法为：
[0013][0014][0015]
式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；
[0016]
对所述相位梯度在进行相位差分，得到二阶差分结果的方法为：
[0017][0018][0019][0020]
式中，与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和；
[0021]
对所述二阶差分结果进行fft2变换的方法为：
[0022][0023]
式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。
[0024]
优选的，所述s2中，
[0025]
对原干涉相位图二阶差分的fft2变化结果进行泊松方程的求解的方法为：
[0026][0027]
式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点fft2变换后结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位；
[0028]
对求解结果进行ifft2变换得到原干涉相位图的预解缠结果的方法为：
[0029][0030]
式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。
[0031]
优选的，所述s3中，通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图的方法为：
[0032][0033]
式中，为原干涉相位图，为预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。
[0034]
优选的，所述s4中，l1范数的表达式为：
[0035][0036]
式中，ψ
noise
为噪声中残留有效信息绝对相位，与分别为残留有效信息距离向与方位向上的缠绕相位梯度。
[0037]
优选的，所述s5中，将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值的方法为：
[0038][0039]
式中，ψ
pu
为最终的解缠相位。
[0040]
本发明还提供了一种针对矿区大梯度形变区域的相位解缠系统，包括：差分模块、预解缠模块、复共轭模块、相位解缠模块和共轭相乘模块；
[0041]
所述差分模块用于对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，对所述相位梯度进行相位差分求得二阶差分结果，对所述二阶差分结果进行fft2变换；
[0042]
所述预解缠模块用于通过对原干涉相位图二阶差分的fft2变换结果进行泊松方程的求解，对求解结果进行二维傅里叶反变换ifft2得到原干涉相位图的预解缠结果；
[0043]
所述复共轭模块用于通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图；
[0044]
所述相位解缠模块用于通过对所述噪声相位图进行l1范数相位解缠，得到噪声相位中残留的有效相位信息；
[0045]
所述共轭相乘模块用于将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值。
[0046]
优选的，所述差分模块中，
[0047]
对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的过程为：
[0048][0049][0050]
式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；
[0051]
对所述相位梯度在进行相位差分，得到二阶差分结果的过程为：
[0052][0053][0054][0055]
式中，与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和；
[0056]
对所述二阶差分结果进行fft2变换的过程为：
[0057][0058]
式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。
[0059]
优选的，所述预解缠模块中，
[0060]
对原干涉相位图二阶差分的fft2变化结果进行泊松方程的求解的过程为：
[0061]
[0062]
式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点fft2变换后结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位；
[0063]
对求解结果进行ifft2变换得到原干涉相位图的预解缠结果的过程为：
[0064][0065]
式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。
[0066]
优选的，所述复共轭模块中，通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图的过程为：
[0067][0068]
式中，为原干涉相位图，为预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。
[0069]
与现有技术相比，本发明的有益效果为：
[0070]
本发明采用一种针对矿区大梯度形变区域的相位解缠方法及系统，通过对原始干涉图进行二阶差分，再对二阶差分结果进行fft2变换，通过泊松方程对变换结果进行求解，在将结果进行反变换，得到预解缠结果；将预解缠结果与原始干涉图进行复共轭相乘，在稀释条纹的同时得到噪声相位图；再对噪声相位图进行l1范数解缠，得到噪声相位图中残留的有效相位信息，最后将提取出的有效相位信息与预解缠结果进行共轭相乘，得到最终的解缠结果。与已有的其它常规相位解缠方法相比，本发明可以从大梯度变化区域也可以获得较其它方法更大面积有效解缠结果，且该方法具有更好的解缠模型鲁棒性，有效的提高了最终insar产品的精度。
附图说明
[0071]
为了更清楚地说明本发明的技术方案，下面对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0072]
图1是本发明实施例的一种针对矿区大梯度形变区域的相位解缠方法流程图；
[0073]
图2为本发明实施例所用的两个模拟干涉相位的参考解缠相位示意图；
[0074]
图3为本发明实施例所用的两个模拟干涉相位数据示意图；
[0075]
图4为本发明实施例所用的两个模拟干涉相位加噪数据示意图；
[0076]
图5为本发明实施例两个模拟干涉相位加噪数据枝切法相位解缠结果及其对应的解缠结果误差示意图；
[0077]
图6为本发明实施例两个模拟干涉相位加噪数据mcf相位解缠结果及其对应的解缠结果误差示意图；
[0078]
图7为本发明实施例两个模拟干涉相位加噪数据snaphu相位解缠结果及其对应的解缠结果误差示意图；
[0079]
图8为本发明实施例两个模拟干涉相位加噪数据基于fft2的l1范数相位解缠方法的相位解缠结果及其对应的解缠结果误差示意图；
[0080]
图9为本发明实施例的实验数据矿区的google地图范围示意图；
[0081]
图10为本发明实施例的实验数据gf-3数据示意图；
[0082]
图11为本发明实施例的三个实验数据干涉相位示意图；
[0083]
图12为本发明实施例的三个实验数据枝切法相位解缠结果；
[0084]
图13为本发明实施例的三个实验数据mcf相位解缠结果；
[0085]
图14为本发明实施例的三个实验数据snaphu相位解缠结果；
[0086]
图15为本发明实施例的三个实验数据基于fft2的l1范数相位解缠方法的相位解缠结果。
具体实施方式
[0087]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0088]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
[0089]
实施例一
[0090]
如图1所示，本发明提出一种针对矿区大梯度形变区域的相位解缠方法，具体包括以下步骤：
[0091]
s1：对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，再对相位梯度进行差分求得二阶差分结果，对二阶差分结果进行fft2变换。二阶差分具体公式如下：
[0092][0093][0094][0095][0096][0097]
式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、
ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和。fft2变换表达式为：
[0098][0099]
式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。
[0100]
s2、通过对干涉的二阶差分fft2变换结果进行泊松方程的求解以及二维傅里叶反变换(ifft2)得到影像的预解缠结果，泊松方程求解表达式如下：
[0101][0102]
式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点二阶差分梯度fft2变换结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位。ifft2变换具体表达式为：
[0103][0104]
式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。
[0105]
s3、通过原影像与s2预解缠结果进行复共轭相乘以稀释干涉条纹，且得到噪声相位图，具体表达式如下：
[0106][0107]
式中，为原干涉相位图，为s2预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。
[0108]
s4、通过对噪声相位图进行l1范数相位解缠得到噪声相位中残留的有效相位信息。理论上说，相邻像素之间的缠绕相位差与绝对相位差是相等的，但由于噪声等因素的影响。缠绕相位差与绝对相位差往往不相等。l1最小范数相位解缠是通过使相邻像素的缠绕相位差与绝对相位差之间的差距达到最小来实现相位解缠。具体表达式如下：
[0109][0110]
式中，ψ
noise
为噪声中残留有效信息绝对相位，与分别为残留有效信息距离向与方位向上的缠绕相位梯度；
[0111]
s5、将s2预解缠结果与s4求得的残留有效相位进行共轭相乘，得到最终的解缠相
位值，具体表达式如下：
[0112][0113]
式中，ψ
pu
为最终的解缠相位。
[0114]
为了验证本发明的技术效果，分别采用枝切法、mcf、snaphu以及本发明的相位解缠方法对相同的干涉图进行解缠处理实验。实验所用模拟参考解缠相位如图2所示，模拟干涉相位数据如图3所示，模拟干涉相位加噪数据如图4所示。模拟干涉相位枝切法解缠结果及其对应的解缠结果误差如图5所示，mcf解缠结果及其对应的解缠结果误差如图6所示，snaphu解缠结果及其对应的解缠结果误差如图7所示，本发明解缠结果及其解缠结果误差如图8所示。为了定量描述相位解缠的质量，我们计算各种相位解缠方法获得误差图的均方根误差，其中，对于模拟数据a，枝切法的解缠结果均方根误差为0.6918rad；mcf方法的解缠结果均方根误差为0.4340rad；snaphu方法的解缠结果均方根误差为0.3379；本发明的基于fft2的l1范数相位解缠方法的解缠结果均方根误差为0.3215rad。对于模拟数据b，枝切法的解缠结果均方根误差为0.7662rad；mcf方法的解缠结果的均方根误差为0.4029rad；snaphu方法的解缠结果均方根误差为0.4308rad；而本发明的基于fft2的l1范数相位解缠方法的解缠结果均方根误差为0.1010rad。结果表明，本发明的解缠结果明显优于枝切法、mcf与snaphu相位解缠方法。同时，实验采用枝切法、mcf、snaphu与本发明的相位解缠方法对山西大同矿区gf-3真实数据进行了解缠处理。实验所用数据的google地图范围如图9所示；实验所用山西大同矿区gf-3数据如图10所示，从图中可以看出，该地区存在着较多的大梯度变形区域，且在变形中心存在着相位不连续现象；实验所用的三个真实数据的干涉相位图如图11所示，从图中可以看出，真实数据中存在着高噪声与大梯度变形；实验数据的枝切法相位解缠结果如图12所示；实验数据mcf相位解缠结果如图13所示；实验数据snaphu相位解缠结果如图14所示；实验数据本发明相位解缠方法的相位解缠结果如图15所示。区域1存在明显的相位不连续现象，从区域1的枝切法、mcf、snaphu相位解缠结果来看，枝切法与snaphu出现明显的解缠误差，mcf出现明显的误差传递现象，而本发明解缠结果较为理想。从区域2的枝切法、mcf与snaphu的相位解缠结果来看，沉陷中心出现了明显的解缠误差，不仅如此，mcf在沉降变形的下部产生了明显的误差传递现象。然而，本发明所提方法不仅能获得形变中心的解缠结果，而且可以避免解缠误差的传递。从区域3的枝切法、mcf与snaphu相位解缠结果来看，形变中心出现了明显的解缠误差。然而，本发明的相位解缠方法仍然能够获得较为理想的解缠结果。综上结果表明，本发明的相位解缠方法能够获得更高精度的相位解缠结果。
[0115]
综上所述，本发明提出了一种针对矿区大梯度形变区域的相位解缠方法。与现有的用于沉降中心大梯度变化的解缠方法不同，该方法不依赖于外部先验数据、形变模型和其他技术。相反，它在算法上能够进一步改善高噪声与大梯度形变的相位解缠问题。首先，对原始的干涉相位图做二阶差分，再将差分结果运用fft2变换转换到频率域，在频率域中去除图像中所存在的噪声；其次，对变换结果进行泊松方程的求解，得到原始干涉图在频率域上的预解缠结果，再对结果进行ifft2变化得到空间域预解缠结果；此时，通过原始干涉图与预解缠结果的复共轭相乘，该步骤在系数干涉条纹的同时获得了原干涉图的噪声相位图；然后，通过对噪声相位图进行l1范数的相位解缠，可以得到噪声相位图中残留的有效相位信息；最后，将提取出的有效相位信息与预解缠结果进行共轭相乘，即可得到最终的相位
解缠结果。通过模拟数据与真实数据实验，结果表明，本发明基于fft2的l1范数相位解缠方法可以有效地解决高噪声与大梯度形变的相位解缠问题。
[0116]
实施例二
[0117]
本发明还提供了一种针对矿区大梯度形变区域的相位解缠系统，包括：差分模块、预解缠模块、复共轭模块、相位解缠模块和共轭相乘模块；
[0118]
所述差分模块用于对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，对相位梯度进行相位差分求得二阶差分结果，对二阶差分结果进行fft2变换；
[0119]
所述预解缠模块用于通过对原干涉相位图二阶差分的fft2变换结果进行泊松方程的求解，对求解结果进行二维傅里叶反变换ifft2得到原干涉相位图的预解缠结果；
[0120]
所述复共轭模块用于通过将原干涉相位图与预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图；
[0121]
所述相位解缠模块用于通过对噪声相位图进行l1范数相位解缠，得到噪声相位中残留的有效相位信息；
[0122]
所述共轭相乘模块用于将预解缠结果与有效相位信息进行共轭相乘，得到最终的解缠相位值。
[0123]
在本实施例中，所述差分模块中，
[0124]
对原干涉相位图进行一阶相位差分，得到相位梯度的过程为：
[0125][0126][0127]
式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；
[0128]
对相位梯度在进行相位差分，得到二阶差分结果的过程为：
[0129][0130][0131][0132]
式中，与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和；
[0133]
对二阶差分结果进行fft2变换的过程为：
[0134][0135]
式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。
[0136]
在本实施例中，所述预解缠模块中，
[0137]
对原干涉相位图二阶差分的fft2变化结果进行泊松方程的求解的过程为：
[0138][0139]
式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点fft2变换后结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位；
[0140]
对求解结果进行ifft2变换得到原干涉相位图的预解缠结果的过程为：
[0141][0142]
式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。
[0143]
在本实施例中，所述复共轭模块中，通过将原干涉相位图与预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图的过程为：
[0144][0145]
式中，为原干涉相位图，为预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。
[0146]
在本实施例中，所述相位解缠模块中，l1范数的表达式为：
[0147][0148]
式中，ψ
noise
为噪声中残留有效信息绝对相位，与分别为残留有效信息距离向与方位向上的缠绕相位梯度。
[0149]
在本实施例中，所述共轭相乘模块中，将预解缠结果与有效相位信息进行共轭相乘，得到最终的解缠相位值的过程为：
[0150][0151]
式中，ψ
pu
为最终的解缠相位。
[0152]
以上所述的实施例仅是对本发明优选方式进行的描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。

技术特征：
1.一种针对矿区大梯度形变区域的相位解缠方法，其特征在于，包括以下步骤：s1：对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，对所述相位梯度进行相位差分求得二阶差分结果，对所述二阶差分结果进行fft2变换；s2：通过对原干涉相位图二阶差分的fft2变换结果进行泊松方程的求解，对求解结果进行二维傅里叶反变换ifft2得到原干涉相位图的预解缠结果；s3：通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图；s4：通过对所述噪声相位图进行l1范数相位解缠，得到噪声相位中残留的有效相位信息；s5：将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值。2.根据权利要求1所述的针对矿区大梯度形变区域的相位解缠方法，其特征在于，所述s1中，对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的方法为：对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的方法为：式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；对所述相位梯度在进行相位差分，得到二阶差分结果的方法为：对所述相位梯度在进行相位差分，得到二阶差分结果的方法为：对所述相位梯度在进行相位差分，得到二阶差分结果的方法为：式中，与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和；对所述二阶差分结果进行fft2变换的方法为：式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。
3.根据权利要求1所述的针对矿区大梯度形变区域的相位解缠方法，其特征在于，所述s2中，对原干涉相位图二阶差分的fft2变化结果进行泊松方程的求解的方法为：式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点fft2变换后结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位；对求解结果进行ifft2变换得到原干涉相位图的预解缠结果的方法为：式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。4.根据权利要求1所述的针对矿区大梯度形变区域的相位解缠方法，其特征在于，所述s3中，通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图的方法为：式中，为原干涉相位图，为预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。5.根据权利要求1所述的针对矿区大梯度形变区域的相位解缠方法，其特征在于，所述s4中，l1范数的表达式为：式中，ψ
noise
为噪声中残留有效信息绝对相位，与分别为残留有效信息距离向与方位向上的缠绕相位梯度。6.根据权利要求5所述的针对矿区大梯度形变区域的相位解缠方法，其特征在于，所述s5中，将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值的方法为：式中，ψ
pu
为最终的解缠相位。7.一种针对矿区大梯度形变区域的相位解缠系统，其特征在于，包括：差分模块、预解缠模块、复共轭模块、相位解缠模块和共轭相乘模块；所述差分模块用于对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度，对所述相位梯度进行相位差分求得二阶差分结果，对所述二阶差分结果进行fft2
变换；所述预解缠模块用于通过对原干涉相位图二阶差分的fft2变换结果进行泊松方程的求解，对求解结果进行二维傅里叶反变换ifft2得到原干涉相位图的预解缠结果；所述复共轭模块用于通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图；所述相位解缠模块用于通过对所述噪声相位图进行l1范数相位解缠，得到噪声相位中残留的有效相位信息；所述共轭相乘模块用于将所述预解缠结果与所述有效相位信息进行共轭相乘，得到最终的解缠相位值。8.根据权利要求7所述的针对矿区大梯度形变区域的相位解缠系统，其特征在于，所述差分模块中，对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的过程为：对矿区大梯度形变区域的原干涉相位图进行一阶相位差分，得到相位梯度的过程为：式中，和分别为(i,j)像素点在距离向与方位向上的相位梯度估计值，分别为(i+1,j)、(i,j)、(i,j+1)像素点的缠绕相位值，ψ
i+1,j
、ψ
i,j
、ψ
i,j+1
分别为(i+1,j)、(i,j)、(i,j+1)像素点的绝对相位值；对所述相位梯度在进行相位差分，得到二阶差分结果的过程为：对所述相位梯度在进行相位差分，得到二阶差分结果的过程为：对所述相位梯度在进行相位差分，得到二阶差分结果的过程为：式中，与分别为干涉相位图距离向与方位向上的二阶差分值，为干涉图(i,j)像素点的二阶差分梯度和；对所述二阶差分结果进行fft2变换的过程为：式中，l为干涉图二阶差分梯度fft2变换结果，fft2(
·
)为二维快速傅里叶变换操作，为干涉图二阶差分梯度图。9.根据权利要求7所述的针对矿区大梯度形变区域的相位解缠系统，其特征在于，所述
预解缠模块中，对原干涉相位图二阶差分的fft2变化结果进行泊松方程的求解的过程为：式中，m与n分别为干涉图的行数与列数，l
i,j
为(i,j)像素点fft2变换后结果，δ
i,j
为(i,j)像素点泊松方程求解后在频率域上的解缠相位；对求解结果进行ifft2变换得到原干涉相位图的预解缠结果的过程为：式中，为预解缠结果，ifft2(
·
)为二维快速傅里叶反变换，δ为频率域预解缠结果。10.根据权利要求7所述的针对矿区大梯度形变区域的相位解缠系统，其特征在于，所述复共轭模块中，通过将所述原干涉相位图与所述预解缠结果进行复共轭相乘以稀释干涉条纹，得到噪声相位图的过程为：式中，为原干涉相位图，为预解缠结果，代表乘法操作，conj(
·
)代表取共轭复数，为噪声相位图。

技术总结
本发明公开了一种针对矿区大梯度形变区域的相位解缠方法及系统，实现在解缠困难区域获取高精度解缠结果。具体包括：基于二维快速傅里叶变换(FFT2)将影像(干涉相位图)转换到频率域，通过对影像的二阶差分结果进行泊松方程的求解以及二维傅里叶反变换(iFFT2)得到影像的预解缠结果；通过原影像与预解缠结果进行复共轭相乘以稀释干涉条纹且得到噪声相位图；通过对噪声相位图进行L1范数相位解缠得到噪声相位中残留的有效相位信息；最后，将预解缠结果与残留相位进行共轭相乘，得到最终的解缠结果。本发明能解决大梯度变化区域的解缠精度低的问题，有效提高解缠精度。有效提高解缠精度。有效提高解缠精度。


技术研发人员：高延东 郑南山 杨化超 张秋昭 李世金 卞和方 张书毕 张帝 阎超
受保护的技术使用者：中国矿业大学
技术研发日：2023.06.30
技术公布日：2023/10/15