CN115808688A - 融合随机介质模型的沉陷区差分干涉影像图相位解缠方法 - Google Patents

Abstract

本发明提供了一种融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，属于图像处理技术领域。针对矿区开采沉陷的特点，从开采沉陷的机理及规律出发，基于随机介质模型推演工作面开采引起的地表变形估值，然后以地表变形的估值来修正相位解缠结果，从而准确获取地表任意点的解缠相位；该方法充分考虑了开采沉陷机理，利用由于地下煤炭开采而引起地表变形的估值来指导相位解缠结果和提高相位解缠精度，所提方法弥补了现有解缠算法在矿区地表快速变形区域无能为力的缺陷。其能够完成对矿区地表快速变形区域干涉图相位的解缠，与流行的最小费用流法相比，具有更为合理的解缠结果。

Description

融合随机介质模型的沉陷区差分干涉影像图相位解缠方法

技术领域

本发明涉及一种融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，属于图像处理技术领域。

背景技术

合成孔径雷达干涉测量(InSAR)技术凭借其地面分辨率高、覆盖范围广以及基于面观测的独特优势，在构建地表数字高程模型和监测地表形变领域已经得到了广泛的应用。然而，干涉图中的相位信息仅是以2π为模的相位值，为了获得绝对相位，必须进行相位解缠。相位解缠的准确性直接影响到所构建DEM以及所监测地表变形的精度，因此，它在SAR干涉处理中是一个极其重要的环节。

现有的相位解缠算法主要有两类：路径跟踪法、最小范数法。前者是通过选择积分路径对相邻像元相位差值进行积分，探测2π相位的整周数，该类算法中最为经典的方法是Goldstein枝切法，此外还有最小生成树法、最小费用流法等；后者是以最小化缠绕相位梯度值与解缠相位梯度值之差为准则，将相位解缠转化为求解数学中最小范数的全局优化问题，典型代表为最小二乘法。后来也有许多学者针对某些特定的问题或实现特定的目标提出了具有不同思路的相位解缠算法，许军毅等人针对大尺度、超宽带P波段SAR干涉相位图提出了基于区域划分的相位解缠思路；黄海峰等人针对地形变化剧烈的区域提出一种快速的多基线、多频域相位解缠算法；刘万丽等人针对矿区高噪声条件的干涉相同图提出了基于卡尔曼滤波的相位解缠方法。此外，Danudirdjo等人根据干涉影像投影缩短效应提出了合成孔径雷达各向异性相位解缠思路。

由于InSAR技术在监测地表形变中具有大范围、高精度、持续动态连续监测的优势，其也被广泛应用于矿区开采沉陷的监测。然而，由煤炭开采引起的地表沉陷不同于其它地表变形，具有变形速率快、沉降量级大的特点，在干涉相位图中表现为相邻像元间相位梯度大，条纹密集度高。相邻像元间相位差的不连续，极大的增加了相位解缠的难度，使得现有相位解缠算法难以实现对矿区地表快速变形区域干涉相位图的相位解缠。

发明内容

为解决现有技术中存在的问题，本发明提供一种以随机介质模型为基础的融合沉陷预计估值的相位解缠方法，该方法从开采沉陷的机理及规律出发，首先基于随机介质模型推演工作面开采引起的地表变形估值，然后以地表变形的估值来修正相位解缠结果，以达到准确获取解缠相位的目的。

为实现上述目标，本发明提出了一种融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，针对矿区开采沉陷的特点，从开采沉陷的机理及规律出发，基于随机介质模型推演工作面开采引起的地表变形估值，然后以地表变形的估值来修正相位解缠结果，从而准确获取地表任意点的解缠相位；

其步骤如下：

S1、对目标矿区地表SAR影像进行差分干涉处理，得到目标矿区地表的含有形变信息的差分干涉条纹图，差分干涉处理步骤包括：配准、重采样、干涉、去地形相位；

