CN106767380A - 一种基于两景sar强度影像的矿区地表大量级三维形变估计方法 - Google Patents

Abstract

本发明公开了一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，利用SAR强度偏移量追踪算法从两景SAR强度影像中生成矿区地表雷达视线向和方位向形变；基于矿区地表水平移动与下沉梯度的比例关系建立东西、南北方向水平移动与下沉之间的函数模型；分别建立矿区地表下沉与视线向和方位向形变的观测方程组；利用最小二乘法估计求解矿区地表下沉；基于求解的地表下沉，使用矿区地表水平移动与下沉梯度的比例关系估计矿区地表在东西、南北方向的水平移动。实现了仅利用两景SAR强度影像估计矿区地表大量级三维形变，有效地克服了传统方法无法从单个InSAR干涉对中获取大量级三维形变的局限。

Description

一种基于两景SAR强度影像的矿区地表大量级三维形变估计 方法

技术领域

本发明涉及一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法。

背景技术

矿区地表在东西、南北和垂直方向的三维形变监测对于提前评估和控制地下开采导致的潜在地质灾害和建构筑物损坏有着重要作用。

合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)技术是一种新兴的雷达遥感技术。该技术基于两景SAR单视复数(single look complex,SLC)影像生成的单个InSAR干涉对，实现了地表雷达视线方向(line-of-sight,LOS)形变的大范围、高精度、低成本、非接触监测。然而，由于InSAR技术仅能从单个InSAR干涉对中获取一维LOS向形变，而非真正的地表三维形变。从而大大地制约了InSAR技术在矿区的应用前景。为了获取矿区地表三维形变，传统的方法至少需要三个来自于不同雷达成像几何学SAR数据的InSAR干涉对，即至少需要六景SLC影像。受限于目前较少的可用SAR卫星数量，该方法在实际应用中很难满足。

发明专利ZL201210440875.4提出了一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，该发明有效地克服了传统方法无法从单个InSAR干涉对中获取矿区地表三维形变的局限，大大降低了传统方法对于数据的苛刻要求，为矿区地质灾害和建构筑物风险评估提供了重要的参考依据。由于InSAR技术是基于两景SLC影像的干涉相位的，所以当干涉相位失相关时，InSAR技术则无法准确地获取地表雷达视线向形变。然而，大量级形变常常导致InSAR干涉相位失相关，从而致使发明专利ZL201210440875.4无法获取矿区大量级三维形变。

发明内容

本发明提供了一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，其目的在于克服传统基于单个InSAR干涉对的方法无法获取矿区地表大量级三维形变的局限，利用矿区地表水平移动与下沉梯度的比例关系从SAR强度偏移量追踪算法得到的矿区地表视线向和方位向形变中估计矿区地表沿着垂直、东西和南北方向的三维形变，大大地拓宽了SAR技术在矿区的应用前景。

一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，包括以下几个步骤：

步骤1：利用SAR强度偏移量追踪算法从两景SAR强度影像中生成待监测矿区地表雷达视线方向d_LOS和方位方向形变d_AZI；

待监测的矿区的形变图的总行数和总列数分别为n和m；

步骤2：基于矿区地表水平移动与下沉梯度的比例关系建立东西、南北方向的水平移动U_E、U_N与下沉W的函数关系；

步骤3：分别根据视线向形变d_LOS和方位向形变d_AZI与地表真实三维形变的投影关系建立矿区地表下沉与利用偏移量追踪算法获取的视线向与方位向形变之间的投影方程；

步骤4：按照步骤3构建的投影方程建立待监测矿区地表各像素点的下沉值与视线向形变的观测方程；

由于矿区边缘受地下开采影响较小，在待监测矿区最后一行和最后一列的水平移动设置为0；

步骤5：求解待监测矿区各像素点的下沉值；

步骤6：基于求解的下沉值并利用步骤2中所述的东西、南北方向的水平移动与下沉梯度的关系估计待监测矿区地表在东西、南北方向的水平移动。

