CN111650587B - 一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 - Google Patents
一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 Download PDFInfo
- Publication number
- CN111650587B CN111650587B CN202010589927.9A CN202010589927A CN111650587B CN 111650587 B CN111650587 B CN 111650587B CN 202010589927 A CN202010589927 A CN 202010589927A CN 111650587 B CN111650587 B CN 111650587B
- Authority
- CN
- China
- Prior art keywords
- mining area
- ground surface
- dynamic
- deformation
- settlement
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000005065 mining Methods 0.000 title claims abstract description 246
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000012544 monitoring process Methods 0.000 claims abstract description 57
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 24
- 238000006073 displacement reaction Methods 0.000 claims abstract description 24
- 238000004364 calculation method Methods 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 claims description 56
- 238000004062 sedimentation Methods 0.000 claims description 23
- 238000004590 computer program Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 12
- 230000002068 genetic effect Effects 0.000 claims description 11
- XOFYZVNMUHMLCC-ZPOLXVRWSA-N prednisone Chemical compound O=C1C=C[C@]2(C)[C@H]3C(=O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 XOFYZVNMUHMLCC-ZPOLXVRWSA-N 0.000 claims description 11
- 230000001131 transforming effect Effects 0.000 claims description 6
- 238000010276 construction Methods 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 238000013507 mapping Methods 0.000 claims description 4
- 230000000007 visual effect Effects 0.000 abstract description 6
- 230000000694 effects Effects 0.000 abstract description 5
- 238000010586 diagram Methods 0.000 description 10
- 238000005516 engineering process Methods 0.000 description 9
- 238000012545 processing Methods 0.000 description 9
- 238000003384 imaging method Methods 0.000 description 5
- 230000007246 mechanism Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 239000011435 rock Substances 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000012905 input function Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000005305 interferometry Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000004075 alteration Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000000126 substance Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/885—Radar or analogous systems specially adapted for specific applications for ground probing
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质，获得待监测矿区地表沿雷达视线向的多时相形变监测值；将多时相一维雷达视线向形变转换为垂直沉降；利用一种新型矿区地表动态沉降模型和与其对应的参数估计智能算法对多时相垂直沉降进行拟合解算，得到矿区地表各点对应的动态沉降模型参数值，既而得到矿区地表动态沉降；基于矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和解算得到的矿区地表动态沉降，对矿区地表动态二维水平位移建模并求解。所提出的矿区地表动态三维形变模型，不仅能够对矿区地下单工作面开采导致的地表形变进行准确预计，而且对于多工作面重复采动导致的地表形变同样具有良好的拟合效果。
Description
技术领域
本发明属于InSAR领域，特别涉及一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质。
背景技术
对矿区地表动态三维形变进行监测和预计，对于掌握矿区形变机理和进行地质灾害防治与评估具有重要意义。合成孔径雷达干涉测量(Interferometric SyntheticAperture Radar，InSAR)技术以其无接触、高精度、高空间分辨率以及不受云雨天气影响等独特优势在矿区地表形变监测中发挥着越来越重要的作用。但是，由于目前星载SAR传感器数量、监测目标优先级等的限制，对于特定矿区而言，在一定时间段内覆盖其地表的SAR数据往往是有限的。这导致InSAR技术得到的矿区地表形变监测数据仅仅是某些特定时刻的稀疏结果，其时间分辨率往往受到较大限制，从而不足以对矿区地表真实动态形变进行有效反应。而且SAR传感器侧视成像的几何机制，也使得其得到的监测数据仅是地表真实三维形变在雷达视线向的一维投影。这同样大大限制了其在矿区地表形变监测中的作用和监测精度。因此，发展一种基于InSAR的矿区地表三维动态形变监测技术迫在眉睫。然而，目前为止，该方面的研究很少有人涉足。
鉴于矿区地表水平移动与对应方向的垂直沉降梯度之间存在线性比例关系，之前提出“一种利用单个InSAR干涉对获取矿区地表三维形变场的方法”(专利号：CN201210440875)。但该方法只能对组成单个干涉对的两景SAR影像期间的地表形变进行反演，而无法用于矿区地表动态三维形变反演。此外，矿区地表形变遵循线性叠加原理。该理论已经广泛应用于概率积分法(一种矿区地表形变拟合预计通用方法)等矿区形变预计理论中。对于单/多工作面开采，其本质上都是地下物质损失后使得周边岩层原始应力平衡状态破坏，从而导致开采工作面上覆岩层及其地表发生移动变形。上述理论在单工作面开采导致的地表单时段或时序形变中虽已有应用，但是在矿区动态三维形变研究以及多工作面重复采动形变研究中的应用仍然很少。
发明内容
本发明的目的是提供一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质，该方法基于矿区开采岩层移动变形机理，利用监测得到的沿SAR传感器雷达视线向的一维多时相形变监测数据，建立其与矿区地表多时相垂直沉降之间的投影转换关系函数模型，并利用一种新型矿区地表动态沉降模型以及与其配套的参数估计智能算法对矿区地表各点进行动态沉降拟合预计，得到矿区地表动态沉降后，利用矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建矿区地表二维水平移动估计方程。
本发明不仅能够对矿区地下单工作面开采导致的地表变形进行准确预计，而且对于多工作面重复采动导致的地表变形同样具有良好的拟合效果。基于所提出的参数估计智能算法，本发明可以对复杂地质采矿条件下的矿区地表动态三维形变进行准确估计；突破已有方法在多工作面重复采动沉陷预计中的不适用性以及无法解算动态三维形变的限制；拓展了InSAR技术在矿区形变监测中应用前景的同时，为矿区形变监测与解译、地质灾害预计、评估和治理等提供了一种有效的手段。
本发明的提供的技术方案如下：
一方面，一种顾及移动规律的矿区地表三维动态形变估计方法，包括如下步骤：
S1：获取待监测矿区的多时相InSAR干涉对，并基于多时相InSAR干涉对得到待监测矿区地表沿雷达视线(Line-of-Sight，LOS)向的多时相DInSAR形变监测结果；
S2：基于DInSAR中所用的DEM坐标系统为参考，将S1得到的所有多时相DInSAR形变监测结果进行地理编码，得到空间坐标系统一致的多时相LOS向形变dLOS；
S3：基于单个DInSAR干涉对LOS向形变dlos与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直形变UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
S4：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降模型；
S5：将S3得到的矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，基于S4构建的地表动态沉降模型，利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
S6：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和S5计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，得到矿区地表动态二维水平位移。
进一步地，所述S4中的矿区多工作面重复采动情况下的地表动态沉降模型为：
其中，U(i,j,t)表示矿区地表待监测点(i,j)位置处在t时刻的地表动态沉降；fl(i,j,t)表示在矿区地表待监测点(i,j)位置处采用的第l个已有的单工作面开采地表动态沉降拟合函数，为Knothe、Weibull、Logistic和Richards函数中的一种；K为对待监测点进行模型拟合所需函数个数，K的最大值为Kmax，Kmax≤5。
进一步地，所述S5中的利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合具体步骤如下：
S51：根据矿区开采具体情况，选取Knothe、Weibull、Logistic和Richards中的一种或几种函数，作为矿区地表待监测点处动态沉降模型的拟合函数库；
比如，选取Logistic函数为模型拟合输入函数，即
S52：令K的初值为1，即K＝1，并设置拟合误差阈值Δ；
此处设置(i,j)点处拟合得到的沉降最大值的0.1倍作为阈值Δ；
S53：根据K值和矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM，利用遗传算法对地表动态沉降模型的待估参数值进行求解，得到矿区地表待监测点处在K下对应的动态沉降输出值；
在利用遗传算法对动态沉降模型参数进行求解的过程中，设置遗传算法的解空间搜索范围、种群大小、迭代次数、停滞代数、交叉变异因子参数，并设置遗传算法的目标函数为：
其中，Upara为地表动态沉降模型的待估参数值；之后，基于上述设置的遗传算法参数在解空间中对模型参数进行搜索，获得待估模型参数值的全局最优解；
S54：计算地表动态沉降模型输出值和观测值之间的拟合误差δ并判断拟合误差是否满足阈值，即δ≤Δ是否成立；
S55：若δ≤Δ不成立，则令K＝K+1并重复S53～S55，直到δ≤Δ成立或K达到阈值Kmax，进入S56；
S56：将当前得到的待监测点的地表动态沉降模型的参数值输出，得到该点动态沉降模型；
利用上述方法对矿区地表所有监测点依次进行模型拟合，即可得到各监测点的最优拟合模型参数值；之后，输入起始采动后的任意时刻的时间t，即可得到该时刻对应的地表动态沉降U(t)。
进一步地，所述S6中的矿区地表沿东西和南北方向的二维水平位移求解模型表达式为：
其中，U(i,j,t)表示矿区地表待监测点(i,j)位置处在t时刻的动态沉降，E(i,j,t)和N(i,j,t)分别表示矿区地表待监测点(i,j)位置处在t时刻的东西向和南北向的位移，i＝1,2,…,Rmax-1，j＝1,2,…,Cmax-1，Rmax和Cmax表示影像的最大行和列数，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率。
进一步地，所述单个DInSAR干涉对LOS向形变dlos与地表真实三维形变U、E、N之间的投影转换关系为：
dlos＝U·cosθ-[E·sin(α-3π/2)+N·cos(α-3π/2)]·sinθ
其中，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角。
进一步地，所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
其中，U(i,j)、E(i,j)和N(i,j)分别表示矿区地表待监测点(i,j)位置处在垂直、东西和南北向的形变分量，Rmax和Cmax表示影像的最大行和列数，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率。
进一步地，所述多时相LOS向形变dLOS与矿区地表多时相垂直沉降UM之间的几何映射模型为：
其中，UM(i,j)表示矿区地表待监测点(i,j)位置处的多时相垂直沉降值，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角。
另一方面，一种顾及移动规律的矿区地表三维动态形变估计装置，其特征在于：包括：
数据获取单元：获取待监测矿区的多时相InSAR干涉对，并基于多时相InSAR干涉对得到待监测矿区地表沿雷达视线向的多时相DInSAR形变监测结果；
坐标转换单元：基于DInSAR中所用的DEM坐标系统为参考，将得到的所有多时相DInSAR形变监测结果进行地理编码，得到空间坐标系统一致的多时相LOS向形变dLOS；
矿区地表多时相垂直沉降转换单元：基于单个DInSAR干涉对LOS向形变dlos与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直形变UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
地表动态沉降模型构建单元：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降模型；
矿区地表动态沉降计算单元：将得到的矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，基于地表动态沉降模型构建单元构建的地表动态沉降模型，利用参数估计智能算法单元对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
矿区地表动态二维水平位移计算单元：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和矿区地表动态沉降计算单元计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，得到矿区地表动态二维水平位移。
进一步地，所述地表动态沉降模型构建单元构建的矿区多工作面重复采动情况下的地表动态沉降模型为：
其中，U(i,j,t)表示矿区地表待监测点(i,j)位置处在t时刻的地表动态沉降；fl(i,j,t)表示在矿区地表待监测点(i,j)位置处采用的第l个已有的单工作面开采地表动态沉降拟合函数，为Knothe、Weibull、Logistic和Richards函数中的一种；K为对待监测点进行模型拟合所需函数个数，K的最大值为Kmax，Kmax≤5。
再一方面，一种计算机存储介质，包括计算机程序，所述计算机程序被处理器执行时实现所述的一种顾及移动规律的矿区地表三维动态形变估计方法。
有益效果
本发明提供了一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质，获得覆盖待监测矿区的所有可用InSAR数据，利用差分InSAR技术得到待监测矿区地表沿雷达视线向的多时相形变监测值；根据SAR成像投影几何以及矿区地表水平移动与对应方向的垂直沉降梯度之间的线性比例关系函数模型，将多时相一维雷达视线向形变转换为垂直沉降；利用一种新型矿区地表动态沉降模型和与其对应的参数估计智能算法对多时相垂直沉降进行拟合求解，得到矿区地表各点的动态沉降；基于矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和解算得到的矿区地表动态沉降，构建矿区地表动态二维水平位移模型并解算得到矿区地表动态二维水平位移。该方法计算简单，实现方便；本发明的技术效果主要体现在以下几点：
第一、提出一种顾及移动规律的新型矿区地表动态沉降模型，该模型在单/多工作面开采地表沉降监测中均具有较好的拟合效果，突破了已有矿区地表动态沉降模型在多工作面重复采动情况下对地表沉降拟合适用性差的局限；
第二、发展了一种参数估计智能算法，通过设定拟合残差阈值，实现矿区地表沉降模型参数的自适应选取和拟合，增强了本发明在实际中的应用程度和推广性；
第三、本发明结合InSAR技术和矿区地表岩层移动变形机理，大大降低了由InSAR一维雷达视线向形变解算矿区地表动态三维形变的解算条件要求，比如对影像数量、成像时间和成像几何等的苛刻要求；
第四、通过对覆盖待监测矿区的所有InSAR数据进行处理，并用于矿区地表动态三维形变解算，在大大提高了SAR数据的利用率的同时，增强了解算得到的矿区地表动态三维形变的可靠性和精度。
附图说明
图1本发明实施例所述方法的流程示意图；
图2模拟的矿区地表动态三维形变，其中(a)-(c)分别表示模拟的沿垂直、东西和南北向的动态三维形变；
图3模拟的雷达视线向形变；
图4基于本发明实施例所述方法拟合的矿区地表动态三维形变，其中(a)-(c)分别表示拟合得到的沿垂直、东西和南北向的动态三维形变；
图5三维形变模拟值与拟合值均方根误差，其中(a)-(c)分别表示垂直、东西和南北向的均方根误差。
具体实施方式
下面将结合实施例对本发明做进一步的说明。
如图1所示，本发明实施例提供了一种顾及移动规律的矿区地表三维动态形变估计方法，包括如下步骤：
S1：将覆盖待监测矿区的所有可用合成孔径雷达干涉测量(InterferometricSynthetic Aperture Radar，InSAR)数据按照数据获取平台类型的差异组成不同数据集，并对各数据集按照设定的时空基线阈值组成多时相InSAR干涉对；
利用合成孔径雷达差分干涉测量(Differential InSAR，DInSAR)技术对步骤S1中的各多时相InSAR干涉对依次进行处理，得到矿区地表沿雷达视线(line-of-sight，LOS)向的多时相DInSAR形变监测结果；
所述DInSAR技术为进行SAR数据差分干涉处理的常用已有技术；
S2：以进行DInSAR技术数据处理时所用的外部数字高程模型(Digital ElevationModel，DEM)坐标系统为参考，将步骤S2中得到的所有多时相DInSAR形变监测结果进行地理编码，即对各多时相DInSAR形变监测结果进行坐标系统框架统一，得到空间坐标系统一致的多时相DInSAR LOS向形变监测结果dLOS；
S3：基于单个DInSAR干涉对得到的LOS向形变与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直沉降UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
所述单个DInSAR干涉对LOS向形变dlos与地表真实三维形变U、E、N之间的投影转换关系为：
dlos＝U·cosθ-[E·sin(α-3π/2)+N·cos(α-3π/2)]·sinθ
其中，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角；
所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
其中，U(i,j)、E(i,j)和N(i,j)分别表示矿区地表待监测点(i,j)位置处在垂直、东西和南北向的形变分量，Rmax和Cmax表示影像的最大行和列数，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率；
所述多时相LOS向形变dLOS与矿区地表多时相垂直沉降UM之间的几何映射模型为：
其中，UM(i,j)表示矿区地表待监测点(i,j)位置处的多时相垂直沉降，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角；
S4：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降预计模型；
所述矿区多工作面重复采动情况下的地表动态沉降模型为：
其中，U(i,j,t)表示矿区地表待监测点(i,j)位置处在t时刻的地表动态沉降；fl(i,j,t)表示在(i,j)位置处采用的第l个已有的单工作面开采地表动态沉降拟合函数，为Knothe、Weibull、Logistic和Richards函数中的一种；K为对待监测点进行模型拟合所需函数个数，K的最大值为Kmax，Kmax≤5；
S5：将矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，根据S5构建的地表动态沉降模型，利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合具体步骤如下：
S51：根据矿区开采具体情况，选取Knothe、Weibull、Logistic和Richards中的一种或几种函数，作为矿区地表动态沉降模型的拟合函数库；
比如，选取Logistic函数为模型拟合输入函数，即
S52：令K的初值为1，即K＝1，并设置拟合误差阈值Δ；
此处建议设置(i,j)点处拟合得到的沉降最大值的0.1倍作为阈值Δ；
S53：根据K值和矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM，利用遗传算法对地表动态沉降模型的待估参数值进行求解，得到(i,j)点处在K下对应的动态沉降输出值；
在利用遗传算法对动态沉降模型参数进行求解的过程中，设置遗传算法的解空间搜索范围(
al(i,j)∈[-50000,50000]，bl(i,j)∈[-0.1,0.1])、种群大小(500)、迭代次数(2000)、停滞代数(5000)，交叉变异因子(0.65)参数，并设置遗传算法的目标函数为：
其中，Upara为地表动态沉降模型的待估参数值；之后，基于上述设置的遗传算法参数在解空间中对模型参数进行搜索，获得待估模型参数值的全局最优解；
S54：计算地表动态沉降模型输出值和观测值之间的拟合误差δ并判断拟合误差是否满足阈值，即δ≤Δ是否成立；
S55：若δ≤Δ不成立，则令K＝K+1并重复S53～S55，直到δ≤Δ成立或K达到阈值Kmax，进入S56；
S56：将当前得到的待监测点的地表动态沉降模型的参数值输出，得到该点动态沉降模型；
利用上述方法对矿区地表所有监测点依次进行模型拟合，即可得到各监测点的最优拟合模型参数值；之后，输入起始采动后的任意时刻的时间t，即可得到该时刻对应的地表动态沉降U(t)。
S7：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和S6计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，并进行求解，得到矿区地表动态二维水平位移；
所述矿区地表沿东西和南北方向的动态二维水平位移求解模型表达式为：
其中，
表示矿区地表待监测点(i,j)位置处在t时刻的动态沉降，E(i,j,t)和N(i,j,t)分别表示矿区地表待监测点(i,j)位置处在t时刻沿东西和南北向的动态二维水平位移，i＝1,2,…,Rmax-1，j＝1,2,…,Cmax-1，Rmax和Cmax表示影像的最大行和列数，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率。
为了更加清楚的说明本发明及其优势，下述将结合具体实施例及其相关的部分图来进一步解释本发明提供的所述方法。以某矿区邻近工作面开采(如图2(a)垂直向形变图中白色矩形WP1和WP2所示)导致的地表三维形变为例，利用概率积分法(一种矿区地表形变拟合预计通用模型)模拟出其垂直、东西和南北方向的动态三维形变(如图2所示)，并以ALOS PALSAR卫星升轨数据的成像几何参数(θ＝38°,α＝-10°)为例，根据雷达视线向形变与地表真实三维形变之间的投影转换关系，将地表真实三维形变投影转换到雷达视线向，得到模拟的InSAR雷达视线向形变监测数据(如图3所示)。利用本发明所述方法对模拟的监测数据进行处理，估计得到矿区地表全盆地动态三维形变模型及其参数值；由估计得到的模型参数计算得到与模拟的SAR影像获取时间一致的地表三维形变结果(如图4所示)。为了评定本发明的拟合精度，分别计算出垂直、东西和南北向形变模拟值与模型拟合值之间的均方根误差，如图5所示。
由图2至图5可以看出，本发明所述方法不仅可以对两工作面外侧的单工作面开采影响区域的地表动态三维形变进行准确估计，而且对两工作面之间的多工作面重复采动下的地表形变同样具有良好的拟合结果。与模拟真值相比，解算得到的沿垂直、东西和南北方向的动态三维形变估计结果的均方根误差普遍较小，沿垂直、东西和南北方向的均方根误差均值分别为0.4cm、2.5cm和2cm，达到了矿区开采沉陷监测预计要求。从而验证了本发明所述方法的可行性和可靠性。
基于上述方法，本发明实施例还提供一种顾及移动规律的矿区地表三维动态形变估计装置，包括：
数据获取单元：获取待监测矿区的多时相InSAR干涉对，并基于多时相InSAR干涉对得到待监测矿区地表沿雷达视线(Line-of-Sight，LOS)向的多时相DInSAR形变监测结果；
坐标转换单元：基于DInSAR中所用的DEM坐标系统为参考，将S1得到的所有多时相DInSAR形变监测结果进行地理编码，得到空间坐标系统一致的多时相LOS向形变dLOS；
矿区地表多时相垂直沉降转换单元：基于单个DInSAR干涉对LOS向形变dlos与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直形变UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
地表动态沉降模型构建单元：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降模型；
矿区地表动态沉降计算单元：将得到的矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，基于地表动态沉降模型构建单元构建的地表动态沉降模型，利用参数估计智能算法单元对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
矿区地表动态二维水平位移计算单元：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和矿区地表动态沉降计算单元计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，得到矿区地表动态二维水平位移。
所述地表动态沉降模型构建单元和参数估计智能算法单元的具体内容参见上述方法中的内容，在此不做重复描述。
应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。
本发明实施例还提供一种计算机存储介质，包括计算机程序，所述计算机程序被处理器执行时实现所述的一种顾及移动规律的矿区地表三维动态形变估计方法；其有益效果参见方法部分的有益效果，在此不再赘述。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详尽的说明，所属领域的普通技术人员应当理解，上述实施例仅仅是对本发明的示意性实现方式的解释，实施例中的细节并不构成对本发明范围的限制，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案的等效变换、简单替换等显而易见的改变，均落在本发明保护范围之内。
Claims (7)
1.一种顾及移动规律的矿区地表三维动态形变估计方法，其特征在于：包括如下步骤：
S1：获取待监测矿区的多时相InSAR干涉对，并基于多时相InSAR干涉对得到待监测矿区地表沿雷达视线向的多时相DInSAR形变监测结果；
S2：基于DInSAR中所用的DEM坐标系统为参考，将S1得到的所有多时相DInSAR形变监测结果进行地理编码，得到空间坐标系统一致的多时相LOS向形变dLOS；
S3：基于单个DInSAR干涉对LOS向形变dlos与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直形变UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
S4：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降模型；
S5：将S3得到的矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，基于S4构建的地表动态沉降模型，利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
S6：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和S5计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，得到矿区地表动态二维水平位移；
所述单个DInSAR干涉对LOS向形变dlos与地表真实三维形变U、E、N之间的投影转换关系为：
dlos＝U·cosθ-[E·sin(α-3π/2)+N·cos(α-3π/2)]·sinθ
其中，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角；
所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
i＝1，2，…，Rmax-1；j＝1，2，…，Cmax-1
其中，U(i，j)、E(i，j)和N(i，j)分别表示矿区地表待监测点(i，j)位置处在垂直、东西和南北向的形变分量，Rmax和Cmax表示影像的最大行和列数，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率；
所述多时相LOS向形变dLOS与矿区地表多时相垂直沉降UM之间的几何映射模型为：
其中，UM(i，j)表示矿区地表待监测点(i，j)位置处的多时相垂直沉降值，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角。
2.根据权利要求1所述方法，其特征在于：所述S4中的矿区多工作面重复采动情况下的地表动态沉降模型为：
其中，U(i，j，t)表示矿区地表待监测点(i，j)位置处在t时刻的地表动态沉降；fl(i，j，t)表示在矿区地表待监测点(i，j)位置处采用的第l个已有的单工作面开采地表动态沉降拟合函数，为Knothe、Weibull、Logistic和Richards函数中的一种；K为对待监测点进行模型拟合所需函数个数，K的最大值为Kmax，Kmax≤5。
3.根据权利要求1所述方法，其特征在于：所述S5中的利用参数估计智能算法对地表各待监测点进行地表动态沉降模型参数拟合具体步骤如下：
S51：根据矿区开采具体情况，选取Knothe、Weibull、Logistic和Richards中的一种或几种函数，作为矿区地表待监测点处动态沉降模型的拟合函数库；
S52：令K的初值为1，即K＝1，并设置拟合误差阈值Δ；
S53：根据K值和矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM，利用遗传算法对地表动态沉降模型的待估参数值进行求解，得到矿区地表待监测点处在K下对应的动态沉降输出值；
S54：计算地表动态沉降模型输出值和观测值之间的拟合误差δ并判断拟合误差是否满足阈值，即δ≤Δ是否成立；
S55：若δ≤Δ不成立，则令K＝K+1并重复S53～S55，直到δ≤Δ成立或K达到阈值Kmax，进入S56；
S56：将当前得到的待监测点的地表动态沉降模型的参数值输出，得到该点动态沉降模型。
4.根据权利要求1所述方法，其特征在于：所述S6中的矿区地表沿东西和南北方向的二维水平位移求解模型表达式为：
其中，U(i，j，t)表示矿区地表待监测点(i，j)位置处在t时刻的动态沉降，E(i，j，t)和N(i，j，t)分别表示矿区地表待监测点(i，j)位置处在t时刻的东西向和南北向的位移，i＝1，2，…，Rmax-1，j＝1，2，…，Cmax-1，Rmax和Cmax表示影像的最大行和列数，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率。
5.一种顾及移动规律的矿区地表三维动态形变估计装置，其特征在于：包括：
数据获取单元：获取待监测矿区的多时相InSAR干涉对，并基于多时相InSAR干涉对得到待监测矿区地表沿雷达视线向的多时相DInSAR形变监测结果；
坐标转换单元：基于DInSAR中所用的DEM坐标系统为参考，将得到的所有多时相DInSAR形变监测结果进行地理编码，得到空间坐标系统一致的多时相LOS向形变dLOS；
矿区地表多时相垂直沉降转换单元：基于单个DInSAR干涉对LOS向形变dlos与地表沿垂直、东西和南北方向的三维形变U、E、N之间的投影转换关系和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型，构建多时相LOS向形变dLOS与多时相垂直形变UM之间的几何映射模型，并将多时相LOS向形变dLOS转换得到矿区地表多时相垂直沉降UM；
地表动态沉降模型构建单元：根据矿区地下开采导致的地表形变的线性叠加原理以及单工作面开采地表动态沉降拟合函数，构建矿区多工作面重复采动情况下的地表动态沉降模型；
矿区地表动态沉降计算单元：将得到的矿区地表多时相垂直沉降UM及其对应的干涉对影像获取时间TM组成拟合数据集，基于地表动态沉降模型构建单元构建的地表动态沉降模型，利用参数估计智能算法单元对地表各待监测点进行地表动态沉降模型参数拟合，得到矿区地表各待监测点对应的动态沉降模型参数值，进而获得矿区地表动态沉降；
矿区地表动态二维水平位移计算单元：根据矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型和矿区地表动态沉降计算单元计算得到的矿区地表动态沉降，构建矿区地表沿东西和南北方向的动态二维水平位移求解模型，得到矿区地表动态二维水平位移；
所述单个DInSAR干涉对LOS向形变dlos与地表真实三维形变U、E、N之间的投影转换关系为：
dlos＝U·cosθ-[E·sin(α-3π/2)+N·cos(α-3π/2)]·sinθ
其中，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角；
所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
i＝1，2，…，Rmax-1；j＝1，2，…，Cmax-1
其中，U(i，j)、E(i，j)和N(i，j)分别表示矿区地表待监测点(i，j)位置处在垂直、东西和南北向的形变分量，Rmax和Cmax表示影像的最大行和列数，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率；
所述多时相LOS向形变dLOS与矿区地表多时相垂直沉降UM之间的几何映射模型为：
其中，UM(i，j)表示矿区地表待监测点(i，j)位置处的多时相垂直沉降值，b和r分别为矿区地表形变的水平移动系数和主要影响半径，ΔE和ΔN表示地理编码后的矿区地表形变监测数据中像元沿东西和南北方向的空间分辨率，θ和α分别表示SAR监测平台的雷达入射角和飞行方位角。
6.根据权利要求5所述的装置，其特征在于，所述地表动态沉降模型构建单元构建的矿区多工作面重复采动情况下的地表动态沉降模型为：
其中，U(i，j，t)表示矿区地表待监测点(i，j)位置处在t时刻的地表动态沉降；fl(i，j，t)表示在矿区地表待监测点(i，j)位置处采用的第l个已有的单工作面开采地表动态沉降拟合函数，为Knothe、Weibull、Logistic和Richards函数中的一种；K为对待监测点进行模型拟合所需函数个数，K的最大值为Kmax，Kmax≤5。
7.一种计算机存储介质，包括计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-4任一项所述的一种顾及移动规律的矿区地表三维动态形变估计方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010589927.9A CN111650587B (zh) | 2020-06-24 | 2020-06-24 | 一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010589927.9A CN111650587B (zh) | 2020-06-24 | 2020-06-24 | 一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111650587A CN111650587A (zh) | 2020-09-11 |
| CN111650587B true CN111650587B (zh) | 2022-03-04 |
Family
ID=72347655
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010589927.9A Active CN111650587B (zh) | 2020-06-24 | 2020-06-24 | 一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111650587B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113253270A (zh) * | 2021-06-11 | 2021-08-13 | 中国测绘科学研究院 | 基于InSAR和Okada模型反演地下采矿参数的方法和系统 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106767380A (zh) * | 2017-01-19 | 2017-05-31 | 中南大学 | 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| ITRM20070399A1 (it) * | 2007-07-19 | 2009-01-20 | Consiglio Nazionale Ricerche | Metodo di elaborazione di dati rilevati mediante radar ad apertura sintetica (synthetic aperture radar - sar) e relativo sistema di telerilevamento. |
| CN102927934B (zh) * | 2012-11-07 | 2015-01-28 | 中南大学 | 一种利用单个InSAR干涉对获取矿区地表三维形变场的方法 |
| KR101490981B1 (ko) * | 2012-12-28 | 2015-02-09 | 서울시립대학교 산학협력단 | 위성레이더 간섭도의 이온왜곡 보정방법 및 그 장치 |
| CN108983229B (zh) * | 2018-05-03 | 2022-04-19 | 电子科技大学 | 基于sar层析技术的高压输电铁塔高度及形变提取方法 |
| CN109738892B (zh) * | 2019-01-24 | 2020-06-30 | 中南大学 | 一种矿区地表高时空分辨率三维形变估计方法 |
| CN110058234A (zh) * | 2019-05-20 | 2019-07-26 | 太原理工大学 | 一种解算矿区地表沉降三维形变的方法 |
-
2020
- 2020-06-24 CN CN202010589927.9A patent/CN111650587B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106767380A (zh) * | 2017-01-19 | 2017-05-31 | 中南大学 | 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111650587A (zh) | 2020-09-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhu et al. | Integration of three dimensional discontinuous deformation analysis (DDA) with binocular photogrammetry for stability analysis of tunnels in blocky rockmass | |
| Dong et al. | Sensitivity analysis of augmented reality-assisted building damage reconnaissance using virtual prototyping | |
| CN111650579B (zh) | 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 | |
| Daout et al. | Strain partitioning and present‐day fault kinematics in NW Tibet from Envisat SAR interferometry | |
| CN102607534B (zh) | 基于运动恢复结构的卫星相对姿态测量方法 | |
| EP3455655B1 (en) | System and method for 3d restoration of complex subsurface models | |
| CN107016649A (zh) | 一种基于局部低秩张量估计的视觉数据补全方法 | |
| Hu et al. | Experiment and application of NATM tunnel deformation monitoring based on 3D laser scanning | |
| Thomas et al. | High resolution (400 m) motion characterization of sea ice using ERS-1 SAR imagery | |
| Wang et al. | PrecipGAN: Merging microwave and infrared data for satellite precipitation estimation using generative adversarial network | |
| CN106767380A (zh) | 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法 | |
| Wang et al. | Crustal density structure, lithosphere flexure mechanism, and isostatic state throughout the Qinling Orogen revealed by in situ dense gravity observations | |
| Beddingfield et al. | Shallow normal fault slopes on Saturnian icy satellites | |
| Yu et al. | Structural state estimation of earthquake-damaged building structures by using UAV photogrammetry and point cloud segmentation | |
| CN111650587B (zh) | 一种顾及移动规律的矿区地表三维动态形变估计方法、装置及存储介质 | |
| Feng et al. | A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica | |
| KR102492075B1 (ko) | 실시간 해양 예측 시스템 (koos-opem) 및 이를 이용한 해양 예측 방법 | |
| Ariza‐López et al. | Spline quasi‐interpolation in the Bernstein basis and its application to digital elevation models | |
| CN116068511B (zh) | 一种基于深度学习的InSAR大尺度系统误差改正方法 | |
| Chen et al. | Improving completeness and accuracy of 3D point clouds by using deep learning for applications of digital twins to civil structures | |
| CN111780660B (zh) | 矿区三维多量级形变优化方法及优化装置 | |
| Watson et al. | An InSAR‐GNSS velocity field for Iran | |
| Barazzetti | Spatio-temporal analysis of georeferenced time-series applied to structural monitoring | |
| Bulyshev et al. | A super-resolution algorithm for enhancement of FLASH LIDAR data | |
| Ghassemi et al. | Salt extrusion kinematics: insights from existing data, morphology and InSAR modelling of the active emergent Anguru diapir in the Zagros fold and thrust belt, Iran |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |