CN105806303B - 融合D-InSAR和模矢法求取概率积分参数的方法 - Google Patents

Abstract

本发明提供一种融合D‑InSAR和模矢法求取概率积分参数的方法，包括：(1)利用概率积分法得到目标像元的预计下沉值、南北方向预计水平移动和东西方向预计水平移动的值，基于三维变形与LOS向变形关系，计算得到目标像元LOS向预计移动变形r'_iLOS的值；(2)利用D‑InSAR技术实测目标像元的LOS向移动变形r_iLOS的值，计算目标像元的移动变形预计残差v_i＝r_iLOS‑r'_iLOS的值，并构造求参误差函数ε(B)＝∑|v_i|；(3)基于ε(B)＝min准则，利用模矢法求取全部概率积分参数。优点为：本发明有效融合了D‑InSAR和模矢法，可求取到全部概率积分参数，并且还具有参数求取精度高的优点。

Description

融合D-InSAR和模矢法求取概率积分参数的方法

技术领域

本发明属于矿山变形监测技术领域，具体涉及一种融合D-InSAR和模矢法求取概率积分参数的方法。

背景技术

开采沉陷预计对建筑物下、铁路下、水体下(简称“三下”)采煤实践具有重要指导意义。概率积分法是我国“三下”采煤规程指定的开采沉陷预计方法，因而概率积分参数的求取也是矿山变形监测数据处理的核心内容。常规而言，开展概率积分参数反演的前提为准确获取采动地表三维变形。D-InSAR(合成孔径雷达差分干涉测量)由于具有全天时、高精度、似面测量和低费用等特点，已经成为当前矿山变形监测的研究热点。然而，当前单视向D-InSAR技术主要用于监测开采沉陷的下沉分量(对下沉敏感)，无法提取到地表三维变形(特别是水平移动)，因而，基于单视向D-InSAR技术无法求取全部概率积分参数，特别是水平移动系数。

解决上述难题，目前采用的主要方法为：建立适用于开采沉陷特点的D-InSAR三维变形监测方法，然后直接进行求参。现有技术出现的解决方案主要有：多视向D-InSAR、D-InSAR+Offset tracking、D-InSAR+MAI、融合GPS和D-InSAR法、地表先验模型+D-InSAR等方法，但上述方法主要用于解决地震、滑坡、火山喷发、冰川移动等特点的移动变形问题，然而上述对象的移动机理不同于矿山开采沉陷，并不能妥善解决矿山开采地表的三维变形监测问题。所以利用常规方法无法求取全部概率积分参数。

可见，如何基于单视向D-InSAR技术，求取全部概率积分参数，特别是水平移动系数，从而实现开采沉陷预计，具有重要意义。现有技术中尚未出现有效的解决方案。

发明内容

针对现有技术存在的缺陷，本发明提供一种融合D-InSAR和模矢法求取概率积分参数的方法，可有效解决上述问题。

本发明采用的技术方案如下：

本发明提供一种融合D-InSAR和模矢法求取概率积分参数的方法，包括以下步骤：

步骤1，根据地质采矿条件，预设工作面开采沉陷预计概率积分参数初值；

步骤2，通过雷达D-InSAR技术获取采动地表变形场影像；从所述采动地表变形场影像中选取n个目标像元；对于每个所述目标像元，均执行步骤2.1-步骤2.3：

步骤2.1，获取任意形状工作面地质采矿条件参数，并对工作面形状进行剖分，同时获取每个剖分单元的参数；另外，提取目标像元中心地理坐标(x,y)；以步骤1的工作面开采沉陷预计概率积分参数初值、剖分单元的参数、目标像元中心地理坐标(x,y)为概率积分法模型的输入，利用所述概率积分法模型预计得到目标像元的预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值；

步骤2.2，将步骤2.1计算得到的目标像元预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值，以及获取到的雷达卫星的入射角θ_i的值和卫星飞行方向方位角α_i的值，代入公式(1)，得到目标像元LOS向预计移动变形r'_iLOS的值；

步骤2.3，利用D-InSAR技术提取目标像元的LOS向实测移动变形r_iLOS的值，采用公式(2)计算得到目标像元的LOS向变形预计残差v_i的值；

v_i＝r_iLOS-r'_iLOS (2)

由此分别计算得到n个目标像元的LOS向变形预计残差v_i的值；

步骤3，构筑求参误差函数ε(B)：

其中，B为全部概率积分参数组成的矩阵；

将所计算得到的n个目标像元的LOS向变形预计残差v_i的值代入公式(3)，计算得到误差函数值；基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数。

