CN115201822B - 一种钻井水溶岩盐矿区采卤量估计方法 - Google Patents

Abstract

本发明提供了一种钻井水溶岩盐矿区采卤量估计方法，包括基于InSAR生成相干系数图和差分干涉图，生成基于高相干点的时序差分干涉相位矩阵；构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型；对动态函数模型中的模型参数进行最优求解；估计InSAR影像干涉时段内的岩盐矿区地表下沉量的时序形变；通过岩盐矿区地表下沉量的时序形变获取总采卤量的估计值等步骤。本发明利用时序InSAR技术反演岩盐矿区任一时刻采卤量的计算方法，可实现采卤量重要组分的定量预计，拓宽了InSAR技术应用领域，可用于指导矿产资源的开采利用，进一步提升矿产资源的回采率，为盐类矿床的合理开发及利用提供数据参考。

Description

一种钻井水溶岩盐矿区采卤量估计方法

技术领域

本发明涉及岩盐矿区开采技术领域，具体涉及一种钻井水溶岩盐矿区采卤量估计方法。

背景技术

由于矿产资源为不可再生资源，所以地下资源的合理开发显得尤为重要。我国是世界上少数拥有极为丰富盐类矿产资源的国家之一，已探明的储量约44500亿吨，90％以上的岩盐矿区开采方式为全方位的钻井热水热溶开采方式。矿区卤水提取过程中回采率过低的问题严重制约了矿产资源的有效利用。因此，研究钻井水溶式开采矿区卤水的采卤量和卤水中重要组分的估计方法可为钻井水溶开采工艺的改进和回采率的提高提供指导，最大限度的开发和利用卤水矿产资源；为盐类矿床的合理开发与采卤量重要组分的提取提供参考；为政府和企业进行岩盐矿产资源的可持续发展，矿山环境的保护，以及矿产资源的经济价值分析提供重要依据。

传统矿区地下溶腔演化规律与卤水重要组分含量反演的方法有有限元数值迭代法、水化模拟及特征研究法、自动检测法等，这些方法能够较好的反应地下卤水资源的变化规律，但需要监测人员频繁出入勘测现场获取地表形变数据，且依赖于开采历史数据的统计，计算极为繁琐。获取地表形变数据的方法主要为水准测量法，全站仪测量和全球卫星导航(GPS/GNSS)等，虽然能满足矿区监测的精度要求，但成本较高，时空分辨率较差。同时，上述方法均无法实现对后期开采引起的地表形变及采卤量的预测。

综上所述，急需一种钻井水溶岩盐矿区采卤量估计方法以解决现有技术中存在的问题。

发明内容

本发明目的在于提供一种钻井水溶岩盐矿区采卤量估计方法，以解决由于采用钻井水溶开采方式导致的形变量和采卤量预测问题。

为实现上述目的，本发明提供了一种钻井水溶岩盐矿区采卤量估计方法，包括以下步骤：

步骤一：基于InSAR生成相干系数图和差分干涉图，生成基于高相干点的时序差分干涉相位矩阵；

步骤二：构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型；

步骤三：对动态函数模型中的模型参数进行最优求解；

步骤四：估计InSAR影像干涉时段内的岩盐矿区地表下沉量的时序形变；

步骤五：通过岩盐矿区地表下沉量的时序形变获取总采卤量的估计值。

优选的，所述步骤一具体包括以下分步骤：

1.1、将多景SAR卫星数据批量格式转换，生成单视复数影像，并进行初始影像的预滤波，减弱影像的噪声相位；

1.2、对所有影像进行基线估计，然后选取一景影像作为所有影像的超级主影像，设置时空阈值，生成所有符合设定的干涉对；利用生成的干涉对生成干涉图，在此过程中附带生成相干系数图和强度图；

1.3、对干涉图进行以下操作：借助精密轨道数据进行去除平地相位、轨道误差；利用DEM数据去除地形相位；然后进行残余相位滤波；生成差分干涉图，对差分干涉图进行相位解缠，生成解缠后的时序差分干涉相位；

1.4、利用强度、振幅离差和相干系数法提取高相干点，并在时序差分干涉相位中提取高相干点的对应值，生成基于高相干点的时序差分干涉相位矩阵。

优选的，所述步骤二中：时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型如表达式6)所示：

其中，

