基于单个雷达成像几何学SAR 影像的矿区三维时序形变监测方法

摘要

本发明公开了一种基于单个雷达成像几何学SAR影像的矿区三维时序形变监测方法，利用单个雷达成像几何学SAR数据生成可用InSAR干涉对；并生成矿区地表多时域差分下沉观测值；建立时间相邻SAR影像期间的下沉速度与多时域下沉观测值之间的观测方程，并使用加权最小二乘法求解下沉速度；使用解出的下沉速度计算矿区地表在SAR影像获取时间的时序下沉；基于矿区地表东西、南北方向水平移动与下沉在对应方向的梯度成比例的关系，使用计算的矿区时序下沉估计出东西、南北方向的时序形变。本发明实现了仅利用单个雷达成像几何学SAR数据监测矿区地表三维时序形变，大大降低了传统方法中需要三个或以上的不同雷达成像几何学SAR数据的苛刻要求。

说明书

技术领域

本发明涉及一种基于单个雷达成像几何学SAR影像的矿区三维时序形变监测方法。

背景技术

矿区地表三维时序形变监测对于提前评估和控制地下开采导致的潜在地质灾害和建构筑物损坏有着重要作用。利用传统大地测量技术(如GPS、水准测量等)监测矿区地表三维时序形变不仅成本高、效率低，且监测范围较小。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，简称InSAR)是一种能够获得地表厘米甚至毫米级形变的遥感技术。其基本原理就是通过对两幅或以上的合成孔径雷达(Synthetic Aperture Radar，简称SAR)卫星影像进行差分干涉处理，从相位差中提取厘米甚至毫米级的雷达视线方向(line-of-sight,LOS)形变值。相对于传统的大地测量技术而言，InSAR具有全天候、大范围覆盖、高精度、低成本等优势。

时序InSAR技术为高级InSAR技术。该技术基于同一雷达成像几何学的多时域SAR影像实现了地表时序形变监测。然而，由于当前的雷达传感器均为斜视，所以时序InSAR技术获取的地表时序形变是沿着雷达LOS方向的，，其为地表真实三维形变按照雷达成像几何学的投影。因此，若想从一维LOS向时序形变中分解出地表三维时序形变，则至少需要三个不同雷达成像几何学的多时域SAR影像。然而，由于当前可用SAR卫星较少，且矿区大形变梯度容易导致InSAR干涉对失相关而不可用。所以，找到覆盖同一个矿区同一时间段内的三个以上的不同雷达成像几何学的多时域SAR数据几乎不可能。这种局限也大大制约了InSAR技术在矿区三维时序形变获取中的应用。

鉴于矿区地表水平移动与下沉梯度之间为比例关系，之前提出了“一种利用单个InSAR干涉对获取矿区地表三维形变场的方法”(专利号：CN201210440875)(简称“SIP”)。其仅利用单个InSAR干涉对即可获取矿区地表三维形变场，大大降低了传统方法对SAR数据的苛刻要求，拓展了InSAR技术在矿区的应用前景。然而，由于该方法仅能获取矿区地表在组成单个InSAR干涉对的两景SAR影像期间的地表三维形变，而不能获取三维时序形变。因此，该方法的应用前景受到了很大的制约，比如无法实现矿区地表动态地质灾害风险评估。其次，该方法中估计的下沉值的误差在传播到水平移动的过程中容易被放大，从而降低其获取的三维形变的精度。

发明内容

本发明提供了一种基于单个雷达成像几何学SAR影像的矿区三维时序形变监测方法，其目的在于克服传统InSAR时序三维形变获取方法至少需要三个不同雷达成像几何学SAR数据的苛刻要求，降低了监测成本、提高了监测精度。

一种基于单个雷达成像几何学SAR影像的矿区三维时序形变监测方法，包括以下几个步骤：

步骤1：利用单个雷达成像几何学SAR影像生成覆盖待监测矿区的可用InSAR干涉对；

