CN111650579A - 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 - Google Patents

一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 Download PDF

Info

Publication number
CN111650579A
CN111650579A CN202010533034.2A CN202010533034A CN111650579A CN 111650579 A CN111650579 A CN 111650579A CN 202010533034 A CN202010533034 A CN 202010533034A CN 111650579 A CN111650579 A CN 111650579A
Authority
CN
China
Prior art keywords
deformation
east
mining area
west
radar
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
CN202010533034.2A
Other languages
English (en)
Other versions
CN111650579B (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.)
Central South University
Original Assignee
Central South University
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 Central South University filed Critical Central South University
Priority to CN202010533034.2A priority Critical patent/CN111650579B/zh
Publication of CN111650579A publication Critical patent/CN111650579A/zh
Application granted granted Critical
Publication of CN111650579B publication Critical patent/CN111650579B/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
    • 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

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)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质，通过获得待监测矿区地表沿多个SAR观测几何的雷达视线向形变；联合不同SAR观测几何监测得到的形变值以及InSAR成像几何关系模型，在忽略南北向形变对雷达视线向形变观测量贡献的前提下，解算矿区地表垂直和东西向的二维形变分量；利用矿区地表水平移动与垂直沉降梯度之间的线性比例关系函数，对垂直向/东西向形变的视线向形变贡献量和岩移参数进行迭代优化，进而解算得到矿区地表三维形变。本方法实现了岩层移动参数的自适应获取和基于InSAR的矿山地表三维形变全自动化估计，克服了现有方法中岩层移动参数难以人工收集并准确获取的局限。

Description