所述利用偏移量追踪算法获取的视线向与方位向形变中的偏移量追踪算法属于现有技术。

进一步地，所述基于矿区地表水平移动与下沉梯度的比例关系建立的东西、南北方向的水平移动U_E、U_N与下沉W的函数关系如下：

其中，i和j为地表点的像素坐标；W(i,j)表示地表点(i,j)处的下沉值，C_E和C_N分别示地表点(i,j)处在东西、南北方向的下沉梯度比例系数，b、H和β分别表示矿区水平移动系数、采深和主要影响角，b、H和β从监测矿区获取；R_E和R_N分别表示视线向和方位向形变图在东西和南北方向的空间分辨率，R_E和R_N从视线向和方位向形变图中获取。

进一步地，所述矿区地表下沉与利用偏移量追踪算法获取的视线向与方位向形变之间的投影方程如下：

和/或

其中，为投影系数矩阵；

A₁₁(i,j)＝-A₁₂(i,j)-A₁₃(i,j)，A₂₁(i,j)＝cosθ-A₂₂(i,j)-A₂₃(i,j)；

A₁₂(i,j)＝sin(α-3π/2)C_E(i,j)，A₂₂(i,j)＝-b·C_E(i,j)·sin(α-3π/2)；

A₁₃(i,j)＝cos(α-3π/2)C_N(i,j)，A₂₃(i,j)＝-b·C_N(i,j)·cos(α-3π/2)；

α和θ表示雷达飞行方位角和入射角，从SAR强度影像的头文件中获得。

用方位向或者视线向监测三维形变，即用任意一个投影方程来估计三维形变时，采用线性方程求解算法直接求解矿区下沉值；

两者联合使用时，采用最小二乘求解算法直接求解矿区下沉值，联合使用时增加了多余观测，使得下沉值的求解结果精度更高。

进一步地，所述待监测矿区地表各像素点的下沉值与视线向形变的观测方程如下：

其中，W_i和W_i+1分别表示待监测的矿区形变图中第i行和第i+1行像素中各像素的下沉值；

W_i＝[W(i,1),W(i,2),…,W(i,m)]^T，W_i+1＝[W(i+1,1),W(i+1,2),…,W(i+1,m)]^T；

d_LOSi和d_AZIi分别表示待监测的矿区形变图中第i行像素中各像素的视线方向和方位方向形变观测值；

d_LOSi＝[d_LOS(i,1),d_LOS(i,2),…,d_LOS(i,m)]^T，d_AZIi＝[d_AZI(i,1),d_AZI(i,2),…,d_AZI(i,m)]^T；

C_i、C_i′、D_i及D_i′表示观测系数矩阵；

C_i′＝diag[A₁₃(i,1) A₁₃(i,2) … A₁₃(i,m-1) 0]_m×m,

D_i′＝diag[A₂₃(i,1) A₂₃(i,2) … A₂₃(i,m-1) 0]_m×m

其中，待监测的矿区形变图中最后一行像素中各像素的视线方向和方位方向形变观测值为D″＝diag(cosθ cosθ … cosθ)_m×m，diag表示对角阵。

进一步地，按照下述公式求解待监测矿区各像素点的下沉值：

利用最小二乘法求解：

其中，E_i和L_i表示中间变量矩阵， 为E_i的转置矩阵， 表示的逆矩阵；

其中，最后一行各像素的下沉值按照公式计算获得：W_n＝(D″)^-1·d_LOSn。

有益效果

本发明公开了一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，利用SAR强度偏移量追踪算法从两景SAR强度影像中生成矿区地表雷达视线向和方位向形变；基于矿区地表水平移动与下沉梯度的比例关系建立东西、南北方向水平移动与下沉之间的函数模型；分别建立矿区地表下沉与视线向和方位向形变的观测方程组；利用最小二乘法估计求解矿区地表下沉；基于求解的地表下沉，使用矿区地表水平移动与下沉梯度的比例关系估计矿区地表在东西、南北方向的水平移动。本发明实现了仅利用两景SAR强度影像估计矿区地表大量级三维形变，有效地克服了传统方法无法从单个InSAR干涉对中获取大量级三维形变的局限。此外，该方法数据制约少，监测结果准确有效，大大拓宽了InSAR技术的应用前景。本方法其对于拓宽SAR技术的应用空间，降低矿区三维形变监测成本有着重要意义。此外，其对于指导矿区安全生产、分析矿区沉降机理、预警矿区地表地质灾害以及生态环境保护也起着重要作用。