优选的，步骤1中，所预设的工作面开采沉陷预计概率积分参数初值包括左拐点偏移距S₁、右拐点偏移距S₂、上拐点偏移距S₃、下拐点偏移距S₄、下沉系数q、主要影响角正切tanβ、最大下沉角θ和水平移动系数b。

优选的，步骤2.1中，对工作面形状进行剖分，同时获取每个剖分单元的参数，具体为：首先获取工作面采深H，然后以边长为H/10的正方形为剖分单元，沿走向和倾向对工作面进行剖分，设任意一个剖分单元为剖分单元j，获取剖分单元j如下的剖分参数：单元面积A_j、单元中心采深H_j和单元中心坐标(x_j,y_j)；

所述获取任意形状工作面地质采矿条件参数包括：工作面尺寸D₁和D₃、采深H、采高m和煤层倾角α。

优选的，步骤2.1中，所述概率积分法模型为：

W_i＝[∑W_j(x,y)·A_j]·W₀ (4)

U_iSN＝[r_j·(-2π/r_j ²)·(x-x_j)·W_j(x,y)·A_j]·b·W₀ (5)

其中：W_j(x,y)＝(1/r_j ²)·exp(-π·(x-x_j)²/r_j ²)·exp(-π·(y-y_j+l_j)²/r_j ²)、l_j＝H_jcotθ、r_j＝H_j/tanβ、W₀＝mqcosα。

优选的，步骤3中，基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数具体为：

步骤3.1，令B₁＝B，并以B₁作为初始基点进行搜索；

其中：

B₁所包含的各个预计参数即为B₁的分量，值为步骤1所设置的值；

步骤3.2，确定B₁的各个分量的探索步长；

即：对于任意分量X_i0(i＝1,2,...,8)，其步长为Δ_i，表示为：

Δ_i＝(0,…,ΔX_i0,…,0)^T (8)

步骤3.3，根据步骤1-步骤3，计算初始基点B₁的误差函数值ε(B₁)，另外计算B₁+Δ₁的误差函数值ε(B₁+Δ₁)；

比较ε(B₁)和ε(B₁+Δ₁)：

若ε(B₁+Δ₁)＜ε(B₁)，则探测成功，以点B₁+Δ₁作为临时矢点，并记作T₁₁，T₁₁的第一个下标1表示建立第一个模矢，第二个下标1表示X₁已被摄动；

若ε(B₁+Δ₁)＞ε(B₁)，则本次探测失败，进行反方向探测，即：试验B₁-Δ₁点，若ε(B₁-Δ₁)＜ε(B₁)，则本次探测成功，以点B₁-Δ₁作为临时矢点，并记作T₁₁；否则，仍以B₁作为临时矢点，并记作T₁₁；上述搜索可用公式(9)进行描述，

步骤3.4，当B₁的第1个分量X₁被摄动时，则用T₁₁代替原来的基点B₁，继续对B₁的第2个分量X₂，进行摄动，得临时矢点T₁₂；依此类推，直到对B₁的第8个分量X₈进行摄动，得到临时矢点T₁₈；

步骤3.5，在得到临时矢点T₁₈后，令T₁₈＝B₂，由初始基点B₁和新基点B₂构成第一个模矢；

步骤3.6，将第一个模矢延长一倍，得到第二个模矢的初始临时矢点T₂₀，即：

T₂₀＝B₁+2(B₂-B₁)＝2B₂-B₁ (10)

步骤3.7，同样，在T₂₀附近进行探索，依次得到临时矢点T₂₁,T₂₂,Λ,T₂₈，并以T₂₈为第三个基点B₃，此时由B₂、B₃构成第二个模矢；

将第二个模矢延长一倍，得到第三个模矢的初始临时矢点T₃₀，即：

T₃₀＝B₂+2(B₃-B₂)＝2B₃-B₂ (11)

步骤3.8，如此不断循环，若对于第i个模矢存在以下关系：

ε(T_i0)＜ε(B_i) (12)

则由T_i0产生不出比B_i更好的点，则应退回到B_i，模矢不再延长，并在B_i附近进行搜索；如果能得出新的下降点，即可引出新的模矢；否则，将步长缩小一半，开展更精细的探查，当步长缩小到足够高的精度时，即可停止迭代，此时得到的概率积分参数为全局最优解。

本发明提供的融合D-InSAR和模矢法求取概率积分参数的方法具有以下优点：

本发明有效融合了D-InSAR和模矢法，可求取全部概率积分参数，并且，还具有求取概率积分参数精度高的优点。

附图说明

图1为本发明提供的融合D-InSAR和模矢法求取概率积分参数的方法的流程示意图；

