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

摘要

本发明提供一种融合D‑InSAR和模矢法求取概率积分参数的方法，包括：(1)利用概率积分法得到目标像元的预计下沉值、南北方向预计水平移动和东西方向预计水平移动的值，基于三维变形与LOS向变形关系，计算得到目标像元LOS向预计移动变形r'

说明书

技术领域

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

背景技术

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

解决上述难题，目前采用的主要方法为：建立适用于开采沉陷特点的 D-InSAR三维变形监测方法，然后直接进行求参。现有技术出现的解决方案主 要有：多视向D-InSAR、D-InSAR+Offsettracking、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的值；

${r^{\prime }}_{iLOS}=-U_{iSN}{\mathrm{sin\theta }}_{i}cos({\alpha }_i-\frac{3}{2}\pi )-U_{iEW}{\mathrm{sin\theta }}_i\sin ({\alpha }_i-\frac{3}{2}\pi )+W_i{\mathrm{cos\theta }}_i---\left(1\right)$

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

v_i＝r_iLOS-r'_iLOS(2)

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

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

$ϵ\left(B\right)=\underset{i=1}{\overset{n}{\Sigma }}\left|v_i\right|---\left(3\right)$

其中，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)

$\begin{array}{c}{U}_{iEW}=[r_j·(-2\pi /{r_j}^2)·(y-y_j+l_j)·W_j(x,y)·A_j\\ +W_j(x,y)·A_j·ctg\theta ]·b·W_0\end{array}---\left(6\right)$

其中：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₁作为初始基点进行搜索；

其中：

$\begin{array}{c}{B}_1={(S_{10},S_{20},S_{30},S_{40},q_0,{\mathrm{tan\beta }}_0,{\theta }_0,b_0)}^T\\ ={(X_{10},X_{20},X_{30},X_{40},X_{50},X_{60},X_{70},X_{80})}^T\end{array}---\left(7\right)$

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)

$\begin{array}{c}{U}_{iEW}=[r_j·(-2\pi /{r_j}^2)·(y-y_j+l_j)·W_j(x,y)·A_j\\ +W_j(x,y)·A_j·ctg\theta ]·b·W_0\end{array}---\left(6\right)$

其中：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的值；

${r^{\prime }}_{iLOS}=-U_{iSN}{\mathrm{sin\theta }}_{i}cos({\alpha }_i-\frac{3}{2}\pi )-U_{iEW}{\mathrm{sin\theta }}_i\sin ({\alpha }_i-\frac{3}{2}\pi )+W_i{\mathrm{cos\theta }}_i---\left(1\right)$

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

v_i＝r_iLOS-r'_iLOS(2)

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

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

$ϵ\left(B\right)=\underset{i=1}{\overset{n}{\Sigma }}\left|v_i\right|---\left(3\right)$

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

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

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

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

其中：

$\begin{array}{c}{B}_1={(S_{10},S_{20},S_{30},S_{40},q_0,{\mathrm{tan\beta }}_0,{\theta }_0,b_0)}^T\\ ={(X_{10},X_{20},X_{30},X_{40},X_{50},X_{60},X_{70},X_{80})}^T\end{array}---\left(7\right)$

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所示。由此证明了本发明求取 的概率积分参数的精确度高的优点。

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