附图说明

图1为本发明的所述方法的流程示意图。

图2为应用本发明所述方法得到的形变图，其中，图(a)为视线向形变图，图(b)为方位向形变图；

图3为应用本发明所述方法获得的模拟的三维形变与实际监测获得的三维形变对比示意图，其中，图(a)为应用本发明所述方法模拟获得的垂直方向形变图，图(b)为图(a)对应的监测获得的形变图，图(c)为AA’剖面的垂直模拟形变值和监测值的对比图，图(d)为应用本发明所述方法模拟获得的东西方向形变图，图(e)为图(d)对应的监测获得的形变图，图(f)为AA’剖面在东西方向上的模拟形变值和监测值的对比图，图(g)为应用本发明所述方法模拟获得的南北方向形变图，图(h)为图(g)对应的监测获得的形变图，图(i)为BB’剖面在南北方向上的模拟形变值和监测值得对比图。

具体实施方式

下面将结合附图1-3对本发明做进一步的说明。

如图1所示，一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，包括以下步骤：

步骤1：实施例利用快速拉格朗日分析软件FLAC3D模拟生成采深为550m，煤层倾角为0°，采厚为5m，采空区尺寸为700m×150m的矿区地表三维形变场，如图3中图(a)、图(d)、图(g)所示。之后将生成的三维形变按照入射角为38°、飞行方位角为350°投影到雷达视线向和方位向(即用式(2)计算)；为了模拟偏移量追踪技术导致的形变误差，实施例在投影后的视线向和方位向形变中加入-0.1到0.1之间服从均匀分布的误差，并以此模拟利用SAR强度偏移量追踪算法从两景SAR强度影像中获取矿区地表沿着视线方向d_LOS和方位方向的形变d_AZI，其结果分别如图2中(a)和(b)所示；

步骤2：基于矿区地表水平移动与下沉梯度的比例关系建立东西、南北方向的水平移动U_E、U_N与下沉W的函数关系

其中，i和j为地表点的像素坐标；W(i,j)表示地表点(i,j)处的下沉值，C_E和C_N分别示地表点(i,j)处在东西、南北方向的下沉梯度比例系数，b、H和tanβ分别表示矿区水平移动系数、采深和主要影响角正切，其中，b＝0.35、H＝550m和tanβ＝2.0；R_E和R_N分别表示视线向和方位向形变图在东西和南北方向的空间分辨率；其中，R_E＝3m和R_N＝3m；

步骤3：视线向d_LOS和方位向形变d_AZI是地表真实三维形变按照雷达成像几何学的投影：

其中，θ和α为雷达的入射角和飞行方位角，其中，θ＝38°，α＝350°。

将式(1)代入式(2)可得：

其中，为投影系数矩阵；

A₁₁(i,j)＝-A₁₂(i,j)-A₁₃(i,j)，A₂₁(i,j)＝cosθ-A₂₂(i,j)-A₂₃(i,j)；

A₁₂(i,j)＝sin(α-3π/2)C_E(i,j)，A₂₂(i,j)＝-b·C_E(i,j)·sin(α-3π/2)；

A₁₃(i,j)＝cos(α-3π/2)C_N(i,j)，A₂₃(i,j)＝-b·C_N(i,j)·cos(α-3π/2)；

步骤4：由于矿区边缘受地下开采影响较小，因此假设在待监测矿区最后一行和最后一列的水平移动为0，即：

其中，n和m为待监测矿区的像素总行数和总列数，在实施例中，n＝301，m＝467。

