CN102927934B - 一种利用单个InSAR干涉对获取矿区地表三维形变场的方法 - Google Patents

Abstract

本发明公开了一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，通过利用InSAR技术获取矿区雷达视线向形变场，对视线向形变场中相干性低于解缠阈值的像素进行空间插值，得到空间连续的形变场；利用矿区工作面分布和主要影响角正切计算各像素点的主要影响半径；获取矿区水平移动系数之后，将矿区地表水平移动转换为下沉值的表达式，并按照雷达成像原理组成方程组；求解方程组，得出地表下沉值的解，然后根据下沉值计算出东西、南北方向倾斜值；最后利用东西、南北方向的倾斜值与水平移动的比例关系计算东西、南北方向的形变场；本发明突破了InSAR求解三维形变场对于数据的苛刻要求及监测费用高等制约，大大的拓宽了InSAR技术在矿区的应用空间。

Description

一种利用单个InSAR干涉对获取矿区地表三维形变场的方法

技术领域

本发明涉及一种利用单个InSAR干涉对获取矿区地表三维形变场的方法。

背景技术

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，简称InSAR)是90年代发展起来的一种新型地表形变监测技术，其基本原理就是通过对两幅或以上的合成孔径雷达(Synthetic Aperture Radar，简称SAR)卫星影像进行差分干涉处理，从相位差中提取厘米甚至毫米级的雷达视线方向形变值。相对于传统的大地测量技术而言，InSAR具有全天候、连续空间覆盖、高精度、低成本等优势。但是，传统的InSAR技术只能监测雷达视线方向的一维形变，由于其只获取了一个方向的形变场，因此不能满足矿区开采控制管理及沉降预计所需的地表三维形变的要求。为了利用InSAR技术获取矿区三维形变场，一些学者提出了GPS与InSAR融合、多平台InSAR资料融合、多孔径InSAR等技术。但是这些方法均存在局限和不足，多平台InSAR资料融合由于目前可用SAR卫星较少，且飞行方向大致相同，同时获得3个以上的同一地区不同平台干涉数据对非常困难。GPS与InSAR结合法要求GPS监测时间分辨率较高，对于建立连续GPS监测网的城市较为适用，对于矿区，由于布设条件以及费用等制约致使矿区很难得到高时间、空间分辨率的GPS监测资料，所以该方法对于矿区而言局限较大。

发明内容

本发明的目的在于提供一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，克服现有的矿区地表三维形变获取方法观测费用高、数据要求苛刻、观测精度不高等缺陷。

一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，包括以下步骤：

步骤1：利用InSAR技术获取待测矿区地理编码后的雷达视线向形变场；所述地理编码，是指将雷达影像坐标系转换到通用横轴墨卡托投影(Universal Transverse Mercatolprojection，简称UTM)坐标系；

步骤2：对视线向形变场中相干性低于解缠时的相干性阈值的像素进行空间插值，得到空间连续的形变场，解缠时的相干性阈值大于0.3；

步骤3：利用对矿区地表监测的水平移动数据、下沉值数据、采深及工作面分布，基于概率积分法模型求取该矿区的水平移动系数和主要影响角正切，之后利用各像素点对应的采深计算各像素点的主要影响半径，主要影响半径等于采深除以主要影响角正切，并按照雷达成像原理和矿区东西、南北方向的水平移动值等于水平移动系数、主要影响半径及该方向上该点倾斜值的乘积的关系，组成矿区地表东西、南北方向的任意点的水平移动值与矿区地表东西、南北方向上对应点的下沉值之间的方程组；

所述雷达成像原理由下式表征：

LOS＝W cosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]；

其中，LOS为雷达视线向的形变值，利用步骤1中的InSAR技术获得，W为下沉值，U_N为南北方向水平移动，U_E为东西方向水平移动；θ为雷达卫星入射角，α_h为卫星飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

