一种钻井水溶盐矿开采沉陷InSAR预计方法与流程

本发明涉及预防盐矿开采沉陷的技术领域，特别地，涉及一种钻井水溶盐矿开采沉陷insar预计方法。

背景技术：

我国岩盐矿产资源丰富，其中岩盐开采90％以上采用钻井水溶法开采，即通过钻井注入溶剂，溶解矿层硐室，生成富含开采矿物的卤水，再将其从钻井中抽出。在开采中因为多井溶腔相连，多进行两井或者多井同时开采。水溶开采法与传统煤矿巷道开采不同，传统煤矿主要采用井工开采定向推进方式，而水溶法则为多方向推进，其采深比煤矿开采深度大。由于卤水溶解的时间和采深较大两方面因素，水溶法开采地表沉降在时间上往往呈现出滞后。随着水溶开釆釆厚加大，形成的采空区导致上覆盐层发生下沉，其对地表的破坏不仅仅是地表塌陷导致的矿区岩体结构、地表地形、构筑物等破坏，还极易引发地下卤水涌出，导致附近土地盐碱化，对环境破坏程度极为惊人。因此，水溶开采岩盐矿山的地面沉降需进行长期监测与预计。

合成孔径雷达差分干涉测量(dinsar：differentialinterferometricsyntheticapertureradar)技术通过利用覆盖矿区的sar影像，将稳定散射特性的高质量点提取出来建立形变模型，解算形变速率参数进而实现时序形变反演，理论上可达到亚毫米级的监测精度。将时序insar技术监测获取的时序形变与静态概率积分模型结合，可反演出概率积分模型预计参数，预计出矿区开采沉陷盆地中任一点的沉降量，目前已在传统煤矿开采沉陷预计中得到广泛应用。然而，概率积分模型实为一种静态沉陷预计模型，仅局限于生成某一特定时刻的矿山开采沉陷盆地，无法实现时序沉降的预计。且其应用多局限于煤矿开采导致的沉陷盆地预计。因此，迫切需要一种适用于水溶开采岩盐矿山的动态沉陷预计方法，以辅助保障岩盐矿山安全生产和健康的生态环境。

技术实现要素：

本发明目的在于提供一种钻井水溶盐矿开采沉陷insar预计方法，以解决目前对于盐矿开采沉陷监控预测不够准确，易产生安全隐患的技术问题。

为实现上述目的，本发明提供了一种钻井水溶盐矿开采沉陷insar预计方法，包括以下步骤：

a、利用时序insar技术反演盐矿区地表时序沉降，获取雷达视线向时序形变量；

b、提取钻井布设位置附近的可靠形变样本点；

c、构建传统静态概率积分法模型，所述模型表示如下式：

上式中，w(x，y)为任一点上矿山开采引发的地表形变量，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θ0为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距；

计算高相干点垂直方向的沉降量，雷达视线向形变量与沉降量的关系为：

wlos＝wvcosθ(3)

式中，wlos为雷达视线向形变量，wv为沉降量，θ为雷达卫星入射角；

d、利用遗传算法估计模型的未知参数：步骤a、步骤c计算出的矿区高相干点沉降量为输入值，根据残差最小原则，建立适应度函数f：f＝||wv-w′v||；式中：wv为可靠样本点沉降集，wv为概率积分预测模型计算的沉降量；利用遗传算法搜索获取满足适应度函数最小的参数集；

e、构建weibull时间函数动态概率积分模型，该模型为下式所示：

式中，w(x，y，tj)为观测点tj时刻矿区沉降量，ti、tj分别为起始时间、终止时间，ve、te分别为矿山开采速度、开采时间，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θ0为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距；

根据公式(8)、(9)已求得的概率积分法模型参数值，计算出各高相干点的概率积分法模型沉降量；根据各时间量与其对应的沉降量，建立关于待求参数c、k的方程组，并计算出c、k结果值，其中c为时间模型系数，k为待求动态沉陷预计参数；

f、将计算得到的q、tanβ、θ0、s1、s2、s3、c、k参数结果值，代入weibull时间函数动态概率积分模型，计算出各点位任意时刻的沉降量，并预计矿区时序形变场。

本发明具有以下有益效果：