根据式(3)和(4)建立矿区各像素下沉值与视线向和方位形变之间的观测方程组，该观测方程组以矩阵形式表示为：

其中，W_i和W_i+1分别表示待监测的矿区形变图中第i行和第i+1行像素中各像素的下沉值；

W_i＝[W(i,1),W(i,2),…,W(i,m)]^T，W_i+1＝[W(i+1,1),W(i+1,2),…,W(i+1,m)]^T；

d_LOSi和d_AZIi分别表示待监测的矿区形变图中第i行像素中各像素的视线方向和方位方向形变观测值；

d_LOSi＝[d_LOS(i,1),d_LOS(i,2),…,d_LOS(i,m)]^T，d_AZIi＝[d_AZI(i,1),d_AZI(i,2),…,d_AZI(i,m)]^T；

C_i、C_i′、D_i及D_i′表示观测系数矩阵；

C_i′＝diag[A₁₃(i,1) A₁₃(i,2) … A₁₃(i,m-1) 0]_m×m,

D_i′＝diag[A₂₃(i,1) A₂₃(i,2) … A₂₃(i,m-1) 0]_m×m

D″＝diag(cosθ cosθ … cosθ)_m×m，diag表示对角阵

步骤5：式(6)中有n未知数W_n和n个线性无关的方程，求解式(6)可得：

W_n＝(D″)^-1·d_LOSn (7)

将解出的W_n代入式(5)后，式(5)中有n未知数W_n-1和2n个方程，所以W_n-1的最小二乘为：

其中，E_i和L_i表示中间变量矩阵， 为E_i的转置矩阵，将解出的W_n-1代入式(5)，即可利用式(8)求解出W_n-2，重复该过程直至解出W₁，如图3(b)所示；

步骤6：基于求解的下沉值利用式(1)估计待监测矿区地表在东西、南北方向的水平移动，如图3(e)和3(h)所示。

通过对比模拟的三维形变与监测的三维形变(如图3所示)发现，两者吻合较好，其在垂直、东西和南北方向的均方根误差分别为0.003、0.03和0.06m。该结果表明本发明是可行和可靠的。为了更进一步验证本发明的应用效果，实施例中在三维形变场中分别选取了AA’剖面的下沉和东西方向水平移动以及BB’剖面在南北方向的水平移动对比分析，结果亦表明：两个剖面上的模拟值和计算值吻合非常好，从而说明本发明是可行的和可靠的。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于两景SAR强度影像的矿区地表大量级三维形变估计方法，其特征在于，包括以下几个步骤：

步骤1：利用SAR强度偏移量追踪算法从两景SAR强度影像中生成待监测矿区地表雷达视线方向d_LOS和方位方向形变d_AZI；

待监测的矿区的形变图的总行数和总列数分别为n和m；

步骤2：基于矿区地表水平移动与下沉梯度的比例关系建立东西、南北方向的水平移动U_E、U_N与下沉W的函数关系；

步骤3：分别根据视线向形变d_LOS和方位向形变d_AZI与地表真实三维形变的投影关系建立矿区地表下沉与利用偏移量追踪算法获取的视线向与方位向形变之间的投影方程；

步骤4：按照步骤3构建的投影方程分别建立待监测矿区地表各像素点的下沉值与视线向形变和方位向形变的观测方程组；

步骤5：求解待监测矿区各像素点的下沉值；

步骤6：基于求解的下沉值并利用步骤2中所述的东西、南北方向的水平移动与下沉梯度的关系估计待监测矿区地表在东西、南北方向的水平移动。

2.根据权利要求1所述的方法，其特征在于，所述基于矿区地表水平移动与下沉梯度的比例关系建立的东西、南北方向的水平移动U_E、U_N与下沉W的函数关系如下：

