融合时空模拟相位和DInSAR相位获取地表沉陷序列的方法

融合时空模拟相位和dinsar相位获取地表沉陷序列的方法
技术领域
1.本发明涉及矿区地表沉陷监测技术领域，尤其涉及一种融合时空模拟相位和dinsar相位获取地表沉陷序列的方法。


背景技术：

2.目前，在煤炭资源采集过程中，通常需要获取准确的矿区开采地表沉陷信息，为后续开采损害预测、控制及沉陷区综合治理等提供依据，其中，相关技术通常采用水准测量技术或dinsar(合成孔径雷达差分干涉测量技术)获取矿区开采地表沉陷信息。
3.然而，相关技术的问题在于，一方面在采用水准测量技术获取矿区开采地表沉陷信息的过程中，受矿区监测范围影响，水准点分布稀疏，观测周期长等影响，导致水准测量技术很难获取大区域面状地表形变信息，另一方面在采用dinsar技术获取矿区开采地表沉陷信息的过程中，受煤矿开采沉陷区域地表复杂、地表植被密集、地表沉陷速度快、开采沉陷范围相对较小等影响，导致dinsar技术在矿区开采地表沉陷监测中的准确性低，因此，亟需一种可以高准确性的获取大区域面状地表形变信息的矿区开采地表沉陷信息监测方法。


技术实现要素：