图2为本发明提供的时段2011.12.25-2012.03.11LOS方向累计地表移动变形量示意图；

图3为本发明提供的求参拟合效果图。

具体实施方式

为了使本发明所解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

发明人研究发现，基于单视向D-InSAR技术虽然无法直接监测开采沉陷三维移动变形，但LOS向移动变形本质为目标点的下沉分量、水平移动分量沿LOS向投影，理论上可根据这种关系模型，结合目标区域地质采矿条件，利用非线性规划理论反演概率积分参数。

鉴于此，发明人开展重点研究，本发明人结合采动变形规律和雷达几何成像特点，采用模矢法构建基于单视向D-InSAR的开采沉陷概率积分参数反演模型，可基于单视向D-InSAR技术反演全部概率积分参数，并且，反演得到的概率积分参数的准确度较高，研究成果对D-InSAR矿山变形监测具有重要参考价值。

具体的，结合图1，本发明提供的融合D-InSAR和模矢法求取概率积分参数的方法，包括以下步骤：

步骤1，根据地质采矿条件，预设工作面开采沉陷预计概率积分参数初值；

本步骤中，所预设的工作面开采沉陷预计概率积分参数初值包括左拐点偏移距S₁、右拐点偏移距S₂、上拐点偏移距S₃、下拐点偏移距S₄、下沉系数q、主要影响角正切tanβ、最大下沉角θ和水平移动系数b。

步骤2，通过雷达D-InSAR技术获取采动地表变形场影像；从所述采动地表变形场影像中选取n个目标像元；对于每个所述目标像元，均执行步骤2.1-步骤2.3：

步骤2.1，获取任意形状工作面地质采矿条件参数，包括：工作面尺寸D₁和D₃、采深H、采高m和煤层倾角α。

对工作面形状进行剖分，同时获取每个剖分单元的参数；此处具体为：首先获取工作面采深H，然后以边长为H/10的正方形为剖分单元，沿走向和倾向对工作面进行剖分，设任意一个剖分单元为剖分单元j，获取剖分单元j如下的剖分参数：单元面积A_j、单元中心采深H_j和单元中心坐标(x_j,y_j)；

另外，提取目标像元中心地理坐标(x,y)；

以步骤1的工作面开采沉陷预计概率积分参数初值、剖分单元的参数、目标像元中心地理坐标(x,y)为概率积分法模型的输入，利用所述概率积分法模型预计得到目标像元的预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值；

本步骤中，所建立的概率积分法模型为：

W_i＝[ΣW_j(x,y)·A_j]·W₀ (4)

U_iSN＝[r_j·(-2π/r_j ²)·(x-x_j)·W_j(x,y)·A_j]·b·W₀ (5)

其中：W_j(x,y)＝(1/r_j ²)·exp(-π·(x-x_j)²/r_j ²)·exp(-π·(y-y_j+l_j)²/r_j ²)、l_j＝H_jcotθ、r_j＝H_j/tanβ、W₀＝mqcosα。

步骤2.2，将步骤2.1计算得到的目标像元预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值，以及获取到的雷达卫星的入射角θ_i的值和卫星飞行方向方位角α_i的值，代入公式(1)，得到目标像元LOS向预计移动变形r'_iLOS的值；

步骤2.3，利用D-InSAR技术提取目标像元的LOS向实测移动变形r_iLOS的值，采用公式(2)计算得到目标像元的LOS向变形预计残差v_i的值；

v_i＝r_iLOS-r'_iLOS (2)

由此分别计算得到n个目标像元的LOS向变形预计残差v_i的值；

步骤3，构筑求参误差函数ε(B)：

其中，B为全部概率积分参数组成的矩阵；

将所计算得到的n个目标像元的LOS向变形预计残差v_i的值代入公式(3)，计算得到误差函数值；基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数。

本步骤中，基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数具体为：

步骤3.1，令B₁＝B，并以B₁作为初始基点进行搜索；

其中：

B₁所包含的各个预计参数即为B₁的分量，值为步骤1所设置的值；

步骤3.2，确定B₁的各个分量的探索步长；

即：对于任意分量X_i0(i＝1,2,...,8)，其步长为Δ_i，表示为：

Δ_i＝(0,…,ΔX_i0,…,0)^T (8)

步骤3.3，根据步骤1-步骤3，计算初始基点B₁的误差函数值ε(B₁)，另外计算B₁+Δ₁的误差函数值ε(B₁+Δ₁)；

比较ε(B₁)和ε(B₁+Δ₁)：