步骤2：利用现有的SIP方法生成待监测矿区多时域差分下沉观测值ΔW：

ΔW＝[ΔW₁,ΔW₂,…,ΔW_m]^T；

其中，m表示可用InSAR干涉对的个数；

所述多时域差分下沉观测值指一系列与可用InSAR干涉对时间同步的下沉值；

步骤3：建立矿区地表各点在两相邻SAR影像时刻间的平均下沉速度V＝[V₁,V₂,…,V_n]^T与多时域差分下沉观测值ΔW＝[ΔW₁,ΔW₂,…,ΔW_m]^T之间的函数关系：B_m×n·V_n×1＝ΔW_m×1；

其中，SAR影像个数为n+1；B为m×n的可用InSAR干涉对辅影像和主影像获取时间差系数矩阵；

对于B的任意第k行，第IM_k个元素前的所有元素均为0，从第IM_k到第IS_k-1个元素，依次为：第IS_k-1个元素以后的所有元素均为0；

其中，IM_k和IS_k分别为生成第k个可用InSAR干涉对的主、辅影像获取时间索引，根据InSAR干涉对的组成情况获得；

本质上，即为在和之间沉降的累计相加；

对于任意的第k个InSAR干涉对，其差分下沉观测值为ΔW_k，主影像的时间为辅影像时间为

用速度表示该观测值为(即为在和之间沉降的累计相加)。在系数矩阵B_m×n中，只有从第IM_k个元素开始才有值，一直到IS_k-1结束，在这一行里，其他元素全部为0：

即其余元素均为0；

在系数矩阵B_m×n中的第k行的具体表达形式如下：

$BV=\left(\begin{array}{cccccccccc}\mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& & & & & \\ .& .& .& .& .& & & & & \\ .& .& .& .& .& & & & & \\ .& .& .& .& .& & & & & \\ 0& 0& \mathrm{...}& t_{{\mathrm{IM}}_k+1}-t_{{\mathrm{IM}}_k}& t_{{\mathrm{IM}}_k+2}-t_{{\mathrm{IM}}_k-1}& \mathrm{...}& t_{{\mathrm{IS}}_k}-t_{{\mathrm{IS}}_k-1}& 0& 0& \mathrm{...}\\ .& & & & & & & & .& .\\ .& & & & & & & & .& .\\ .& & & & & & & & .& .\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\\ .\\ .\\ .\\ v_n\end{array}\right).$

步骤4：求解地表各点在时间相邻SAR影像时刻的平均下沉速度V：

V＝[B^T·B]^-1·[B^T·ΔW]；

步骤5：估计矿区地表时序下沉W，并利用时序下沉计算获得东西、南北方向的时序水平移动，完成待监测矿区三维时序形变监测；

$W=[W_1,W_2,...,W_n],W_k=\underset{l=1}{\overset{k}{\Sigma }}(t_l-t_{l-1})·V_l,(k=1,2,...,n);$

其中，W_k表示在t_k时刻矿区地表的时序下沉值，t_l表示时间，V_l表示t_l-1至t_l期间的平均下沉速度；

$E_k(i,j)=\frac{b·H(i,j)}{tan\beta }·\frac{W_k(i,j+1)-W_k(i,j)}{R_E}$

$N_k(i,j)=\frac{b·H(i,j)}{tan\beta }·\frac{W_k(i+1,j)-W_k(i,j)}{R_N}$

其中，E(i,j)和N(i,j)分别表示像素坐标为(i,j)的点在t_k时刻的东西、南北方向的水平移动，k＝1,2,…,n；b为待监测矿区水平移动系数，H为开采采深，tanβ为主要影响角正切，R_E和R_N分别为LOS向形变图在东西和南北方向的空间分辨率。

所述求解地表各点在时间相邻SAR影像时刻的平均速度V时采用加权最小二乘法求解：

V＝[B^T·Q·B]^-1·[B^T·Q·ΔW]

其中，Q地表各点多时域差分下沉值的加权值，