4.本发明旨在至少在一定程度上解决相关技术中的技术问题之一。为此，本发明的第一个目的在于提出一种融合时空模拟相位和dinsar相位获取地表沉陷序列的方法，能够根据融合相位获取矿区开采监测区域的地表沉陷序列，提高相位解缠的准确性，以提高地表沉陷序列的分辨率，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
5.为了达到上述目的，本发明实施例提出的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法，包括：获取矿区开采监测区域的概率积分法模型参数和时间函数模型参数，并根据所述概率积分法模型参数和所述时间函数模型参数，构建矿区开采沉陷时空模型；获取所述矿区开采监测区域的dinsar相位，并根据所述矿区开采沉陷时空模型获取所述矿区开采监测区域的模拟相位；根据所述dinsar相位和所述模拟相位获取融合相位，并对所述融合相位进行解缠，以获取所述矿区开采监测区域的地表沉陷序列。
6.根据本发明实施例提出的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法，获取矿区开采监测区域的概率积分法模型参数和时间函数模型参数，并根据概率积分法模型参数和时间函数模型参数，构建矿区开采沉陷时空模型，进而，获取矿区开采监测区域的dinsar相位，并根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位，以及，根据dinsar相位和模拟相位获取融合相位，并对融合相位进行解缠，以获取矿区开采监测区域的地表沉陷序列。由此，根据融合相位获取矿区开采监测区域的地表沉陷序列，提高相位解缠的准确性，以提高地表沉陷序列的分辨率，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
7.另外，根据本发明上述实施例的融合时空模拟相位和dinsar相位获取地表沉陷序
列的方法，还可以具有如下的附加技术特征：
8.根据本发明的一个实施例，所述获取矿区开采监测区域的概率积分法模型参数，包括：获取水准监测时序数据，并构建所述概率积分法模型参数的反演目标函数；根据所述水准监测时序数据和所述反演目标函数，利用多种预设算法，获取多组所述概率积分法模型参数；利用加权组合预测模型对多组所述概率积分法模型参数进行加权融合，并将加权融合后的组合预测参数作为所述矿区开采监测区域的概率积分法模型参数。
9.根据本发明的一个实施例，所述获取矿区开采监测区域的时间函数模型参数，包括：确定所述水准监测时序数据的监测点时序启动时刻，并根据所述监测点时序启动时刻调整所述水准监测时序数据；根据所述监测点时序启动时刻和调整后的水准监测时序数据，利用迭代求函数极值lm算法，得到所述矿区开采监测区域的时间函数模型参数。
10.根据本发明的一个实施例，所述获取所述矿区开采监测区域的dinsar相位，包括：获取所述矿区开采监测区域的历史雷达监测数据；按照dinsar成像时间生成所述历史雷达监测数据对应的干涉像对；根据所述干涉像对得到所述矿区开采监测区域的dinsar相位。
11.根据本发明的一个实施例，所述根据所述干涉像对得到所述矿区开采监测区域的dinsar相位，包括：对所述干涉像对依次进行影像裁剪、精密轨道数据加载、去平地相位、去地形相位、干涉、相位滤波、地理编码，得到所述矿区开采监测区域的dinsar相位。
12.根据本发明的一个实施例，根据所述矿区开采沉陷时空模型获取所述矿区开采监测区域的模拟相位，包括：根据所述矿区开采沉陷时空模型计算与所述dinsar成像时间相对应的地表垂直沉降量；将所述地表垂直沉降量转换为雷达视线向的斜距变化量；将所述斜距变化量进行相位转换，得到所述矿区开采监测区域的模拟相位。
13.根据本发明的一个实施例，所述根据所述dinsar相位和所述模拟相位获取融合相位，包括：对所述dinsar相位和所述模拟相位进行空间重采样；计算空间重采样后的dinsar相位和模拟相位的相位差；利用卡尔曼滤波算法对所述相位差进行滤波处理，得到滤波相位；对所述模拟相位和所述滤波相位进行求和，得到所述融合相位。
14.根据本发明的一个实施例，在获取所述融合相位之后，所述方法还包括：将所述融合相位修正至区间[-π,π]。
[0015]
根据本发明的一个实施例，所述对所述融合相位进行解缠，包括：利用非连续路径可靠性排序的快速二维相位展开算法对所述融合相位进行解缠。
[0016]
根据本发明的一个实施例，在对所述融合相位进行解缠之后，所述方法还包括：对解缠结果进行loess平滑处理，得到平滑后的所述矿区开采监测区域的地表沉陷序列。
[0017]
本发明附加的方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0018]
图1是根据本发明实施例的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法的流程示意图；
[0019]
图2是根据本发明实施例的获取矿区开采监测区域的概率积分法模型参数的流程示意图；
[0020]
图3是根据本发明的一个具体实施例的获取矿区开采监测区域的概率积分法模型
参数的流程示意图；
[0021]
图4是根据本发明实施例的获取矿区开采监测区域的时间函数模型参数的流程示意图；
[0022]
图5是根据本发明的一个具体实施例的获取矿区开采监测区域的时间函数模型参数的流程示意图；
[0023]
图6是根据本发明实施例的获取矿区开采监测区域的dinsar相位的流程示意图；
[0024]
图7是根据本发明实施例的根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位的流程示意图；
[0025]
图8是根据本发明的一个具体实施例的根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位的流程示意图；
[0026]
图9是根据本发明实施例的根据dinsar相位和模拟相位获取融合相位的流程示意图；
[0027]
图10是根据本发明的一个具体实施例的根据dinsar相位和模拟相位获取融合相位的流程示意图。
具体实施方式
[0028]
下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。
[0029]
下面参考附图描述本发明实施例的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法。
[0030]
图1是根据本发明实施例的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法的流程示意图。
[0031]
如图1所示，融合时空模拟相位和dinsar相位获取地表沉陷序列的方法，包括以下步骤：
[0032]
s101，获取矿区开采监测区域的概率积分法模型参数和时间函数模型参数，并根据概率积分法模型参数和时间函数模型参数，构建矿区开采沉陷时空模型。
[0033]
具体而言，在本发明的一些实施例中，可以通过获取得到的矿区开采监测区域的概率积分法模型参数par
pim
，构建对应矿区开采监测区域的概率积分法模型(pim模型)，以及，通过获取得到的矿区开采监测区域的时间函数模型参数，构建对应矿区开采监测区域的时间函数模型(例如，gompertz模型)，进而，根据概率积分法模型和时间函数模型构建对应矿区开采监测区域的矿区开采沉陷时空模型(例如，gompertz-pim模型)。
[0034]
s102，获取矿区开采监测区域的dinsar相位，并根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位。
[0035]
具体而言，在本发明的一些实施例中，可以基于雷达监测数据获取矿区开采监测区域的dinsar(合成孔径雷达差分干涉)相位以及，可以基于水准监测时序数据，根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位
[0036]
s103，根据dinsar相位和模拟相位获取融合相位，并对融合相位进行解缠，以获取
矿区开采监测区域的地表沉陷序列。
[0037]
应理解的是，在本发明的上述实施例中，首先，获取矿区开采监测区域的概率积分法模型参数和时间函数模型参数，并根据概率积分法模型参数和时间函数模型参数，构建矿区开采沉陷时空模型，进而，获取矿区开采监测区域的dinsar相位，并根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位，然后，根据dinsar相位和模拟相位获取融合相位，并对融合相位进行解缠，以获取矿区开采监测区域的地表沉陷序列。由此，根据融合相位获取矿区开采监测区域的地表沉陷序列，提高相位解缠的准确性，以提高地表沉陷序列的分辨率，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
[0038]
进一步地，如图2所示，获取矿区开采监测区域的概率积分法模型参数，包括：
[0039]
s201，获取水准监测时序数据，并构建概率积分法模型参数的反演目标函数。
[0040]
具体而言，在本发明的一些实施例中，可以获取水准监测时序数据，并构建概率积分法模型参数的反演目标函数，从而利用反演目标函数和水准监测时序数据，对概率积分法模型参数进行参数反演。
[0041]
s202，根据水准监测时序数据和反演目标函数，利用多种预设算法，获取多组概率积分法模型参数。
[0042]
具体而言，在本发明的一些实施例中，可以根据水准监测时序数据和反演目标函数，利用多种预设算法(例如，遗传算法ga、粒子群算法pso、差分进化算法de、鲸鱼优化算法woa)，获取多组概率积分法模型参数。
[0043]
s203，利用加权组合预测模型对多组概率积分法模型参数进行加权融合，并将加权融合后的组合预测参数作为矿区开采监测区域的概率积分法模型参数。
[0044]
需要说明的是，由于在概率积分法模型参数反演过程中，受参数相关性和预设算法的不同，存在多次反演结果不同的问题，即得到多组概率积分法模型参数，因此，在本发明的一些实施例中，还可以利用加权组合预测模型(owcm模型)对多组概率积分法模型参数进行加权融合，并将加权融合后的组合预测参数作为矿区开采监测区域的概率积分法模型参数，从而提高构建对应矿区开采监测区域的概率积分法模型的精确度和可靠性。
[0045]
下面结合附图3与本发明的具体实施例，对获取矿区开采监测区域的概率积分法模型参数的具体过程进行相应的说明。
[0046]
具体地，获取矿区开采监测区域的概率积分法模型参数，包括以下步骤：
[0047]
s21，根据pim原理，计算矿区开采地表上任意点(x,y)的沉降量。
[0048]
具体地，可以通过以下计算公式得到矿区开采地表上任意点(x,y)的沉降量：
[0049][0050]
[0051]
其中，w0为最大下沉值，w0(x)为走向主断面上投影点的下沉值、w0(y)倾向主断面上投影点的下沉值，u为概率积分的开采单元，m为开采厚度，h0为平均开采深度，h1为下山开采深度，、h2为上山开采深度，q为下沉系数，tanβ为主要影响角正切值，s4为左边界的拐点偏移距，s3为右边界的拐点偏移距，s2为上山拐点偏移距，s1为下山拐点偏移距，θ为开采影响传播角，d1为走向工作面长度，d3倾向工作面长度。
[0052]
需要说明的是，上山指的是在开采水平之上，沿着煤层倾向向上，掘进的巷道是上山(即斜向上的巷道)，下山指的是在开采水平之下，沿着煤层倾向向下，掘进的巷道是下山(就是斜向下的巷道)，走向指的是岩层面与水平面的交线，走向线两端所指方向是岩层走向，其中，走向线有无数条平行线，但走向只有两个，且相差180
°
，倾向指的是倾向线在水平面投影的方向，其中，层面上与走向线垂直并沿斜面向下所引的直线叫倾向线，拐点偏移距指的是实际开采边界与计算边界之间的距离。
[0053]
基于上述公式可知，需要进行反演的矿区开采监测区域的概率积分法模型参数可以包括：q、tanβ、θ、s4、s3、s2、s1。
[0054]
因此，在本发明的一个具体实施例中，可以构建概率积分法参数反演目标函数f
vv
，其中，f
vv
表达式如下：
[0055][0056]
其中，n为参与参数反演的水准监测时序数据监测点个数，w为pim计算值，为水准监测时序数据值。
[0057]
s22，利用多种预设算法(例如，遗传算法、粒子群算法、差分进化算法、鲸鱼优化算法)，循环迭代使f
vv
达到最小。
[0058]
具体地，基于矿区开采监测区域的概率积分法模型参数实际值而由于受参数相关性及预设算法选择的不同，以预设算法分别为上述遗传算法ga、粒子群算法pso、差分进化算法de、鲸鱼优化算法woa为例，可以得到不同的矿区开采监测区域的概率积分法模型参数：(par1,par2,par3,par4)，
[0059]
因此，在本发明的一个具体实施例中，还利用加权组合预测模型对多组概率积分法模型参数进行加权融合。
[0060]
具体地，加权组合预测模型中的权重为w＝(w1,w2,w3,w4)
t
，利用第j种预设算法的矿区开采监测区域的概率积分法模型参数预测值parj＝[qj,bj,tanβj,θj,(s4,s3,s2,s1)j]，其中j＝1,2,3,4。
[0061]
进一步地，可以根据矿区开采监测区域的概率积分法模型参数实际值和预测值计算拟合偏差，其中，计算公式如下：
[0062][0063]
其中，e
t
为拟合偏差、为矿区开采监测区域的概率积分法模型参数实际值，parj为利用第j种预设算法的矿区开采监测区域的概率积分法模型参数预测值。
[0064]
另外，加权组合预测模型的预测误差平方和为：
[0065][0066][0067]
目标函数表达式为：
[0068][0069]
r＝(1,1,1
…
,1)
t
[0070]
利用拉格朗日乘子算法对目标函数进行求解，引入拉格朗日乘数λ，q可以表示为：
[0071]
q＝w
t
ew 2λ(r
t
w-1)
[0072]
由得：
[0073]rt
w＝1
[0074]rt
λe-1
r＝1
[0075][0076]
因此，最优权重向量为：
[0077][0078]
目标函数的最小值为：
[0079][0080]
构造加权组合预测模型：
[0081][0082]
s23，利用加权组合预测模型对多组概率积分法模型参数进行加权融合，得到加权融合后的矿区开采监测区域的概率积分法模型参数，并根据矿区开采监测区域的概率积分法模型参数，构建对应矿区开采监测区域的概率积分法模型。
[0083]
应理解的是，在本发明的上述实施例中，可以利用多种预设算法(例如，遗传算法ga、粒子群算法pso、差分进化算法de、鲸鱼优化算法woa)，反演多组矿区开采监测区域的概率积分法模型参数，并利用加权组合预测模型对多组概率积分法模型参数进行加权融合，得到加权融合后的矿区开采监测区域的概率积分法模型参数par
pim
，从而可以利用概率积分法模型参数par
pim
构建对应矿区开采监测区域的概率积分法模型。
[0084]
进一步地，如图4所示，获取矿区开采监测区域的时间函数模型参数，包括：
[0085]
s301，确定水准监测时序数据的监测点时序启动时刻，并根据监测点时序启动时刻调整水准监测时序数据。
[0086]
具体地，在本发明的一些实施例中，可以将工作面开采至对应监测点正下方时的时刻确定为水准监测时序数据的监测点时序启动时刻，并根据监测点时序启动时刻调整水准监测时序数据，从而确保每个监测点的水准监测时序数据的一致性，有利于提高构建时间函数模型构的精确度与可靠性。
[0087]
s302，根据监测点时序启动时刻和调整后的水准监测时序数据，利用迭代求函数极值lm算法，得到矿区开采监测区域的时间函数模型参数。
[0088]
具体而言，在本发明的一些实施例中，可以利用迭代求函数极值lm算法，结合监测点时序启动时刻和调整后的水准监测时序数据进行迭代计算，从而得到矿区开采监测区域的时间函数模型参数。
[0089]
下面结合附图5和本发明的具体实施例，对获取矿区开采监测区域的时间函数模型参数的具体过程进行相应的说明。
[0090]
具体地，获取矿区开采监测区域的时间函数模型参数，包括以下步骤：
[0091]
s31，确定任意监测点时间函数。
[0092]
其中，任意监测点gompertz时间函数表达式为：
[0093][0094]
其中，w
max
为最大沉降量，k为增长系数，b为拐点，其中，拐点具体用于描述矿区开采地表动态形变，分别受到开始沉降时间与沉陷速度影响。
[0095]
进一步地，在本发明的一个具体实施例，将工作面开采至对应监测点正下方时的时刻确定为水准监测时序数据的监测点时序启动时刻，并根据监测点时序启动时刻调整水准监测时序数据。
[0096]
s32，以水准监测时序数据的监测点时序启动时刻作为gompertz时间函数g(t)的启动时刻，并利用迭代求函数极值lm算法结合调整后的水准监测时序数据，对gompertz时间函数进行循环运算，得到矿区开采监测区域的时间函数模型参数(k，b)，从而根据时间函数模型参数(k，b)构建对应矿区开采监测区域的时间函数模型。
[0097]
应理解的是，在本发明的上述实施例中，可以利用水准监测时序数据，结合迭代求函数极值(lm，levenberg-marquardt)算法，反演矿区开采监测区域的时间函数模型参数(k，b)，从而可以利用时间函数模型参数(k，b)构建对应矿区开采监测区域的时间函数模型。
[0098]
进一步地，在前述利用概率积分法模型参数par
pim
构建对应矿区开采监测区域的概率积分法模型w(x,y)，以及前述利用时间函数模型参数(k，b)构建对应矿区开采监测区域的时间函数模型g(t)之后，对根据概率积分法模型和时间函数模型构建对应矿区开采监测区域的矿区开采沉陷时空模型的具体过程进行说明。
[0099]
具体地，可以通过计算平均每天的工作面推进距离[l1,l2,
…
,ln]，并将开采工作面以天为单位划分n块，其中，n的大小由工作面开采速度决定，此时，可以根据概率积分法模型得到每块矿区开采监测区域的沉降量[w1,w2,
…
,wn]，以及，可以根据时间函数模型获取每块矿区开采监测区域的动态沉降量和地表动态沉降量：wi(t)＝g(t)wi，其中，wi(t)为
每块矿区开采监测区域的动态沉降量；其中，w(t)为每块矿区开采监测区域地表动态沉降量。
[0100]
进一步地，如图6所示，获取矿区开采监测区域的dinsar相位，包括：
[0101]
s401，获取矿区开采监测区域的历史雷达监测数据。
[0102]
具体而言，在本发明的一些实施例中，矿区开采监测区域的历史雷达监测数据可以为c波段合成孔径雷达卫星预设时间段内监测矿区开采监测区域得到的雷达监测数据。
[0103]
s402，按照dinsar成像时间生成历史雷达监测数据对应的干涉像对。
[0104]
具体而言，在本发明的一些实施例中，可以按照dinsar成像时间对历史雷达监测数据进行差分处理，进而生成历史雷达监测数据对应的干涉像对。
[0105]
s403，根据干涉像对得到矿区开采监测区域的dinsar相位。
[0106]
具体而言，在本发明的一些实施例中，可以对干涉像对进行相应的图像处理，进而得到矿区开采监测区域的dinsar相位。
[0107]
下面结合本发明的具体实施例，对获取矿区开采监测区域的dinsar相位的具体过程进行相应的说明。
[0108]
具体地，insar可以通过传感器两次回波信号的相位差和卫星轨道的参数获取相位信息，其中，相位信息可表达为：
[0109][0110]
其中，为insar相位信息，为平地相位，为地形相位，为地表形变相位，为电离层延迟和对流层延迟引起的大气相位，以及为噪声相位。
[0111]
进一步地，往返斜距分别对应斜距相位ψ1和ψ2，可以分别表示为：
[0112][0113][0114]
其中，n1，n2分别为雷达沿斜距r1和r2往返传播的相位整周数，和分别为雷达沿斜距r1和r2往返传播的非整周相位，λ为雷达波长。
[0115]
更进一步地，在本发明的一个具体实施例中，在去除大气相位平地相位和噪声相位之后，地表形变相位表达式为：
[0116][0117]
其中，r
1-r2为受地形起伏影响下的相位变化，δr为dinsar监测地表形变沿雷达视线方向的斜距变化量，此时，可以通过外部数字高程模型(dem)模拟地形相位并将项移除，则地表形变相位表达式可以进一步地更新为：
[0118][0119]
进而，在本发明的一个具体实施例中，可以利用矿区开采监测区域的历史雷达监测数据，结合上述地表形变相位表达式获得干涉像对，并利用snap软件对干涉像对进
行相应的图像处理，进而得到矿区开采监测区域的dinsar相位
[0120]
应理解的是，在本发明的上述实施例中，可以从历史雷达监测数据中获取矿区开采监测区域的历史雷达监测数据，并按照dinsar成像时间生成历史雷达监测数据对应的干涉像对，从而可以根据干涉像对得到矿区开采监测区域的dinsar相位
[0121]
进一步地，根据干涉像对得到矿区开采监测区域的dinsar相位，包括：对干涉像对依次进行影像裁剪、精密轨道数据加载、去平地相位、去地形相位、干涉、相位滤波、地理编码，得到矿区开采监测区域的dinsar相位。
[0122]
应理解的是，在本发明的上述实施例中，可以利用snap软件对干涉像对依次进行影像裁剪、精密轨道数据加载、去平地相位、去地形相位、干涉、相位滤波、地理编码等数据处理，从而得到矿区开采监测区域的dinsar相位。
[0123]
进一步地，如图7所示，根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位，包括：
[0124]
s501，根据矿区开采沉陷时空模型计算与dinsar成像时间相对应的地表垂直沉降量。
[0125]
具体而言，在本发明的一些实施例中，可以利用构建好的矿区开采沉陷时空模型，计算与dinsar成像时间相对应的地表垂直沉降量。
[0126]
s502，将地表垂直沉降量转换为雷达视线向的斜距变化量。
[0127]
具体而言，在本发明的一些实施例中，可以将地表垂直沉降量转换为雷达视线向的斜距变化量，以便于根据雷达视线向的斜距变化量得到矿区开采监测区域的模拟相位
[0128]
s503，将斜距变化量进行相位转换，得到矿区开采监测区域的模拟相位。
[0129]
具体而言，在本发明的一些实施例中，可以对雷达视线向的斜距变化量进行相位转换，从而得到矿区开采监测区域的模拟相位
[0130]
下面结合附图8与本发明的具体实施例，对得到矿区开采监测区域的模拟相位的具体过程进行相应的说明。
[0131]
s51，通过构建的gompertz-pim矿区开采沉陷时空模型计算与dinsar对应的地表垂直沉降量。
[0132]
s52，将地表垂直沉降量转换为雷达视线向los的斜距变化量。
[0133]
其中，雷达视线向los的斜距变化量的计算公式如下：
[0134]
δr(t)＝w(t)cosθ；
[0135]
其中，θ为雷达入射角，δr(t)为雷达视线向los的斜距变化量，w(t)为地表垂直沉降量。
[0136]
s53，将雷达视线向的斜距变化量转换为矿区开采监测区域的模拟相位值。
[0137]
其中，矿区开采监测区域的模拟相位值的计算公式如下：
[0138][0139]
其中，为矿区开采监测区域的模拟相位值。
[0140]
基于上述公式，由地表动态沉降量w(t)可生成雷达视线向的模拟相位(即矿区开采监测区域的模拟相位)，计算公式如下：
[0141][0142]
应理解的是，在本发明的上述实施例中，利用构建的矿区开采沉陷时空模型计算与dinsar成像时间相对应的地表垂直沉降量w(t)，从而根据相位生成原理，以w(t)计算得到矿区开采监测区域的模拟相位
[0143]
进一步地，如图9所示，根据dinsar相位和模拟相位获取融合相位，包括：
[0144]
s601，对dinsar相位和模拟相位进行空间重采样。
[0145]
具体而言，在本发明的一些实施例中，可以对dinsar相位和模拟相位进行空间重采样，以使dinsar相位和模拟相位的空间分辨率一致，从而后续将dinsar相位和模拟相位进行相位融合。
[0146]
s602，计算空间重采样后的dinsar相位和模拟相位的相位差。
[0147]
具体而言，在本发明的一些实施例中，由于空间重采样后的dinsar相位和模拟相位的空间分辨率一致，因此，可以将空间重采样后的dinsar相位和模拟相位的差值作为dinsar相位和模拟相位的相位差。
[0148]
s603，利用卡尔曼滤波算法对相位差进行滤波处理，得到滤波相位。
[0149]
具体而言，在本发明的一些实施例中，可以利用卡尔曼滤波算法对空间重采样后的dinsar相位和模拟相位的相位差进行滤波处理，从而得到滤波相位。
[0150]
s604，对模拟相位和滤波相位进行求和，得到融合相位。
[0151]
可选地，在本发明的另一些实施例中，还可以采用其他融合算法(例如，赫尔默特方差分量估计、神经网络等)对dinsar相位和模拟相位进行融合处理，得到融合相位。
[0152]
具体而言，在本发明的一些实施例中，在得到滤波相位之后，可以将模拟相位和滤波相位相加，从而得到融合相位。
[0153]
下面结合附图10与本发明的具体实施例对得到融合相位的具体过程进行相应的说明。
[0154]
s61，将dinsar相位空间重采样至与模拟相位的空间分辨率一致。
[0155]
可以理解的是，在本发明的一些实施例中，通过将dinsar相位空间重采样至与模拟相位的空间分辨率一致，从而为后续相位融合做准备。
[0156]
s62，计算模拟相位与dinsar相位的相位差。
[0157]
其中，模拟相位与dinsar相位的相位差的计算公式如下：
[0158][0159]
其中，为模拟相位，为dinsar相位，为相位差。
[0160]
s63，利用卡尔曼滤波算法对相位差滤波处理得到滤波相位。
[0161]
具体地，在本发明的一个具体实施例中，通过卡尔曼滤波算法以的1/3标准差为观测噪音方差r对所有像元的dinsar相位与模拟相位进行滤波处理，并执行如下计算：
[0162]
(1)计算t 1时刻的状态预报值：
[0163]
(2)计算预报值的协方差矩阵：p(t 1|t)＝φp(t|t)φ
t
γqγ
t
；
[0164]
(3)计算t 1时刻的增益矩阵：k(t 1)＝p(t 1|t)h
t
[hp(t 1|t)h
t
r]-1
；
[0165]
(4)计算t 1时刻的状态估值：(4)计算t 1时刻的状态估值：
[0166]
(5)计算状态估值得协方差矩阵：p(t 1|t 1)＝[i
n-k(t 1)h]p(t 1|t)；
[0167]
由此，通过上述循环运算，获得滤波相位
[0168]
s64，将模拟相位与滤波相位相加得到融合相位。
[0169]
其中，融合相位的计算公式如下：
[0170][0171]
其中，为滤波相位，为融合相位。
[0172]
应理解的是，在本发明的上述实施例中，可以利用卡尔曼滤波算法融合模拟相位与dinsar相位得到融合相位从而根据融合相位获取矿区开采监测区域的地表沉陷序列，提高相位解缠的准确性，以提高地表沉陷序列的分辨率，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
[0173]
进一步地，在获取融合相位之后，方法还包括：将融合相位修正至区间[-π,π]。
[0174]
应理解的是，在本发明的上述实施例中，可以将融合相位修正至区间[-π,π]，其中，计算公式如下：
[0175][0176]
进一步地，对融合相位进行解缠，包括：利用非连续路径可靠性排序的快速二维相位展开算法对融合相位进行解缠。
[0177]
可选地，在本发明的另一些实施例中，可以利用lp最小范数法、最小费用流法等对融合相位进行解缠。
[0178]
应理解的是，在本发明的上述实施例中，可以利用非连续路径可靠性排序的快速二维相位展开算法，对融合相位进行解缠，从而得到矿区开采监测区域的地表沉陷序列。
[0179]
进一步地，在对融合相位进行解缠之后，方法还包括：对解缠结果进行loess平滑处理，得到平滑后的矿区开采监测区域的地表沉陷序列。
[0180]
可选地，在本发明的另一些实施例中，可以利用bisquare算法、adjacent-averaging算法等对解缠结果进行平滑处理。
[0181]
应理解的是，在本发明的上述实施例中，可以对解缠结果进行loess平滑处理，得到平滑后的矿区开采监测区域的地表沉陷序列，以提高融合相位解缠结果的准确性，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
[0182]
综上，根据本发明实施例提出的融合时空模拟相位和dinsar相位获取地表沉陷序列的方法，获取矿区开采监测区域的概率积分法模型参数和时间函数模型参数，并根据概率积分法模型参数和时间函数模型参数，构建矿区开采沉陷时空模型，进而，获取矿区开采监测区域的dinsar相位，并根据矿区开采沉陷时空模型获取矿区开采监测区域的模拟相位，以及，根据dinsar相位和模拟相位获取融合相位，并对融合相位进行解缠，以获取矿区
开采监测区域的地表沉陷序列。由此，根据融合相位获取矿区开采监测区域的地表沉陷序列，提高相位解缠的准确性，以提高地表沉陷序列的分辨率，从而，在提高矿区开采监测效率和监测结果可靠性的同时，降低矿区开采监测成本。
[0183]
需要说明的是，在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读介质中，以供指令执行系统、装置或设备(如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统)使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，"计算机可读介质"可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。计算机可读介质的更具体的示例(非穷尽性列表)包括以下：具有一个或多个布线的电连接部(电子装置)，便携式计算机盘盒(磁装置)，随机存取存储器(ram)，只读存储器(rom)，可擦除可编辑只读存储器(eprom或闪速存储器)，光纤装置，以及便携式光盘只读存储器(cdrom)。另外，计算机可读介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以例如通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。
[0184]
应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(pga)，现场可编程门阵列(fpga)等。
[0185]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
[0186]
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“顺时针”、“逆时针”、“轴向”、“径向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
[0187]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。
[0188]
在本发明中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内
部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。
[0189]
在本发明中，除非另有明确的规定和限定，第一特征在第二特征“上”或“下”可以是第一和第二特征直接接触，或第一和第二特征通过中间媒介间接接触。而且，第一特征在第二特征“之上”、“上方”和“上面”可是第一特征在第二特征正上方或斜上方，或仅仅表示第一特征水平高度高于第二特征。第一特征在第二特征“之下”、“下方”和“下面”可以是第一特征在第二特征正下方或斜下方，或仅仅表示第一特征水平高度小于第二特征。
[0190]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。