CN111323776A - 一种矿区形变的监测方法 - Google Patents

一种矿区形变的监测方法 Download PDF

Info

Publication number
CN111323776A
CN111323776A CN202010122691.8A CN202010122691A CN111323776A CN 111323776 A CN111323776 A CN 111323776A CN 202010122691 A CN202010122691 A CN 202010122691A CN 111323776 A CN111323776 A CN 111323776A
Authority
CN
China
Prior art keywords
deformation
expression
phase
mining
time
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.)
Granted
Application number
CN202010122691.8A
Other languages
English (en)
Other versions
CN111323776B (zh
Inventor
邢学敏
刘祥彬
朱懿凯
朱凌杰
鲍亮
杨东
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hunan Tianxianghe Information Technology Co ltd
Original Assignee
Changsha University of Science and Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Changsha University of Science and Technology filed Critical Changsha University of Science and Technology
Priority to CN202010122691.8A priority Critical patent/CN111323776B/zh
Publication of CN111323776A publication Critical patent/CN111323776A/zh
Application granted granted Critical
Publication of CN111323776B publication Critical patent/CN111323776B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/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
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric 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

本发明提供一种矿区形变的监测方法，结合地表移动坐标时间函数，在传统概率积分静态模型中融入时间参数，利用雷达视线向形变与地表下沉的函数关系构建一种时序InSAR概率积分动态形变模型；利用遗传算法搜索得出模型中的概率积分预计参数初始值；使用单纯型搜索法进行二次优化搜索，获取模型未知参数的最优估值；将获取的概率积分参数代入时序InSAR概率积分动态形变模型中，实现矿区形变估计及开采沉陷变形预计。本发明克服了纯经验InSAR矿区形变模型中未考虑矿区岩石物理特性的缺陷，可更合理解译形变结果；在传统静态概率积分模型中融入时间因素，实现矿区的动态变形预计，拓宽InSAR技术在矿区开采沉陷预计中的应用。

Description

一种矿区形变的监测方法
技术领域
本发明涉及矿区地表形变监测及变形预计领域，特别地，涉及一种顾及开采沉陷预计参数的时序InSAR的矿区形变的监测方法。
背景技术
传统的矿区地表沉降监测方法，如全站仪/棱镜法、水准测量法、全球导航卫星系统等，虽然能够满足矿区监测精度要求。但这些方法的成本高、时空分辨率较差，对矿区地表整体沉降观测仍存在不足，而且需要监测人员频繁出入现场勘测等特性严重制约了其应用。差分合成孔径雷达干涉测量(Differential Interferometric Synthetic ApertureRadar，简称D-InSAR)技术作为一种新兴的空间对地观测技术，主要用于监测雷达视线方向厘米级或更微小的地球表面形变，以其监测范围广、空间分辨率高、非接触式测量等优点，大大弥补了传统测量手段的不足，为地下开采沉降监测提供了新的契机。然而，时空失相关和大气延迟影响限制了D-InSAR技术的监测精度。为克服D-InSAR技术缺陷，时间序列InSAR技术被学者们提出，该项技术主要是以永久散射体(Permanent Scatterer，简称PS)、小基线集(Small Baseline Subset，简称SBAS)、时域相关点干涉(Temporarily CoherentPoint InSAR，简称TCP-InSAR)等技术为代表，通过对时间序列上多幅SAR影像进行相位分析，利用具有稳定散射特性的高质量点，实现地表形变提取，使其在矿区地表变形监测中更具潜力。
时间序列InSAR数据处理过程中，建立形变模型是至关重要的环节，即建立高相干点形变分量与时间、形变参数之间的函数关系。准确可靠的形变模型可以提高形变估计精度、控制残余相位在小于整周π范围内，而且可以为后续形变结果的解译提供参考。然而，目前InSAR技术采用的多为纯经验数学模型，很少顾及矿区开采沉陷变形机制。矿区开采是一个复杂的非线性过程，其呈现的地表移动随时间的变化呈现复杂的非线性特性，用单一的经验数学函数建模这一物理过程，显然不能真实反应其变形规律，严重影响变形监测的精度及变形预测的准确性，也不利于后期的形变分析解译。概率积分法(ProbabilityIntegral Method，简称PIM)是开采沉陷理论中较为成熟有效的矿区变形预计方法之一，成功解决了地表移动空间问题和上覆岩层移动预计等问题。概率积分预计模型中的预计参数，包括即下沉系数q、影响角正切tanβ、拐点偏距s、水平移动系数b、开采影响传播角θ等，是反应矿区地下沉陷移动变形特征的重要参数。虽然目前有大量学者利用InSAR形变估计的结果反演概率积分变形参数，但InSAR形变估计中仍多采用线性模型建模视线向形变与时间之间的函数过程，这不但会影响InSAR形变估计的精度，更会直接影响后续开采沉陷预计参数的反演。
发明内容
本发明目的在于提供一种顾及概率积分开采沉陷预计参数的时序InSAR的矿区形变监监测方法，以解决现有InSAR形变模型中没有顾及监测矿区物理特性的缺陷，同时也解决了传统利用纯经验InSAR数学模型获取的矿区形变精度不高的缺陷。
本发明采用的技术方案如下：
首先利用时间序列InSAR技术获取监测矿区解缠后的差分干涉相位图，并选取高相干点目标；结合地表移动坐标时间函数，在传统概率积分静态下沉模型中融入时间参数，利用雷达视线向形变与地表下沉的函数关系构建了一种时序InSAR概率积分动态形变模型；其次利用遗传算法搜索得出模型中未知的概率积分预计参数的初始值；再使用单纯型搜索法进行二次优化搜索，进而获取模型未知概率积分预计参数的最优估值；最后利用估计出的模型参数实现矿区时序地表形变的反演及SAR影像时段后续开采沉陷变形预计，详细方案如下：
一种矿区形变的监测方法，包括以下步骤：
步骤一、基于高相干点的时序InSAR差分干涉图生成，具体是：利用时间序列InSAR技术获取监测矿区解缠后的差分干涉相位图；
步骤二、构建InSAR干涉相位与概率积分预计参数的时间函数模型；
步骤三、概率积分预计参数估计，具体是：先利用遗传算法估计概率积分预计参数的初始值，再利用单纯形搜索法对遗传算法得到的概率积分预计参数进行优化；
步骤四、基于高相干点概率积分预计参数的时序形变估计，具体是：估计InSAR 影像干涉时段内的矿区时序形变。
以上技术方案中优选的，步骤一包括如下步骤：
步骤1.1、对N+1幅SAR数据进行干涉组合、影像配准及重采样，生成干涉图和相干图；
步骤1.2、对步骤1.1中干涉图进行去轨道、平地相位和地形相位处理，生成差分干涉图；
步骤1.3、对步骤1.2中的差分干涉图进行相位解缠，得到解缠后的差分干涉图；
步骤1.4、利用步骤1.1中生成的相干图和步骤1.3获得的解缠后的差分干涉图进行高相干点提取，生成基于高相干点目标的时序差分干涉相位矩阵。
以上技术方案中优选的，步骤二包括如下步骤：
步骤2.1、结合地表移动的坐标-时间函数，将时间参数融入概率积分静态模型，构建地表下沉与概率积分参数间的时间函数模型如表达式1)：
Figure RE-GDA0002467328550000031
式中：t为开采时间，W(x，y，t)为工作面推进时间为t时矿山地表任一点(x，y)上开采引发的形变量，w0为该地质采矿条件下的最大沉降值，w0＝mq·cosα，m为开采厚度，q为下沉系数，α为矿层倾角；表达式1)中在t时刻地表下沉W(x，t)和W(y，t)表示如下：
当矿区地下开采工作面倾向充分开采，走向主断面为有限开采时，造成的地表下沉W(x，t)的公式如下：
当x≤nH时，为表达式2)：
Figure RE-GDA0002467328550000032
当nH＜x＜1.3r时，为表达式3)：
Figure RE-GDA0002467328550000033
当1.3r≤x时，为表达式4)：
Figure RE-GDA0002467328550000034
其中：x为地表点的计算横坐标值，
Figure RE-GDA0002467328550000035
为静态下沉公式，erf 为概率积分函数，其形式为
Figure RE-GDA0002467328550000036
u为积分参数，
Figure RE-GDA0002467328550000037
为影响半径，H为开采深度，tanβ为影响角正切值，nH是地表开始移动时工作面的推进距离，n为起动距系数，ω静态超前影响角，v是工作面在匀速开采条件下的推进速度，系数
Figure RE-GDA0002467328550000038
系数
Figure RE-GDA0002467328550000039
当矿区地下开采工作面走向充分开采，倾向有限开采时，造成地表下沉W(y，t)的公式如表达式5)：
Figure RE-GDA0002467328550000041
其中，y为矿区地表点的计算纵坐标，
Figure RE-GDA0002467328550000042
分别为下、上山影响半径，
Figure RE-GDA0002467328550000043
分别为下、上山的采深，
Figure RE-GDA0002467328550000044
为倾向工作面计算长度，D1为工作面倾向斜长，s1、s2分别为下、上山拐点偏移距，θ0为开采影响传播角，θ0＝90°-kα，式中k为小于1的常数，一般取值在0.5-0.8之间；
在传统静态概率积分模型基础上融入坐标-时间参数后，表达式1)中的矿区地质参数GP＝[m，α，H，D1，ω]为已知参数，可根据矿区地质条件和开采条件确定；概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]]为未知参数，需通过遗传算法估计得到，其中下沉系数q的取值范围为0.01-1，走向、倾向下山、倾向上山影响角正切tanβ、tanβ1、tanβ2的取值范围均为1-3.8，下、上山拐点偏移距s1、s2的取值范围均为0.05H-0.3H，开采影响传播角θ0＝90°-kα，α为矿层倾角，其中k的取值范围为0.5-0.8，起动距系数 n的取值范围为1/7-1/2；解缠后的差分干涉相位中高程改正值Δh的取值范围为-20-20m。
步骤2.2、构建InSAR视线向形变与概率积分预计参数的时间函数模型，如下：忽略地表水平移动时，将雷达视线向形变表示为表达式6)：
dLOS＝dLOS(tB)-dLOS(tA)＝[W(x，y，tB)-W(x，y，tA)]·cosθ 6)；
式中：dLOS(tB)和dLOS(tA)分别是在时间tB和tA相对于开始时间t0＝0沿雷达视线方向的形变，dLOS(t0)≡0；W(x，y，tB)和W(x，y，tA)分别在tB和tA相对于起始时间t0＝0的矿区地表动态沉降量，W(t0＝0)≡0；θ为雷达入射角；
步骤2.3、构建InSAR差分干涉相位与概率积分预计参数之间的函数模型，包含以下步骤：
步骤2.3.1、由步骤一获取的时序差分干涉相位矩阵建立如表达式7)的时间序列函数模型：
Figure RE-GDA0002467328550000045
式中：δφ是对于在tA和tB时刻获取两幅SAR图像上解缠后生成的差分干涉相位，为步骤1.4中生成，λ、θ和B⊥分别代表雷达波长、雷达入射角和垂直基线，r为相干目标与雷达卫星所在位置之间的距离，Δh为高程改正值，残差相位δφresidual＝δφNoise+ δφatm+δφNon，其中δφNoise表示噪声相位，δφatm为大气延迟相位，δφNon为高通形变分量部分；
步骤2.3.2、将表达式6)代入表达式7)中，得到InSAR相位与概率积分法预计参数的函数模型如表达式8)：
Figure RE-GDA0002467328550000051
以上技术方案中优选的，步骤三包括以下步骤：
步骤3.1、利用遗传算法估计概率积分预计参数的初始值，包括以下步骤：
步骤3.1.1、根据残差最小原则，建立算法的适应度函数如表达式9)：
fm＝||δφ′-δφ|| 9)；
式中：δφ′为由表达式8)的概率积分预计参数计算的雷达差分干涉相位，δφ为时序 InSAR计算的差分干涉相位；
步骤3.1.2、根据实验需要确定种群大小、迭代次数、种群交叉概率及变异概率；根据步骤二中概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]和高程改正值Δh的取值范围随机生成初始种群个体；
步骤3.1.3、按照步骤3.1.1中表达式9)计算目标函数的个体适应度；
步骤3.1.4、按照设定的迭代终止条件，判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，然后在重复步骤3.1.3；若满足，则输出种群中适应度函数值最小的个体；
步骤3.2、利用单纯形搜索法进行解的优化，包括以下步骤：
步骤3.2.1、遗传算法得到的概率积分预计参数值UP＝[q，tanβ，s1，s2，k，n]，将其作为单纯形搜索法初始值，选择7个顶点构成初始单纯形；
步骤3.2.2、计算单纯形各顶点的函数值fmin和它们的最大点xmax、次大点xsec和最小点xmin；
步骤3.2.3、采用表达式10)计算除最大点xmax外的点的重心xgra和反射点xref：
Figure RE-GDA0002467328550000052
其中：ξ为反射因子，其值通过试验来选择，通常取值为1；
步骤3.2.4、若fmin(xref)≤fmin(xsec)，但fmin(xref)≥fmin(xmin)，进行反射，令反射点替换最大点，组成新的单纯形，然后从步骤3.2.2开始；
若fmin(xref)＜fmin(xmin)，则沿反射方向延伸，然后从反射点和延伸点这两点中选一个函数值较小的替代最大值，并以新单纯形从步骤3.2.2开始；
若fmin(xref)＞fmin(xsec)，但fmin(xref)＞fmin(xmax)，说明此时反射点改进不大，则对反射点进行收缩，收缩后的点xsys依据表达式11)进行试算：
xsys＝xgra+γ(xref-xgra) 11)；
式中：γ为收缩因子，常取值为0.5；
若fmin(xsys)＜fmin(xmax)，则用xsys代替xmax，并以新单纯形从步骤3.2.2开始重复；
若fmin(xref)＞fmin(xmax)，此时新点收缩在xmaxxgra之间，则收缩点xshr依据表达式12)进行计算：
xshr＝xgra+γ(xgra-xmax) 12)；
若fmin(xshr)＜fmin(xmax)，则用xshr代替xmax，并以新单纯形从步骤3.2.2开始重复；满足收敛条件fm≤0.0005，结束；
步骤3.3、将步骤3.2获取的概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]代入表达式1)中，计算相干点在坐标(x，y)上的低通形变分量dLP。
以上技术方案中优选的，步骤四包括以下步骤：
步骤4.1、先通过使用均值滤波对表达式7)中的残余相位δφresidual进行空间维低通滤波，去除噪声，然后使用三角滤波对去除噪声的残余相位进行时间维高通滤波，得到时序大气延迟相位；最后将残余相位中减去大气延迟相位，得到高通形变差分干涉相位δφNon如表达式13)：
Figure RE-GDA0002467328550000061
则所有干涉图高通形变相位表示为表达式14)所示的矩阵：
Figure RE-GDA0002467328550000062
式中：tB，i和tA，i是第m个干涉对对应主影像tB和从影像tA的索引，vk为线性形变速率；
步骤4.2、求解出表达式14)中速度矢量vk，对每个时段形变速率在时间域上进行积分得到每个时段的高通形变相位值，最后通过乘以相位形变转换系数(λ/4π)得到高通形变分量dHP；
步骤4.3、将步骤3.3输出的低通形变分量dLP与步骤4.2输出的高通形变分量dHP累加，获取相干点视线向总形变量；
步骤4.4、将步骤4.3生成的相干点视线向总形变量进行地理编码，生成矿区地表垂直向时序形变。
应用本发明的技术方案，具有以下有益效果：
本发明利用时间序列InSAR技术，在传统静态概率积分模型基础上，结合地表移动变形的坐标-时间函数(the Coordinate-Time Functions，简称CT)，构建InSAR视线向形变与概率积分预计参数的时间函数(简称InSAR-CTPIM)，改进原有InSAR矿区变形监测的纯经验数学模型，不但可极大提高InSAR技术在矿区长期变形监测精度，更可使得形变估计结果解译更为合理。此外，利用InSAR相位模型方程组可直接估计出的概率积分预计参数，可为矿区开采沉陷预计提供一种新的预计参数估计方法，辅助保障矿山安全生产和保护生态环境。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为InSAR-CTPIM模型构建及预计参数解算流程图；
图2为模拟出的10个解缠后的差分干涉相位图；
图3为模拟出的时序形变结果图；
图4为概率积分预计参数计算得到的时序形变结果图；
图5为高相干点动态模型预测值与真实值对比图；
图6为InSAR-CTPIM模型预测值残差中误差图。
图3、图4及图5中均示意了第24天、第120天、第228天、第468天、第528 天及第576天的相关形变结果。
具体实施方式
以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。
实施例：
一种矿区形变监监测方法，详见图1，具体实施步骤如下：
步骤一、利用时序InSAR技术获取待监测矿区解缠后的差分干涉图，并选取高相干点，详情如下：
步骤1.1、将收集的多幅SAR卫星数据批量进行格式转换，生成可由处理软件读取的单视复数影像，并进行初始影像的预滤波，减弱影像的噪声相位。
步骤1.2、对所有影像进行基线估计。处理软件会自动选取一景影像作为影像的主影像；设置时空基线域值，根据设定的阈值进行影像的干涉组合，生成符合阈值条件的高质量干涉对，这一步骤需要根据显示出的干涉组合时空基线分布和干涉图相应质量，反复调整参数和系数，控制最终干涉组合对数，并确定最终的超级主影像。将所有的副影像配准并重采样到超级主影像，配准过程中采用多项式拟合。本实施例中采用模拟试验，模拟生成10个干涉对，后续处理采用SBAS小基线集技术生成差分干涉图。
步骤1.3、干涉图生成：这一步骤主要是利用之前组建好的干涉组合，依次进行干涉图的生成。生成干涉图的同时，各干涉对会生成对应的相干系数图。生成干涉图后借助精密轨道数据去除干涉图序列中的平地相位、轨道误差；利用外部DEM数据去除干涉图中的地形相位，并进行残余相位滤波，去除影像中的噪声相位。
相位解缠：将去除了平地、轨道、地形及噪声相位后的残余相位进行相位解缠，主要采用最小费用流法解缠，进而生成解缠后的时间序列差分干涉相位，详见图2。
步骤1.4、提取测试矿区的高相干目标点，并组建高相干目标点的时间序列差分相位矩阵。SBAS技术是利用测试区域内的高质量点进行建模分析，进而提取出形变分量，反演矿区地表形变的过程。利用强度、振幅离差和相干系数法提取高相干目标点 (本实施例中高相干点位置利用程序的随机模拟器生成)。利用提取出的矿区高相干点序列，在解缠后的差分干涉相位序列中提取出对应的相位值，组建成时间序列高相干点差分干涉矩阵(本实施例中模拟生成的干涉图如图2)。
步骤二、构建InSAR干涉相位与概率积分预计参数的时序函数模型，详情如下：
步骤2.1、结合地表移动的坐标-时间函数，将时间参数融入概率积分静态模型，构建地表下沉与概率积分参数间的时间函数模型如表达式1)：
Figure RE-GDA0002467328550000091
式中：W(x，y，t)为工作面推进时间为t时矿山地表任一点(x，y)上开采引发的形变量， w0为该地质采矿条件下的最大沉降值，w0＝ma·cosα，m为开采厚度，q为下沉系数，α为矿层倾角；表达式1)中在t时刻地表下沉W(x，t)和W(y，t)表示如下：
当矿区地下开采工作面倾向充分开采，走向主断面为有限开采时，造成的地表下沉W(x，t)的公式如下：
当x≤nH时，为表达式2)：
Figure RE-GDA0002467328550000092
当nH＜x＜1.3r时，为表达式3)：
Figure RE-GDA0002467328550000093
当1.3r≤x时，为表达式4)：
Figure RE-GDA0002467328550000094
其中：x为地表点的计算横坐标值，
Figure RE-GDA0002467328550000095
为静态下沉公式，erf 为概率积分函数，其形式为
Figure RE-GDA0002467328550000096
u为积分参数，
Figure RE-GDA0002467328550000097
为影响半径，H为开采深度，tanβ为影响角正切值，nH是地表开始移动时工作面的推进距离，n为启动系数，ω静态超前影响角，v是工作面在匀速开采条件下的推进速度，系数
Figure RE-GDA0002467328550000098
系数
Figure RE-GDA0002467328550000099
当矿区地下开采工作面走向充分开采，倾向有限开采时，造成地表下沉W(y，t)的公式如表达式5)：
Figure RE-GDA00024673285500000910
其中，y为矿区地表点的计算纵坐标，
Figure RE-GDA00024673285500000911
分别为下、上山影响半径，
Figure RE-GDA0002467328550000101
分别为下、上山的采深，
Figure RE-GDA0002467328550000102
为倾向工作面计算长度，D1为工作面倾向斜长，s1、s2分别为下、上山拐点偏移距，θ0为开采影响传播角，θ0＝90°-kα，式中k为小于1的常数，一般取值在0.5-0.8之间。
在传统静态概率积分模型基础上融入坐标-时间参数后，表达式1)中的矿区地质参数GP＝[m，α，H，D1，ω]为已知参数，可根据矿区地质条件和开采条件确定；概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]为未知参数，需通过遗传算法估计得到，其中下沉系数q的取值范围为0.01-1，走向、倾向下山、倾向上山影响角正切tanβ、tanβ1、tanβ2的取值范围均为1-3.8，下、上山拐点偏移距s1、s2的取值范围均为0.05H-0.3H，开采影响传播角θ0＝90°-kα，α为矿层倾角，其中k的取值范围为0.5-0.8，起动距系数 n的取值范围为1/7-1/2；解缠后的差分干涉相位中高程改正值Δh的取值范围为-20-20m。
步骤2.2、构建InSAR视线向形变与概率积分预计参数的时间函数模型，如下：忽略地表水平移动时，将雷达视线向形变表示为表达式6)：
dLOS＝dLOS(tB)-dLOS(tA)＝[W(x，y，tB)-W(x，y，tA)]·cosθ 6)；
式中：dLOS(tB)和dLOS(tA)分别是在时间tB和tA相对于开始时间t0＝0沿雷达视线方向的形变，dLOS(t0)≡0；W(x，y，tB)和W(x，y，tA)分别在tB和tA相对于起始时间t0＝0的矿区地表动态沉降量，W(t0＝0)≡0；θ为雷达入射角；
步骤2.3、构建InSAR差分干涉相位与概率积分预计参数之间的函数模型，包含以下步骤：
步骤2.3.1、由步骤一获取的时序差分干涉相位矩阵建立如表达式7)的时间序列函数模型：
Figure RE-GDA0002467328550000103
式中：δφ是对于在tA和tB时刻获取两幅SAR图像上解缠后生成的差分干涉相位，为步骤1.4中生成，λ、θ和B⊥分别代表雷达波长、雷达入射角和垂直基线，r为相干目标与雷达卫星所在位置之间的距离，Δh为高程改正值，残差相位δφresidual＝δφNoise+ δφatm+δφNon，其中δφNoise表示噪声相位，δφatm为大气延迟相位，δφNon为高通形变分量部分；
步骤2.3.2、将表达式6)代入表达式7)中，得到InSAR相位与概率积分法预计参数的函数模型如表达式8)：
Figure RE-GDA0002467328550000111
步骤三、先利用遗传算法估计模型的概率积分预计参数初始值，再利用单纯形搜索法对遗传算法得到的概率积分预计参数进行优化，详情如下：
步骤3.1、遗传算法估计模型的概率积分的预计参数初始值，主要包括以下步骤：
遗传算法(Genetic Algorithm，简称GA)是一类借鉴生物界的进化规律(适者生存，优胜劣汰遗传机制)演化而来的随机化搜索方法，具有全局搜寻能力和能自动获取和指导优化的搜索空间。本实施例利用遗传算法估计具体包含如下步骤：
步骤3.1.1、根据残差最小原则，建立算法的适应度函数如表达式9)：
fm＝||δφ′-δφ|| 9)；
式中：δφ′为由表达式8)的概率积分预计参数计算的雷达差分干涉相位，δφ为时序 InSAR计算的差分干涉相位；
步骤3.1.2、根据实验需要确定种群大小、迭代次数、种群交叉概率及变异概率；根据步骤2.3中概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]和高程改正值Δh的取值范围随机生成初始种群个体；
步骤3.1.3、按照步骤3.1.1中表达式9)计算目标函数的个体适应度；
步骤3.1.4、按照设定的迭代终止条件，判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，然后在重复步骤3.1.3；若满足，则输出种群中适应度函数值最小的个体；
步骤3.2、利用单纯形搜索法对遗传算法得到的概率积分预计参数初始值进行优化，具体是：
单纯形搜索法(Simplex Search Method)是一种无约束最优化的直接方法，不必计算目标函数梯度的算法。在n维空间构成一个具有n+1个顶点的单纯形，计算目标函数fmin在这些顶点的函数值并比较大小，从而判断出fmin的变化趋势，确定搜索方向，找出新点，并以fmin较小的点替换原单纯形中具有最大函数值的点，构成新单纯形；重复上述过程，找出极小值。具体包括如下步骤：
步骤3.2.1、遗传算法得到的概率积分预计参数值UP＝[q，tanβ，s1，s2，k，n]，将其作为单纯形搜索法初始值，选择7个顶点构成初始单纯形；
步骤3.2.2、计算单纯形各顶点的函数值fmin和它们的最大点xmax、次大点xsec和最小点xmin；
步骤3.2.3、采用表达式10)计算除最大点xmax外的点的重心xgra和反射点xref：
Figure RE-GDA0002467328550000121
其中：ξ为反射因子，其值通过试验来选择，通常取值为1；
步骤3.2.4、若fmin(xref)≤fmin(xsec)，但fmin(xref)≥fmin(xmin)，进行反射，令反射点替换最大点，组成新的单纯形，然后从步骤3.2.2开始；
若fmin(xref)＜fmin(xmin)，则沿反射方向延伸，然后从反射点和延伸点这两点中选一个函数值较小的替代最大值，并以新单纯形从步骤3.2.2开始；
若fmin(xref)＞fmin(xsec)，但fmin(xref)＜fmin(xmax)，说明此时反射点改进不大，则对反射点进行收缩，收缩后的点xsys依据表达式11)进行试算：
xsys＝xgra+γ(xref-xgra) 11)；
式中：γ为收缩因子，常取值为0.5；
若fmin(xsys)＜fmin(xmax)，则用xsys代替xmax，并以新单纯形从步骤3.2.2开始重复；
若fmin(xref)＞fmin(xmax)，此时新点收缩在xmaxxgra之间，则收缩点xshr依据表达式12)进行计算：
xshr＝xgra+γ(xgra-xmax) 12)；
若fmin(xshr)＜fmin(xmax)，则用xshr代替xmax，并以新单纯形从步骤3.2.2开始重复；满足收敛条件fm≤0.0005，结束。
此处目标阈值设为0.0005。
步骤3.3、将步骤3.2获取的概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]代入表达式1)中，计算相干点在坐标(x，y)上的低通形变分量dLP；
通过步骤三的计算，本实施列中获得的概率积分预计参数q＝1.53、tanβ＝ tanβ1＝tanβ2＝0.531，s1＝30.55m，s2＝16.93m，k＝0.501，n＝0.167。
步骤四、基于相干点概率积分预计参数的矿区时序形变估计，具体包括如下步骤：
步骤4.1、先通过使用均值滤波对表达式7)中的残余相位δφresidual进行空间维低通滤波，去除噪声，然后使用三角滤波对去除噪声的残余相位进行时间维高通滤波，得到时序大气延迟相位；最后将残余相位中减去大气延迟相位，可得到高通形变差分干涉相位δφNon如表达式13)：
Figure RE-GDA0002467328550000131
则所有干涉图高通形变相位表示为表达式14)所示的矩阵：
Figure RE-GDA0002467328550000132
式中：tB，i和tA，i是第m个干涉对对应主影像tB和从影像tA的索引，vk为线性形变速率；
步骤4.2、求解出式14)中速度矢量vk，对每个时段形变速率在时间域上进行积分可得到每个时段的高通形变相位值，最后通过乘以相位形变转换系数(λ/4π)可得到高通形变分量dHP；
步骤4.3、将步骤3.2输出的低通形变分量dLP与步骤4.2输出的高通形变分量dHP累加，获取相干点视线向总形变量；
步骤4.4、将步骤4.3生成的视线向形变进行地理编码，生成矿区地表垂直向时序形变，详见图4。
通过对矿区测试区域的地质水文资料、开采状态及SAR影像数据调查，本实施例中设定的概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]，详见表1：
表1概率积分预计参数真值和估计值
参数 q tanβ s<sub>1</sub>(m) s<sub>2</sub>(m) k n
真值 1.56 0.504 30.31 16.08 0.5244 0.167
估计值 1.53 0.531 30.55 16.93 0.501 0.167
在范围[-20,20]内选用高斯函数模拟高程改正值Δh，随机生成300个高相干点目标。模拟实验中采用的卫星参数均依据搭载了C波段的Sentinel-1A卫星生成的宽幅干涉模式和升轨的SAR影像参数设定。根据采集到的真实卫星数据干涉组合后，提取出其中10个干涉效果好的干涉对，将其时空基线参数作为本实施例中的时空基线参数。模拟差分干涉相位场过程中，分别加入方差为0，均方根误差为0rad-0.65rad的噪声。图2 为模拟出的10个解缠后的差分干涉相位图(噪声水平为0.65rad)；图3为根据模拟真实参数值计算得到的时序形变结果图；图4为利用概率积分预计参数估计值计算得到的时序形变结果图；图5为高相干点动态模型预测值与真实值对比图；表1为概率参数的最终估值与模拟真值之间的对比结果；图6为InSAR-CTPIM模型预测值残差中误差图。由图可见，在加入一定的噪声前提下，图3的模拟形变场与图4的估计值计算获取的形变场在时空分布上保持一致。表1、图5与图6的定量分析显示，本发明估计出的概率积分参数与真实值吻合仍然较好。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种矿区形变的监测方法，其特征在于，包括以下步骤：
步骤一、基于高相干点的时序InSAR差分干涉图生成，具体是：利用时间序列InSAR技术获取监测矿区解缠后的差分干涉相位图；
步骤二、构建InSAR干涉相位与概率积分预计参数的时间函数模型；
步骤三、概率积分预计参数估计，具体是：先利用遗传算法估计概率积分预计参数的初始值，再利用单纯形搜索法对遗传算法得到的概率积分预计参数进行优化；
步骤四、基于高相干点概率积分预计参数的时序形变估计，具体是：估计InSAR影像干涉时段内的矿区时序形变。
2.根据权利要求1所述矿区形变的监测方法，其特征在于，步骤一包括如下步骤：
步骤1.1、对M+1幅SAR数据进行干涉组合、影像配准及重采样，生成干涉图和相干图；
步骤1.2、对步骤1.1中干涉图进行去轨道、平地相位和地形相位处理，生成差分干涉图；
步骤1.3、对步骤1.2中的差分干涉图进行相位解缠，得到解缠后的差分干涉图；
步骤1.4、利用步骤1.1中生成的相干图和步骤1.3获得的解缠后的差分干涉图进行高相干点提取，生成基于高相干点目标的时序差分干涉相位矩阵。
3.根据权利要求1所述的矿区形变的监测方法，其特征在于，步骤二具体是：
步骤2.1、结合地表移动的坐标-时间函数，将时间参数融入概率积分静态模型，构建地表下沉与概率积分参数间的时间函数模型如表达式1)：
Figure FDA0002393465170000011
式中：t为开采时间，W(x，y，t)为工作面推进时间为t时矿山地表任一点(x，y)上开采引发的形变量，w0为该地质采矿条件下的最大沉降值，w0＝mq·cosα，m为开采厚度，q为下沉系数，α为矿层倾角；表达式1)中在t时刻地表下沉W(x，t)和W(y，t)表示如下：
当矿区地下开采工作面倾向充分开采，走向主断面为有限开采时，造成的地表下沉W(x，t)的公式如下：
当x≤nH时，为表达式2)：
Figure FDA0002393465170000021
当nH＜x＜1.3r时，为表达式3)：
Figure FDA0002393465170000022
当1.3r≤x时，为表达式4)：
Figure FDA0002393465170000023
其中：x为地表点的计算横坐标值，
Figure FDA0002393465170000024
为静态下沉公式，erf为概率积分函数，其形式为
Figure FDA0002393465170000025
u为积分参数，
Figure FDA0002393465170000026
为影响半径，H为开采深度，tanβ为影响角正切值，nH是地表开始移动时工作面的推进距离，n为起动距系数，ω静态超前影响角，v是工作面在匀速开采条件下的推进速度，系数
Figure FDA0002393465170000027
系数
Figure FDA0002393465170000028
当矿区地下开采工作面走向充分开采，倾向有限开采时，地表下沉W(y，t)的公式如表达式5)：
Figure FDA0002393465170000029
其中，y为矿区地表点的计算纵坐标；
Figure FDA00023934651700000210
分别为下、上山影响半径，β1、β2分别为倾向下山、上山影响角，
Figure FDA00023934651700000211
分别为下、上山的采深；
Figure FDA00023934651700000212
为倾向工作面计算长度，D1为工作面倾向斜长，s1、s2分别为下、上山拐点偏移距，θ0为开采影响传播角，θ0＝90°-kα，式中k为小于1的常数；
步骤2.2、构建InSAR视线向形变与概率积分预计参数的时间函数模型，如下：忽略地表水平移动时，将雷达视线向形变表示为表达式6)：
dLOS＝dLOS(tB)-dLOS(tA)＝[W(x，y，tB)-W(x，y，tA)]·cosθ 6)；
式中：dLOS(tB)和dLOS(tA)分别是在时间tB和tA相对于开始时间t0＝0沿雷达视线方向的形变，dLOS(t0)≡0；W(x，y，tB)和W(x，y，tA)分别在tB和tA相对于起始时间t0＝0的矿区地表动态沉降量，W(t0＝0)≡0；θ为雷达入射角；
步骤2.3、构建InSAR差分干涉相位与概率积分预计参数之间的函数模型，包含以下步骤：
步骤2.3.1、由步骤一获取的时序差分干涉相位矩阵建立如表达式7)的时间序列函数模型：
Figure FDA0002393465170000031
式中：δφ是对于在tA和tB时刻获取两幅SAR图像上解缠后生成的差分干涉相位，为步骤1.4中生成，λ、θ和B⊥分别代表雷达波长、雷达入射角和垂直基线，r为相干目标与雷达卫星所在位置之间的距离，Δh为高程改正值，残差相位δφresidual＝δφNoise+δφatm+δφNon，其中δφNoise表示噪声相位，δφatm为大气延迟相位，δφNon为高通形变分量部分；
步骤2.3.2、将表达式6)代入表达式7)中，得到InSAR相位与概率积分法预计参数的函数模型如表达式8)：
Figure FDA0002393465170000032
4.根据权利要求3所述的矿区形变的监测方法，其特征在于，在传统静态概率积分模型基础上融入坐标-时间参数后，表达式1)中的矿区地质参数GP＝[m，α，H，D1，ω]为已知参数，根据矿区地质条件和开采条件确定；概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]为未知参数，需通过遗传算法估计得到，其中下沉系数q的取值范围为0.01-1，走向、倾向下山、倾向上山影响角正切tanβ、tanβ1、tanβ2的取值范围均为1-3.8，下、上山拐点偏移距s1、s2的取值范围均为0.05H-0.3H，开采影响传播角θ0＝90°-kα，α为矿层倾角，其中k的取值范围为0.5-0.8，启动距系数n的取值范围为1/7-1/2；解缠后的差分干涉相位中高程改正值Δh的取值范围为-20-20m。
5.根据权利要求3所述矿区形变的监测方法，其特征在于，步骤三包括以下步骤：
步骤3.1、利用遗传算法估计概率积分预计参数的初始值，包括以下步骤：
步骤3.1.1、根据残差最小原则，建立算法的适应度函数如表达式9)：
fm＝||δφ′-δφ|| 9)；
式中：δφ＇为由表达式8)的概率积分预计参数计算的雷达差分干涉相位，δφ为时序InSAR计算的差分干涉相位；
步骤3.1.2、根据实验需要确定种群大小、迭代次数、种群交叉概率及变异概率；根据步骤二中概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]和高程改正值Δh的取值范围随机生成初始种群个体；
步骤3.1.3、按照步骤3.1.1中表达式9)计算目标函数的个体适应度；
步骤3.1.4、按照设定的迭代终止条件，判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，然后在重复步骤3.1.3；若满足，则输出种群中适应度函数值最小的个体；
步骤3.2、利用单纯形搜索法进行解的优化，包括以下步骤：
步骤3.2.1、遗传算法得到的概率积分预计参数值UP＝[q，tanβ，s1，s2，k，n]，将其作为单纯形搜索法初始值，选择7个顶点构成初始单纯形；
步骤3.2.2、计算单纯形各顶点的函数值fmin和它们的最大点xmax、次大点xsec和最小点xmin；
步骤3.2.3、采用表达式10)计算除最大点xmax外的点的重心xgra和反射点xref：
Figure FDA0002393465170000041
其中：ξ为反射因子，其值通过试验来选择，通常取值为1；
步骤3.2.4、若fmin(xref)≤fmin(xsec)，但fmin(xref)≥fmin(xmin)，进行反射，用反射点替换最大点，组成新的单纯形，然后从步骤3.2.2开始；
若fmin(xref)＜fmin(xmin)，则沿反射方向延伸，从反射点和延伸点这两点中选一个函数值较小的替代最大值，并以新单纯形从步骤3.2.2开始重复；
若fmin(xref)＞fmin(xsec)，但fmin(xref)＜fmin(xmax)，说明此时反射点改进不大，则对反射点进行收缩，收缩后的点xsys依据表达式11)进行试算：
xsys＝xgra+γ(xref-xgra) 11)；
式中：γ为收缩因子，常取值为0.5；
若fmin(xsys)＜fmin(xmax)，则用xsys代替xmax，并以新单纯形从步骤3.2.2开始重复；
若fmin(xref)＞fmin(xmax)，此时新点收缩在xmaxxgra之间，则收缩点xshr依据表达式12)进行计算：
xshr＝xgra+γ(xgra-xmax) 12)；
若fmin(xshr)＜fmin(xmax)，则用xshr代替xmax，并以新单纯形从步骤3.2.2开始重复；满足收敛条件fm≤0.0005，结束；
步骤3.3、将步骤3.2获取的概率积分预计参数UP＝[q，tanβ，s1，s2，k，n]代入表达式1)中，计算相干点在坐标(x，y)上的低通形变分量dLP。
6.根据权利要求5所述矿区形变的监测方法，其特征在于，步骤四包括以下步骤：
步骤4.1、先通过使用均值滤波对表达式7)中的残余相位δφresidual进行空间维低通滤波和去除噪声，然后使用三角滤波对去除噪声的残余相位进行时间维高通滤波，得到时序大气延迟相位；最后将残余相位中减去大气延迟相位，得到高通形变差分干涉相位δφNon如表达式13)：
Figure FDA0002393465170000051
则所有干涉图高通形变相位表示为表达式14)所示的矩阵：
Figure FDA0002393465170000052
式中：tB，i和tA，i是第m个干涉对对应主影像tB和从影像tA的索引，vk为线性形变速率；
步骤4.2、求解出表达式14)中速度矢量vk，对每个时段形变速率在时间域上进行积分得到每个时段的高通形变相位值，最后通过乘以相位形变转换系数(λ/4π)得到高通形变分量dHP；
步骤4.3、将步骤3.3输出的低通形变分量dLP与步骤4.2输出的高通形变分量dHP累加，获取相干点视线向总形变量；
步骤4.4、将步骤4.3生成的相干点视线向总形变量进行地理编码，生成矿区地表垂直向时序形变。
CN202010122691.8A 2020-02-27 2020-02-27 一种矿区形变的监测方法 Active CN111323776B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010122691.8A CN111323776B (zh) 2020-02-27 2020-02-27 一种矿区形变的监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010122691.8A CN111323776B (zh) 2020-02-27 2020-02-27 一种矿区形变的监测方法