是第m个差分干涉图中第i个高相干点在t_A时刻和t_B时刻两景SAR图像上解缠后生成的差分干涉相位，△t＝t_B-t_A，λ为雷达波长，D为扩散系数，C_s为溶液饱和浓度，C_x为溶液未达到饱和状态时的浓度，δ为饱和溶液层厚度，θ为雷达入射角，B_m为第m个干涉对垂直基线的长度，R_p为相干目标与雷达卫星所在位置之间的距离，δH_i为第i个高相干点对应的高程改正值，w_i ^m为第m个差分干涉图中第i个高相干点对应的残差相位。

优选的，所述步骤二中，t时刻对应的由水溶开采导致的沉降量S_sol(t)通过表达式3)计算：

S_sol(t)＝(D(C_s-C_x)/C_sδ)·t 3)。

优选的，所述步骤三中，根据最小二乘平差理论，利用附有不等式约束条件的动态优化迭代算法模型对模型参数X＝[D，C_x，δH_i]进行最优求解；构建的动态优化迭代算法模型如表达式7)所示：

其中，φ(X)为关于X的平差准则，V为m×1的改正数矩阵，V^T为V的转置矩阵，P为最小二乘原理中的权阵，W₁、W₂和W₃分别为C_x、D、δH_i对应的限值。

优选的，所述步骤三中，将动态优化迭代算法模型转换为步骤二中动态函数模型的虚拟观测方程，对InSAR观测相位矩阵B增加扰动，此时X的正则化解如表达式8)所示：

X＝[I-[B^TB+μI]^-1μI]^-1[B^TB+μI]^-1B^TL8)；

其中，I为单位矩阵，B^T为B的转置矩阵，μ为正则参数，L为m×1的观测值。

优选的，所述步骤四中，将步骤三中求解得到的模型参数代入动态函数模型中，并通过表达式4)计算雷达视线向形变d_Los：

d_Los＝d_Los(t_B)-d_Los(t_A)＝[S_sol(t_B)-S_sol(t_A)]·cosθ4)；

其中，d_Los(t_B)和d_Los(t_A)分别是t_A时刻和t_B时刻相对于开始时间t ₀＝0沿雷达视线方向的形变；S_sol(t_A)和S_sol(t_B)分别为t_A时刻和t_B时刻对应的由水溶开采导致的沉降量。

优选的，所述步骤四中，根据大气延迟相位和噪声相位在时空域上的特性，对动态函数模型中的w_i ^m进行空间维低通滤波和时间维高通滤波，从而提取残余相位高通形变分量d_HP，将雷达视线向形变d_Los和高通形变分量d_HP累加，实现岩盐矿区中任意一点的地表下沉量S_i的时序形变反演；对S_i进行地理编码，生成岩盐矿区地表垂直向时序形变，从而获取地表变形规律。

优选的，所述步骤五中，t时刻的总采卤量Q通过表达式9)进行计算：

其中，k为一景影像中高相干点的总数，S_i(t)为t时刻对应的地表下沉量，F₀为由雷达获取的高相干点的像素在实际地理坐标中大小；

卤水中盐类重要组分含量T_y通过表达式10)计算：

T_y＝Q×g10)；

其中，g为卤水中主要组分平均含量。

优选的，所述步骤五中，溶蚀通道直径d通过表达式13)计算：

其中，L_S为井组间距，H为开采矿层厚度，G为矿石体重，E为矿石平均品位；

溶腔发展过程中任意一井组连接方式的溶采总面积F通过表达式11)进行计算：

应用本发明的技术方案，具有以下有益效果：

(1)本发明中，通过将InSAR技术与水溶动力学参数结合，构建动态函数模型，对动态函数模型中的模型参数进行最优求解后，可通过动态函数模型估计InSAR影像干涉时段内的岩盐矿区地表下沉的时序形变，并可通过岩盐矿区地表下沉的时序形变获取总采卤量的估计值。其中，InSAR具有监测范围广、空间分辨率高、非接触式测量等优点，大大弥补了传统测量手段的不足，为地表形变的研究与采卤量估计提供了更为有效的监测手段；本申请利用时序InSAR技术反演岩盐矿区任一时刻采卤量的计算方法，可实现采卤量重要组分的定量预计，拓宽了InSAR技术应用领域，可用于指导矿产资源的开采利用，进一步提升矿产资源的回采率，为盐类矿床的合理开发及利用提供数据参考。