Q＝diag[1000·γ₁³,1000·γ₂³,…,1000·γ_m³]，diag表示主对角矩阵；γ为各干涉对在该点的相干性。

所述SIP方法是一种利用单个InSAR干涉对获取矿区地表三维形变场的方法，基于矿区地表水平移动与下沉梯度的比例关系，利用单个InSAR干涉对获取的一维雷达视线向形变获取了矿区地表在该时间段内沿着垂直、东西和南北方向的三维差分形变。

有益效果

本发明公开了一种基于单个雷达成像几何学SAR影像的矿区三维时序形变监测方法，利用单个雷达成像几何学SAR数据生成可用InSAR干涉对；利用现有的基于单个InSAR干涉对的矿区地表三维形变获取方法处理每一个可用InSAR干涉对，从而生成矿区地表多时域下沉观测值；建立时间相邻SAR影像期间的下沉速度与多时域下沉观测值之间的观测方程，并使用加权最小二乘法求解下沉速度；使用解出的下沉速度计算矿区地表在SAR影像获取时间的时序下沉；基于矿区地表东西、南北方向水平移动与下沉在对应方向的梯度成比例的关系，使用计算的矿区时序下沉估计出东西、南北方向的时序形变。本发明实现了仅利用单个雷达成像几何学SAR数据监测矿区地表三维时序形变，大大降低了传统方法中需要三个或以上的不同雷达成像几何学SAR数据的苛刻要求，该方法构思巧妙，数据制约少，监测结果准确有效，大大拓宽了InSAR技术的应用前景。通过加权最小二乘和大量的观测值，有效地抑制了InSAR视线向形变误差传播到时序下沉中。此外，本发明利用误差抑制后的时序下沉计算东西、南北方向的时序水平移动，其比直接使用SPI方法估计获得的结果的精度高。因而能有效地降低了水平移动中的误差传播，提高了三维时序形变的精度和可靠性。

因此，本方法其对于拓宽InSAR应用空间，降低矿区三维时序形变监测成本，提高InSAR矿区三维形变监测精度有着重要意义。此外，其对于指导矿区安全生产、预警矿区地表地质灾害以及生态环境保护也起着重要作用。

附图说明

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

图2表示走向长壁工作面模拟示意图；

图3表示垂直Ws、东西Es和南北Ns方向的三维时序形变的模拟示意图，其中Ws₀＝0，Es₀＝0，Ns₀＝0；

图4表示InSAR监测LOS向形变模拟示意图，其中，d_LOS(t_i-t_j)表示该干涉对由获取时间为t_i和t_j的两景SAR影像生成；

图5为采用本发明所述方法估计的垂直W、东西E和南北N方向三维时序形变示意图，其中，W₀＝0，E₀＝0，N₀＝0。

具体实施方式

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

步骤1：利用单个雷达成像几何学SAR影像生成可用InSAR干涉对；

假设有n+1幅覆盖待监测矿区的单个雷达成像几何学SAR影像，其获取时间分别为(t₀,t₁,…,t_n)(其中t₀＜t₁＜…＜t_n)。

设定InSAR干涉对的时间基线和空间基线阈值，并利用SAR影像生成时间基线和空间基线均小于对应阈值的可用InSAR干涉对。假定可用InSAR干涉对的数量m，且覆盖整个时域过程(即从时间t₀到t_n期间均有InSAR干涉对覆盖)。令IM和IS分别表示组成m个可用InSAR干涉对的时间索引，即