S2、利用现有常规相位解缠方法：枝切法、最小费用流法，对目标矿区地表的差分干涉条纹图进行解缠处理，得到目标矿区地表的初始解缠相位图；

S3、根据井下开采信息及目标矿区地质采矿情况，计算矿区地下资源开采所引起地表变形沉陷的估值；以目标矿区地表SAR影像分辨率为基准，对与SAR影像分辨率不一致的地表变形估值进行重采样使两者分辨率一致，重采样采用双线性内插法；地表变形估值的坐标系统一到SAR影像坐标系下，并生成地表变形图，之后计算地表变形图中相邻的地表变形估值点间的梯度；

S4、基于相位干涉的基本原理，以λ*u/2为基准，对地表变形图中的区域进行分类，将形变梯度小于λ*u/2的区域标定为准真实形变区，其余标定为待修订形变区，其中λ为SAR影像传感器的雷达波长，u为SAR影像传感器的影像分辨率；

S5、利用SAR的雷达参数，将完成分类的两类区域的地表变形信息分别转化为相位信息，并分别标记为准真实相位区和待修订相位区；

S6、将步骤S5得到的准真实相位区与步骤S2得到的初始解缠相位图两者之间坐标范围相同的区域定义为公共部分，采用最小二乘法对公共部分进行相位匹配处理，将准真实相位修正到初始解缠相位，得到修正多项式f(i，j)；

S7、利用修正多项式f(i，j)对待修订相位区的所有相位信息进行叠加，从而计算出修正后相邻像元间相位整周数；

S8、将步骤S7中获得的相位整周数与步骤S2得到的初始解缠相位进行融合，即可得到沉陷区差分干涉影像图的最终解缠相位。

进一步，相位解缠中，为得到实际的地表形变量，必须清除2π整数倍的相位模糊性，恢复到真实相位，因此真实相位与缠绕相位之间满足以下关系：

即从缠绕相位

中寻求最佳2π周期整数倍n(i,j)从而实现相位解缠；

式中，

为缠绕相位值；φ(i,j)为真实相位值；n(i,j)为整数；干涉相位图的像素为M×N，像元中直接观测到的相位值处于(-π,π)之间。

进一步，步骤3中地表变形值的具体计算过程为：利用随机介质模型，通过设定预计拐点，拐点内的地表变形值直接使用随机介质模型进行计算，即公式(2)，拐点以外地表下沉值不能直接使用随机介质模型的计算结果，需要对模型修正后使用，修正后的结果即为公式(3)，采用分段函数的方式，预计工作面开采后地表变形值，进而根据变形预计值探测相邻像元间的相位整周数：

拐点以内地表下沉值：

随机介质模型满足a)开采引起的各方向的移动与方向无关；b)地下大工作面开采引起的地表变形是多个小工作面开采引起地表变形的线性叠加；二维平面内单元开采引起地表点x的下沉值W_e(x)表示为：

式中，r＝H₀/tanβ为常数，称为主要影响半径，H₀是平均开采深度，tanβ是主要影响角正切，式中x表示地面任意点与拐点间的水平距离；

拐点以外地表下沉值：

应用随机介质模型进行开采沉陷预计时，拐点以外的部分存在收敛过快的缺陷，影响下沉曲线收敛性的主要因素是常数r，通过改变r值，即可实现对下沉曲线的控制；拐点以外r的取值与拐点内r的取值之比近似为一常数k，故拐点以外，在二维平面内单元开采引起地表点的下沉值表示为：

式中，k为常数，定义为沉陷曲线收敛系数；

三维平面任意点地表下沉值的表达：

地表点A(x,y)的下沉是由于煤层沿走向方向和倾斜方向两个方向开采的结果，在三维条件下，某采区开采后引起地表点A(x,y)的下沉值W(x,y)_A表示为：

式中，f(x,y)为空间概率密度函数，F表示区域范围，即x∈[0,l],y∈[0,L]的区域；