(2)本发明中，通过建立由水溶开采导致的沉降量与水溶动力学参数和时间之间的函数关系，再通过雷达视线向形变与由水溶开采导致的沉降量之间的对应关系，构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型，可以建立雷达参数与水溶动力学参数之间的关系，进一步实现采卤量的估计和地表下沉量的时序反演。

(3)本发明中，根据InSAR干涉相位方程组的特征，利用附有不等式约束的最小二乘算法实现模型高效、稳健的参数估计，将动态优化迭代算法模型中的不等式约束条件转换成动态函数模型的虚拟观测方程，克服原方程的秩亏问题；对InSAR观测相位矩阵B中对角线上的各元素增加微小扰动，做B^TB+μI的处理，这样的处理可克服原法方程系数阵的病态问题，并能降低参数对初值的依赖性，可加速算法的收敛，最终获取模型参数的最优估值。

(4)本发明中，通过将雷达视线向形变d_Los和高通形变分量d_HP累加，可以实现岩盐矿区中任意一点的地表下沉量S_i的时序形变反演；对S_i进行地理编码，生成岩盐矿区地表垂直向时序形变，从而获取地表变形规律。

(5)本发明中，通过高相干点对应的地表下沉量结合由雷达获取的高相干点的像素在实际地理坐标中大小即可实现水溶岩盐矿区内采卤量的估计，并进一步实现卤水中盐类重要组分含量的定量反演。

(6)本申请中，通过总采卤量和卤水中盐类重要组分含量可获取溶蚀通道直径，进而获取地下溶蚀通道直径随时间的演化过程，合理探究地下溶腔的演化与采卤量重要组分含量变化之间的规律，可以指导矿产资源的开发利用与回采率的提升。

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

附图说明

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

图1是本申请实施例中一种钻井水溶岩盐矿区采卤量估计方法的流程图；

图2是通过模拟实验获取的水溶岩盐矿区时序形变结果图；

图3是本申请实施例中通过最优求解的模型参数估计值计算得到的水溶岩盐矿区时序形变结果图；

图4是高相干点上通过模拟实验获取的模拟值与通过本申请方法获取的预测值之间的对比图；

图5是高相干点上通过本申请方法获得的预测值相对于通过模拟实验获得的时序形变量模拟值的偏差的占比图；

图6是通过本申请方法反演的采卤量重要组分含量与溶蚀通道直径的变化规律的示意图。

具体实施方式

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

实施例：

参见图1至图6，一种钻井水溶岩盐矿区采卤量估计方法，本实施例应用于钻井水溶岩盐矿区中总采卤量的估计。

一种钻井水溶岩盐矿区采卤量估计方法，具体为一种基于合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)技术的钻井水溶岩盐矿区采卤量估计方法，如图1所示，包括以下步骤：

步骤一：基于InSAR生成相干系数图和差分干涉图，生成基于高相干点的时序差分干涉相位矩阵；

1.1、将多景SAR(合成孔径雷达)卫星数据(即N+1景SAR影像)批量格式转换，生成N+1景可由处理软件读取的单视复数影像(SLC)，并进行初始影像的预滤波，减弱影像的噪声相位；在格式转换的过程中，每景数据都会生成对应的参数文件；在GAMMA数据处理过程中，哨兵数据需要经过前期数据的精确配准。

1.2、对所有影像进行基线估计，然后选取一景影像作为所有影像的超级主影像，设置时空阈值，生成所有符合设定的干涉对；在此基础上，所有影像都会配准重采样到超级主影像；在配准的过程中，采用多项式拟合的方法；利用生成的干涉对生成干涉图，在此过程中附带生成相干系数图和强度图；

1.3、对干涉图进行以下操作：借助精密轨道数据进行去除平地相位、轨道误差；利用外部30m分辨率的DEM(数字高程模型)数据去除地形相位；然后进行残余相位滤波，去除影像中的噪声；生成差分干涉图，对差分干涉图进行相位解缠(部分干涉对进行3D解缠，无法进行3D解缠的利用最小费用流法进行解缠)，生成解缠后的时序差分干涉相位；

1.4、利用强度、振幅离差和相干系数法提取高相干点，并在时序差分干涉相位中提取高相干点的对应值，生成基于高相干点的时序差分干涉相位矩阵。

步骤二：构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型(即物理模型)；

步骤2.1、结合盐类矿床水溶开采沉陷机制，建立单位时间饱和卤水量与水溶开采动力学参数之间的函数关系。