$\{\begin{array}{c}IM=[{\mathrm{IM}}_1,{\mathrm{IM}}_2,...,{\mathrm{IM}}_m]\\ IS=[{\mathrm{IS}}_1,{\mathrm{IS}}_2,...,{\mathrm{IS}}_m]\end{array},---\left(1\right)$

其中，且比如第k个可用InSAR干涉对由时间为分别为t₄和t₅的两幅主、辅SAR影像生成，则IM_k＝4且IS_k＝5。

步骤2：生成多时域差分下沉观测值；

随后，利用现有的方法“一种利用单个InSAR干涉对获取矿区地表三维形变场的方法”(以下简称SIP方法)处理每一个可用InSAR干涉对，得到m个沿着垂直方向的多时域差分下沉观测值ΔW＝[ΔW₁,ΔW₂,…,ΔW_m]^T。

步骤3：建立矿区地表各点在两相邻SAR影像时刻间的平均下沉速度V＝[V₁,V₂,…,V_n]^T与多时域差分下沉观测值ΔW＝[ΔW₁,ΔW₂,…,ΔW_m]^T之间的观测方程；

令像素坐标为(i,j)的地表点在时间相邻SAR影像期间的平均下沉速度为V(i,j)＝[V₁(i,j),V₂(i,j),…,V_n(i,j)]^T，即

$V_k(i,j)=\frac{W_k(i,j)-W_{k-1}(i,j)}{t_k-t_{k-1}},(k=1,2,\mathrm{...},n),---\left(2\right)$

其中，W_k(i,j)和W_k-1(i,j)表示点(i,j)在t_k和t_k-1时刻相对于t₀时刻(即W₀(i,j)＝0)的时序下沉值。因此，对于在点(i,j)在SAR影像获取时刻的时序沉降值可表示为：

$W_k(i,j)=\underset{l=1}{\overset{k}{\Sigma }}(t_l-t_{l-1})·V_l(i,j),(k=1,2,...,n).---\left(3\right)$

从式(3)可以看出，若想获取矿区地表点(i,j)处的时序沉降值，其在时间相邻SAR影像期间的平均速度V(i,j)＝[V₁(i,j),V₂(i,j),…,V_n(i,j)]^T首先需要被估计。也就是说，对于矿区地表任意一点，其有n个未知数需要估计。理论上，若能够组建至少n个关于未知数的线性无关方程组，即可解出这n个未知数，从而估计出矿区该点在SAR影像获取时刻的下沉值。在矿区地表点(i,j)上，第k个SIP方法获取的下沉观测值ΔW_k(i,j)可以表示为：

$\begin{array}{c}{\mathrm{\Delta W}}_k(i,j)=W_{{\mathrm{IM}}_k}(i,j)-W_{{\mathrm{IS}}_k}(i,j)\\ =\underset{l={\mathrm{IS}}_k+1}{\overset{{\mathrm{IM}}_k}{\Sigma }}(t_l-t_{l-1})·V_l(i,j)\end{array},(k=1,2,\mathrm{...},m).---\left(4\right)$

从式(4)中可以看出，对于每个SIP方法获取的差分下沉值ΔW_k(i,j)(k＝1,2,…,m)均可建立如式(4)的线性方程。若可用InSAR干涉对覆盖整个时域过程(即m≥n)，即可解出点(i,j)的n个未知数。为了便于表示，将式(4)表示为矩阵形式：

B_m×n·V(i,j)_n×1＝ΔW(i,j)_m×1(5)

其中B为m×n的系数矩阵，其形式取决于可用InSAR干涉对的组成，对于任意第k行的形式为：其余元素均为0。比如ΔW₁(i,j)＝W₁(i,j)-W₀(i,j)、ΔW₂(i,j)＝W₃(i,j)-W₀(i,j)，则式(5)可表示为：

${\left(\begin{array}{ccccc}{t}_1-t_0& 0& 0& 0& \mathrm{...}\\ t_1-t_0& t_2-t_1& t_3-t_2& 0& \mathrm{...}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\end{array}\right)}_{m\times n}·{\left(\begin{array}{c}{V}_1(i,j)\\ V_2(i,j)\\ .\\ .\\ .\\ V_N(i,j)\end{array}\right)}_{n\times 1}={\left(\begin{array}{c}{\mathrm{\Delta W}}_1(i,j)\\ {\mathrm{\Delta W}}_2(i,j)\\ .\\ .\\ .\\ {\mathrm{\Delta W}}_M(i,j)\end{array}\right)}_{m\times 1}---\left(6\right)$

假设有n+1景SAR影像，其获取时间按照先后顺序排序为t＝[t₁,t₂,…,t_n]，生成了m个可用InSAR干涉对。令IM＝[IM₁,IM₂,…,IM_m]为InSAR干涉对的主影像时间索引，IS＝[IS₁,IS₂,…,IS_m]为辅影像时间索引。所述时间索引为影像获取时间的下标。

1)假设干涉对1是时间为t₁和t₂的SAR影像生成，则ΔW₁为t₁和t₂之间的差分下沉观测值，其等于ΔW₁＝v₁·(t₂-t₁)；