一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装 置及介质
技术领域
本发明属于InSAR领域，特别涉及一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质。
背景技术
对矿区地表进行全盆地的三维形变监测，对于理解矿区形变机理和进行灾害防治与评估具有重要意义。合成孔径雷达干涉测量(Interferometric Synthetic ApertureRadar，InSAR)技术由于其无接触、高精度、高空间分辨率以及不受云雨天气影响等优势在矿区地表形变监测中发挥着越来越重要的作用。但是，InSAR技术监测得到的矿区地表形变是地表真实三维形变在雷达视线向的一维投影。这大大限制了其在矿区地表形变监测中的作用和监测精度。目前利用InSAR技术获取矿区地表三维形变的方法主要分为两类：1)利用成像几何结构具有显著差异的多轨道SAR数据联合监测法；2)基于矿区地表形变先验模型的单轨InSAR监测技术。
由于目前SAR卫星大多是沿极轨飞行，因此第一类方法对南北向形变的敏感性不高。此外，由于SAR卫星数量有限，且矿区地表形变具有高度非线性的特征，这也导致第一类方法的数据可行性和监测精度难以保证，使得该类方法在矿区地表三维形变监测中的可行性不高。鉴于矿区地表水平位移与垂直沉降梯度之间存在线性比例关系，之前提出“一种利用单个InSAR干涉对获取矿区地表三维形变场的方法”(专利号：CN201210440875)。该方法利用上述线性比例关系，从一个InSAR干涉对中获取了矿区地表三维形变场，从而大大降低了传统InSAR技术监测矿区地表三维形变中对监测数据的苛刻要求，拓展了InSAR技术在矿区地表形变监测中的应用前景。但是，该方法在进行三维形变估计过程中，需要已知待监测矿区的岩层移动参数(以下简称岩移参数)。这些岩移参数在矿山一般是保密的，通常很难获取；而且由于矿山环境的复杂性，所获得的岩移参数往往与真实情况相比存在一定的误差。这些问题阻碍了该方法的实际应用。
发明内容
本发明的目的是提供一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质，该方法基于矿区地表形变先验模型，利用多平台/多轨道SAR卫星监测得到的雷达视线向形变建立矿区地表真实三维形变与多SAR几何观测结果之间的投影转换关系模型，并利用迭代求解的方法对岩移参数进行迭代优化，从而实现岩移参数的自适应获取和矿区地表三维形变的全自动化估计，克服现有方法中矿山岩移参数难以人工收集并准确获取的问题，实现真正意义上的基于InSAR的矿区地表三维形变监测。
本发明提供的技术方案如下：
一方面，一种岩移参数自适应获取的InSAR矿区三维形变估计方法，包括如下步骤：
S1：获取矿区地表沿不同雷达成像几何的视线向形变；
S2：将S1中得到的视线向形变进行地理编码转换到DEM的参考地理坐标系统下，得到参考坐标系统统一的多个雷达成像几何的视线向形变；
S3：根据SAR卫星成像几何关系，将坐标转换后的多个雷达成像几何的视线向形变在忽略南北向形变贡献量后进行解算，得到垂直向、东西向形变分量的初始解；
S4：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
求解岩移参数的初始值时，只需要利用矿区地表东西向水平移动与东西向的沉降梯度之间的线性比例关系函数；
在忽略矿区地表南北向形变对雷达视线向形变贡献量的假设条件下，利用雷达成像几何不同的多平台/多轨道SAR数据解算得到地表垂直和东西向形变，并以此为初值，利用稳健估计对岩移参数进行初步估计，并应用于接下来的迭代优化计算。
S5：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由坐标转换后的多个雷达成像几何的视线向形变到矿区地表三维形变的解算方程，将获得的岩移参数代入解算方程，获得矿区地表三维形变；
S5中建立解算方程时，需要利用矿区地表东西向、垂直向水平移动与东西向、垂直向的沉降梯度之间的线性比例关系函数；
S6：将S5中解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
S7：结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复S5～S7进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1，进入S8；
S8：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
进一步地，采用基于M-估计的稳健估计方法，并通过迭代重加权的方式对岩移参数的初始值c0进行迭代求解，即：
Figure BDA0002536071790000021
其中，dE表示东西向形变分量初始解的观测矩阵，
Figure BDA0002536071790000031
为第i次迭代求解过程中得到的重加权矩阵，重加权矩阵中的每个元素表示参与对c0进行解算的每个差分InSAR监测点的权重，GE为东西向的沉降梯度矩阵，
Figure BDA0002536071790000032
为GE的转置；稳健估计中
Figure BDA0002536071790000033
的阈值设置为10-5，
Figure BDA0002536071790000034
和
Figure BDA0002536071790000035
分别表示第i+1次和第i次迭代求解得到的岩移参数的初始值。
在岩移参数的初始解的稳健估计中，观测矩阵中存在误差的监测点处的权重将接近于0，从而削除其对解算结果的影响；
由于矿区差分InSAR结果易受解缠误差的影响，从而在接下来的参数估计中造成较大的偏差，稳健估计类型较多，这里设定为M-估计(一种已有的用于参数稳健估计的经典算子)；
进一步地，所述将S5中解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除是指按照以下公式计算获得：
Figure BDA0002536071790000036
其中，
Figure BDA0002536071790000037
是指通过第i次迭代后获得的垂直和东西向形变的视线向形变贡献量，即忽略了南北向形变贡献量的雷达视线向形变量，
Figure BDA0002536071790000038
表示第i次迭代解算得到的南北向形变分量，λ表示雷达入射角，α为卫星飞行方位角。
由于垂直和东西向形变分量的初始值是在忽略南北向形变贡献量dN的情况下解算得到的，此处得到了南北向形变分量的解算值
Figure BDA0002536071790000039
进一步地，S3中在忽略南北向形变贡献量的情况下，由坐标转换后的k个雷达成像几何的视线向形变组成的观测矩阵
Figure BDA00025360717900000310
解算得到垂直向和东西向形变分量的初始解dU和dE的计算公式为：
Figure BDA00025360717900000311
其中，
Figure BDA00025360717900000312
表示第j个SAR成像几何的观测值，λj表示第j个SAR成像几何的雷达入射角，αj为第j个SAR成像几何的卫星飞行方位角，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量。
进一步地，S3中所述SAR卫星成像几何关系式为：
Figure BDA0002536071790000041
其中，
Figure BDA0002536071790000042
表示第j个SAR成像几何的观测值，λj表示第j个SAR成像几何的雷达入射角，αj为第j个SAR成像几何的卫星飞行方位角，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量；dU和dE分别表示垂直向和东西向形变分量，dN表示南北向形变贡献量。
进一步地，S4中所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
Figure BDA0002536071790000043
式中，p＝1,…,M，q＝1,…,N，M和N分别表示InSAR监测结果中影像的行列数；dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，其中，GE(p,q)＝[dU(p+1,q)-dU(p,q)]RE，GN(p,q)＝[dU(p,q)-dU(p,q+1)]RN，dU(p,q)为地表(p,q)点处的垂直沉降解算值；RE和RN为地理编码后的矿区地表形变场沿东西和南北向的像元分辨率。
进一步地，S5中所述由坐标转换后的雷达视线向形变dLOS到矿区三维形变的解算方程为：
Figure BDA0002536071790000044
其中，
Figure BDA0002536071790000045
dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，Bj为第j轨道InSAR监测结果的系数矩阵，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量。
另一方面，一种岩移参数自适应获取的InSAR矿区三维形变估计装置，包括：
数据获取单元：用于获取矿区地表沿不同雷达成像几何的视线向形变；
坐标转换单元：用于将数据获取单元得到的视线向形变进行地理编码转换到DEM的参考地理坐标系统下，得到参考坐标系统统一的多个雷达成像几何的视线向形变；
形变分量初始解算单元：根据SAR卫星成像几何关系，将坐标转换后的多个雷达成像几何的视线向形变在忽略南北向形变贡献量后进行解算，得到垂直向、东西向形变分量的初始解；
岩移参数初始解算单元：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
地表三维形变解算单元：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由坐标转换后的多个雷达成像几何的视线向形变到矿区地表三维形变的解算方程，将获得的岩移参数代入解算方程，获得矿区地表三维形变；
迭代计算单元：地表三维形变解算单元解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
并结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复调用地表三维形变解算单元和迭代计算单元进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1时，调用矿区三维形变估计单元；
矿区三维形变估计单元：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
进一步地，所述岩移参数初始解算单元采用基于M-估计的稳健估计方法，并通过迭代重加权的方式对岩移参数的初始值c0进行迭代求解，即：
Figure BDA0002536071790000051
其中，dE表示东西向形变分量初始解的观测矩阵，
Figure BDA0002536071790000052
为第i次迭代求解过程中得到的重加权矩阵，重加权矩阵中的每个元素表示参与对c0进行解算的每个差分InSAR监测点的权重，GE为东西向的沉降梯度矩阵，
Figure BDA0002536071790000053
为GE的转置；稳健估计中
Figure BDA0002536071790000054
的阈值设置为10-5，
Figure BDA0002536071790000055
和
Figure BDA0002536071790000056
分别表示第i+1次和第i次迭代求解得到的岩移参数的初始值。
再一方面，一种计算机存储介质，包括计算机程序，所述计算机程序被处理器执行时实现所述的一种岩移参数自适应获取的InSAR矿区三维形变估计方法。
有益效果
本发明提供了一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质，利用多平台/多轨道DInSAR技术从多观测几何InSAR干涉对中获得矿区地表沿不同雷达视线向的形变场；然后，联合不同SAR观测几何监测得到的形变值以及InSAR成像几何关系模型，在忽略南北向形变对雷达视线向形变观测量贡献的前提下，解算得到矿区地表沿垂直向和东西向的形变分量；之后，利用矿区地表水平移动与垂直沉降梯度之间的线性比例关系函数模型以及上述解算得到的矿区地表垂直和东西向形变分量，对垂直向/东西向形变的雷达视线向形变贡献量和岩层移动参数(简称岩移参数)进行迭代优化，进而求解得到矿区地表三维形变。该方法计算简单，实现方便；
本发明的技术效果主要体现在以下几点：
第一、相对于传统的利用成像几何具有显著差异的多几何SAR数据联合监测技术，本方法结合矿区地表形变先验模型，在实现矿区三维形变监测的同时，显著提高了多轨道InSAR解算地表三维形变的监测精度，丰富了InSAR技术在矿区地表三维形变监测领域的功效；
第二、相对于已有的基于单个InSAR干涉对获取矿区地表三维形变的方法，本方法结合多平台/多轨道观测数据，在增加观测值的同时，可显著增强模型参数解算结果的鲁棒性，使得监测结果更加稳定、可靠；
第三、相对于已有的基于单个InSAR干涉对获取矿区地表三维形变的方法，本方法结合矿区地表形变先验模型，对多几何InSAR数据进行联合解算，通过对地表三维形变和岩移参数的迭代求解，克服了已有方法对矿山岩层移动参数的依赖性，从而真正意义上实现了基于InSAR的矿山地表三维形变全自动获取，大大拓展了InSAR技术在矿区地表三维形变监测预计领域的应用前景。
附图说明
图1本发明实例所述方法的流程示意图；
图2模拟的矿区地表三维形变示意图；
图3基于本发明实例所述方法计算得到的矿区地表三维形变示意图。
具体实施方式
下面将结合实施例对本发明做进一步的说明。
如图1所示，本发明提供了一种岩移参数自适应获取的InSAR矿区三维形变估计方法，包括如下步骤：
S1：利用多平台/多轨道合成孔径雷达差分干涉测量(DifferentialInterferometric Synthetic Aperture Radar，DInSAR)技术监测得到矿区地表沿不同雷达成像几何的视线向形变dlos；其中，
Figure BDA0002536071790000071
k为多平台/多轨道DInSAR观测值的数量；
S2：以进行DInSAR技术数据处理时所用的外部数字高程模型(Digital ElevationModel，DEM)坐标系统为参考，将步骤S1中得到的dlos进行地理编码，得到坐标系统统一后的雷达视线向形变dLOS；
所述地理编码，是指将不同成像几何InSAR数据获取的各形变场的坐标系转换到DEM的参考地理坐标系统，从而实现不同SAR成像几何形变场数据的坐标系统统一；
S3：根据SAR卫星成像几何关系，将坐标系统统一后的雷达视线向形变dLOS在忽略南北向形变贡献量dN的情况下按照下述公式解算得到垂直向和东西向形变分量的初始解dU和dE：
Figure BDA0002536071790000072
所述SAR卫星成像几何关系式为：
Figure BDA0002536071790000073
其中λi表示雷达入射角，αi为卫星飞行方位角；
S4：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
需要说明的是，由于矿区差分InSAR结果易受解缠误差的影响，从而在接下来的参数估计中造成较大的偏差，因此，本实例采用一种稳健估计方法对初代岩移参数进行求解，即：
Figure BDA0002536071790000074
其中，
Figure BDA0002536071790000075
为第i次迭代求解过程中得到的重加权矩阵，重加权矩阵中的每个元素表示参与对c0进行解算的每个差分InSAR监测点的权重，GE为东西向的沉降梯度，
Figure BDA0002536071790000076
为GE的转置；稳健估计中
Figure BDA0002536071790000077
的阈值设置为10-5，
Figure BDA0002536071790000078
和
Figure BDA0002536071790000079
分别表示第i+1次和第i次迭代求解得到的岩移参数的初始解；
在岩移参数的初始解的稳健估计中，观测矩阵中存在误差的监测点处的权重将接近于0，从而削除其对解算结果的影响；
由于矿区差分InSAR结果易受解缠误差的影响，从而在接下来的参数估计中造成较大的偏差，稳健估计类型较多，这里设定为M-估计(一种已有的用于参数稳健估计的经典算子)；
所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数为：
Figure BDA0002536071790000081
式中，p＝1,…,M，q＝1,…,N，M和N分别表示InSAR监测结果影像的行列数；dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，GE(p,q)＝[dU(p+1,q)-dU(p,q)]RE和GN(p,q)＝[dU(p,q)-dU(p,q+1)]/RN，其中，dU(p,q)为地表(p,q)点处的垂直沉降解算值；RE和RN为地理编码后的矿区地表形变场沿东西和南北向的像元分辨率。
S5：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由dLOS到矿区三维形变的解算方程，将计算得到的c0带入该方程，获得矿区地表初始三维形变场；
所述由坐标转换后的雷达视线向形变dLOS到矿区三维形变的解算方程为：
Figure BDA0002536071790000082
其中，
Figure BDA0002536071790000083
dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，Bj为第j轨道InSAR监测结果的系数矩阵，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量。
S6：将S5中解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
S7：结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复S5～S7进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1，进入S8；
S8：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
为了更加清楚的说明本发明及其优势，下述将结合部分图来进一步解释本发明提供的所述方法。本节以某矿区真实地表形变的岩移参数(c＝75)为参考，模拟了在该岩移参数条件下的矿区地下开采导致的地表三维形变(如图2所示)，并以ALOS PALSAR卫星升轨(λA＝38°，αA＝-10°)和降轨(λD＝38°，αD＝190°)SAR数据的成像几何参数为例，根据几何投影转换关系模型，将模拟的地表三维形变分别投影转换到升/降轨SAR数据的雷达视线向，组成dLOS。之后，利用本发明所述方法对模拟的升/降轨数据形变监测场进行处理，反演得到矿区地表三维形变场(如图3所示)。
从图2和图3中可以看出，本发明所述方法在无需外部输入岩移参数的情况下，依靠沿雷达视线向和飞行方位向的监测数据直接反演出了矿区地表三维形变。为了更直观的反应本发明所述方法的可靠性，我们计算了沿垂直、东西和南北方向的三维形变估值与模拟真值之间的均方根误差，分别为0.04cm、0.056cm和0.13cm。从而验证了本发明所述方法的可行性和可靠性。
基于上述方法，本发明实施例还提供一种岩移参数自适应获取的InSAR矿区三维形变估计装置，包括：
数据获取单元：用于获取矿区地表沿不同雷达成像几何的视线向形变；
坐标转换单元：用于将数据获取单元得到的视线向形变进行地理编码转换到DEM的参考地理坐标系统下，得到参考坐标系统统一的多个雷达成像几何的视线向形变；
形变分量初始解算单元：根据SAR卫星成像几何关系，将坐标转换后的多个雷达成像几何的视线向形变在忽略南北向形变贡献量后进行解算，得到垂直向、东西向形变分量的初始解；
岩移参数初始解算单元：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
地表三维形变解算单元：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由坐标转换后的多个雷达成像几何的视线向形变到矿区地表三维形变的解算方程，将获得的岩移参数代入解算方程，获得矿区地表三维形变；
迭代计算单元：地表三维形变解算单元解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
并结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复调用地表三维形变解算单元和迭代计算单元进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1时，调用矿区三维形变估计单元；
矿区三维形变估计单元：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中，可以采用硬件或软件的形式来实现。
本发明实施例还提供一种计算机存储介质，包括计算机程序，所述计算机程序被处理器执行时实现所述的一种岩移参数自适应获取的InSAR矿区三维形变估计方法；其有益效果参见方法部分的有益效果，在此不再赘述。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详尽的说明，所属领域的普通技术人员应当理解，上述实施例仅仅是对本发明的示意性实现方式的解释，实施例中的细节并不构成对本发明范围的限制，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案的等效变换、简单替换等显而易见的改变，均落在本发明保护范围之内。