为建立雷达视线向形变d_Los与开采水溶动力学参数之间的关联，首先需要将垂直向沉降量与开采期间总采卤量Q建立关联，如表达式1)所示：

Q＝Q₀×F₀×t＝S_sol×F₀1)；

其中，t为开采时间，可根据两景SAR影像的过境时间来进行计算，Q为影像干涉组合期间总采卤量，Q₀为单位时间内单位面积的采卤量，F₀为由雷达获取的高相干点的像素在实际地理坐标中大小，S_sol为由水溶开采导致的沉降量(垂直向)。本实施例中，采用Sentinel-1A卫星数据的的高相干点的像素在实际地理坐标中所代表的面积F₀为66.47m²。

由扩散溶解动力学方程可知，Q₀可由表达式2)进行求解：

Q₀＝D(C_s-C_x)/C_sδ2)；

则由表达式1)和表达式2)可建立由水溶开采导致的沉降量S_sol(t)与水溶动力学参数和时间之间的函数关系，如表达式3)所示：

S_sol(t)＝(D(C_S-C_x)/C_sδ)·t3)；

其中，S_sol(t)为t时刻对应的由水溶开采导致的沉降量，D为扩散系数，属于需要通过求解获取最优值的未知参数；C_s为溶液饱和浓度；C_x为t时刻对应的溶液未达到饱和状态时的浓度，在不同的影像时刻量值不同，也属于需要通过求解获取最优值的未知参数；δ为饱和溶液层厚度，取决于溶剂的运动速度，可视为常数，可根据研究区域地质情况查阅相关资料获取。

步骤2.2、构建InSAR雷达视线向形变d _Los与水溶动力学参数、时间之间的动态函数 模型，其中雷达视线向形变d _Los是基于SBAS-InSAR(短基线集时间序列InSAR分析)技术获取的，是利用测试区域内的高质量点进行建模分析，进而提取出形变分量，反演矿区地表形变的过程。

在水溶开采同一时刻同一溶腔内，基本的开采条件相同的情况下，主要受溶液在此时刻的浓度和溶解面与理想水平面夹角的大小影响。根据现有技术的试验资料已经充分表明，夹角为180度时(即上溶)的溶解速率是夹角为90度(侧溶)时的溶解速率的2倍，而夹角为0度(底溶)时溶解速率最小，接近于0。溶腔引起的地表形变主要沿垂直方向进行，当忽略水平移动时，雷达视线向形变d_Los如表达式4)所示：

d_Los＝d_Los(t_B)-d_Los(t_A)＝[S_sol(t_B)-S_sol(t_A)]·cosθ4)；

其中，d_Los(t_B)和d_Los(t_A)分别是t_A时刻和t_B时刻相对于开始时间t₀＝0沿雷达视线方向的形变，d_Los(t₀)＝0；S_sol(t_A)和S_sol(t_B)分别为t_A时刻和t_B时刻对应的由水溶开采导致的沉降量，S_sol(t₀)＝0，θ为雷达入射角。

步骤2.3、构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型：对于每一个差分干涉相位，其与雷达视线向形变之间存在如表达式5)所示的关系：

其中，

是第m个差分干涉图中第i个高相干点在t_A时刻和t_B时刻两景SAR图像上解缠后生成的差分干涉相位，λ为雷达波长，B_m为第m个干涉对垂直基线的长度，R_p为相干目标与雷达卫星所在位置之间的距离，δH_i为第i个高相干点对应的高程改正值，w_i ^m为第m个差分干涉图中第i个高相干点对应的残差相位，残差相位包括噪声相位、大气延迟相位和高通形变分量。

将公式3)和公式4)代入公式5)中，可建立差分干涉相位与水溶动力学参数和时间之间的动态函数模型，如表达式6)所示：

其中，△t＝t_B-t_A；动态函数模型中的矿区地质参数GP＝[C_s，δ]为已知参数，根据矿区地质条件和开采条件确定，本实施例中，C_s为48.8g/100gH₂O，δ为3m。

水溶动力学参数UP＝[D，C_x]和雷达参数δH_i为未知参数，因此需对模型参数X＝[D，C_x，δH_i]进行最优求解。

步骤三：对动态函数模型中的模型参数进行最优求解；具体是根据InSAR干涉相位方程组的特征，利用附有不等式约束的最小二乘算法实现模型高效、稳健的参数估计，详情如下：