设矩形采区沿走向开采长度为l，沿倾向开采宽度为L，开采尺寸分别为x∈[0,l],y∈[0,L]，考虑到x,y方向概率的独立性，当以采区左下角为坐标原点时，对随机模型进行分段函数的表述，地表变形值可进一步表示为：

式中，r为煤层走向主要影响半径，r₁₍₂₎分别为煤层倾向在上山和下山方向上的主要影响半径；W_max＝m·q·cosα为地表最大下沉值，m为煤层开采厚度，q为下沉系数，α为煤层倾角，则地面点A(x,y)在方位角为β方向上的沉降梯度函数为下沉的一阶导数，某一方向的梯度计为W’(x,y)_β。

进一步，区域分类的原则是比较某一方向的梯度W’(x,y)_β与λ*u/2的大小关系，若前者小则划定为准真实形变区；否则，划定为待修订形变区；

利用公式(7)将完成分类的两类区域的地表变形信息分别转化为相位信息，并分别标记为准真实相位区和待修订相位区；

利用公式f(i,j)＝a*i²+b*j²+c*ij+d*i+e*j+f计算修正多项，式中，a、b、c、d、e、f为待求参数，(i,j)为平面坐标；确定f(i,j)的表达式，就是求取a、b、c、d、e、f这6个参数；该6参数的确定方法是：在公共部分内通过选取相同的点位，利用最小二乘法进行求解；

利用公式：计算修正相位φ′_defo(i,j)＝φ_defo(i,j)+f(i，j)，及相位整周数，

进一步，基于相位修正后的的解缠相位表达：

根据D-InSAR技术的基本原理，地表垂直沉降W与沿雷达视线向的形变相位

之间的关系式为：

式中，λ为雷达卫星波长；θ为视线角；

为沿雷达视线向的形变相位，是从干涉相位图中去除参考椭球相位、地形相位、大气延迟相位、噪声相位等之后的剩余部分，则步骤S5中将地表变形转换为相位的公式为：

因此，根据步骤S7，表下沉预计值探测的相位整周数n(x,y)：

n(x,y)＝Int[φ′_defo/(2π)] (8)

式中Int[·]表示取整运算，φ′_defo表示修正后的形变相位；

将式(8)代入(1)即可得步骤S8中地面任意点的解缠相位：

利用相位整周数n(x，y)对初始相位解缠结果进行修正。

有益效果：

本发明针对现有常规相位解缠方法难以准确实现矿区差分干涉相位解缠的问题，针对矿区开采沉陷的特点，提出的方法从开采沉陷的机理及规律出发，首先基于随机介质模型推演工作面开采引起的地表变形估值，然后以地表变形的估值来修正相位解缠结果，以达到准确获取解缠相位的目的。

本发明充分考虑了开采沉陷机理，利用由于地下煤炭开采而引起地表变形的估值来指导相位解缠结果和提高相位解缠精度，所提方法弥补了现有解缠算法在矿区地表快速变形区域无能为力的缺陷。提出的基于沉陷预计的相位解缠算法能够完成对矿区地表快速变形区域干涉图相位的解缠，与流行的最小费用流法相比，具有更为合理的解缠结果。

附图说明

图1为本发明的相位解缠示意图，上图为绝对相位；下图为缠绕相位；

图2为本发明融合随机介质模型的沉陷区差分干涉影像图相位解缠方法流程示意图：

图3为二维空间坐标系图；

图4为三维空间坐标系图；

图5为实施例一示意图，图中(a)：模拟的开采工作面及随机介质模型预计的地表下沉等值线；(b)：模拟的地表变形的真实相位；(c)：模拟的地表变形的缠绕相位；(d)：本方法探测的相位整周数；(e)：本文方法解缠后的解缠相位；(f)：真实相位与解缠相位差值的直方图；

图6为实施例二示意图，图中(a)：开采工作面及随机介质模型预计的地表下沉等值线；(b)：地面三维激光扫描仪获取的地下开采引起的地表变形；(c)：地下开采引起地表变形的真实相位；(d)：地下开采引起地表变形的缠绕相位；(e)：本方法探测的相位整周数；(f)：本方法解缠后的解缠相位；