本发明提供了一种钻井水溶盐矿开采沉陷insar预计方法，首先利用时序insar技术获取盐矿钻井分布线上可靠样本点的沉降序列；再反演静态概率积分模型预计参数和weibull时间函数模型的时间参数，构建动态开采沉陷预计模型；实现盐矿地表任意点在任意时刻的开采沉陷动态预计。不仅克服了传统概率积分法模型多应用于煤矿区的局限性，且根据步骤5与公式8、公式9，将weibull时间函数模型与传统概率积分法模型结合，构建出weibull时间函数动态概率积分模型，通过该模型可预测出矿区任意时刻沉降场，克服了传统概率积分法模型局限于静态下沉盆地开采沉陷监测的缺陷，拓宽了insar技术在矿区沉陷预计中的应用空间，有利于保障钻井水溶岩盐矿山的生态环境及生产安全。

除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。

附图说明

构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：

图1是本发明优选实施例的流程图；

图2是本发明优选实施例的矿山年平均沉降速率模拟图，图中白色点为高相干点；

图3是本发明优选实施例的模拟矿山沉降量数据图(参考时间为2016年2月10日)；

图4是本发明优选实施例的静态概率积分法模型计算所得沉降场图(参考时间为2016年2月10日)；

图5是本发明优选实施例的weibull时间函数动态概率积分法模型计算所得沉降场图(参考时间为2016年2月10日)；

图6是本发明优选实施例的高相干点动态模型预测值与真实值对比图；

图7是本发明优选实施例的weibull时间函数动态概率积分法模型预测值残差中误差图。

具体实施方式

以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。

参见图1，本方法操作步骤如下：

步骤1：利用时序insar技术反演盐矿区地表时序沉降。该过程主要包括：选取超级主影像、配准、干涉图生成、去平地效应、去轨道、滤波、相位解缠、提取高相干点与奇异值分解法解算形变速率，并通过时间维积分获取雷达视线向时序形变量。

步骤2：提取钻井布设位置附近的可靠形变样本点。由于受到insar技术相位梯度的限制，insar获取的大部分大变形相干点的形变结果并不可靠，无法用于后续参数的反演。在此综合考虑样本点的散射特性,样本点的地理分布特征及变形结果的精度来构建样本点集可靠性的综合评价指标，以这一指标值作为可靠样本点的选取标准，提取分布于钻井位置附近的可靠样本点集。

步骤3：构建传统静态概率积分法模型。模型可表示如下：

上式中，w(x，y)为任一点上矿山开采引发的地表形变量，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θ0为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距。

计算高相干点垂直方向的沉降量，雷达视线向形变量与沉降量的关系为：

wlos＝wvcosθ(3)式中，wlos为雷达视线向形变量，wv为沉降量，θ为雷达卫星入射角。

步骤4：利用遗传算法估计模型的未知参数。

概率积分法模型的未知参数包括q、tanβ、θ0、s1、s2、s3其余为已知参数。步骤1、步骤3计算出的矿区高相干点沉降量为输入值，根据残差最小原则，建立适应度函数

f：f＝||wv-w′v||；式中：wv为可靠样本点沉降集，w′v为概率积分预测模型计算的沉降量。利用遗传算法搜索获取满足适应度函数最小的参数集。

步骤5：构建weibull时间函数动态概率积分模型，并解算动态预计参数。

weibull时间序列函数关系式为：

式中：t为时间间隔，c为时间模型系数，k为待求动态沉陷预计参数。

weibull函数与沉降量的函数关系为：

式中：w(t)为该观测点t时刻的下沉量，wm为其漏斗最大下沉量(mm)。

结合式(1)、式(2)、式(4)、式(5)可得weibull时间函数动态概率积分法模型为：

式中，w(x，y，tj)为观测点tj时刻矿区沉降量，ti、tj分别为起始时间、终止时间，ve、te分别为矿山开采速度、开采时间，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θo为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距。

根据已求得的概率积分法模型参数值，可计算出各高相干点的概率积分法模型沉降量。根据各时间量与其对应的沉降量，建立关于待求参数c、k的方程组，并计算出c、k结果值。再利用残差最小原理对c、k结果值进行优化，得到最优c、k值。

步骤6：weibull时间函数动态概率积分模型预测后期矿山形变。将计算得到的q、tanβ、θ0、s1、s2、s3、c、k参数结果值，代入weibull时间函数动态概率积分模型，计算出各点位任意时刻的沉降量，并预计矿区时序形变场。