步骤3.1：将动态函数模型表达式6)写成V＝f(X)-L的形式，f(X)＝AX是关于模型参数X＝[D，C_x，δH_i]的函数，A为m×i的系数矩阵，L为m×1的观测值，V为m×1的改正数矩阵，此方程中参数C_x在不同的SAR影像干涉对中都不同，因此当方程个数为N时，总未知参数个数为N+2个，为典型的秩亏方程，由于C_x的浓度范围小于或等于溶液饱和浓度C_s，可利用附有不等式约束的动态优化迭代算法模型进行求解，根据最小二乘平差理论，可构建如表达式7)所示的动态优化迭代算法模型：

其中，φ(X)为关于X的平差准则，V为m×1的改正数矩阵，V^T为V的转置矩阵，P为最小二乘原理中的权阵，W₁、W₂和W₃分别为C_x、D、δH_i对应的限值。本实施例中，未饱和溶液的浓度C_x的取值范围为[0～48.8]g/100gH₂O，扩散系数D的取值范围为[0～1.87]×10^-9mm²/s，解缠后的差分干涉相位中高程改正值δH_i的取值范围为[-20～20]m，因此本实施例的动态优化迭代算法模型中，W₁＝48.8，W₂＝1.87，W₃＝20。

步骤3.2、将表达式7)的不等式约束条件转换成原待求解动态函数模型(即表达式6))的虚拟观测方程，克服原方程的秩亏问题。

步骤3.3、补充了虚拟观测方程后，可利用线性最小二乘平差法估计未知参数，但由于方程仍然呈现严重病态，可能导致程序运行时不收敛，且高度依赖未知参数初始值的设定。对InSAR观测相位矩阵B中对角线上的各元素增加微小扰动，做B^TB+μI的处理，这样的处理可克服原法方程系数阵的病态问题，并能降低参数对初值的依赖性，可加速算法的收敛，此时X的正则化解如表达式8)所示：

X＝[I-[B^TB+μI]^-1μI]^-1[B^TB+μI]^-1B^TL8)；

其中，I为单位矩阵，B^T为B的转置矩阵；μ为正则参数，可利用矩阵B^TB的最大特征值和最小特征值进行计算。

步骤四：估计InSAR影像干涉时段内的岩盐矿区地表下沉的时序形变；

步骤4.1、将步骤三中计算得出的模型参数X＝[D，C_x，δH_i]代入动态函数模型的表达式6)中，并结合表达式4)和表达式5)计算高相干点在坐标(x，y)上的雷达视线向形变d_Los，计算所得的雷达视线向形变即为基于高相干点的低通形变分量。

步骤4.2、根据大气延迟相位和噪声相位在时空域上的特性，对动态函数模型中的w_i ^m进行空间维低通滤波和时间维高通滤波，从而提取残余相位高通形变分量d_HP。

步骤4.3、将雷达视线向形变d_Los(低通形变分量)和高通形变分量d_HP累加，实现岩盐矿区中任意一点(即第i个高相干点)的地表下沉量S_i的时序形变反演；对S_i进行地理编码，生成岩盐矿区地表垂直向时序形变，从而获取地表变形规律。

步骤五：通过岩盐矿区地表下沉的时序形变获取总采卤量的估计值。

步骤5.1、当忽略水平移动时，t时刻的总采卤量Q通过表达式9)进行计算：

其中，k为一景影像中高相干点的总数，S _i(t)为第i个高相干点在t时刻对应的地表下 沉量，F ₀为由雷达获取的高相干点的像素在实际地理坐标中大小。通过表达式9)即可实 现水溶岩盐矿区内采卤量的估计。

卤水中盐类重要组分含量T _y通过表达式10)进行定量反演：

T _y＝Q×g 10)；

其中，g为卤水中主要组分平均含量。本实施例中，卤水中主要组分平均含量g为0.52t/m³。

步骤5.2、计算溶蚀通道直径d，进而获取地下溶蚀通道直径d随时间的演化过程，合理探究地下溶腔的演化与采卤量重要组分含量变化之间的规律，如图6所示，可以指导矿产资源的开发利用与回采率的提升，具体是：

溶腔形状到开采后期一般近似取作两端为半圆柱、中部为长方柱的长槽状，当忽略水平移动时，溶腔发展过程中任意一井组连接方式的溶采总面积F通过表达式11)进行计算：