Claims (10)

1.一种岩移参数自适应获取的InSAR矿区三维形变估计方法，其特征在于：包括如下步骤：
S1：获取矿区地表沿不同雷达成像几何的视线向形变；
S2：将S1中得到的视线向形变进行地理编码转换到DEM的参考地理坐标系统下，得到参考坐标系统统一的多个雷达成像几何的视线向形变；
S3：根据SAR卫星成像几何关系，将坐标转换后的多个雷达成像几何的视线向形变在忽略南北向形变贡献量后进行解算，得到垂直向、东西向形变分量的初始解；
S4：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
S5：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由坐标转换后的多个雷达成像几何的视线向形变到矿区地表三维形变的解算方程，将获得的岩移参数代入解算方程，获得矿区地表三维形变；
S6：将S5中解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
S7：结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复S5～S7进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1，进入S8；
S8：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
2.根据权利要求1所述的方法，其特征在于，采用基于M-估计的稳健估计方法，并通过迭代重加权的方式对岩移参数的初始值c0进行迭代求解，即：
Figure FDA0002536071780000011
其中，dE表示东西向形变分量初始解的观测矩阵，
Figure FDA0002536071780000012
为第i次迭代求解过程中得到的重加权矩阵，重加权矩阵中的每个元素表示参与对c0进行解算的每个差分InSAR监测点的权重，GE为东西向的沉降梯度矩阵，
Figure FDA0002536071780000013
为GE的转置；稳健估计中
Figure FDA0002536071780000014
的阈值设置为10-5，
Figure FDA0002536071780000015
和
Figure FDA0002536071780000016
分别表示第i+1次和第i次迭代求解得到的岩移参数的初始值。
3.根据权利要求1所述的方法，其特征在于，所述将S5中解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除是指按照以下公式计算获得：
Figure FDA0002536071780000017
其中，
Figure FDA0002536071780000021
是指通过第i次迭代后获得的垂直和东西向形变的视线向形变贡献量，即忽略了南北向形变贡献量的雷达视线向形变量，
Figure FDA0002536071780000022
表示第i次迭代解算得到的南北向形变分量，λ表示雷达入射角，α为卫星飞行方位角。
4.根据权利要求1所述方法，其特征在于：S3中在忽略南北向形变贡献量的情况下，由坐标转换后的k个雷达成像几何的视线向形变组成的观测矩阵
Figure FDA0002536071780000023
解算得到垂直向和东西向形变分量的初始解dU和dE的计算公式为：
Figure FDA0002536071780000024
其中，
Figure FDA0002536071780000025
表示第j个SAR成像几何的观测值，λj表示第j个SAR成像几何的雷达入射角，αj为第j个SAR成像几何的卫星飞行方位角，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量。
5.根据权利要求1所述的方法，其特征在于，S3中所述SAR卫星成像几何关系式为：
Figure FDA0002536071780000026
其中，
Figure FDA0002536071780000027
表示第j个SAR成像几何的观测值，λj表示第j个SAR成像几何的雷达入射角，αj为第j个SAR成像几何的卫星飞行方位角，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量；dU和dE分别表示垂直向和东西向形变分量，dN表示南北向形变贡献量。
6.根据权利要求1所述的方法，其特征在于，S4中所述矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数模型为：
Figure FDA0002536071780000028
式中，p＝1,…,M，q＝1,…,N，M和N分别表示InSAR监测结果中影像的行列数；dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，其中，GE(p,q)＝[dU(p+1,q)-dU(p,q)]/RE，GN(p,q)＝[dU(p,q)-dU(p,q+1)]/RN，dU(p,q)为地表(p,q)点处的垂直沉降解算值；RE和RN为地理编码后的矿区地表形变场沿东西和南北向的像元分辨率。
7.根据权利要求1所述方法，其特征在于：S5中所述由坐标转换后的雷达视线向形变dLOS到矿区三维形变的解算方程为：
Figure FDA0002536071780000031
其中，
Figure FDA0002536071780000032
dE(p,q)和dN(p,q)分别表示影像中(p,q)像素对应的矿区地表位置坐标处分别沿东西和南北向的形变分量，c为待计算的岩移参数；GE(p,q)和GN(p,q)表示解算得到的矿区地表(p,q)点处在东西向和南北向的沉降梯度，Bj为第j轨道InSAR监测结果的系数矩阵，j＝1,…,k，k表示参与解算的多个雷达成像几何的视线向形变观测值的数量。
8.一种岩移参数自适应获取的InSAR矿区三维形变估计装置，其特征在于：包括：
数据获取单元：用于获取矿区地表沿不同雷达成像几何的视线向形变；；
坐标转换单元：用于将数据获取单元得到的视线向形变进行地理编码转换到DEM的参考地理坐标系统下，得到参考坐标系统统一的多个雷达成像几何的视线向形变；
形变分量初始解算单元：根据SAR卫星成像几何关系，将坐标转换后的多个雷达成像几何的视线向形变在忽略南北向形变贡献量后进行解算，得到垂直向、东西向形变分量的初始解；
岩移参数初始解算单元：利用垂直向形变分量的初始解，计算地表沿东西向的沉降梯度，并结合东西向形变分量的初始解和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，获取岩移参数的初始值；
地表三维形变解算单元：联立SAR卫星成像几何关系式和地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，建立由坐标转换后的多个雷达成像几何的视线向形变到矿区地表三维形变的解算方程，将获得的岩移参数代入解算方程，获得矿区地表三维形变；
迭代计算单元：地表三维形变解算单元解算得到的南北向形变分量从坐标转换后的多个雷达成像几何的视线向形变中剔除后，进行重新解算，获取新的垂直向、东西向形变分量，并更新东西向的沉降梯度；
并结合新解算得到的东西向形变分量和矿区地表水平移动与对应方向的沉降梯度之间的线性比例关系函数，对岩移参数进行优化重解；重复调用地表三维形变解算单元和迭代计算单元进行迭代计算，直到岩移参数的当前值与前一次迭代得到的计算值之差小于0.1时，调用矿区三维形变估计单元；
矿区三维形变估计单元：利用岩移参数的当前值代入到由坐标转换后的多个雷达成像几何的视线向形变到矿区三维形变的解算方程中，实现岩移参数自适应获取的矿区三维形变估计。
9.根据权利要求8所述的装置，其特征在于，所述岩移参数初始解算单元采用基于M-估计的稳健估计方法，并通过迭代重加权的方式对岩移参数的初始值c0进行迭代求解，即：
Figure FDA0002536071780000041
其中，dE表示东西向形变分量初始解的观测矩阵，
Figure FDA0002536071780000042
为第i次迭代求解过程中得到的重加权矩阵，重加权矩阵中的每个元素表示参与对c0进行解算的每个差分InSAR监测点的权重，GE为东西向的沉降梯度矩阵，
Figure FDA0002536071780000043
为GE的转置；稳健估计中
Figure FDA0002536071780000044
的阈值设置为10-5，
Figure FDA0002536071780000045
和
Figure FDA0002536071780000046
分别表示第i+1次和第i次迭代求解得到的岩移参数的初始值。
10.一种计算机存储介质，包括计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7任一项所述的一种岩移参数自适应获取的InSAR矿区三维形变估计方法。
CN202010533034.2A 2020-06-12 2020-06-12 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质 Active CN111650579B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010533034.2A CN111650579B (zh) 2020-06-12 2020-06-12 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010533034.2A CN111650579B (zh) 2020-06-12 2020-06-12 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质