若ε(B₁+Δ₁)＜ε(B₁)，则探测成功，以点B₁+Δ₁作为临时矢点，并记作T₁₁，T₁₁的第一个下标1表示建立第一个模矢，第二个下标1表示X₁已被摄动；

若ε(B₁+Δ₁)＞ε(B₁)，则本次探测失败，进行反方向探测，即：试验B₁-Δ₁点，若ε(B₁-Δ₁)＜ε(B₁)，则本次探测成功，以点B₁-Δ₁作为临时矢点，并记作T₁₁；否则，仍以B₁作为临时矢点，并记作T₁₁；上述搜索可用公式(9)进行描述，

步骤3.4，当B₁的第1个分量X₁被摄动时，则用T₁₁代替原来的基点B₁，继续对B₁的第2个分量X₂，进行摄动，得临时矢点T₁₂；依此类推，直到对B₁的第8个分量X₈进行摄动，得到临时矢点T₁₈；

步骤3.5，在得到临时矢点T₁₈后，令T₁₈＝B₂，由初始基点B₁和新基点B₂构成第一个模矢；

步骤3.6，将第一个模矢延长一倍，得到第二个模矢的初始临时矢点T₂₀，即：

T₂₀＝B₁+2(B₂-B₁)＝2B₂-B₁ (10)

步骤3.7，同样，在T₂₀附近进行探索，依次得到临时矢点T₂₁,T₂₂,Λ,T₂₈，并以T₂₈为第三个基点B₃，此时由B₂、B₃构成第二个模矢；

将第二个模矢延长一倍，得到第三个模矢的初始临时矢点T₃₀，即：

T₃₀＝B₂+2(B₃-B₂)＝2B₃-B₂ (11)

步骤3.8，如此不断循环，若对于第i个模矢存在以下关系：

ε(T_i0)＜ε(B_i) (12)

则由T_i0产生不出比B_i更好的点，则应退回到B_i，模矢不再延长，并在B_i附近进行搜索；如果能得出新的下降点，即可引出新的模矢；否则，将步长缩小一半，开展更精细的探查，当步长缩小到足够高的精度时，即可停止迭代，此时得到的概率积分参数为全局最优解。

由此可见，本发明克服了单视向D-InSAR技术由于不能准确获取采动地表水平移动，导致无法求取与水平移动相关的水平移动参数的难题，本发明利用LOS向移动本质为目标点的下沉分量、水平移动分量沿LOS向投影之和的关系，基于模矢法提出了一种融合D-InSAR和模矢法求取概率积分参数的方法，可基于单视向D-InSAR技术反演全部概率积分参数，并且，反演得到的概率积分参数的精确度较高。

利用本发明提出的融合D-InSAR和模矢法求取概率积分参数的方法，对兖州南屯煤矿9310工作面的开采沉陷规律积分参数求取进行了应用实验。图2为获取的2011.12.25-2012.03.11时段LOS方向累计地表移动变形量，根据本发明方法求取2012.3.11日的工作面开采沉陷的动态概率积分参数为：q＝0.170、b＝0.12、tanβ＝2.12、θ＝88°、S₁＝-11、S₂＝10、S₃＝-90、S₄＝33，拟合误差约在-40～40mm之间，拟合中误差m＝±14.67mm，拟合效果图如图3所示。由此证明了本发明求取的概率积分参数的精确度高的优点。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (1)

1.一种融合D-InSAR和模矢法求取概率积分参数的方法，其特征在于，包括以下步骤：

步骤1，根据地质采矿条件，预设工作面开采沉陷预计概率积分参数初值；

步骤2，通过雷达D-InSAR技术获取采动地表变形场影像；从所述采动地表变形场影像中选取n个目标像元；对于每个所述目标像元，均执行步骤2.1-步骤2.3：

步骤2.1，获取任意形状工作面地质采矿条件参数，并对工作面形状进行剖分，同时获取每个剖分单元的参数；另外，提取目标像元中心地理坐标(x,y)；以步骤1的工作面开采沉陷预计概率积分参数初值、剖分单元的参数、目标像元中心地理坐标(x,y)为概率积分法模型的输入，利用所述概率积分法模型预计得到目标像元的预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值；

步骤2.2，将步骤2.1计算得到的目标像元预计下沉值W_i、南北方向预计水平移动U_iSN和东西方向预计水平移动U_iEW的值，以及获取到的雷达卫星的入射角θ_i的值和卫星飞行方向方位角α_i的值，代入公式(1)，得到目标像元LOS向预计移动变形r'_iLOS的值；

步骤2.3，利用D-InSAR技术提取目标像元的LOS向实测移动变形r_iLOS的值，采用公式(2)计算得到目标像元的LOS向变形预计残差v_i的值；