其中，F表示溶腔发展过程中任意一井组连接方式的溶采总面积，d为溶蚀通道直径，L_S为井组间距，可通过研究区域矿业公司提供的设计资料中获取。本实施例中，井组间距L_S为56.4m。

根据水溶动力学原理，卤水中盐类重要组分含量T_y还可通过表达式12)表示：

其中，V为溶腔体积，且V＝FH，H为开采矿层厚度，G为矿石体重，E为矿石平均品位本实施例中，开采矿层厚度H为3m，矿石体重G为2.165t/m³，矿石平均品位E为80％。

则通过表达式9)至表达式12)可获得溶蚀通道直径d的计算表达式13)：

设置模拟实验与本实施例进行对照，模拟实验中采用的卫星参数均依据C波段的Sentinel-1A卫星生成的宽幅干涉模式和升轨的SAR影像参数设定。其中D和C_x均利用二维高斯函数模型进行模拟，高程改正值δH_i则利用二维高斯随机函数模拟器模拟，取值范围为[-20～20]m，根据采集到的真实卫星数据干涉组合后，提取出其中48个干涉效果好的干涉对，并将其时空基线参数作为本实施例中的时空基线参数。

模拟实验在模拟差分干涉相位生成过程中，分别加入(0～0.65)rad的噪声，图2为模拟实验中根据模拟真实参数值计算得到的时序形变结果图，分别示意了第24天、第120天、第288天、第432天、第528天、第612天、第708天及第792天的与地表沉降量相关形变结果。

图3为本实施例中通过动态函数模型进行时序形变反演得到的地表下沉量S_i在第24天、第120天、第288天、第432天、第528天、第612天、第708天及第792天的形变结果，由图2和图3可知：即使在较大噪声的干扰下，图3中通过模型参数最优估计值计算得到的时序形变结果仍能与图2中模拟获取的矿区时序形变结果在时空分布上保持良好的一致性，为采卤量估计和卤水中盐类重要组分含量的定量反演提供了准确的分析基础。

图4为500个高相干点通过模拟实验获取的模拟值与通过本申请方法反演的预测值(地表下沉量S_i)之间的对比图；从第24天、第120天、第288天、第432天、第528天、第612天、第708天及第792天的数据统计来看，预测值中差值分布在[-4,4]mm之间的高相干点分别占100％、99.2％、98.4％、96％、94％、88％、86％、83％，其中的最大差值17.2mm占最大形变量223mm的7.7％，满足精度要求。

图5为高相干点上通过本申请方法获得的预测值相对于通过模拟实验获得的时序形变量模拟值的偏差的占比图，从图中可以看出，52.3％的高相干点的形变量偏差分布于[0，2]mm之间，2.4％的高相干点形变量偏差分布于[9，10]mm之间。由本申请中的模型参数求解得到的预测形变值(即预测值)与模拟值的均方根误差(RMSE)为±6.2mm，从而验证了本申请的动态函数模型与附有不等式约束的最小二乘算法实现模型参数最优求解的可靠性。

本申请利用时间序列InSAR技术，将水溶开采动力学原理引入时序InSAR形变建模，建立了一种适用于描述盐类矿床水溶开采沉陷动力学原理的时序InSAR物理形变模型(即动态函数模型)；利用附有不等式约束的平差方法解算模型参数(包括水溶动力学参数和雷达参数)，进而反演岩盐矿区沉陷盆地垂直向的时序形变；建立矿区总采卤量与矿区形变之间时序函数，实现矿区采卤量估计，可实现卤水中盐类重要组分的定量预计，拓宽了InSAR技术应用领域，有利于合理探究地表变形规律，指导矿产资源的开采利用，进一步提升矿产资源的回采率，为盐类矿床的合理开发及利用提供数据参考。

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

Claims (1)

1.一种钻井水溶岩盐矿区采卤量估计方法，其特征在于，包括以下步骤：

步骤一：基于InSAR生成相干系数图和差分干涉图，生成基于高相干点的时序差分干涉相位矩阵；

所述步骤一具体包括以下分步骤：

1.1、将多景SAR卫星数据批量格式转换，生成单视复数影像，并进行初始影像的预滤波，减弱影像的噪声相位；

1.2、对所有影像进行基线估计，然后选取一景影像作为所有影像的超级主影像，设置时空阈值，生成所有符合设定的干涉对；利用生成的干涉对生成干涉图，在此过程中附带生成相干系数图和强度图；