图7为实施例三示意图，图中(a)：开采工作面及随机介质模型预计的地表下沉等值线；(b)：地下开采引起地表变形的缠绕相位；(cb)：最小费用流法相位解缠方法及本文所提方法探测的相位整周数；(d)：本方法解缠后的解缠相位；

图8为解缠相位重新缠绕后与原始缠绕相位差值的统计直方图。

具体实施方式

下面结合附图对本发明的实施例做进一步说明：

解缠方法

在一幅M×N的干涉相位图中，像元中直接观测到的相位值处于(-π,π)之间，称之为缠绕相位。为得到实际的地表形变量，必须清除2π整数倍的相位模糊性，恢复到真实相位，使真实相位与缠绕相位之间满足以下关系：

上式中，

为缠绕相位值；φ(i,j)为真实相位值；n(i,j)为整数。相位解缠就是从缠绕相位

中寻求最佳2π周期整数倍n(i,j)的过程，一维缠绕和解缠相位关系如下图1所示。

由图1可知，从理论上讲，如果相邻像元相位差的绝对值小于π，通过对相位差值进行积分，就能完成相位解缠。然而，在真实的SAR影像干涉图中，雷达阴影、相位噪声、地表快速形变等因素会造成相位数据的不连续，使得对相位差的积分变得极为困难。特别的，当利用D-InSAR技术对矿区地表形变进行监测时，由于地下煤炭资源采出后地表短时间内即会出现大级别的沉降，导致干涉相位图相邻像元间出现显著的不连续现象，在相位解缠处理过程中，造成对相位差值积分的不准确，影响解缠效果。

解缠策略

与由于地震、滑坡、冰川移动等造成地表变形不可预知性的特点相比，由地下开采而引起的地表变形符合一定的规律，并且地表任意点的形变量值及其空间分布状态是可以预计的。因此，我们就可利用两相邻像元间预计结果的差值探测相位跳变的整周数，以此来达到相位解缠的目的。

图2给出了基于该思路的流程图，具体流程说明如下：

1)利用SAR影像进行差分干涉处理(含配准、重采样、干涉、去地形相位等处理步骤)，得到目标区的含有形变信息的差分干涉条纹图。

2)利用现有常规相位解缠方法(如枝切法、最小费用流法等)对干涉相位图进行解缠处理，得到目标区的初始解缠相位图。

3)根据井下开采信息及目标区地质采矿情况，计算地下资源开采所引起地表变形的估值；以SAR影像分辨率为基准，对地表变形计算值进行重采样处理；将其统一到雷达坐标系统下，并计算相邻点间的形变梯度。

4)基于相位干涉的基本原理，以λ*u/2为基准(λ和u分别为SAR影像传感器的雷达波长和影像分辨率)，对第3步得到的形变图进行分类。分类原则为：形变梯度小于λ*u/2的区域标定为准真实形变区，其余为待修订形变区。

5)利用雷达传感器参数，基于公式(7)，将第4步中已分类的两类区域的地表变形分别转化为相位信息，并分别记为准真实相位区和待修订相位区。

6)将第5步得到的准真实相位区与第2步得到的初始解缠相位图的公共部分，采用最小二乘法进行相位匹配处理。将准真实相位修正到初始解缠相位，并得到修正多项式f(i，j)。

7)利用修正多项式f(i，j)对待修订相位区的相位信息进行修正，利用公式(8)计算修正后的相位整周数。

8)将第7步中获得的相位整周数与第2步得到的初始解缠相位进行融合，即可得到最终解缠相位。

要实现上述方法，有两个问题需要解决，第一个是地表变形值的预计，第二个是解缠相位的表达。

开采引起地表变形预计