$\left\{\begin{array}{c}{U}_E(i,j)=\frac{b\·H(i,j)}{\tan \mathrm{\β}}\frac{W(i,j+1)-W(i,j)}{R_E}=C_E\[W(i,j+1)-W(i,j)\]\\ U_N(i,j)=\frac{b\·H(i,j)}{\tan \mathrm{\β}}\frac{W(i+1,j)-W(i,j)}{R_N}=C_N\[W(i+1,j)-W(i,j)\]\end{array}\right.$

其中，i和j为地表点的像素坐标；W(i,j)表示地表点(i,j)处的下沉值，C_E和C_N分别示地表点(i,j)处在东西、南北方向的下沉梯度比例系数，H和β分别表示矿区水平移动系数、采深和主要影响角，b、H和β从监测矿区获取；R_E和R_N分别表示视线向和方位向形变图在东西和南北方向的空间分辨率，R_E和R_N从视线向和方位向形变图中获取。

3.根据权利要求2所述的方法，其特征在于，所述矿区地表下沉与利用偏移量追踪算法获取的视线向与方位向形变之间的投影方程如下：

$d_{AZI}(i,j)=A_1(i,j)\left[\begin{array}{c}W(i,j)\\ W(i,j+1)\\ W(i+1,j)\end{array}\right]$

和/或

$d_{LOS}(i,j)=A_2(i,j)\left[\begin{array}{c}W(i,j)\\ W(i,j+1)\\ W(i+1,j)\end{array}\right]$

其中，为投影系数矩阵；

A₁₁(i,j)＝-A₁₂(i,j)-A₁₃(i,j)，A₂₁(i,j)＝cosθ-A₂₂(i,j)-A₂₃(i,j)；

A₁₂(i,j)＝sin(α-3π/2)C_E(i,j)，A₂₂(i,j)＝-b·C_E(i,j)·sin(α-3π/2)；

A₁₃(i,j)＝cos(α-3π/2)C_N(i,j)，A₂₃(i,j)＝-b·C_N(i,j)·cos(α-3π/2)；

α和θ表示雷达飞行方位角和入射角，从SAR强度影像的头文件中获得。

4.根据权利要求3所述的方法，其特征在于，所述待监测矿区地表各像素点的下沉值与视线向形变的观测方程组如下：

${\left[\begin{array}{cc}{C}_i& C_i^{\′}\\ D_i& D_i^{\′}\end{array}\right]}_{2m\×2m}{\left[\begin{array}{c}{W}_i\\ W_{i+1}\end{array}\right]}_{2m\×1}={\left[\begin{array}{c}{d}_{LOSi}\\ d_{AZIi}\end{array}\right]}_{2m\×1}$

其中，W_i和W_i+1分别表示待监测的矿区形变图中第i行和第i+1行像素中各像素的下沉值；

W_i＝[W(i,1),W(i,2),…,W(i,m)]^T，W_i+1＝[W(i+1,1),W(i+1,2),…,W(i+1,m)]^T；

d_LOSi和d_AZIi分别表示待监测的矿区形变图中第i行像素中各像素的视线方向和方位方向形变观测值；

d_LOSi＝[d_LOS(i,1),d_LOS(i,2),…,d_LOS(i,m)]^T，d_AZIi＝[d_AZI(i,1),d_AZI(i,2),…,d_AZI(i,m)]^T；

C_i、C_i′、D_i及D_i′表示观测系数矩阵；

C_i′＝diag[A₁₃(i,1) A₁₃(i,2) … A₁₃(i,m-1) 0]_m×m,

D_i′＝diag[A₂₃(i,1) A₂₃(i,2) … A₂₃(i,m-1) 0]_m×m

其中，待监测的矿区形变图中最后一行像素中各像素的视线方向和方位方向形变观测值为D″＝diag(cosθ cosθ … cosθ)_m×m，diag表示对角阵。

5.根据权利要求4所述的方法，其特征在于，按照下述公式求解待监测矿区各像素点的下沉值：

利用最小二乘法求解：

其中，E_i和L_i表示中间变量矩阵， 为E_i的转置矩阵， 表示的逆矩阵；

其中，最后一行各像素的下沉值按照公式计算获得：W_n＝(D″)^-1·d_LOSn。