本发明属于矿山变形监测技术领域，具体涉及一种融合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)为概率积分法模型的输入，利用所述概率积分法模型预计得到目标像元的预计下沉值Wi、南北方向预计水平移动UiSN和东西方向预计水平移动UiEW的值；
步骤2.2，将步骤2.1计算得到的目标像元预计下沉值Wi、南北方向预计水平移动UiSN和东西方向预计水平移动UiEW的值，以及获取到的雷达卫星的入射角θi的值和卫星飞行方向方位角αi的值，代入公式(1)，得到目标像元LOS向预计移动变形r'iLOS的值；
步骤2.3，利用D-InSAR技术提取目标像元的LOS向实测移动变形riLOS的值，采用公式(2)计算得到目标像元的LOS向变形预计残差vi的值；
vi＝riLOS-r'iLOS (2)
由此分别计算得到n个目标像元的LOS向变形预计残差vi的值；
步骤3，构筑求参误差函数ε(B)：
其中，B为全部概率积分参数组成的矩阵；
将所计算得到的n个目标像元的LOS向变形预计残差vi的值代入公式(3)，计算得到误差函数值；基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数。
优选的，步骤1中，所预设的工作面开采沉陷预计概率积分参数初值包括左拐点偏移距S1、右拐点偏移距S2、上拐点偏移距S3、下拐点偏移距S4、下沉系数q、主要影响角正切tanβ、最大下沉角θ和水平移动系数b。
优选的，步骤2.1中，对工作面形状进行剖分，同时获取每个剖分单元的参数，具体为：首先获取工作面采深H，然后以边长为H/10的正方形为剖分单元，沿走向和倾向对工作面进行剖分，设任意一个剖分单元为剖分单元j，获取剖分单元j如下的剖分参数：单元面积Aj、单元中心采深Hj和单元中心坐标(xj,yj)；
所述获取任意形状工作面地质采矿条件参数包括：工作面尺寸D1和D3、采深H、采高m和煤层倾角α。
优选的，步骤2.1中，所述概率积分法模型为：
Wi＝[∑Wj(x,y)·Aj]·W0 (4)
UiSN＝[rj·(-2π/rj2)·(x-xj)·Wj(x,y)·Aj]·b·W0 (5)
其中：Wj(x,y)＝(1/rj2)·exp(-π·(x-xj)2/rj2)·exp(-π·(y-yj+lj)2/rj2)、lj＝Hjcotθ、rj＝Hj/tanβ、W0＝mqcosα。
优选的，步骤3中，基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数具体为：
步骤3.1，令B1＝B，并以B1作为初始基点进行搜索；
其中：
B1所包含的各个预计参数即为B1的分量，值为步骤1所设置的值；
步骤3.2，确定B1的各个分量的探索步长；
即：对于任意分量Xi0(i＝1,2,...,8)，其步长为Δi，表示为：
Δi＝(0,…,ΔXi0,…,0)T (8)
步骤3.3，根据步骤1-步骤3，计算初始基点B1的误差函数值ε(B1)，另外计算B1+Δ1的误差函数值ε(B1+Δ1)；
比较ε(B1)和ε(B1+Δ1)：
若ε(B1+Δ1)＜ε(B1)，则探测成功，以点B1+Δ1作为临时矢点，并记作T11，T11的第一个下标1表示建立第一个模矢，第二个下标1表示X1已被摄动；
若ε(B1+Δ1)＞ε(B1)，则本次探测失败，进行反方向探测，即：试验B1-Δ1点，若ε(B1-Δ1)＜ε(B1)，则本次探测成功，以点B1-Δ1作为临时矢点，并记作T11；否则，仍以B1作为临时矢点，并记作T11；上述搜索可用公式(9)进行描述，
步骤3.4，当B1的第1个分量X1被摄动时，则用T11代替原来的基点B1，继续对B1的第2个分量X2，进行摄动，得临时矢点T12；依此类推，直到对B1的第8个分量X8进行摄动，得到临时矢点T18；
步骤3.5，在得到临时矢点T18后，令T18＝B2，由初始基点B1和新基点B2构成第一个模矢；
步骤3.6，将第一个模矢延长一倍，得到第二个模矢的初始临时矢点T20，即：
T20＝B1+2(B2-B1)＝2B2-B1 (10)
步骤3.7，同样，在T20附近进行探索，依次得到临时矢点T21,T22,Λ,T28，并以T28为第三个基点B3，此时由B2、B3构成第二个模矢；
将第二个模矢延长一倍，得到第三个模矢的初始临时矢点T30，即：
T30＝B2+2(B3-B2)＝2B3-B2 (11)
步骤3.8，如此不断循环，若对于第i个模矢存在以下关系：
ε(Ti0)＜ε(Bi) (12)
则由Ti0产生不出比Bi更好的点，则应退回到Bi，模矢不再延长，并在Bi附近进行搜索；如果能得出新的下降点，即可引出新的模矢；否则，将步长缩小一半，开展更精细的探查，当步长缩小到足够高的精度时，即可停止迭代，此时得到的概率积分参数为全局最优解。
本发明提供的融合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，根据地质采矿条件，预设工作面开采沉陷预计概率积分参数初值；
本步骤中，所预设的工作面开采沉陷预计概率积分参数初值包括左拐点偏移距S1、右拐点偏移距S2、上拐点偏移距S3、下拐点偏移距S4、下沉系数q、主要影响角正切tanβ、最大下沉角θ和水平移动系数b。
步骤2，通过雷达D-InSAR技术获取采动地表变形场影像；从所述采动地表变形场影像中选取n个目标像元；对于每个所述目标像元，均执行步骤2.1-步骤2.3：
步骤2.1，获取任意形状工作面地质采矿条件参数，包括：工作面尺寸D1和D3、采深H、采高m和煤层倾角α。
对工作面形状进行剖分，同时获取每个剖分单元的参数；此处具体为：首先获取工作面采深H，然后以边长为H/10的正方形为剖分单元，沿走向和倾向对工作面进行剖分，设任意一个剖分单元为剖分单元j，获取剖分单元j如下的剖分参数：单元面积Aj、单元中心采深Hj和单元中心坐标(xj,yj)；
另外，提取目标像元中心地理坐标(x,y)；
以步骤1的工作面开采沉陷预计概率积分参数初值、剖分单元的参数、目标像元中心地理坐标(x,y)为概率积分法模型的输入，利用所述概率积分法模型预计得到目标像元的预计下沉值Wi、南北方向预计水平移动UiSN和东西方向预计水平移动UiEW的值；
本步骤中，所建立的概率积分法模型为：
Wi＝[ΣWj(x,y)·Aj]·W0 (4)
UiSN＝[rj·(-2π/rj2)·(x-xj)·Wj(x,y)·Aj]·b·W0 (5)
其中：Wj(x,y)＝(1/rj2)·exp(-π·(x-xj)2/rj2)·exp(-π·(y-yj+lj)2/rj2)、lj＝Hjcotθ、rj＝Hj/tanβ、W0＝mqcosα。
步骤2.2，将步骤2.1计算得到的目标像元预计下沉值Wi、南北方向预计水平移动UiSN和东西方向预计水平移动UiEW的值，以及获取到的雷达卫星的入射角θi的值和卫星飞行方向方位角αi的值，代入公式(1)，得到目标像元LOS向预计移动变形r'iLOS的值；
步骤2.3，利用D-InSAR技术提取目标像元的LOS向实测移动变形riLOS的值，采用公式(2)计算得到目标像元的LOS向变形预计残差vi的值；
vi＝riLOS-r'iLOS (2)
由此分别计算得到n个目标像元的LOS向变形预计残差vi的值；
步骤3，构筑求参误差函数ε(B)：
其中，B为全部概率积分参数组成的矩阵；
将所计算得到的n个目标像元的LOS向变形预计残差vi的值代入公式(3)，计算得到误差函数值；基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数。
本步骤中，基于ε(B)＝min准则，利用模矢法搜索，获取全局最优概率积分预计参数具体为：
步骤3.1，令B1＝B，并以B1作为初始基点进行搜索；
其中：
B1所包含的各个预计参数即为B1的分量，值为步骤1所设置的值；
步骤3.2，确定B1的各个分量的探索步长；
即：对于任意分量Xi0(i＝1,2,...,8)，其步长为Δi，表示为：
Δi＝(0,…,ΔXi0,…,0)T (8)
步骤3.3，根据步骤1-步骤3，计算初始基点B1的误差函数值ε(B1)，另外计算B1+Δ1的误差函数值ε(B1+Δ1)；
比较ε(B1)和ε(B1+Δ1)：
若ε(B1+Δ1)＜ε(B1)，则探测成功，以点B1+Δ1作为临时矢点，并记作T11，T11的第一个下标1表示建立第一个模矢，第二个下标1表示X1已被摄动；
若ε(B1+Δ1)＞ε(B1)，则本次探测失败，进行反方向探测，即：试验B1-Δ1点，若ε(B1-Δ1)＜ε(B1)，则本次探测成功，以点B1-Δ1作为临时矢点，并记作T11；否则，仍以B1作为临时矢点，并记作T11；上述搜索可用公式(9)进行描述，
步骤3.4，当B1的第1个分量X1被摄动时，则用T11代替原来的基点B1，继续对B1的第2个分量X2，进行摄动，得临时矢点T12；依此类推，直到对B1的第8个分量X8进行摄动，得到临时矢点T18；
步骤3.5，在得到临时矢点T18后，令T18＝B2，由初始基点B1和新基点B2构成第一个模矢；
步骤3.6，将第一个模矢延长一倍，得到第二个模矢的初始临时矢点T20，即：
T20＝B1+2(B2-B1)＝2B2-B1 (10)
步骤3.7，同样，在T20附近进行探索，依次得到临时矢点T21,T22,Λ,T28，并以T28为第三个基点B3，此时由B2、B3构成第二个模矢；
将第二个模矢延长一倍，得到第三个模矢的初始临时矢点T30，即：
T30＝B2+2(B3-B2)＝2B3-B2 (11)
步骤3.8，如此不断循环，若对于第i个模矢存在以下关系：
ε(Ti0)＜ε(Bi) (12)
则由Ti0产生不出比Bi更好的点，则应退回到Bi，模矢不再延长，并在Bi附近进行搜索；如果能得出新的下降点，即可引出新的模矢；否则，将步长缩小一半，开展更精细的探查，当步长缩小到足够高的精度时，即可停止迭代，此时得到的概率积分参数为全局最优解。
由此可见，本发明克服了单视向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°、S1＝-11、S2＝10、S3＝-90、S4＝33，拟合误差约在-40～40mm之间，拟合中误差m＝±14.67mm，拟合效果图如图3所示。由此证明了本发明求取的概率积分参数的精确度高的优点。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。