参见图2至图7所示的实施例，选取sentinel-1a、c波段、升轨模式参数，共模拟生成9幅差分干涉图(3/5/2016-28/2/2017)，进行模拟实验得年平均沉降速率见图2，图中白色点为高相干点。在模拟矿区累积沉降量时，加入均值为3、标准差为0.2的噪声数据。设置下沉系数q取值范围为：[0.01，1]；走向主要影响角正切tanβ为[1，3.8]；走向左拐点偏移距为s3为[0.05h，0.h]；倾向下山、倾向上山拐点偏移距分别为s1、s2，取值为[0.05h，0.3h]；开采影响传播角θ0＝90°-kα；其中α为矿层倾角，取值为[0～45°]；k为取值范围为[0.5，0.8]。随机生成200个点作为研究区内高相干点。图3为模拟的矿区时序沉降场。选取前4景反演概率积分法模型与weibull时间函数中的待求参数值，图5用于检验weibull时间函数动态概率积分法模型获取的预计沉降场质量；图4为利用概率积分模型获取的静态沉降场；图5为weibull时间函数动态概率积分法模型获取的沉降场；图6为高相干点真值与weibull时间函数动态概率积分法模型计算值的均方根误差。由图6可见，在加入一定的噪声前提下，本发明获取的预测沉降量与模拟真值吻合较好。

本发明的具体实施步骤为：

步骤1：时序insar技术反演盐矿区地表时序沉降。

以时间序列insar技术中小基线集技术(sbas-insar：smallbaselinesubset)为例，首先选取超级主影像，设置合适的时空基线阈值，进行干涉基线组合，并对任一干涉对进行配准、干涉图生成、去平地效应、去轨道、滤波和相位解缠(本实施例中模拟生成9景解缠后的差分干涉图)。通过设置相干性与振幅离差指数阈值提取矿区的高相干点(本实施例模拟生成200个高相干点)。利用高相干点相位值构建相位，形变速率的函数模型，并利用奇异值分解法解算出各时段形变速率，再通过时间维积分获取雷达视线向时序形变量(本实施例中利用模拟生成的沉降速率场和真实sar数据的时间基线、轨道参数共同模拟生成时序沉降场，见图3)。

步骤2：提取钻井布设位置附近的可靠形变样本点。

水溶盐矿钻井沿主断面的走向和倾向布设。由于受到insar技术相位梯度的限制，insar获取的大部分大变形相干点的形变结果并不可靠，无法用于后续参数的反演。在此综合考虑样本点的散射特性，样本点的地理分布特征及变形结果的精度来构建样本点集的可靠性综合评价指标，以这一指标值作为可靠样本点的选取标准，提取分布于钻井位置附近的可靠样本点集。具体步骤如下：

步骤a：判断insar可检测的最大最小变形临界值，根据矿区内开展的大地测量变形先验信息，划分形变图上的大变形区域，这一区域内的高相干点目标与小变形区域的相干点目标进行分类处理。大变形区域内高相干目标的insar形变采用偏移量估计方法进行获取，这样可保留大变形样本点的形变信息，小变形区域的样本点则采用时间序列方法进行获取；

步骤b：在可检测的变形区域范围内，分别计算各相干点的点位稳定性参数，包括：相干系数γ，幅度时间标准差σamp，相位时间标准差σpha，相干系数时间标准差σγ，由这些参数共同构建点位稳定性综合评定指标，分别以下式依次计算各相干点的σpoi值：

步骤c：计算相干点的变形稳定性指标σdef。在时间序列insar雷达相位方程组变形速率求解过程中，应用lamdba方法实现变形速率参数的求解后，可依据最小二乘原理计算出各高相干点变形结果的中误差，在此用将其作为σdef值；

步骤d：由于样本点集在地理位置上的分布为随机，在此结合钻井水溶盐矿的开采特点及钻井的分布特征设置样本点的地理分布指标σdis，原则为：与矿区走向和倾向分布越接近的样本点σdis取值越接近1，此类分布的样本点可作为开采沉陷的一级观测线，并利用移动角法来计算观测线的长度；对于样本点集中超出观测线长度限制的点剔除；在提取出的一级观测线两端选取稳定区域，落在此区的样本点可选为不受采动影响的稳定参考点，σdis取为0；每一端一般选取不少于3个点，另保证稳定点与一级观测线上的监测点间距在50m以内；

步骤e：依据点位稳定性综合评定指标σpoi值，点位变形稳定性指标σdef值和点的地理分布指标σdis，以σpoi值作初选，初选时对步骤a中的大变形样本点集和小变形区域分开进行处理，再使用σdef值二次筛选；最后根据σdis取值筛选样本点集，保证样本点沿主断面附近分布。经过这样的筛选，最终可构建出概率积分参数反演算法的输入样本集。