Publications (2)

Publication Number Publication Date
CN111323776A true CN111323776A (zh) 2020-06-23
CN111323776B CN111323776B (zh) 2021-04-13

Family

ID=71165362

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010122691.8A Active CN111323776B (zh) 2020-02-27 2020-02-27 一种矿区形变的监测方法

Country Status (1)

Country Link
CN (1) CN111323776B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111859786A (zh) * 2020-07-03 2020-10-30 安徽理工大学 一种改进动态预计模型约束的全尺度梯度开采沉陷D-InSAR三维预测方法
CN112184902A (zh) * 2020-09-21 2021-01-05 东华理工大学 一种面向越界开采识别的地下开采面反演方法
CN112857312A (zh) * 2021-03-31 2021-05-28 中铁上海设计院集团有限公司 依据沉降速率的不同时序差分干涉地面沉降测量融合方法
CN113064170A (zh) * 2021-03-29 2021-07-02 长安大学 一种基于时序InSAR技术的膨胀土区域地表形变监测方法
CN113091600A (zh) * 2021-04-06 2021-07-09 长沙理工大学 一种利用时序InSAR技术监测软土地基形变的监测方法
CN114325693A (zh) * 2021-12-19 2022-04-12 南京市测绘勘察研究院股份有限公司 一种基于InSAR时序形变结果的采空区中心形变预测方法
CN115164828A (zh) * 2022-05-05 2022-10-11 四川省地质工程勘察院集团有限公司 一种地表沉降预测方法、系统、装置及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091675A (zh) * 2013-01-11 2013-05-08 中南大学 一种基于InSAR技术的矿区开采监测方法
CN105926569A (zh) * 2016-04-28 2016-09-07 河北地质大学 一种基于沉降监测数据的煤矿老采空区场地稳定性定量评价方法
CN108802727A (zh) * 2018-04-13 2018-11-13 长沙理工大学 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法
CN109471107A (zh) * 2018-11-27 2019-03-15 长沙理工大学 基于PS-InSAR技术的桥梁永久形变分析方法
EP3540462A1 (en) * 2018-03-14 2019-09-18 Elta Systems Ltd. Coherence change detection techniques

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091675A (zh) * 2013-01-11 2013-05-08 中南大学 一种基于InSAR技术的矿区开采监测方法
CN105926569A (zh) * 2016-04-28 2016-09-07 河北地质大学 一种基于沉降监测数据的煤矿老采空区场地稳定性定量评价方法
EP3540462A1 (en) * 2018-03-14 2019-09-18 Elta Systems Ltd. Coherence change detection techniques
CN108802727A (zh) * 2018-04-13 2018-11-13 长沙理工大学 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法
CN109471107A (zh) * 2018-11-27 2019-03-15 长沙理工大学 基于PS-InSAR技术的桥梁永久形变分析方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
ZE FA YANG等: "InSAR-Based Model Parameter Estimation of Probability Integral Method and Its Application for Predicting Mining-Induced Horizontal and Vertical Displacements", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
ZEFA YANG等: "An InSAR-Based Temporal Probability Integral Method and its Application for Predicting Mining-Induced Dynamic Deformations and Assessing Progressive Damage to Surface Buildings", 《IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING》 *
成晓倩等: "联合ＤＩｎＳＡＲ和ＰＩＭ 技术的沉陷特征模拟和时序分析", 《中国矿业大学学报》 *
汪磊等: "融合概率积分模型与D－InSAR 的开采沉陷预计", 《金属矿山》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111859786A (zh) * 2020-07-03 2020-10-30 安徽理工大学 一种改进动态预计模型约束的全尺度梯度开采沉陷D-InSAR三维预测方法
CN112184902A (zh) * 2020-09-21 2021-01-05 东华理工大学 一种面向越界开采识别的地下开采面反演方法
CN113064170A (zh) * 2021-03-29 2021-07-02 长安大学 一种基于时序InSAR技术的膨胀土区域地表形变监测方法
CN112857312A (zh) * 2021-03-31 2021-05-28 中铁上海设计院集团有限公司 依据沉降速率的不同时序差分干涉地面沉降测量融合方法
CN113091600A (zh) * 2021-04-06 2021-07-09 长沙理工大学 一种利用时序InSAR技术监测软土地基形变的监测方法
CN114325693A (zh) * 2021-12-19 2022-04-12 南京市测绘勘察研究院股份有限公司 一种基于InSAR时序形变结果的采空区中心形变预测方法
CN115164828A (zh) * 2022-05-05 2022-10-11 四川省地质工程勘察院集团有限公司 一种地表沉降预测方法、系统、装置及存储介质