矿山开采沉陷预计是开采沉陷学研究的核心问题，在众多沉陷预计方法中，波兰学者李特威尼申(J.Litwinisyn)提出的随机介质模型是最为成熟的方法之一，该方法已在中国得到了广泛应用，并且大量的生产实践已经证实，该方法所预计拐点以内变形值的准确度能达到90％以上(拐点位置如图3所示)；拐点以外的部分，通过模型修正以后，地表变形结果的预计值也可达到90％。因此，应用随机介质模型，采用分段函数的方式，预计工作面开采后地表变形值，进而根据变形预计值探测相邻像元间的相位整周数在理论上是可信的。

1)拐点以内地表下沉值的表达

随机介质模型有两个基本假设条件：a)开采引起的各方向的移动与方向无关；b)地下大工作面开采引起的地表变形是多个小工作面开采引起地表变形的线性叠加。在这两个假定条件下，二维平面内单元开采引起地表点x的下沉值W_e(x)可以表示为：

式中，r＝H₀/tanβ为常数，称为主要影响半径。其中，H₀是平均开采深度；tanβ是主要影响角正切。

2)拐点以外地表下沉值的表达

应用随机介质模型进行开采沉陷预计时，拐点以外的部分存在收敛过快的缺陷。从式(2)可以看出影响下沉曲线收敛性的主要因素是r，只要r发生变化，边界收敛性就会发生变化。因此，通过改变r值，即可实现对下沉曲线的控制。研究发现，拐点以外r的取值与拐点内r的取值之比近似为一常数k。故拐点以外，在二维平面内单元开采引起地表点的下沉值可以表示为：

式中，k为常数，可定义为沉陷曲线收敛系数。

3)三维平面任意点地表下沉值的表达

地表移动是个三维问题，即地表点A(x,y)的下沉是由于煤层沿走向方向(x轴方向)和倾斜方向(y轴方向)两个方向开采的结果。在三维条件下，某采区开采后引起地表点A(x,y)的下沉值W(x,y)_A可表示为：

式中，f(x,y)为空间概率密度函数。

如图4所示，假设一矩形采区沿走向开采长度为l，沿倾向开采宽度为L，开采尺寸分别为x∈[0,l],y∈[0,L]，考虑到x,y方向概率的独立性，当以采区左下角为坐标原点时方程式(4)可以表示为：

式中，r为煤层走向主要影响半径，r₁₍₂₎分别为煤层倾向在上山和下山方向上的主要影响半径；W_max＝m·q·cosα为地表最大下沉值，m为煤层开采厚度，q为下沉系数，α为煤层倾角。

则地面点A(x,y)在方位角为β方向上的沉降梯度函数为下沉的一阶导数，计为W’(x,y)_β。

基于预计结果的解缠相位表达

根据D-InSAR技术的基本原理，地表垂直沉降W与沿雷达视线向的形变相位

之间的关系可以表示为：

式中，λ为雷达卫星波长；θ为视线角；

为沿雷达视线向的形变相位，是从干涉相位图中去除参考椭球相位、地形相位、大气延迟相位、噪声相位等之后的剩余部分。则：

因此，根据地表下沉预计值探测的相位整周数n(x,y)可以表示为：

式中Int[·]表示取整运算，其余各参数含义与前述相同。

将式(8)代入(1)即可得地面任意点的解缠相位，如下方程式(9)所示：

式中，W(x,y)为分段函数，具体表达式为(5)。

实施例一、

以表1中的地质采矿条件模拟了一个开采工作面，其中工作面的开采深度和开采尺寸以及应用随机介质模型预计开采后的地表下沉等值线如图5(a)所示；图5(b)和图5(c)分别是以表2中的卫星参数模拟RADARSAT-2卫星利用差分干涉测量监测的地表形变真实相位和缠绕相位(即，差分干涉处理中需要解缠的缠绕相位)。

表1模拟开采的地质采矿条件

表2模拟SAR卫星的基本参数

缠绕相位估计