v_i＝r_iLOS-r'_iLOS (2)

由此分别计算得到n个目标像元的LOS向变形预计残差v_i的值；

步骤3，构筑求参误差函数ε(B)：

其中，B为全部概率积分参数组成的矩阵；

将所计算得到的n个目标像元的LOS向变形预计残差v_i的值代入公式(3)，计算得到误差函数值；基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数；

其中，步骤1中，所预设的工作面开采沉陷预计概率积分参数初值包括左拐点偏移距S₁、右拐点偏移距S₂、上拐点偏移距S₃、下拐点偏移距S₄、下沉系数q、主要影响角正切tanβ、最大下沉角θ和水平移动系数b；

其中，步骤2.1中，对工作面形状进行剖分，同时获取每个剖分单元的参数，具体为：首先获取工作面采深H，然后以边长为H/10的正方形为剖分单元，沿走向和倾向对工作面进行剖分，设任意一个剖分单元为剖分单元j，获取剖分单元j如下的剖分参数：单元面积A_j、单元中心采深H_j和单元中心坐标(x_j,y_j)；

所述获取任意形状工作面地质采矿条件参数包括：工作面尺寸D₁和D₃、采深H、采高m和煤层倾角α；

其中，步骤2.1中，所述概率积分法模型为：

W_i＝[∑W_j(x,y)·A_j]·W₀ (4)

U_iSN＝[r_j·(-2π/r_j ²)·(x-x_j)·W_j(x,y)·A_j]·b·W₀ (5)

其中：W_j(x,y)＝(1/r_j ²)·exp(-π·(x-x_j)²/r_j ²)·exp(-π·(y-y_j+l_j)²/r_j ²)、

l_j＝H_jcotθ、r_j＝H_j/tanβ、W₀＝mqcosα；

其中，步骤3中，基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数具体为：

步骤3.1，令B₁＝B，并以B₁作为初始基点进行搜索；

其中：

B₁所包含的各个预计参数即为B₁的分量，值为步骤1所设置的值；

步骤3.2，确定B₁的各个分量的探索步长；

即：对于任意分量X_i0(i＝1,2,...,8)，其步长为△_i，表示为：

△_i＝(0,…,△X_i0,…,0)^T (8)

步骤3.3，根据步骤1-步骤3，计算初始基点B₁的误差函数值ε(B₁)，另外计算B₁+△₁的误差函数值ε(B₁+△₁)；

比较ε(B₁)和ε(B₁+△₁)：

若ε(B₁+△₁)<ε(B₁)，则探测成功，以点B₁+△₁作为临时矢点，并记作T₁₁，T₁₁的第一个下标1表示建立第一个模矢，第二个下标1表示X₁已被摄动；

若ε(B₁+△₁)>ε(B₁)，则本次探测失败，进行反方向探测，即：试验B₁-△₁点，若ε(B₁-△₁)<ε(B₁)，则本次探测成功，以点B₁-△₁作为临时矢点，并记作T₁₁；否则，仍以B₁作为临时矢点，并记作T₁₁；上述搜索可用公式(9)进行描述，

步骤3.4，当B₁的第1个分量X₁被摄动时，则用T₁₁代替原来的基点B₁，继续对B₁的第2个分量X₂，进行摄动，得临时矢点T₁₂；依此类推，直到对B₁的第8个分量X₈进行摄动，得到临时矢点T₁₈；

步骤3.5，在得到临时矢点T₁₈后，令T₁₈＝B₂，由初始基点B₁和新基点B₂构成第一个模矢；

步骤3.6，将第一个模矢延长一倍，得到第二个模矢的初始临时矢点T₂₀，即：

T₂₀＝B₁+2(B₂-B₁)＝2B₂-B₁ (10)

步骤3.7，同样，在T₂₀附近进行探索，依次得到临时矢点T₂₁,T₂₂,…,T₂₈，并以T₂₈为第三个基点B₃，此时由B₂、B₃构成第二个模矢；

将第二个模矢延长一倍，得到第三个模矢的初始临时矢点T₃₀，即：

T₃₀＝B₂+2(B₃-B₂)＝2B₃-B₂ (11)

步骤3.8，如此不断循环，若对于第i个模矢存在以下关系：

ε(T_i0)<ε(B_i) (12)

则由T_i0产生不出比B_i更好的点，则应退回到B_i，模矢不再延长，并在B_i附近进行搜索；如果能得出新的下降点，即可引出新的模矢；否则，将步长缩小一半，开展更精细的探查，当步长缩小到足够高的精度时，即可停止迭代，此时得到的概率积分参数为全局最优解。