Also Published As

Publication number Publication date
CN111323776B (zh) 2021-04-13

Similar Documents

Publication Publication Date Title
CN111323776B (zh) 一种矿区形变的监测方法
CN109918781B (zh) 一种钻井水溶盐矿开采沉陷InSAR预计方法
CN108663017A (zh) 一种监测城市地铁沿线地表沉降的方法
Fan et al. A model for extracting large deformation mining subsidence using D-InSAR technique and probability integral method
Baker et al. The transition from complex craters to multi‐ring basins on the Moon: Quantitative geometric properties from Lunar Reconnaissance Orbiter Lunar Orbiter Laser Altimeter (LOLA) data
Zhang et al. Monitoring of urban subsidence with SAR interferometric point target analysis: A case study in Suzhou, China
Hani et al. A method for computation of surface roughness of digital elevation model terrains via multiscale analysis
Fedi et al. 2D continuous wavelet transform of potential fields due to extended source distributions
Ai et al. High-precision ice-flow velocities from ground observations on Dalk Glacier, Antarctica
Hu et al. Methods for monitoring fast and large gradient subsidence in coal mining areas using SAR images: A review
Liu et al. [Retracted] Interpolation Parameters in Inverse Distance‐Weighted Interpolation Algorithm on DEM Interpolation Error
Ma et al. Decision-making fusion of InSAR technology and offset tracking to study the deformation of large gradients in mining areas-Xuemiaotan mine as an example
Zhang et al. Predictable condition analysis and prediction method of SBAS-InSAR coal mining subsidence
Ochoa et al. Analysis and correction of digital elevation models for plain areas
Puniach et al. Determination of the coefficient of proportionality between horizontal displacement and tilt change using UAV photogrammetry
Yang et al. Resolving 3-D Mining Displacements From Multi-Track InSAR by Incorporating With a Prior Model: The Dynamic Changes and Adaptive Estimation of the Model Parameters
Waelder An application of the fuzzy theory in surface interpolation and surface deformation analysis
Siska et al. Predicting ordinary kriging errors caused by surface roughness and dissectivity
Li et al. A Global-Scale DEM Elevation Correction Model Using ICESat-2 Laser Altimetry Data
CN111859786B (zh) 一种全尺度梯度开采沉陷D-InSAR三维预测方法
Zhang et al. A novel approach to investigating 3d fracture Connectivity in Ultrahigh steep rock slopes
Becek et al. Identifying land subsidence using global digital elevation models
Wieczorek EVALUATION OF DEFORMATIONS IN THE URBAN AREA OF OLSZTYN USING SENTINEL-1 SAR INTERFEROMETRY.
Schiavone et al. Near‐station topographic masses correction for high‐accuracy gravimetric prospecting
Beshr et al. Using modified inverse distance weight and principal component analysis for spatial interpolation of foundation settlement based on geodetic observations

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
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20220804

Address after: Room 03 and room 04, 10 / F, business office building, area B, Tianxin Software Industrial Park, 66-68 Xinling Road, Tianxin District, Changsha City, Hunan Province, 410000

Patentee after: Hunan tianxianghe Information Technology Co.,Ltd.

Address before: Changsha University of science and technology

Patentee before: CHANGSHA University OF SCIENCE AND TECHNOLOGY