在现实情况下，开采沉陷预计时所采用的参数并不能完全真实的表示地质采矿条件，这是利用随机介质模型预计地下煤层开采引起地表变形不完全准确的原因。因此，在应用推导公式进行相位整周数探测时，采用表3中的另一套参数，以此来反应与实际开采条件不符合的现实。

表3相位解缠所用参数

利用式(8)、式(9)和表3中的参数分别进行相位整周数的探测和最终解缠相位的计算，其结果如图5中的(d)和(e)所示；图中(a)：模拟的开采工作面及随机介质模型预计的地表下沉等值线；(b)：模拟的地表变形的真实相位；(c)：模拟的地表变形的缠绕相位；(d)：本方法探测的相位整周数；(e)：本文方法解缠后的解缠相位；(f)：真实相位与解缠相位差值的直方图；

结果评价

将图5(e)中的解缠相位与图5(b)中的真实相位进行差值运算，图5(f)为相位差的统计直方图。从图5(f)中可知相位差为负数的像元数远大于为正数的像元数，其原因是探测相位整周数时所采用的下沉系数q偏小，使得探测到的相位整周数小于真实的整周数，从而导致解缠相位小于真实相位的概率增大；相位差的平均差和标准差分别为-0.0058rad和1.7405rad，两者都不足π，说明本文所提出的针对矿区差分干涉影像的相位解缠方法在无噪声条件下是可行的。

实施例二、万年矿13266工作面

万年矿位于河北省邯郸市峰峰矿区，其中13266工作面在2009.12.4-2010.3.11时间段内开采范围及应用随机介质模型预计开采后的地表下沉等值线如下图6(a)所示(预计参数如下表4所示)；图6(b)为工作面开采期间采用三维激光扫描仪获取的地表变形空间分布情况；以表2中的卫星参数模拟的地表形变的真实相位和缠绕相位(即，差分干涉处理中需要解缠的缠绕相位)分别如图6(c)和图6(d)所示。

表4相位解缠所用参数

利用式(8)、式(9)和表4中的参数分别进行相位整周数的探测和最终解缠相位的计算，其结果如图6(e)和图6(f)所示。

将解缠后的相位与真实相位相减，得到相位差的平均值为-1.56rad，小于π；如果将解缠后的相位不符值当作噪声，则计算的信噪比为58.56。两项指标均说明本文所提出的方法在解缠地下开采引起地表变形的缠绕相位时具有较高的有效性。

实施例三、

研究区域为九龙矿15223工作面，同样属于河北省邯郸市峰峰矿区，工作面在2014.1.3-2014.3.6时间段内开采范围及应用随机介质模型预计开采后的地表下沉等值线如下图7(a)所示，预计参数如表5所示：

表5相位解缠所用参数

该实验用于进行差分干涉测量的影像数据为RadarSat-2雷达卫星产品，影像的采集日期及相关参数如下表6所示。

表6 RadarSat-2影像参数信息

缠绕相位估计：

实验过程中，应用GAMMA软件，采用常规的“二轨法”进行SAR影像的差分干涉处理，得到的差分干涉相位如图7中(b)所示。

为获得地表变形引起的真实相位，我们采用两种方法进行了缠绕相位的估计：最小费用流法和本文所提的基于沉陷预计的相位解缠方法，其结果分别如图7中(c)和(d)。

结果评价：

从图7(c)可知，由于干涉相位图中高的条纹密集度以及相位噪声，使得最小费用流法无法实现研究区域内所有像元的相位解缠；而本文所提出的基于沉陷预计的相位解缠方法不仅充分考虑了开采沉陷的机理，而且还克服了常规相位解缠方法在矿区地表快速变形区域无能为力的缺陷。

将图7(d)中的解缠相位进行重新缠绕并与图7(b)中的干涉进行差值运算，相位差的统计直方图如图8所示。从图8中可看出解缠相位重新缠绕后与原始缠绕相位差符合正态分布，相位差的平均差和标准差分别为-0.0020rad和0.9011rad，两者均不足π，进一步表明该方法在较好的完成相位解缠的同时又保持了与原有缠绕相位的高度一致性。

针对矿区差分干涉相位图难以实现相位解缠的问题，从开采沉陷的机理出发，利用随机介质模型的基本理论，提出了基于沉陷预计估值的相位解缠方法。根据地下开采引起地表变形的预计值，推导出了地表任意点在干涉相位图中的解缠相位。

本方法分别利用无噪声条件下的模拟数据、含噪声的实测三维激光扫描数据和真实的InSAR影像对该方法的可行性、有效性和适用性进行了验证。理论分析和实例验证表明，该方法不仅能够实现对矿区地表快速变形区域差分干涉影像的相位解缠，而且具有很好的解缠效果，从而进一步提高了D-InSAR技术在矿区的应用。

Claims (7)

1.一种融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于：针对矿区开采沉陷的特点，从开采沉陷的机理及规律出发，基于随机介质模型推演工作面开采引起的地表变形估值，然后以地表变形的估值来修正相位解缠结果，从而准确获取地表任意点的解缠相位；

其步骤如下：

S1、对目标矿区地表SAR影像进行差分干涉处理，得到目标矿区地表的含有形变信息的差分干涉条纹图，差分干涉处理步骤包括：配准、重采样、干涉、去地形相位；

S2、利用现有常规相位解缠方法：枝切法、最小费用流法，对目标矿区地表的差分干涉条纹图进行解缠处理，得到目标矿区地表的初始解缠相位图；

S3、根据井下开采信息及目标矿区地质采矿情况，计算矿区地下资源开采所引起地表变形沉陷的估值；以目标矿区地表SAR影像分辨率为基准，对与SAR影像分辨率不一致的地表变形估值进行重采样使两者分辨率一致，重采样采用双线性内插法；地表变形估值的坐标系统一到SAR影像坐标系下，并生成地表变形图，之后计算地表变形图中相邻的地表变形估值点间的梯度；

S4、基于相位干涉的基本原理，以λ*u/2为基准，对地表变形图中的区域进行分类，将形变梯度小于λ*u/2的区域标定为准真实形变区，其余标定为待修订形变区，其中λ为SAR影像传感器的雷达波长，u为SAR影像传感器的影像分辨率；

S5、利用SAR的雷达参数，将完成分类的两类区域的地表变形信息分别转化为相位信息，并分别标记为准真实相位区和待修订相位区；

S6、将步骤S5得到的准真实相位区与步骤S2得到的初始解缠相位图两者之间坐标范围相同的区域定义为公共部分，采用最小二乘法对公共部分进行相位匹配处理，将准真实相位修正到初始解缠相位，得到修正多项式f(i，j)；

S7、利用修正多项式f(i，j)对待修订相位区的所有相位信息进行叠加，从而计算出修正后相邻像元间相位整周数；

S8、将步骤S7中获得的相位整周数与步骤S2得到的初始解缠相位进行融合，即可得到沉陷区差分干涉影像图的最终解缠相位。

2.根据权利要求1所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，相位解缠中，为得到实际的地表形变量，必须清除2π整数倍的相位模糊性，恢复到真实相位，因此真实相位与缠绕相位之间满足以下关系：

即从缠绕相位

中寻求最佳2π周期整数倍n(i,j)从而实现相位解缠；

式中，

为缠绕相位值；φ(i,j)为真实相位值；n(i,j)为整数；干涉相位图的像素为M×N，像元中直接观测到的相位值处于(-π,π)之间。

3.根据权利要求2所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，步骤3中地表变形值的具体计算过程为：利用随机介质模型，通过设定预计拐点，拐点内的地表变形值直接使用随机介质模型进行计算，即公式(2)，拐点以外地表下沉值不能直接使用随机介质模型的计算结果，需要对模型修正后使用，修正后的结果即为公式(3)，采用分段函数的方式，预计工作面开采后地表变形值，进而根据变形预计值探测相邻像元间的相位整周数：

拐点以内地表下沉值：

随机介质模型满足a)开采引起的各方向的移动与方向无关；b)地下大工作面开采引起的地表变形是多个小工作面开采引起地表变形的线性叠加；二维平面内单元开采引起地表点x的下沉值W_e(x)表示为：

式中，r＝H₀/tanβ为常数，称为主要影响半径，H₀是平均开采深度，tanβ是主要影响角正切，式中x表示地面任意点与拐点间的水平距离；

拐点以外地表下沉值：

应用随机介质模型进行开采沉陷预计时，拐点以外的部分存在收敛过快的缺陷，影响下沉曲线收敛性的主要因素是常数r，通过改变r值，即可实现对下沉曲线的控制；拐点以外r的取值与拐点内r的取值之比近似为一常数k，故拐点以外，在二维平面内单元开采引起地表点的下沉值表示为：

式中，k为常数，定义为沉陷曲线收敛系数；

三维平面任意点地表下沉值的表达：

地表点A(x,y)的下沉是由于煤层沿走向方向和倾斜方向两个方向开采的结果，在三维条件下，某采区开采后引起地表点A(x,y)的下沉值W(x,y)_A表示为：

式中，f(x,y)为空间概率密度函数，F表示区域范围，即x∈[0,l],y∈[0,L]的区域；

设矩形采区沿走向开采长度为l，沿倾向开采宽度为L，开采尺寸分别为x∈[0,l],y∈[0,L]，考虑到x,y方向概率的独立性，当以采区左下角为坐标原点时，对随机模型进行分段函数的表述，地表变形值可进一步表示为：

式中，r为煤层走向主要影响半径，r₁₍₂₎分别为煤层倾向在上山和下山方向上的主要影响半径；W_max＝m·q·cosα为地表最大下沉值，m为煤层开采厚度，q为下沉系数，α为煤层倾角，则地面点A(x,y)在方位角为β方向上的沉降梯度函数为下沉的一阶导数，某一方向的梯度计为W’(x,y)_β。

4.根据权利要求3所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，区域分类的原则是比较某一方向的梯度W’(x,y)_β与λ*u/2的大小关系，若前者小则划定为准真实形变区；否则，划定为待修订形变区；

利用公式(7)将完成分类的两类区域的地表变形信息分别转化为相位信息，并分别标记为准真实相位区和待修订相位区；

利用公式f(i,j)＝a*i²+b*j²+c*ij+d*i+e*j+f计算修正多项，式中，a、b、c、d、e、f为待求参数，(i,j)为平面坐标；确定f(i,j)的表达式，就是求取a、b、c、d、e、f这6个参数；该6参数的确定方法是：在公共部分内通过选取相同的点位，利用最小二乘法进行求解；

利用公式：计算修正相位φ’_defo(i,j)＝φ_defo(i,j)+f(i，j)，及相位整周数，；

为沿雷达视线向的形变相位。

5.根据权利要求4所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，基于相位修正后的的解缠相位表达：

根据D-InSAR技术的基本原理，地表垂直沉降W与沿雷达视线向的形变相位

之间的关系式为：

式中，λ为雷达卫星波长；θ为视线角；

为沿雷达视线向的形变相位，是从干涉相位图中去除参考椭球相位、地形相位、大气延迟相位、噪声相位等之后的剩余部分，则步骤S5中将地表变形转换为相位的公式为：

6.根据权利要求5所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，根据步骤S7，表下沉预计值探测的相位整周数n(x,y)：

n(x,y)＝Int[φ’_defo/(2π)] (8)

式中Int[·]表示取整运算，φ’_defo表示修正后的形变相位。

7.根据权利要求6所述融合随机介质模型的沉陷区差分干涉影像图相位解缠方法，其特征在于，将式(8)代入(1)即可得步骤S8中地面任意点的解缠相位：

利用相位整周数n(x，y)对初始相位解缠结果进行修正。

Images