Publications (2)

Publication Number Publication Date
CN111650579A true CN111650579A (zh) 2020-09-11
CN111650579B CN111650579B (zh) 2022-09-30

Family

ID=72347665

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010533034.2A Active CN111650579B (zh) 2020-06-12 2020-06-12 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质

Country Status (1)

Country Link
CN (1) CN111650579B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114236541A (zh) * 2021-12-08 2022-03-25 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114777633A (zh) * 2022-04-02 2022-07-22 山西省煤炭地质勘查研究院有限公司 一种关闭矿区阶段性形变监测及分析方法
CN115343710A (zh) * 2022-07-11 2022-11-15 中国地质大学（武汉） 基于多平台时序InSAR的数据融合及冻融监测方法及装置
CN115951354A (zh) * 2023-02-14 2023-04-11 中国铁道科学研究院集团有限公司铁道建筑研究所 一种融合升降轨的D-InSAR形变监测方法
CN114777633B (zh) * 2022-04-02 2024-05-14 山西省煤炭地质勘查研究院有限公司 一种关闭矿区阶段性形变监测及分析方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040090360A1 (en) * 2002-10-24 2004-05-13 The Regents Of The University Of California Using dynamic interferometric synthetic aperature radar (InSAR) to image fast-moving surface waves
US20090284758A1 (en) * 2008-05-19 2009-11-19 Toyonaka Kenkyusho Co., Ltd. Displacement measuring method, displacement measuring apparatus and displacement measuring target
CN106526590A (zh) * 2016-11-04 2017-03-22 山东科技大学 一种融合多源sar影像工矿区三维地表形变监测及解算方法
CN106767380A (zh) * 2017-01-19 2017-05-31 中南大学 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法
CN107144213A (zh) * 2017-06-29 2017-09-08 中南大学 基于sar强度影像的矿区大量级三维时序形变估计方法及装置
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
CN110058234A (zh) * 2019-05-20 2019-07-26 太原理工大学 一种解算矿区地表沉降三维形变的方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040090360A1 (en) * 2002-10-24 2004-05-13 The Regents Of The University Of California Using dynamic interferometric synthetic aperature radar (InSAR) to image fast-moving surface waves
US20090284758A1 (en) * 2008-05-19 2009-11-19 Toyonaka Kenkyusho Co., Ltd. Displacement measuring method, displacement measuring apparatus and displacement measuring target
CN106526590A (zh) * 2016-11-04 2017-03-22 山东科技大学 一种融合多源sar影像工矿区三维地表形变监测及解算方法
CN106767380A (zh) * 2017-01-19 2017-05-31 中南大学 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法
CN107144213A (zh) * 2017-06-29 2017-09-08 中南大学 基于sar强度影像的矿区大量级三维时序形变估计方法及装置
CN109738892A (zh) * 2019-01-24 2019-05-10 中南大学 一种矿区地表高时空分辨率三维形变估计方法
CN110058234A (zh) * 2019-05-20 2019-07-26 太原理工大学 一种解算矿区地表沉降三维形变的方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
ZEFA YANG 等: "An InSAR-Based Temporal Probability Integral", 《IEEE JOURNAL OF SELECTED TOPICS IN APPLIED EARTH OBSERVATIONS AND REMOTE SENSING》 *
ZEFA YANG等: "Time-Series 3-D Mining-Induced Large", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
吴立新等: "基于D-InSAR的煤矿区开采沉陷遥感监测技术分析", 《地理与地理信息科学》 *
王永哲等: "融合多平台DInSAR数据解算拉奎拉地震三维同震形变场", 《武汉大学学报(信息科学版)》 *
胡俊等: "DInSAR监测地表三维形变的方法", 《工程勘察》 *
胡俊等: "融合升降轨SAR干涉相位和幅度信息揭示地表三维形变场的研究", 《中国科学:地球科学》 *
郭红: "基于INSAR的矿区地表三维形变解算方法应用研究", 《内蒙古煤炭经济》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114236541A (zh) * 2021-12-08 2022-03-25 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114236541B (zh) * 2021-12-08 2023-05-16 电子科技大学 基于Sentinel-1卫星SAR图像的大面积地表三维形变计算方法
CN114777633A (zh) * 2022-04-02 2022-07-22 山西省煤炭地质勘查研究院有限公司 一种关闭矿区阶段性形变监测及分析方法
CN114777633B (zh) * 2022-04-02 2024-05-14 山西省煤炭地质勘查研究院有限公司 一种关闭矿区阶段性形变监测及分析方法
CN115343710A (zh) * 2022-07-11 2022-11-15 中国地质大学（武汉） 基于多平台时序InSAR的数据融合及冻融监测方法及装置
CN115343710B (zh) * 2022-07-11 2024-04-16 中国地质大学（武汉） 基于多平台时序InSAR的数据融合及冻融监测方法及装置
CN115951354A (zh) * 2023-02-14 2023-04-11 中国铁道科学研究院集团有限公司铁道建筑研究所 一种融合升降轨的D-InSAR形变监测方法