2)假设干涉对2是时间为t₁和t₃的SAR影像生成，则ΔW₂为t₁和t₃之间的差分下沉观测值，其等于ΔW₁＝v₁·(t₂-t₁)+v₂·(t₃-t₂)(两个时间段内的沉降相加)；

3)假设干涉对3是时间为t₂和t₄的SAR影像生成，则ΔW₃为t₂和t₄之间的差分下沉观测值，其等于ΔW₃＝v₂·(t₃-t₂)+v₃·(t₄-t₃)(两个时间段内的沉降相加)；

由于第1-3个InSAR干涉对分别有时间为t₁-t₂，t₁-t₃和t₂-t₄的SAR影像生成，所以就有IM₁＝1，IM₂＝1，IM₃＝2，而IS₁＝2，IS₂＝3，IS₃＝4；

把上面三个假设写成矩阵的形式即为：B·V＝ΔW

$\left(\begin{array}{ccc}(t_2-t_1)& 0& 0\\ (t_2-t_1)& (t_3-t_2)& 0\\ 0& (t_3-t_2)& (t_4-t_3)\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right)=\left(\begin{array}{c}\Delta W_1\\ {\mathrm{\Delta W}}_2\\ {\mathrm{\Delta W}}_3\end{array}\right)$

步骤4：求解地表各点在时间相邻SAR影像时刻的平均速度；

式(6)有m个观测方程和n个未知数，若可用InSAR干涉对覆盖了整个时域过程(即m≥n)，则可利用最小二乘法解出n个未知数。然而，由于不同InSAR干涉对的噪声水平不同，所以SIP方法获取的下沉观测值ΔW_k(k＝1,2,…,M)的精度也是不同的。

为了尽可能抑制SIP方法获取的差分下沉值ΔW_k中误差传播到估计的未知数中，本发明选用加权最小二乘法估计未知数V(i,j)。

考虑到相干性γ是评价InSAR干涉对获取的LOS向形变精度的重要参考指标，本发明设计了基于相干性的加权函数Q，即：

Q_k(i,j)＝[10·γ_k(i,j)]³，(k＝1,2,…,m)>

其中γ_k(i,j)为第k个InSAR干涉对在点(i,j)处的相干性，Q_k(i,j)为ΔW_k的加权值。基于该加权函数，即解出未知数的加权最小二乘解：

V(i,j)＝[B^TQ(i,j)B]^-1·[B^TQ(i,j)ΔW(i,j)]>

式中，Q(i,j)是对角线元素为Q₁(i,j),Q₂(i,j),…,Q_M(i,j)的主对角矩阵。

步骤5：估计矿区地表时序下沉；

在解出点(i,j)处的下沉速度V(i,j)后，即可利用公式(4)计算出该点的时序沉降值。对于矿区其它地表点重复以上步骤，即可估计出整个待监测矿区在SAR影像获取时刻的地表时序下沉值W＝[W₁,W₂,…,W_n]。

步骤6：计算东西、南北方向的时序水平移动；

鉴于矿区地表水平在东西E、南北方向的水平移动N与对应方向的下沉梯度之间是存在比例关系；即点(i,j)在t_k(k＝1,2,…,n)时刻的东西E(i,j)、南北方向的水平移动N(i,j)按以下公式表示：