步骤3：建立静态概率积分预计模型。

上式中，w(x，y)为任一点上矿山开采引发的地表形变量，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θ0为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距；。其中，参数m、h、h1、h2、α、d3、d1可根据矿山实际开采情况计算得到，参数q、tanβ、θ0、s1、s2、s3为待求参数。

建立sbas-insar雷达方向形变量与概率积分法模型关系，由步骤1计算得到矿区高相干点的雷达视线向形变量。因水溶法开采矿山区域主要由采场内水压力支撑地表，导致其地表形变主要以沉降为主，地表水平移动微小。进而雷达视线向形变量与沉降量的关系可表示为：

wlos＝wvcosθ(3)

式中，wlos为雷达视线向形变量，wv为沉降量，θ为雷达卫星入射角；

w(x，y，ti)＝wv(x，y，ti)(10)

根据式(8)、式(9)可使计算出的sbas-insar测量的雷达方向形变量，与概率积分法模型联系起来。

步骤4：利用遗传算法求解静态概率积分法模型待求参数。具体步骤如下：

步骤a：设置种群个体个数为n，迭代次数为iter，分别设置染色体节点q、tanβ、θ0、s1、s2、s3的取值范围为下沉系数q，取值为：[0.01，1]；走向主要影响角正切tanβ为[1，3.8]；走向左拐点偏移距为s3为[0.05h，0.h]；倾向下山、倾向上山拐点偏移距分别为s1、s2，取值为[0.05h，0.3h]；开采影响传播角θ0＝90°-kα；其中α为矿层倾角，取值为[0～45°]；k为取值范围为[0.5，0.8]。根据其取值范围，以随机赋值方式，对染色体进行初始化。

步骤b：建立适应度函数f，如下式：

f＝||wv-w′v||(6)

式中wv为步骤1-2获取的可靠样本点沉降集，wv为式(1)模型计算的沉降量。计算每个初始化的个体染色体适应度f，设置阈值，并判断f是否达到阈值要求。若不满足，则对种群个体进行选择、交叉、变异等操作，产生新的种群个体，再次进行适应度判断，以此循环迭代，直至获得满足要求的种群个体，即为参数集(q、tanβ、s1、s2、s3、θ0)；

步骤c：优化遗传算法所求概率积分法模型参数值。

对于n个参数集，以适应度函数值作权重，计算其加权平均值，作为概率积分法模型参数值初始值goaf；根据goaf计算对应的研究区沉降值。根据模拟实验真值数据，计算沉降场a的残差中误差；设置残差阈值，若残差中误差不满足阈值要求，则依据残差值越大，权值越小的原则调整权值，并以新权值计算新的模型参数值goaf′，开始迭代，直到满足阈值条件，停止迭代。

步骤5：解算动态预计参数，并构建weibull时间函数动态概率积分模型。

weibull函数与沉降量的函数关系为：

式中，wm为其漏斗最大下沉量(mm)，t为时间间隔，c为时间模型系数，k为待求动态沉陷预计参数。

根据计算出的概率积分法模型参数q、tanβ、θ0、s1、s2、s3结果值，代入概率积分法模型中，可计算出各点的最大下沉量wm。将n个高相干点的时间量t，与对应的概率积分法模型计算所得沉降量w(t)作为已知量，建立n个方程求解出m个待求参数c、k。运用残差最小原理对m个待求参数进行优化，得到最优参数值。

结合传统概率积分模型与weibull时间函数关系，可构建weibull时间函数动态概率积分模型为：

式中，w(x，y，tj)为观测点tj时刻矿区沉降量，ti、tj分别为起始时间、终止时间，ve、te分别为矿山开采速度、开采时间，x、y为任意点坐标，m为开采厚度，q为下沉系数，α为矿层倾角，u为积分参数，d1为工作面倾向斜长，d3为工作面走向斜长，θ0为开采影响角，h、h1、h2分别为走向、倾向下山、倾向上山开采深度，tanβ、tanβ1、tanβ2分别为走向、倾向下山、倾向上山主要影响角正切，s1、s2分别为下山拐点偏移距、上山拐点偏移距，s3、s4分别为走向左、右拐点偏移距；

步骤6、将求得的概率积分法模型参数q、tanβ、θ0、s1、s2、s3与动态预计参数c、k代入weibull时间函数动态概率积分模型(式(8)、式(9))，根据研究区各点的坐标值x、y与时间量t，可计算出各点位任意时刻的沉降量，并预计矿区时序形变场。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。