CN112505699A - 融合InSAR和PSO反演地下采空区位置参数的方法 - Google Patents

融合InSAR和PSO反演地下采空区位置参数的方法 Download PDF

Info

Publication number
CN112505699A
CN112505699A CN202011351515.8A CN202011351515A CN112505699A CN 112505699 A CN112505699 A CN 112505699A CN 202011351515 A CN202011351515 A CN 202011351515A CN 112505699 A CN112505699 A CN 112505699A
Authority
CN
China
Prior art keywords
deformation
insar
los
goaf
individual
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
CN202011351515.8A
Other languages
English (en)
Other versions
CN112505699B (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.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
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 China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN202011351515.8A priority Critical patent/CN112505699B/zh
Publication of CN112505699A publication Critical patent/CN112505699A/zh
Application granted granted Critical
Publication of CN112505699B publication Critical patent/CN112505699B/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
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Theoretical Computer Science (AREA)
  • Biophysics (AREA)
  • Mathematical Physics (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Computational Linguistics (AREA)
  • Biomedical Technology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Alarm Systems (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种融合InSAR和PSO反演地下采空区位置参数的方法，适用于矿产地质监测与防护领域。利用InSAR方法获取矿区雷达视线向地表形变；确定采空区参数的取值范围；初始化粒子群中每个个体的位置与速度信息；将粒子群的位置坐标带入到概率积分模型获取地表竖向下沉和东西向、南北向水平移动，根据雷达成像几何将变形转换为LOS向形变；将转换后的LOS向地表变形与InSAR解算的LOS向地表形变带入适应度函数计算适应度值，得到每个粒子最优解和群体的全局最优解并记录；求取迭代中种群最优解对应的适应度值。其操作简单，反演地下采空区位置参数的效果好，可为地下煤炭资源非法开采、煤火区位置反演提供技术支持，同时也拓展了InSAR技术的应用空间。

Description

融合InSAR和PSO反演地下采空区位置参数的方法
技术领域
本发明涉及一种反演地下采空区位置参数的方法，尤其适用于矿区开采沉陷监测使用的一种融合InSAR和PSO反演地下采空区位置参数的方法，属于地理监测领域。
技术背景
我国的地域幅员辽阔，拥有各种自然资源，其中煤炭资源是现在我国主要依赖的能源。正因如此，在我国很多地区存在对煤炭资源进行私挖乱采的情况，这一方面造成国家资源的浪费，损害国有资产；另一方面地下非法开采存在很大的安全隐患，不仅威胁开采人员的生命财产安全，而且还会引发含水层破坏、山体滑坡、地表塌陷和建构筑物损毁等矿山地质环境灾害。因此，为规范煤炭开采的行为，使我国的煤炭资源可以有序开采，提高监测、检查效率和准确度，同时避免地下非法开采形成较大的地质灾害，有必要对地下煤炭资源开采位置及范围进行反演，获取采空区的煤层走向方位角α、煤层走向长度l、煤层倾向长度L、煤层开采深度H、煤层开采厚度m的精确定位。
目前，矿区地下开采监测和采空区的定位方法主要包括：人为巡查、地球物理勘探(如：电法勘探、地震勘探、重力勘探、放射性勘探等)及InSAR技术监测等。人为巡查费时、费力，缺乏具体的调查位置，并且无法做到实时监测。地球物理勘探技术实现复杂，花费较高，且探测区域太小，灵动性不够，效率低，在开采深度较大时，探测精度较低。InSAR技术有着全天候、全天时、覆盖范围广、成本低、精度高等优点，将该技术应用于矿区开采监测可以降低监测成本，由于其覆盖范围较大的特点，可以实现大范围的监测，提高了监管效率，但该技术只能获取地表沉降，难以获取地下开采位置参数。因此，为解决现有地下采空区位置反演方法存在的问题，本发明提出了一种融合InSAR和PSO算法反演地下采空区位置参数的方法，其思想是利用InSAR解算的地表形变作为真值，通过粒子群算法(PSO)不断优化更新地下采空区位置参数，结合开采沉陷理论模型得到多组地表形变作为观测值，比较地表形变真值与观测值间的误差，取最小误差时的采空区位置参数作为结果。
发明内容
针对上述技术问题，本发明提供一种方案简单，使用效果好，反演精度高的融合InSAR和PSO反演地下采空区位置参数的方法。
为实现上述技术目的，本发明的融合InSAR和PSO反演地下采空区位置参数的方法，首先利用InSAR合成孔径雷达干涉测量技术获取矿区雷达视线向地表形变信息，记为LOS；结合开采沉陷理论确定采空区的参数取值范围，并根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》确定研究区域概率积分法预计参数；然后通过以上信息获取地表沉陷形变区域的地表形变点的数量，根据地表形变点的密度决定参与采空区反演的形变点数量，将参与采空区反演的形变点数量作为参数输入到经典的PSO粒子群优化算法中，并将待反演的采空区位置参数作为个体的位置信息初始化粒子群中每个个体的位置与速度信息；将粒子群的位置坐标带入到概率积分模型获取地表竖向下沉、东西向和南北向水平移动并转换为雷达视线向形变los；将los与InSAR解算的LOS代入适应度函数计算适应度值，判断其是否满足迭代的停止条件，如果满足把该次迭代的种群最优解作为地下采空区真实参数输出，不满足迭代停止条件则对种群中每个个体的速度、位置进行更新，继续迭代，直到满足迭代停止条件为止，并输出此时的种群最优解作为地下采空区位置参数。
具体步骤如下：
S1利用InSAR技术解算被测矿区地表各点的视线向形变场，记为LOS；
S2根据开采沉陷理论中所述的拐点位置、边界角、地表形变对称特征，结合由LOS形变场分布给出采空区位置参数取值范围，包括：被测区的开采深度H、工作面走向及倾向长度l和L、采厚m、煤层走向坐标方位角α的取值范围，并根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》确定研究区域概率积分法预计参数；
S3由InSAR技术最终获取的地表沉陷形变区域的形变点的数量，并根据所获取地表形变点的密度决定最终参与采空区反演的形变点数量，使用经典粒子群算法的参数，初始化粒子群中每个个体的位置与速度信息；根据上述最终得到参与反演的地表点数量确定粒子群算法种群大小N；初始化粒子群中每个个体的速度与位置信息，其中个体的位置信息就是待反演的采空区位置参数，包括：开采深度、工作面走向及倾向长度、采厚、煤层走向坐标方位角。
S4将种群中每个个体的位置信息代入开采沉陷理论模型概率积分法中，计算得到地表各点位置处的三维形变场，包括竖向下沉W和南北、东西向水平移动UN、UE；结合InSAR三维形变分解模型将W、UN、UE合成地表各点雷达视线向形变场，记为：los；
所述InSAR三维形变分解模型由下式表征：
los＝Wcosθ-sinθ[UNcos(αh-3π/2)+UEsin(αh-3π/2)]
其中，los为由开采沉陷计算的三维形变合成的雷达视线向形变场；θ为雷达卫星入射角，αh为卫星航向角，W为竖向下沉，UN为南北向水平移动、UE为东西向水平移动；
S5将由InSAR技术解算的LOS视为真值，将概率积分模型计算地表三维形变后合成的los作为观测值，把LOS和los代入到适应度函数中，求得适应度函数值v，v越小则反演结果越优，记录每次迭代中同一个体的最优值和整个种群的最优值；
S6对粒子群算法的种群中每个个体的速度、位置进行更新；
粒子群每个个体新的速度公式如下：
Figure BDA0002800838780000031
粒子群每个个体位置依据速度进行转移的公式如下；
Figure BDA0002800838780000032
式中，
Figure BDA0002800838780000033
是第t次迭代第i个粒子的j维速度；
Figure BDA0002800838780000034
是第t次迭代第i个粒子的j维位置；pbest是个体历史最优解；gbest是种群中的最优状态；c1、c2是种群的学习因子；w是速度在更新时上次迭代时的速度占这次迭代的速度的权重；r1和r2是介于[0，1]之间的随机数；
S7寻求一个满足要求的种群最优值作为地下采空区真实参数，判断种群最优解gbest对应的适应度值是否满足迭代停止的预设阀值，不满足条件则重复步骤S4到步骤S6，直到满足迭代停止的预设阀值，并输出此时的种群最优解作为地下采空区位置参数。
步骤S3中的粒子群中每个个体的位置的初始化使用混沌方式产生，其具体的生成方式为随机生成一个随机数矩阵，确定一个种群第一个个体的初始化位置，种群的其余个体利用前一个个体的随机数矩阵和Logistic映射产生初始位置时的随机数矩阵，依次执行下去直到产生所有个体的初始位置；
计算公式如下：
Figure BDA0002800838780000035
式中，up，low是被测区的开采深度H、工作面走向l、倾向长度L、采厚m、煤层走向坐标方位角α的取值范围的上下边界；rand(k)是第k个个体产生初始位置时对应的一个1×5的介于[0，1]之间的随机数矩阵；pos(k)时第k个个体的初始化位置；z是一个介于(2，4]之间的常数。
步骤S4中的概率积分法计算竖向下沉W和南北、东西向水平移动UN、UE的计算方法为：
W＝W0(x)W0(y)/W0
Figure BDA0002800838780000041
Figure BDA0002800838780000042
Figure BDA0002800838780000043
Figure BDA0002800838780000044
其中：
W0＝mqcosα
Figure BDA0002800838780000045
Figure BDA0002800838780000046
投影长度及主要影响半径的计算方法为：
Figure BDA0002800838780000047
l＝D1-s1-s3
Figure BDA0002800838780000048
其中，m是煤层厚度；q是下沉系数；α是煤层倾角；x、y是地表任意点的矿区坐标系坐标；r、r1、r2是走向、倾向上山、倾向下山主要影响半径；H是煤层采深；tanβ是主要影响角正切；l、L分别为煤层的走向和倾向投影到水平面上的计算长度；D1，D2是工作面走向及倾向长度；s1，s2，s3，s4是左上右下拐点偏移距；θ是开采影响传播角；b是水平移动系数；erf是误差函数；
Figure BDA0002800838780000049
为煤层走向方位角。
步骤S4所述的获取地表形变点的密度决定最终参与采空区反演的形变点数量，InSAR获取的观测点的数据量和密度很大，最终用于反演数据点的密度和传统测量的观测站设置的密度尽量保持一致，当InSAR获取的数据量密度不能达到下表要求时，直接按照InSAR获取的数据量密度进行反演；最终参与采空区位置反演的地表形变点密度见表1；
Figure BDA0002800838780000051
步骤S5所述的适应度函数v的表达式为：
Figure BDA0002800838780000052
式中，los(k)和LOS(k)分别为地表第k个观测点反演三维形变后由概率积分法模型合成的视线向形变和InSAR观测得到的时序视线向形变；n为最终参与反演的地表点个数。
所述步骤S6中的w为动态权重，在迭代开始到结束的过程中w是遵循线性变化的其变化的规则为：
w(t)＝wmax-(wmax-wmin)·t/T
其中，wmin和wmax是动态权重的最小值和最大值；w(t)是第t次迭代过程中的权重取值；T是总迭代数。
有益效果：
本方法根据实际矿区沉陷变形情况和真实地表沉陷评估所需材料，在不损失反演参数精度的前提下将原有的8个参数，精炼为现有的6个有益参数；利用PSO算法寻优使实际采空区参数的反演精度提升更高，同时该算法的使用提升了反演的效率，降低了反演的复杂度，减少了反演的时间，反演精度见本发明后续案例；将InSAR数据引入到参数反演中，增加了初始数据量，提高了反演结果的可信度，节约了传统测量数据采集的人力物力。
附图说明
图1为本发明一种融合InSAR和PSO算法反演地下采空区位置参数的方法流程图；
图2为本发明的矿区沉降和工作位置相对关系图；
图3为本发明的种群个体最优迭代平均值进化曲线图；
图4为本发明种群最优迭代进化曲线图。
具体实施方式
以下将结合具体实施过程对本发明做进一步详细说明：
如图1所示，本发明的融合InSAR和PSO反演地下采空区位置参数的方法，首先利用InSAR合成孔径雷达干涉测量技术获取矿区雷达视线向地表形变信息，记为LOS；结合开采沉陷理论确定采空区的参数取值范围，并根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》确定研究区域概率积分法预计参数；然后通过以上信息获取地表沉陷形变区域的地表形变点的数量，根据地表形变点的密度决定参与采空区反演的形变点数量，将参与采空区反演的形变点数量作为参数输入到经典的PSO粒子群优化算法中，并将待反演的采空区位置参数作为个体的位置信息初始化粒子群中每个个体的位置与速度信息；将粒子群的位置坐标带入到概率积分模型获取地表竖向下沉、东西向和南北向水平移动并转换为雷达视线向形变los；将los与InSAR解算的LOS代入适应度函数计算适应度值，判断其是否满足迭代的停止条件，如果满足把该次迭代的种群最优解作为地下采空区真实参数输出，不满足迭代停止条件则对种群中每个个体的速度、位置进行更新，继续迭代，直到满足迭代停止条件为止，并输出此时的种群最优解作为地下采空区位置参数。
具体的某地倾斜煤层，该煤层的走向长度D1＝750m；倾向长度D2＝301m；煤层的走向方位角
Figure BDA0002800838780000061
工作面中心坐标Xc＝3525.4685m；Yc＝3594.9769m；煤层倾角α＝15°；平均开采深度为H＝600m；煤层开采厚度m＝1800mm。
步骤1：利用InSAR方法获取矿区雷达视线向(LOS)形变，形变示意图如图2所示，图中等值线为开采沉陷量，矩形线条为开采工作面范围，从图可以看出工作面开采位置与地表沉陷间的关系，可依据开采沉陷理论估计地下工作面位置参数的取值范围；数据导出结果为一个n×3的矩阵，第一列为坐标x，第二列为坐标y第三列为该点对应的LOS形变值，n代表使用遥感方法观测到了多少个这样的点；
步骤2：结合开采沉陷理论和形变特征给定采空区的参数取值范围；
对该煤层的采空区情况进行预测时在，根据下沉盆地的实际下沉情况，在范围(90°，150°)内随机取煤层走向方位角的值，在范围(500m，900m)内随机取煤层走向的计算长度，在范围(200m，300m)内随机取煤层倾向的计算长度。同时可以在范围(400m，800m)内取随机取平均高程的值，在范围(1000mm，4000mm)内随机取开采厚度的值。根据规则煤层的实际下沉情况，矿区的开采沉陷的等高线位置是呈现对称分布的原则，直接得到矿区中心坐标，根据当地的煤层地质情况，可以直接获取煤层倾角的大小，矿区位置的示意图如图2。根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》得到该矿区概率积分法预计参数如下：下沉系数q＝0.61；水平移动系数b＝0.3；主要影响角正切tanβ＝2；开采影响传播角θ0＝87.2°，走向左右拐点偏移距s1＝s3＝20m，倾向上拐点偏移距s2＝20m，倾向下拐点偏移距s4＝26m；
步骤3：初始化粒子群中每个个体的位置与速度信息；
鉴于煤层深度在(400m，800m)内，模拟数据的地表点选取，遵循相邻点距离控制在至少35m，最终参与运算的地面预计的个数为2337个点，为了可以让粒子群的种群数目可以均匀的平均分配到整个种群的基础上，对整个粒子种群进行检索，对粒子群的数目N选取为100；为了让种群在确定的种群数的情况下靠近最优解的能力达到最好，同时保证种群的个体之间不丧失信息共享机制故这里学习因子取c1＝c2＝1.5。
粒子群中每个个体的位置的始化使用混沌的方式产生，其生成方式如下：
Figure BDA0002800838780000071
如果粒子群每个个体的位置
Figure BDA0002800838780000072
则up＝[900，300，150，800，4000]，low＝[500，200，90，400，1000]，位置rand(k)是第k个个体产生初始位置时对应的一个1×5的介于[0，1]之间，这里z取4。
先随机生成一个随机数矩阵，先确定一个种群第一个个体的初始化位置，种群的其余个体利用前一个个体的随机数矩阵和Logistic映射产生生成初始位置时的随机数矩阵，依次执行下去直到产生所有个体的初始位置，得到种群的初始位置x，其矩阵大小为100×5。种群的初始速度v定为一个100×5的0矩阵。
步骤4：结合粒子群种群信息和概率积分模型合成雷达视线向形变；
初始位置矩阵x中，第一个粒子的位置如果为
Figure BDA0002800838780000073
根据x1中的位置煤层厚度m、假设中的煤层倾角α＝15°、步骤2中下沉系数q＝0.61可以求：
W0＝mqcosα
根据步骤1，可以知道2337个点的x，y，根据粒子位置D1，D2，H可以求得：
l＝D1-s1-s3
Figure BDA0002800838780000081
Figure BDA0002800838780000082
Figure BDA0002800838780000083
Figure BDA0002800838780000084
根据步骤1，2的其他参数信息最终可以求得：
W(x，y)＝W0(x)W0(y)/W0
Figure BDA0002800838780000085
Figure BDA0002800838780000086
Figure BDA0002800838780000087
Figure BDA0002800838780000088
最终可以合成2337个地表点，由InSAR三维形变分解模型由下式表征的los形变，其中αh＝349°：
los＝Wcosθ-sinθ[UNcos(αh-3π/2)+UEsin(αh-3π/2)]
粒子群的100个个体都按照上面的方法分别计算得到各自2337个地表点的los形变。
步骤5：求取种群每个个体的适应度函数值v，比较v得到个体与种群最优解并记录；
Figure BDA0002800838780000089
如果在v是对单个粒子而言在每次迭代中对于目标函数的最优解叫做pbest，如果在v是对粒子群而言在每次迭代的对于目标函数的最优解叫做gbest：
gbest＝min(pbest)
步骤6：对种群中每个个体的速度进行更新、位置进行转移；
粒子群每个个体的新的速度公式如下：
Figure BDA0002800838780000091
粒子群每个个体位置依据速度进行转移的公式如下；
Figure BDA0002800838780000092
Figure BDA0002800838780000093
是第t次迭代第i个粒子的j维速度；
Figure BDA0002800838780000094
是第t次迭代第i个粒子的j维位置；pbest个体历史最优解；gbest种群中的最优状态；c1、c2是种群的学习因子；w是速度在更新时上次迭代时的速度占这次迭代的速度的权重。
表2：反演结果对照表
Figure BDA0002800838780000095
步骤7：寻求一个满足要求的种群最优值作为地下采空区真实参数，判断种群最优解gbest是否满足迭代的停止条件，不满足条件重复步骤4到6，迭代过程得个体最优解计算得适应度值平均值变化如图3，把每次迭代的个体最优值取平均值，并绘图展示，反映了算法的种群整体的迭代收敛情况；迭代过程得最优解计算得适应度值变化如图4，把每次迭代的种群体最优值绘图展示，反映了算法的最优个体的迭代收敛情况。得满足条件迭代停止输出结果，这里的迭代误差设置为单个点5mm，预计精度见表2。
表2：反演结果对照表
Figure BDA0002800838780000096

Claims (7)

1.一种融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于：首先利用InSAR合成孔径雷达干涉测量技术获取矿区雷达视线向地表形变信息，记为LOS；结合开采沉陷理论确定采空区的参数取值范围，并根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》确定研究区域概率积分法预计参数；然后通过以上信息获取地表沉陷形变区域的地表形变点的数量，根据地表形变点的密度决定参与采空区反演的形变点数量，将参与采空区反演的形变点数量作为参数输入到经典的PSO粒子群优化算法中，并将待反演的采空区位置参数作为个体的位置信息初始化粒子群中每个个体的位置与速度信息；将粒子群的位置坐标带入到概率积分模型获取地表竖向下沉、东西向和南北向水平移动并转换为雷达视线向形变los；将los与InSAR解算的LOS代入适应度函数计算适应度值，判断其是否满足迭代的停止条件，如果满足把该次迭代的种群最优解作为地下采空区真实参数输出，不满足迭代停止条件则对种群中每个个体的速度、位置进行更新，继续迭代，直到满足迭代停止条件为止，并输出此时的种群最优解作为地下采空区位置参数。
2.根据权利要求1所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于具体步骤如下：
S1利用InSAR技术解算被测矿区地表各点的视线向形变场，记为LOS；
S2根据开采沉陷理论中所述的拐点位置、边界角、地表形变对称特征，结合由LOS形变场分布给出采空区位置参数取值范围，包括：被测区的开采深度H、工作面走向及倾向长度l和L、采厚m、煤层走向坐标方位角α的取值范围，并根据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》确定研究区域概率积分法预计参数；
S3由InSAR技术最终获取的地表沉陷形变区域的形变点的数量，并根据所获取地表形变点的密度决定最终参与采空区反演的形变点数量，使用经典粒子群算法的参数，初始化粒子群中每个个体的位置与速度信息；根据上述最终得到参与反演的地表点数量确定粒子群算法种群大小N；初始化粒子群中每个个体的速度与位置信息，其中个体的位置信息就是待反演的采空区位置参数，包括：开采深度、工作面走向及倾向长度、采厚、煤层走向坐标方位角。
S4将种群中每个个体的位置信息代入开采沉陷理论模型概率积分法中，计算得到地表各点位置处的三维形变场，包括竖向下沉W和南北、东西向水平移动UN、UE；结合InSAR三维形变分解模型将W、UN、UE合成地表各点雷达视线向形变场，记为：los；
所述InSAR三维形变分解模型由下式表征：
los＝Wcosθ-sinθ[UNcos(αh-3π/2)+UEsin(αh-3π/2)]
其中，los为由开采沉陷计算的三维形变合成的雷达视线向形变场；θ为雷达卫星入射角，αh为卫星航向角，W为竖向下沉，UN为南北向水平移动、UE为东西向水平移动；
S5将由InSAR技术解算的LOS视为真值，将概率积分模型计算地表三维形变后合成的los作为观测值，把LOS和los代入到适应度函数中，求得适应度函数值v，v越小则反演结果越优，记录每次迭代中同一个体的最优值和整个种群的最优值；
S6对粒子群算法的种群中每个个体的速度、位置进行更新；
粒子群每个个体新的速度公式如下：
Figure FDA0002800838770000021
粒子群每个个体位置依据速度进行转移的公式如下；
Figure FDA0002800838770000022
式中，
Figure FDA0002800838770000023
是第t次迭代第i个粒子的j维速度；
Figure FDA0002800838770000024
是第t次迭代第i个粒子的j维位置；pbest是个体历史最优解；gbest是种群中的最优状态；c1、c2是种群的学习因子；w是速度在更新时上次迭代时的速度占这次迭代的速度的权重；r1和r2是介于[0，1]之间的随机数；
S7寻求一个满足要求的种群最优值作为地下采空区真实参数，判断种群最优解gbest对应的适应度值是否满足迭代停止的预设阀值，不满足条件则重复步骤S4到步骤S6，直到满足迭代停止的预设阀值，并输出此时的种群最优解作为地下采空区位置参数。
3.根据权利要求2所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于：步骤S3中的粒子群中每个个体的位置的初始化使用混沌方式产生，其具体的生成方式为随机生成一个随机数矩阵，确定一个种群第一个个体的初始化位置，种群的其余个体利用前一个个体的随机数矩阵和Logistic映射产生初始位置时的随机数矩阵，依次执行下去直到产生所有个体的初始位置；
计算公式如下：
Figure FDA0002800838770000031
式中，up，low是被测区的开采深度H、工作面走向l、倾向长度L、采厚m、煤层走向坐标方位角α的取值范围的上下边界；rand(k)是第k个个体产生初始位置时对应的一个1×5的介于[0，1]之间的随机数矩阵；pos(k)时第k个个体的初始化位置；z是一个介于(2，4]之间的常数。
4.根据权利要求2所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于：步骤S4中的概率积分法计算竖向下沉W和南北、东西向水平移动UN、UE的计算方法为：
W＝W0(x)W0(y)/W0
Figure FDA0002800838770000032
Figure FDA0002800838770000033
Figure FDA0002800838770000034
Figure FDA0002800838770000035
其中：
W0＝mq cosα
Figure FDA0002800838770000036
Figure FDA0002800838770000037
投影长度及主要影响半径的计算方法为：
Figure FDA0002800838770000038
l＝D1-s1-s3
Figure FDA0002800838770000041
其中，m是煤层厚度；q是下沉系数；α是煤层倾角；x、y是地表任意点的矿区坐标系坐标；r、r1、r2是走向、倾向上山、倾向下山主要影响半径；H是煤层采深；tanβ是主要影响角正切；l、L分别为煤层的走向和倾向投影到水平面上的计算长度；D1，D2是工作面走向及倾向长度；s1，s2，s3，s4是左上右下拐点偏移距；θ是开采影响传播角；b是水平移动系数；erf是误差函数；
Figure FDA0002800838770000042
为煤层走向方位角。
5.根据权利要求2所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于：步骤S4所述的获取地表形变点的密度决定最终参与采空区反演的形变点数量，InSAR获取的观测点的数据量和密度很大，最终用于反演数据点的密度和传统测量的观测站设置的密度尽量保持一致，当InSAR获取的数据量密度不能达到下表要求时，直接按照InSAR获取的数据量密度进行反演；最终参与采空区位置反演的地表形变点密度见表1；
Figure FDA0002800838770000043
6.根据权利要求2所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于步骤S5所述的适应度函数v的表达式为：
Figure FDA0002800838770000044
式中，los(k)和LOS(k)分别为地表第k个观测点反演三维形变后由概率积分法模型合成的视线向形变和InSAR观测得到的时序视线向形变；n为最终参与反演的地表点个数。
7.根据权利要求2所述的融合InSAR和PSO反演地下采空区位置参数的方法，其特征在于：所述步骤S6中的w为动态权重，在迭代开始到结束的过程中w是遵循线性变化的其变化的规则为：
w(t)＝wmax-(wmax-wmin)·t/T
其中，wmin和wmax是动态权重的最小值和最大值；w(t)是第t次迭代过程中的权重取值；T是总迭代数。
CN202011351515.8A 2020-11-26 2020-11-26 融合InSAR和PSO反演地下采空区位置参数的方法 Active CN112505699B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011351515.8A CN112505699B (zh) 2020-11-26 2020-11-26 融合InSAR和PSO反演地下采空区位置参数的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011351515.8A CN112505699B (zh) 2020-11-26 2020-11-26 融合InSAR和PSO反演地下采空区位置参数的方法

Publications (2)

Publication Number Publication Date
CN112505699A true CN112505699A (zh) 2021-03-16
CN112505699B CN112505699B (zh) 2022-08-23

Family

ID=74966471

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011351515.8A Active CN112505699B (zh) 2020-11-26 2020-11-26 融合InSAR和PSO反演地下采空区位置参数的方法

Country Status (1)

Country Link
CN (1) CN112505699B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113281748A (zh) * 2021-05-24 2021-08-20 西南石油大学 一种地表形变监测方法
CN114089335A (zh) * 2021-11-16 2022-02-25 安徽理工大学 一种基于单轨道InSAR的山区开采沉陷三维变形提取方法
CN114325693A (zh) * 2021-12-19 2022-04-12 南京市测绘勘察研究院股份有限公司 一种基于InSAR时序形变结果的采空区中心形变预测方法
CN116520451A (zh) * 2023-04-24 2023-08-01 四川阳光上元能源技术有限公司 一种基于量子探测技术的地下空区探测方法
CN116976751A (zh) * 2023-08-28 2023-10-31 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) 一种基于时空数据融合的煤矿采空区勘察系统

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20100061087A (ko) * 2008-11-28 2010-06-07 영남대학교 산학협력단 Ｐｓｏ를 이용한 스텝주파수 레이더 영상의 요동 보상 방법
CN101770038A (zh) * 2010-01-22 2010-07-07 中国科学院武汉岩土力学研究所 矿山微震源智能定位方法
CN103091675A (zh) * 2013-01-11 2013-05-08 中南大学 一种基于InSAR技术的矿区开采监测方法
KR101397934B1 (ko) * 2014-03-05 2014-05-27 국방과학연구소 펄스 신호를 이용한 입자 기반의 표적 위치 추정 방법
CN103885046A (zh) * 2012-12-20 2014-06-25 河南省电力勘测设计院 基于GPS的InSAR大气延迟改正方法
CN106226764A (zh) * 2016-07-29 2016-12-14 安徽理工大学 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法
CN107271998A (zh) * 2017-07-01 2017-10-20 东华理工大学 一种集成D‑InSAR和GIS技术的地下非法采矿识别方法及系统
CN110991048A (zh) * 2019-12-04 2020-04-10 中国矿业大学 一种关闭井工矿地表沉陷预测方法
CN111076704A (zh) * 2019-12-23 2020-04-28 煤炭科学技术研究院有限公司 一种利用insar精确解算采煤沉陷区地表下沉量的方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20100061087A (ko) * 2008-11-28 2010-06-07 영남대학교 산학협력단 Ｐｓｏ를 이용한 스텝주파수 레이더 영상의 요동 보상 방법
CN101770038A (zh) * 2010-01-22 2010-07-07 中国科学院武汉岩土力学研究所 矿山微震源智能定位方法
CN103885046A (zh) * 2012-12-20 2014-06-25 河南省电力勘测设计院 基于GPS的InSAR大气延迟改正方法
CN103091675A (zh) * 2013-01-11 2013-05-08 中南大学 一种基于InSAR技术的矿区开采监测方法
KR101397934B1 (ko) * 2014-03-05 2014-05-27 국방과학연구소 펄스 신호를 이용한 입자 기반의 표적 위치 추정 방법
CN106226764A (zh) * 2016-07-29 2016-12-14 安徽理工大学 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法
CN107271998A (zh) * 2017-07-01 2017-10-20 东华理工大学 一种集成D‑InSAR和GIS技术的地下非法采矿识别方法及系统
CN110991048A (zh) * 2019-12-04 2020-04-10 中国矿业大学 一种关闭井工矿地表沉陷预测方法
CN111076704A (zh) * 2019-12-23 2020-04-28 煤炭科学技术研究院有限公司 一种利用insar精确解算采煤沉陷区地表下沉量的方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
范洪冬 等: "基于多轨道SAR的老采空区地表三维形变监测", 《采矿与安全工程学报》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113281748A (zh) * 2021-05-24 2021-08-20 西南石油大学 一种地表形变监测方法
CN113281748B (zh) * 2021-05-24 2022-02-11 西南石油大学 一种地表形变监测方法
CN114089335A (zh) * 2021-11-16 2022-02-25 安徽理工大学 一种基于单轨道InSAR的山区开采沉陷三维变形提取方法
CN114325693A (zh) * 2021-12-19 2022-04-12 南京市测绘勘察研究院股份有限公司 一种基于InSAR时序形变结果的采空区中心形变预测方法
CN116520451A (zh) * 2023-04-24 2023-08-01 四川阳光上元能源技术有限公司 一种基于量子探测技术的地下空区探测方法
CN116520451B (zh) * 2023-04-24 2024-03-12 四川阳光上元能源技术有限公司 一种基于量子探测技术的地下空区探测方法
CN116976751A (zh) * 2023-08-28 2023-10-31 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) 一种基于时空数据融合的煤矿采空区勘察系统
CN116976751B (zh) * 2023-08-28 2024-02-23 山东省鲁南地质工程勘察院(山东省地质矿产勘查开发局第二地质大队) 一种基于时空数据融合的煤矿采空区勘察系统

Also Published As

Publication number Publication date
CN112505699B (zh) 2022-08-23

Similar Documents

Publication Publication Date Title
CN112505699B (zh) 融合InSAR和PSO反演地下采空区位置参数的方法
CN103091675B (zh) 一种基于InSAR技术的矿区开采监测方法
Jiang et al. Observe the temporal evolution of deep tunnel's 3D deformation by 3D laser scanning in the Jinchuan No. 2 Mine
CN105806303B (zh) 融合D-InSAR和模矢法求取概率积分参数的方法
CN105678399A (zh) 一种区域矿产资源量估算分析方法及其系统
CN111076704A (zh) 一种利用insar精确解算采煤沉陷区地表下沉量的方法
CN115115180A (zh) 一种基于多参量分析的矿区地表塌陷风险识别与预测方法
Li et al. Study of probability integration method parameter inversion by the genetic algorithm
CN111551932A (zh) 准确获取开采影响边界及确定建筑损害等级的方法
Wang et al. Automatic identification of rock discontinuity and stability analysis of tunnel rock blocks using terrestrial laser scanning
Wang et al. D-InSAR monitoring method of mining subsidence based on Boltzmann and its application in building mining damage assessment
CN113091598A (zh) 一种InSAR划定采空区建筑场地稳定性等级范围的方法
Wang et al. Research on 3D laser scanning monitoring method for mining subsidence based on the auxiliary for probability integral method
An et al. Ground subsidence monitoring in based on UAV-LiDAR technology: a case study of a mine in the Ordos, China
Gao et al. The development mechanism and control technology visualization of the vault cracks in the ancient underground cavern of Longyou
Wang et al. Deriving mining-induced 3-D deformations at any moment and assessing building damage by integrating single InSAR interferogram and gompertz probability integral model (SII-GPIM)
CN106886584B (zh) 基于城市多种地理数据基础上的地下空间开发利用现状估算方法
CN108776854A (zh) 大型露天矿山边坡稳定性等精度评价方法
Pan et al. Assessment method of slope excavation quality based on point cloud data
Wei et al. Fusing minimal unit probability integration method and optimized quantum annealing for spatial location of coal goafs
Meier Geological characterisation of an underground research facility in the Bedretto tunnel
Ba et al. Development status of digital detection technology for unfavorable geological structures in deep tunnels
CN114608527A (zh) 一种基于煤矿开采地表沉陷特征的D-InSAR和UAV摄影测量融合方法
Tan et al. Prediction method for the deformation of deep foundation pit based on neural network algorithm optimized by particle swarm
Wickramathilaka et al. Calculation of road traffic noise, development of data, and spatial interpolations for traffic noise visualization in three-dimensional space

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