$E_k(i,j)=\frac{b·H(i,j)}{tan\beta }·\frac{W_k(i,j+1)-W_k(i,j)}{R_E}---\left(9\right)$

$N_k(i,j)=\frac{b·H(i,j)}{tan\beta }·\frac{W_k(i+1,j)-W_k(i,j)}{R_N}---\left(10\right)$

式中，b为待监测矿区水平移动系数，H为开采采深，tanβ为主要影响角正切，R_E和R_N分别为LOS向形变图在东西和南北方向的空间分辨率。

至此，矿区地表在SAR影像获取时刻沿着东西E＝[E₁,E₂,…,E_n]、南北N＝[N₁,N₂,…,N_n]和垂直三个方向W＝[W₁,W₂,…,W_n]的三维时序形变。

本发明使用模拟数据进行可行性和精度论证，假设有一走向长壁开采工作面，其开采深度为700m，倾向宽度为150m，倾角0°，至西向东以平均速度3.3m/天推进(如图2所示，其中虚线表示各SAR影像获取时刻工作面的推进位置)。假设覆盖该矿区的可用SAR数据为6景升轨的ALOS PALSAR数据(单一雷达成像几何学)，其时间间隔均为46天，即(t₀,t₁,…,t₅)＝(0,46,92,138,184,230天)。采用快速拉格朗日分析软件FLAC^3D(Fast>0,Ws₁,Ws₂,…,Ws₅]，Es＝[Es₀,Es₁,Es₂,…,Es₅]，Ns＝[Ns₀,Ns₁,Ns₂,…,Ns₅]。由于在t₀＝0时刻，地下开采刚开始，所以Ws₀＝0，Es₀＝0，Ns₀＝0。通过分析三维时序形变后发现在该地质采矿条件下的走向和倾向水平移动系数为b＝0.3，主要影响角正切为tanβ＝1.6。

将6景SAR影像中任意两景组成InSAR干涉对，并设置时间基线阈值为180天。剔除时间基线大于时间基线阈值的InSAR干涉对，得到12对可用InSAR干涉对(干涉对构成如图4所示)。随后，按照可用InSAR干涉对的主辅影像时间，并将对应的FLAC^3D模拟的三维形变投影到LOS方向(如图4所示)，从而生成12个模拟的InSAR监测的LOS向形变。比如，第一个可用InSAR干涉对有第一景t₀和第三景t₂SAR影像组成，基于FLAC^3D模拟的三维时序形变可计算出0天到92天之间地表的差分三维形变，即ΔW₁＝Ws₁-Ws₃，ΔE₁＝Es₁-Es₃和ΔN₁＝Ns₁-Ns₃。将三维差分形变投影至LOS方向即可获得模拟的InSAR干涉对监测的LOS向形变d_LOS(t₀-t₂)，即：

$d_{LOS}(t_0-t_2)={\mathrm{\Delta Ws}}_{1}cos\theta -sin\theta \left[{\mathrm{\Delta Ns}}_{1}cos\left({\alpha }_h-\frac{3\pi }{2}\right)+{\mathrm{\Delta Es}}_1\sin \left({\alpha }_h-\frac{3\pi }{2}\right)\right]$

式中θ为雷达入射角，α_h为雷达飞行方位角，在本实施例中θ＝38.7°，α_h＝349.7°。

假设LOS向干涉图中各点的相干性相等，且均为0.6，即γ＝0.6。利用本发明的方法处理12个模拟的LOS向形变图，从而得到在6景SAR影像获取时刻矿区地表的三维时序形变W＝[W₀,W₁,…,W₅]，Es＝[E₀,E₁,…,E₅]，Ns＝[N₀,N₁,…,N₅]，其结果如图5所示。通过对比本发明所述方法获取的三维时序形变和模拟的三维时序形变后得出：两者吻合较好，且在垂直、东西和南北方向的中误差分别为0.006、0.007和0.009m。该结果表明本发明所述方法是可行的，且结果较为可靠。

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