Also Published As

Publication number Publication date
CN111650579B (zh) 2022-09-30

Similar Documents

Publication Publication Date Title
US10762645B2 (en) Stereo visual odometry method based on image gradient joint optimization
CN110738252B (zh) 空间自相关的机器学习卫星降水数据降尺度方法、系统
CN110058237B (zh) 面向高分辨率SAR影像的InSAR点云融合及三维形变监测方法
CN110058236B (zh) 一种面向三维地表形变估计的InSAR和GNSS定权方法
CN111650579B (zh) 一种岩移参数自适应获取的InSAR矿区三维形变估计方法、装置及介质
Li et al. 3-D shoreline extraction from IKONOS satellite imagery
CN108983232B (zh) 一种基于邻轨数据的InSAR二维地表形变监测方法
CN102607534B (zh) 基于运动恢复结构的卫星相对姿态测量方法
CN105677942A (zh) 一种重复轨道星载自然场景sar复图像数据快速仿真方法
CN109738892A (zh) 一种矿区地表高时空分辨率三维形变估计方法
CN103218780B (zh) 基于逆rd定位模型的无控星载sar图像正射校正方法
CN111724465B (zh) 基于平面约束优选虚拟控制点的卫星影像平差方法及装置
CN103454636B (zh) 基于多像素协方差矩阵的差分干涉相位估计方法
CN111077525B (zh) 融合sar与光学偏移量技术的地表三维形变计算方法及系统
CN109061641A (zh) 一种基于序贯平差的InSAR时序地表形变监测方法
CN111398959A (zh) 基于地表应力应变模型的InSAR时序地表形变监测方法
CN107144213A (zh) 基于sar强度影像的矿区大量级三维时序形变估计方法及装置
CN107316280A (zh) 离岛卫星影像rpc模型高精度几何定位方法
CN111666896A (zh) 一种基于线性融合模型的遥感影像时空融合方法
CN115469308A (zh) 多轨道InSAR震间形变速率场拼接方法、装置、设备及介质
CN115629384A (zh) 一种时序InSAR误差的改正方法及相关设备
CN104537614B (zh) 一种环境一号卫星ccd影像正射校正方法
CN115795402B (zh) 一种基于变分法的多源降水数据融合方法和系统
CN116068511B (zh) 一种基于深度学习的InSAR大尺度系统误差改正方法
CN109579796B (zh) 一种投影后影像的区域网平差方法

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