步骤4：求步骤3的方程组的解，得出地表各像素点的下沉值；再求出各像素点在东西、南北方向上下沉值的梯度大小；各像素点在东西、南北方向上下沉值的梯度大小即为各像素点在东西、南北水平方向上的倾斜值，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

步骤5：利用各像素点东西、南北方向的水平移动值等于水平移动系数、主要影响半径及该点在东西、南北方向上的倾斜值的乘积关系，计算东西、南北方向的水平移动值，最终得出各像素点的三维形变值，包括各像素点的东西、南北方向的水平移动值及下沉值，即得到整个矿区的三维形变场。

所述步骤3方程组中的矿区地表水平移动值与下沉值之间的转换关系如下：

$\left\{\begin{array}{c}{U}_E(i,j)=b\·r(i,j)\·T_E(i,j)=b\·r(i,j)\·[W(i,j+1)-W(i,j)]/\mathrm{\ΔE}\\ U_N(i,j)=b\·r(i,j)\·T_N(i,j)=b\·r(i,j)\·[W(i+1,j)-W(i,j)]/\mathrm{\ΔN}\end{array}\right.,(i,j=\mathrm{1,2},\·\·\·,n-1)$

(i，j)为地理编码后的矿区地表区域内的任一像素点，U_E(i，j)、U_N(i，j)分别为矿区地表东西、南北方向任意像素点的水平移动值，b为矿区水平移动系数，T_E(i,j)、T_N(i,j)分别为东西、南北方向上该像素点的倾斜值，r(i，j)为该像素点的主要影响半径；△E、△N分别为地理编码后像素在东西、南北方向的分辨率，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得。

所述水平移动系数和主要影响角正切在不能通过本矿区实测数据求取时，选用临近地质采矿条件类似的矿区水平移动系数和主要影响角正切，所述临近地质矿区是指地质矿区的采煤方法和顶板管理方法相同，煤矿上面的岩石力学性质、岩层分布、开采厚度和深度有70％以上的相同。

所述矿区煤层的倾角小于16°。

有益效果

本发明一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，通过利用InSAR技术获取矿区雷达视线向形变场，对视线向形变场中相干性低于解缠阈值的像素进行空间插值，得到空间连续的形变场；利用矿区工作面分布计算各像素点的主要影响半径；获取矿区水平移动系数之后，将矿区地表水平移动转换为下沉值，并按照雷达成像原理组成方程组；求解方程组，得出满足精度要求的地表下沉值的解，然后根据下沉值计算出东西、南北方向倾斜值；最后利用东西、南北方向的倾斜值与水平移动的比例关系计算东西、南北方向的形变场。突破了InSAR求解三维形变场对于数据的苛刻要求，监测费用高等制约，大大的提升了InSAR技术在矿区的应用前景，为矿区大范围、高精度、低成本三维形变监测奠定基础。

另外，由于三维形变场精度除了受到InSAR技术本身误差外，其他的误差源主要是主要影响角正切和水平移动系数的参数误差，因此，本发明中使用实测数据拟合出满足矿区实际情况的主要影响角正切和水平移动系数进行计算，获取的地表三维形变场精度较传统三维形变场监测方法得到了较大的提高。

附图说明

图1为地表单元下沉盆地及水平移动曲线图，图中曲线(a)表示单元开采导致的地表下沉盆地；曲线(b)表示单元开采时地表水平移动曲线；曲线(c)表示地表；曲线(d)表示地下开采单元；曲线(e)表示单元开采时倾斜曲线；

图2为本发明的数据处理流程图；

图3为模拟形变场加入真实SAR数据差分干涉处理得到的视线向形变图；

图4为某矿区的三维形变图，其中图(a)、图(d)、图(g)分别表示矿区地表模拟的东西、南北、垂直的三维形变场；图(b)、图(e)、图(h)分别表示计算的矿区地表东西、南北、垂直三维形变图；图(c)、图(f)、图(i)分别表示模拟与计算形变场在东西、南北、垂直三方向的差值图；

图5为某矿区形变区的剖面图，其中图(a)、图(b)分别表示东西方向形变场中AA’、BB’两个剖面的模拟形变值与理论形变值比较图；图(c)、图(d)分别表示南北方向形变场中CC’、DD’两个剖面的的模拟形变值与理论形变值比较图；图(e)、图(f)分别表示垂直方向形变场中EE’、FF’两个剖面的的模拟形变值与理论形变值比较图。

具体实施方式

为了使本技术领域的人员能够更好的理解本发明的方法，下面将结合附图，对本发明的实施方案进行清楚、详细的描述。

对于地下水平或近水平煤层单元开采时，地表的下沉和水平移动在任意方向上存在一定的内在关系，本实施例以煤层走向x方向为例，其坐标关系如图1所示，即：

$U_e\left(x\right)=B\frac{{\mathrm{dW}}_e\left(x\right)}{\mathrm{dx}}=b\·r\·\frac{{\mathrm{dW}}_e\left(x\right)}{\mathrm{dx}}=b\·r\·T_e\left(x\right)=-b\·r\·\frac{2\mathrm{\πx}}{r^3}{e}^{-\mathrm{\π}\frac{x^2}{r^2}}---\left(1\right)$

式中，b为水平移动系数，r为该盆地主要影响半径，其值等于该点采深H与主要影响角正切tanβ的比值，即r＝H/tanβ，如图1(c)所示W_e(x)为单元开采导致的地表坐标为x处的下沉值，下沉值曲线如图1(a)所示，U_e(x)为该位置的水平移动值，水平移动值的曲线如图1(b)所示，如图1(e)所示为T_e(x)为倾斜值，其大小等于下沉值的导数。

从式(1)中可以看出，在矿区水平或近水平开采时，矿区地表水平移动与下沉值导数之间存在着比例关系，本发明将利用此关系将地表东西、南北方向的水平移动值转换为下沉值的表达式，然后利用合成孔径雷达成像原理将三维形变场与视线向形变场建立方程组，从而解出地表下沉值，最后根据式(1)的关系算出水平移动场。

实施例利用快速拉格朗日分析软件FLAC3D模拟生成采深为700m，松散层200m，煤层倾角为0°，采空区尺寸为1100m×400m的矿区地表三维形变场，如图4中图(a)、图(d)、图(g)所示。将三维形变场加入真实SAR数据进行差分干涉处理，最后得出地理编码后矿区地表视线向形变场。

本发明具体实施步骤如下，如图2所示：

(1)基于InSAR技术利用GAMMA软件得到地理编码后的矿区雷达视线向形变场，即LOS，如图3所示；

(2)由于形变场中部分点的相干性低于解缠阈值而未能参与解缠，导致这些点上没有形变信息，为了获取空间的连续形变场，本发明利用反距离加权法对视线向形变场中没有形变值的像素进行内插处理，其处理结果如图3所示；

相干性是指两幅影像的相似程度；差分干涉图得到的相位差只是该点相位的对2π取模后的值，不是真实值，因此将恢复相位差真实值的过程叫相位解缠。

(3)利用对矿区地表监测的水平移动数据、下沉值数据、采深及采空区分布，基于概率积分法模型求取该矿区的水平移动系数b和主要影响角正切tanβ，根据主要影响角正切和采深信息，按照公式r＝H/tanβ，r为主要影响半径，H为采深，tanβ为主要影响角正切，计算出各像素点的主要影响半径；其中，b＝0.313；tanβ＝2.03；H＝700m；

(4)假设待监测的沉降区域视线向形变影像范围为n行×n列，本实施例中为100行×100列，任意一个像素的下沉值为W(i,j)(i,j＝1,2,…,n)，分别令东西、南北两方向任意像素的倾斜值为T_E(i,j)(i,j＝1,2,…,n)、T_N(i,j)(i,j＝1,2,…,n)，根据开采沉陷倾斜值定义可知：

$\left\{\begin{array}{c}{T}_E(i,j)=[E(i,j+1)-W(i,j)]/\mathrm{\ΔE}\\ T_N(i,j)=[W(i+1,j)-W(i,j)]/\mathrm{\ΔN}\end{array}\right.,(i,j=\mathrm{1,2},\·\·\·,n-1)---\left(2\right)$

式中ΔE、ΔN分别表示地理编码后像素在东西、南北方向的分辨率；本实例中，△E＝15.74547m/像素、△N＝14.9746m/像素，其值从步骤(1)中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得。

从式(2)中，计算的东西南北方向的倾斜值只有n-1行×n-1列，为了便于矩阵计算，且矿区下沉盆地边缘地区下沉值较小，因此可以认为第n行和第n列像素中视线向形变值只由垂直方向形变值贡献，即东西、南北倾斜矩阵中的第n行和第n列全部为0。

将式(2)代入式(1)得出东西、南北方向水平移动形变值，其表达式为：

$\left\{\begin{array}{c}{U}_E(i,j)=b\·r(i,j)\·T_E(i,j)=b\·r(i,j)\·[W(i,j+1)-W(i,j)]/\mathrm{\ΔE}\\ U_N(i,j)=b\·r(i,j)\·T_N(i,j)=b\·r(i,j)\·[W(i+1,j)-W(i,j)]/\mathrm{\ΔN}\end{array}\right.,(i,j=\mathrm{1,2},\·\·\·,n-1)---\left(3\right)$

(i，j)为地理编码后的矿区地表区域内的任一像素点，U_E(i，j)、U_N(i，j)分别为矿区地表东西、南北方向任意像素点的水平移动值，b为矿区水平移动系数，T_E(i,j)、T_N(i,j)分别为东西、南北方向上该像素点的倾斜值，r(i，j)为该像素点的主要影响半径，依据步骤(3)中的r＝H/tanβ进行计算，△E、△N分别为地理编码后像素在东西、南北方向的分辨率，△E＝15.74547m/像素、△N＝14.9746m/像素。

从式(3)可以看出，对于任意像素(i，j)的水平移动值，均可转换为该像素与周围像素下沉值的表达式；根据合成孔径雷达干涉测量原理可知，雷达视线向形变值LOS由南北方向移动形变值U_N、东西方向移动形变值U_E以及垂直方向下沉值W共同决定，其关系为：

LOS＝W cosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]  (4)

式中，θ为雷达入射角，α_h飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得，本实例中两个参数值分别为38.735°、-10.181°。

将式(3)代入式(4)可得：

$\begin{array}{c}\mathrm{LOS}(i,j)=W(i,j)\cos \mathrm{\θ}-\sin \mathrm{\θ}[U_N(i,j)\cos ({\mathrm{\α}}_h-3\mathrm{\π}/2)+U_E(i,j)\sin ({\mathrm{\α}}_h-3\mathrm{\π}/2)]\\ =[\cos \mathrm{\θ}+\sin \mathrm{\θ}\cos ({\mathrm{\α}}_h-\frac{3\mathrm{\π}}{2})b\·r(i,j)/\mathrm{\ΔN}+\sin \mathrm{\θ}\sin ({\mathrm{\α}}_h-\frac{3\mathrm{\π}}{2})b\·r(i,j)/\mathrm{\ΔE}]\·\\ W(i,j)+[\sin \mathrm{\θ}\cos ({\mathrm{\α}}_h-\frac{3\mathrm{\π}}{2})b\·r(i,j)/\mathrm{\ΔN}=W(i+1,j)+\\ [\sin \mathrm{\θ}\sin ({\mathrm{\α}}_h-\frac{3\mathrm{\π}}{2})b\·r(i,j)/\mathrm{\ΔE}]W(i,j+1)\\ =C1\·W(i,j)+C2\·W(i+1,j)+C3\·W(i,j+1)\end{array}---\left(5\right)$

(i,j＝1,2,…,n-1)

式中，C1、C2、C3分别为W(i，j)、W(i+1，j)、W(i，j+1)的系数，其值为：C1＝3.9566、C2＝-2.6935、C3＝-0.4828。

对于待监测区域的第n行和第n列，由于地下开采导致的地表形变边缘地区形变值较小，所以该区域的水平移动对视线向形变贡献很小，可以忽略水平移动值的影响，因此下沉值和雷达视线向关系表达式为：

$\left\{\begin{array}{c}\mathrm{LOS}(i,n)=W(i,n)\cos \mathrm{\θ}=C4\·W(i,n)\\ \mathrm{LOS}(n,s)=W(n,i)\cos \mathrm{\θ}=C4\·W(n,i)\end{array}\right.,i=\mathrm{1,2},\·\·\·n---\left(6\right)$

式中，C4为W(i，n)和W(n，i)的系数，其值为：C4＝0.7802。

从式(5)和(6)中可以看出，对于n×n的视线向形变场，任意像素都有如式(5)或式(6)的表达式，像素点未处于形变区域最后一行或最后一列采用式(5)，像素点位于形变区域最后一行或最后一列采用式(6)，形变区域总计有未知数下沉值n²个，方程总数也为n²个，因此，三维形变场的估算转换为求解一线性方程组，其方程组为：

$B_{n^2\×n^2}\·W_{n^2\×1}^{\′}={\mathrm{LOS}}_{n^2\×1}---\left(7\right)$

式中，W′、LOS分别为n²行×1列的下沉值矩阵和视线向形变值矩阵，矩阵中第i·j行表示研究区域像素坐标为(i,j)的点的下沉值和视线向形变值；B为n²行×n²列的方程系数矩阵，当像素点未处于最后一行和最后一列时，即(i,j＝1,2,…,n-1)，系数阵B中第i·j行、第i·j列、第i·(j+1)行、第i·j列，第i·(j+n)行、第i·j列的值分别为C1、C2、C3，即B(i·j,i·j)＝C1、B[i·(j+1),i·j]＝C2、B[i·(j+n),i·j]＝C3，其它系数均为0。对于处于最后一行或最后一列的像素(i,j)(i＝n,j＝1,…norj＝n,i＝1,…n)，系数阵B(i·j,i·j)＝C4，其它系数均为0。即方程组矩阵为

$W^{\′}={\left[\begin{array}{c}W\left(\mathrm{1,1}\right)\\ W\left(\mathrm{1,2}\right)\\ .\\ .\\ .\\ W(1,n)\\ W\left(\mathrm{2,1}\right)\\ W\left(\mathrm{2,2}\right)\\ .\\ .\\ .\\ .\\ .\\ .\\ W(n,n-1)\\ W(n,n)\end{array}\right]}_{n^2\×1},\mathrm{LOS}={\left[\begin{array}{c}\mathrm{LOS}\left(\mathrm{1,1}\right)\\ \mathrm{LOS}\left(\mathrm{1,2}\right)\\ .\\ .\\ .\\ \mathrm{LOS}(1,n)\\ \mathrm{LOS}\left(\mathrm{2,1}\right)\\ \mathrm{LOS}\left(\mathrm{2,2}\right)\\ .\\ .\\ .\\ .\\ .\\ .\\ \mathrm{LOS}(n,n-1)\\ \mathrm{LOS}(n,n)\end{array}\right]}_{n^2\×1}$

从式(7)方程组中可以看出，对于方程组中下沉值的求解，直接利用求解线性方程组的方法即可解出所有的未知数W，得出矿区地表下沉场，如图4中图(h)所示。

(5)利用式(2)计算东西、南北方向的倾斜值；

(6)利用式(3)计算出南北、东西方向的水平移动场，如图4中图(b)、图(e)所示，完成三维形变场的获取。

为了定量验证本发明的应用效果，实施例中在三维形变场中分别选取了两个剖面对比分析，即AA’、BB’、CC’、DD’、EE’、FF六个剖面，分别对模拟的三维形变场剖面和计算三维形变场剖面做剖面图，如图5所示，模拟的三维形变场剖面由快速拉格朗日分析软件FLAC3D模拟得到，六个剖面的模拟值和计算值吻合非常好，其均方根误差很小，从而说明本发明是可行的和可靠的。

Claims (4)

1.一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，其特征在于，包括以下步骤：

步骤1：利用InSAR技术获取待测矿区地理编码后的雷达视线向形变场；

所述地理编码，是指将雷达影像坐标系转换到通用横轴墨卡托投影(Universal TransverseMercatol projection，简称UTM)坐标系；

步骤2：对视线向形变场中相干性低于解缠时的相干性阈值像素进行空间插值，得到空间连续的形变场，解缠时的相干性阈值大于0.3；

步骤3：利用对矿区地表监测的水平移动数据、下沉值数据、采深及工作面分布，基于概率积分法模型求取该矿区的水平移动系数和主要影响角正切，之后利用各像素点对应的采深计算各像素点的主要影响半径，主要影响半径等于采深除以主要影响角正切，并按照雷达成像原理和矿区东西、南北方向的水平移动值等于水平移动系数、主要影响半径及该方向上该点倾斜值的乘积的关系，组成矿区地表东西、南北方向的任意点的水平移动值与矿区地表东西、南北方向上对应点的下沉值之间的方程组；

所述雷达成像原理由下式表征：

LOS＝W cosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]；

其中，LOS为雷达视线向的形变值，利用步骤1中的InSAR技术获得，W为下沉值，U_N为南北方向水平移动，U_E为东西方向水平移动；θ为雷达卫星入射角，α_h为卫星飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

步骤4：忽略水平移动对视线向形变的影响，求步骤3的方程组的解，得出地表各像素点的下沉值；再求出各像素点在东西、南北方向上下沉值的梯度大小，各像素点在东西、南北方向上下沉值的梯度大小即为各像素点在东西、南北水平方向上的倾斜值；

步骤5：利用各像素点东西、南北方向的水平移动值等于水平移动系数、主要影响半径及该点在东西、南北方向上的倾斜值的乘积关系，计算东西、南北方向的水平移动值，最终得出各像素点的三维形变值，包括各像素点的东西、南北方向的水平移动值及下沉值，即得到整个矿区的三维形变场。

2.根据权利要求1所述的利用单个InSAR干涉对获取矿区地表三维形变场的方法，其特征在于，所述步骤3方程组中的矿区地表水平移动值与下沉值之间的转换关系如下：

$\left\{\begin{array}{c}{U}_E(i,j)=b\·r(i,j)\·T_E(i,j)=b\·r(i,j)\·[W(i,j+1)-W(i,j)]/\mathrm{\ΔE}\\ U_N(i,j)=b\·r(i,j)\·T_N(i,j)=b\·r(i,j)\·[W(i+1,j)-W(i,j)]/\mathrm{\ΔN}\end{array}\right.,(i,j=\mathrm{1,2},\·\·\·,n-1)$

(i，j)为地理编码后的矿区地表区域内的任一像素点，U_E(i，j)、U_N(i，j)分别为矿区地表东西、南北方向任意像素点的水平移动值，b为矿区水平移动系数，T_E(i,j)、T_N(i,j)分别为东西、南北方向上该像素点的倾斜值，r(i，j)为该像素点的主要影响半径；△E、△N分别为地理编码后像素在东西、南北方向的分辨率，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得。

3.根据权利要求1或2任一项所述的利用单个InSAR干涉对获取矿区地表三维形变场的方法，其特征在于，所述水平移动系数和主要影响角正切，选用临近地质采矿条件类似的矿区水平移动系数和主要影响角正切，所述临近地质矿区是指地质矿区的采煤方法和顶板管理方法相同，工作面上面的岩石力学性质、岩层分布、开采厚度和深度有70％以上的相同。

4.根据权利要求3所述的利用单个InSAR干涉对获取矿区地表三维形变场的方法，其特征在于，所述矿区煤层的倾角小于16°。