1.3、对干涉图进行以下操作：借助精密轨道数据进行去除平地相位、轨道误差；利用DEM数据去除地形相位；然后进行残余相位滤波；生成差分干涉图，对差分干涉图进行相位解缠，生成解缠后的时序差分干涉相位；

1.4、利用强度、振幅离差和相干系数法提取高相干点，并在时序差分干涉相位中提取高相干点的对应值，生成基于高相干点的时序差分干涉相位矩阵；

步骤二：构建时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型；

所述步骤二中：时序差分干涉相位与水溶动力学参数和时间之间的动态函数模型如表达式6)所示：

其中，

是第m个差分干涉图中第i个高相干点在t_A时刻和t_B时刻两景SAR图像上解缠后生成的差分干涉相位，△t＝t_B-t_A，λ为雷达波长，D为扩散系数，C_s为溶液饱和浓度，C_x为溶液未达到饱和状态时的浓度，δ为饱和溶液层厚度，θ为雷达入射角，B_m为第m个干涉对垂直基线的长度，R_p为相干目标与雷达卫星所在位置之间的距离，δH_i为第i个高相干点对应的高程改正值，w_i ^m为第m个差分干涉图中第i个高相干点对应的残差相位；

所述步骤二中，t时刻对应的由水溶开采导致的沉降量S_sol(t)通过表达式3)计算：

S_sol(t)＝(D(C_s-C_x)/C_sδ)·t3)；

步骤三：对动态函数模型中的模型参数进行最优求解；

所述步骤三中，根据最小二乘平差理论，利用附有不等式约束条件的动态优化迭代算法模型对模型参数X＝[D，C_x，δH_i]进行最优求解；构建的动态优化迭代算法模型如表达式7)所示：

其中，φ(X)为关于X的平差准则，V为m×1的改正数矩阵，V^T为V的转置矩阵，P为最小二乘原理中的权阵，W₁、W₂和W₃分别为C_x、D、δH_i对应的限值；

所述步骤三中，将动态优化迭代算法模型转换为步骤二中动态函数模型的虚拟观测方程，对InSAR观测相位矩阵B增加扰动，此时X的正则化解如表达式8)所示：

X＝[I-[B^TB+μI]^-1μI]^-1[B^TB+μI]^-1B^TL8)；

其中，I为单位矩阵，B^T为B的转置矩阵，μ为正则参数，L为m×1的观测值；

步骤四：估计InSAR影像干涉时段内的岩盐矿区地表下沉量的时序形变；

所述步骤四中，将步骤三中求解得到的模型参数代入动态函数模型中，并通过表达式4)计算雷达视线向形变d_Los：

d_Los＝d_Los(t_B)-d_Los(t_A)＝[S_sol(t_B)-S_sol(t_A)]·cosθ4)；

其中，d_Los(t_B)和d_Los(t_A)分别是t_A时刻和t_B时刻相对于开始时间t₀＝0沿雷达视线方向的形变；S_sol(t_A)和S_sol(t_B)分别为t_A时刻和t_B时刻对应的由水溶开采导致的沉降量；

所述步骤四中，根据大气延迟相位和噪声相位在时空域上的特性，对动态函数模型中的w_l ^m进行空间维低通滤波和时间维高通滤波，从而提取残余相位高通形变分量d_HP，将雷达视线向形变d_Los和高通形变分量d_HP累加，实现岩盐矿区中任意一点的地表下沉量S_i的时序形变反演；对S_i进行地理编码，生成岩盐矿区地表垂直向时序形变，从而获取地表变形规律；

步骤五：通过岩盐矿区地表下沉量的时序形变获取总采卤量的估计值；

所述步骤五中，t时刻的总采卤量Q通过表达式9)进行计算：

其中，k为一景影像中高相干点的总数，S_i(t)为t时刻对应的地表下沉量，F₀为由雷达获取的高相干点的像素在实际地理坐标中大小；

卤水中盐类重要组分含量T_y通过表达式10)计算：

T_y＝Q×g10)；

其中，g为卤水中主要组分平均含量；

所述步骤五中，溶蚀通道直径d通过表达式13)计算：

其中，L_S为井组间距，H为开采矿层厚度，G为矿石体重，E为矿石平均品位；

溶腔发展过程中任意一井组连接方式的溶采总面积F通过表达式11)进行计算：

Images