CN103576149A - 一种基于幅度信息的地基干涉雷达三维形变提取方法 - Google Patents

Abstract

本发明公开了一种基于幅度信息的地基干涉雷达三维形变提取方法，该方法在监测区域安置两个测站，采用图像配准技术对每个测站采集的相邻两幅图像进行精密配准，可获得每个像元在四个不同方向的偏移量，其中稳定区域像元的偏移量由仪器重新安装导致，通过建立偏移量与三维形变间的变换关系，可求解三维变形量。本发明中雷达系统在非连续监测模式下工作，系统安装简便快捷；本发明通过对多幅图像进行时间滤波预处理，消弱了斑点噪声的影响；本发明基于幅度信息求解三维形变，不受大气扰动影响，且不需要进行相位解缠；本发明可获得监测区域三维形变，克服了单一视线向变形无法反映监测区域真实变形的缺点。

Description

一种基于幅度信息的地基干涉雷达三维形变提取方法

技术领域

本发明属于雷达干涉测量及变形监测领域，涉及一种基于幅度信息的地基干涉雷达三维形变提取方法。 

背景技术

合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar,InSAR）技术开辟了遥感技术用于监测地表形变的先河。利用合成孔径雷达差分干涉测量技术可提取厘米甚至毫米级的地表形变信息。地基干涉雷达是一种基于微波干涉技术的创新雷达技术，其思想来源于星载InSAR技术。该技术基于微波探测主动成像方式获取监测区域二维影像，通过合成孔径技术和步进频率技术实现雷达影像方位向和距离向的高空间分辨率，克服了星载SAR影像受时空失相干严重和时空分辨率低的缺点，通过干涉测量技术可实现毫米级微变形监测，已经广泛用于大型桥梁、滑坡、冰川和大坝的变形监测。 

使用地基雷达干涉测量技术进行变形监测，仅能获得视线向的一维变形信息，单一方向的位移场无法反映监测区域的真实形变；雷达电磁波在不均匀的对流层中传播时，会产生折射，传播路径和方向发生改变，电磁波存在延迟，即大气扰动，监测精度受到严重影响；干涉测量技术要求图像具有强相干性，系统重复安置会降低图像相干性，所以雷达系统要一直安放在现场，系统的使用效率较低。因此，使用干涉测量技术监测变形存在局限性。 

因此，需要一种新的地基干涉雷达三维形变提取方法以解决上述问题。 

发明内容

发明目的：本发明针对现有技术存在的局限，提供一种基于幅度信息的地基干涉雷达三维形变提取方法。 

技术方案：为解决上述技术问题，本发明的基于幅度信息的地基干涉雷达三维形变提取方法采用如下技术方案： 

一种基于幅度信息的地基干涉雷达三维形变提取方法，包括以下步骤： 

1）、选择监测区域，并在所述监测区域内设置观测墩A和观测墩B，求解所述观测墩A和观测墩B之间的夹角； 

2）、在所述监测区域内均匀设置若干角反射器，所述观测墩A上设置有第一地基干涉雷 达，非连续监测模式下，在t₁和t₂时刻分别采集SAR复图像

和

幅； 

3）、对同一时刻采集的多幅图像进行时间滤波预处理，得到非相干平均图像A_incoh，其中，非相干平均图像A_incoh通过下式计算得到： 

$A_{\mathrm{incoh}}(x,y)=\frac{\underset{k=0}{\overset{N_A}{\mathrm{\Σ}}}{A}_k(x,y)}{M_A}$

式中，(x,y)为像元坐标，A_k(x,y)为第k幅图像像元(x,y)幅度值，A_incoh(x,y)为处理后像元(x,y)幅度值，对观测墩A在t₁和t₂时刻采集的多幅图像进行预处理后，得到主图像M_A和辅图像S_A； 

4）、在主图像M_A上选取Q个像元，在这些像元上进行配准，得到Q对偏移量，Q对偏移量中一部分位于稳定区域，其余位于变形区域，每一对偏移量包括距离向和方位向两个方向的位移(R_ro,R_ao)，主图像M_A和辅图像S_A的配准转换模型可用下式表示： 

$\left\{\begin{array}{c}{x}_S=x_M+R_{\mathrm{ro}}\\ y_S=y_M+R_{\mathrm{ao}}\end{array}\right.$

式中，(x_M,y_M)为主图像上的像元坐标，(x_S,y_S)为辅图像上的像元坐标； 

5）、利用稳定区域的偏移量

及其像元坐标(x,y)构建如下双二次多项式模型： 

$\left\{\begin{array}{c}{R}_{\mathrm{ro}}^0=a_0+a_{1}x+a_{2}y+a_3\mathrm{xy}+a_4{x}^2+a_5{y}^2\\ R_{\mathrm{ao}}^0=b_0+b_{1}x+b_{2}y+b_3\mathrm{xy}+b_4{x}^2+b_5{y}^2\end{array}\right.$

式中，

为距离向偏移量，

为方位向偏移量，a_i(i＝0,1,…,5)和b_i(i＝0,1,…,5)为多项式系数，采用最小二乘准则进行求解系数；利用双二次多项式模型改正变形区域像元偏移量,得到完全由变形引起的距离向和方位向位移量

6）、以所述角反射器为公共点，通过所述公共点，将幅度图像投影到外部DEM上，结合观测墩A的高程，求得选取的像元与观测墩A之间的高程差，利用反正弦函数可求得第一地基干涉雷达入射角θ₁； 

7）、所述观测墩B上设置有第二地基干涉雷达，非连续监测模式下，在t₁和t₂时刻,观测墩B上分别采集SAR复图像

和

幅，采用与步骤3）-6）相同的方法，得到完全由变形引起的距离向和方位向位移量

和第二地基干涉雷达入射角θ₂；以所述角反射器为公共点建立观测墩A和观测墩B之间的图像转换关系,把观测墩A和观测墩B中的像元一一对应，每个像元有四个不同方向的位移量：距离向

和

方位向

和  $R_{\mathrm{ao}}^B;$

8）、设(r_x,r_y,r_z)为像元的三维形变量，步骤7）中四个不同方向的位移量：距离向

和 

方位向

和

根据雷达扫描成像几何关系及观测墩A和观测墩B之间的几何关系，构造下面的误差方程： 

δ＝K·r-R 

式中， $R={[R_{\mathrm{ro}}^A,R_{\mathrm{ro}}^B,R_{\mathrm{ao}}^A,R_{\mathrm{ao}}^B]}^T,$ r＝[r_x,r_y,r_z]^T， $\mathrm{\δ}={[{\mathrm{\δ}}_{\mathrm{ro}}^A,{\mathrm{\δ}}_{\mathrm{ro}}^B,{\mathrm{\δ}}_{\mathrm{ao}}^A,{\mathrm{\δ}}_{\mathrm{ao}}^B]}^T$ 为相应的观测误差，K为投影矩阵： 

$K=\left[\begin{array}{ccc}\cos (\frac{\mathrm{\π}}{2}-{\mathrm{\α}}_1)& \cos {\mathrm{\α}}_1\sin {\mathrm{\θ}}_1& \cos {\mathrm{\θ}}_1\\ -\cos (\frac{\mathrm{\π}}{2}+{\mathrm{\α}}_2-\mathrm{\β})\sin {\mathrm{\θ}}_2& \cos (\mathrm{\β}-{\mathrm{\α}}_2)\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\\ 1& 0& 0\\ \cos \left(\mathrm{\θ}\right)& \cos (\frac{\mathrm{\π}}{2}-\mathrm{\β})& 0\end{array}\right]$

其中，θ₁为第一地基干涉雷达入射角，θ₂为第二地基干涉雷达入射角，α₁为像元斜距在观测墩A的第一坐标系中xy平面内的投影与y轴的夹角，α₂为像元斜距在观测墩B的第二坐标系中xy平面内的投影与y轴的夹角，β为第一坐标系的x轴与第二坐标系的x轴之间的夹角，β与观测墩A和观测墩B之间的夹角相同， 

根据最小二乘法δ^TPδ＝min，权阵P为单位矩阵，通过下式计算得到三维形变： 

r＝(K^TPK)^-1K^TPR。 

更进一步的，步骤1）中，所述监测区域为长期缓慢变形区域，所述观测墩A和观 测墩B的上表面均位于同一水平面上，所述观测墩A和观测墩B之间的夹角大于90度。 

更进一步的，步骤4）中所述Q个像元包括角反射器、裸露的岩石或人工构筑物。 

更进一步的，步骤8）中观测墩A的第一坐标系的x轴为所述第一地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面；观测墩B的第二坐标系的x轴为所述第二地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面。 

更进一步的，步骤2）和7）中每次采集的SAR复图像为15-20幅。 

更进一步的，步骤2）中每次观测结束后均拆下所述第一地基干涉雷达；步骤7）中每次观测结束后均拆下所述第二地基干涉雷达。每次观测结束后，都将地基干涉雷达拆下，这样可以利用有一个地基干涉雷达，实现本发明的监测目的，可以有效节省成本。 

有益效果：本发明中雷达系统在非连续监测模式下工作，系统安装简便快捷；本发明通过对多幅图像进行时间滤波预处理，消弱了斑点噪声的影响；本发明基于幅度信息求解三维形变，不受大气扰动影响，且不需要进行相位解缠；本发明可获得监测区域三维形变，克服了单一视线向变形无法反映监测区域真实变形的缺点。 

附图说明

图1是本发明测站坐标系与目标点几何关系图； 

图2是本发明基于幅度信息的偏移量法流程图； 

图3是本发明三维形变求解流程图。 

具体实施方式

根据下述实施例，可以更好地理解本发明。然而，本领域的技术人员容易理解，实施例所描述的内容仅用于说明本发明，而不应当也不会限制权利要求书中所详细描述的本发明。 

请参阅图1、图2和图3所示，本发明基于幅度信息的地基干涉雷达三维形变提取方法，包括以下步骤： 

1）、选择监测区域，并在监测区域内设置观测墩A和观测墩B，求解观测墩A和观测墩B之间的夹角；其中，监测区域为长期缓慢变形区域，观测墩A和观测墩B的上表面均位于同一水平面上，观测墩A和观测墩B之间的夹角大于90度，利用全站仪测量观测墩A和观测墩B的两端坐标，计算出该夹角，如图1所示。 

2）、在监测区域内均匀设置若干角反射器，观测墩A上设置有第一地基干涉雷达，非连续监测模式下，在t₁和t₂时刻分别采集SAR复图像

和

幅；每次观测结束后，均拆下第一地基干涉雷达； 

3）、对同一时刻采集的多幅图像进行时间滤波预处理，得到非相干平均图像A_incoh，其中，非相干平均图像A_incoh通过下式计算得到： 

$A_{\mathrm{incoh}}(x,y)=\frac{\underset{k=0}{\overset{N_A}{\mathrm{\Σ}}}{A}_k(x,y)}{M_A}$

式中，(x,y)为像元坐标，A_k(x,y)为第k幅图像像元(x,y)幅度值，A_incoh(x,y)为处理后像元(x,y)幅度值，对观测墩A在t₁和t₂时刻采集的多幅图像进行预处理后，得到主图像M_A和辅图像S_A； 

4）、在主图像M_A上选取Q个像元，包括角反射器、裸露的岩石或人工构筑物等。在这些像元上进行配准，得到Q对偏移量，Q对偏移量中一部分位于稳定区域，其余位于变形区域，每一对偏移量包括距离向和方位向两个方向的位移(R_ro,R_ao)，主图像M_A和辅图像S_A的配准转换模型可用下式表示： 

$\left\{\begin{array}{c}{x}_S=x_M+R_{\mathrm{ro}}\\ y_S=y_M+R_{\mathrm{ao}}\end{array}\right.$

式中，(x_M,y_M)为主图像上的像元坐标，(x_S,y_S)为辅图像上的像元坐标； 

5）、利用稳定区域的偏移量及其像元坐标(x,y)构建如下双二次多项式模型： 

$\left\{\begin{array}{c}{R}_{\mathrm{ro}}^0=a_0+a_{1}x+a_{2}y+a_3\mathrm{xy}+a_4{x}^2+a_5{y}^2\\ R_{\mathrm{ao}}^0=b_0+b_{1}x+b_{2}y+b_3\mathrm{xy}+b_4{x}^2+b_5{y}^2\end{array}\right.$

式中，为距离向偏移量，

为方位向偏移量，a_i(i＝0,1,…,5)和b_i(i＝0,1,…,5)为多项式系数，采用最小二乘准则进行求解系数；利用双二次多项式模型改正变形区域像元偏移量,可消除其中由仪器重新安装导致的偏移分量,得到完全由变形引起的距离向 和方位向位移量

6）、以角反射器为公共点，通过公共点，将幅度图像投影到外部DEM上，结合观测墩A的高程，可求得选取的像元与观测墩A之间的高程差，由于像元斜距与高程差位于同一三角形，利用反正弦函数可求得第一地基干涉雷达入射角θ₁； 

7）、观测墩B上设置有第二地基干涉雷达，非连续监测模式下，在t₁和t₂时刻,观测墩B上分别采集SAR复图像

和

幅，每次观测完毕，拆下第二地基干涉雷达。每次采集的SAR复图像为15-20幅。采用与步骤3）-6）相同的方法，得到完全由变形引起的距离向和方位向位移量

和第二地基干涉雷达入射角θ₂；以角反射器为公共点建立观测墩A和观测墩B之间的图像转换关系,把观测墩A和观测墩B中的像元一一对应，每个像元有四个不同方向的位移量：距离向

和

方位向

和

8）、设(r_x,r_y,r_z)为像元的三维形变量，其中，观测墩A的第一坐标系的x轴为第一地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面；观测墩B的第二坐标系的x轴为第二地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面。步骤7）中四个不同方向的位移量：距离向

和

方位向

和

根据雷达扫描成像几何关系及观测墩A和观测墩B之间的几何关系（如图1所示），构造下面的误差方程： 

δ＝K·r-R 

式中， $R={[R_{\mathrm{ro}}^A,R_{\mathrm{ro}}^B,R_{\mathrm{ao}}^A,R_{\mathrm{ao}}^B]}^T,$ r＝[r_x,r_y,r_z]^T， $\mathrm{\δ}={[{\mathrm{\δ}}_{\mathrm{ro}}^A,{\mathrm{\δ}}_{\mathrm{ro}}^B,{\mathrm{\δ}}_{\mathrm{ao}}^A,{\mathrm{\δ}}_{\mathrm{ao}}^B]}^T$ 为相应的观测误差，K为投影矩阵： 

$K=\left[\begin{array}{ccc}\cos (\frac{\mathrm{\π}}{2}-{\mathrm{\α}}_1)& \cos {\mathrm{\α}}_1\sin {\mathrm{\θ}}_1& \cos {\mathrm{\θ}}_1\\ -\cos (\frac{\mathrm{\π}}{2}+{\mathrm{\α}}_2-\mathrm{\β})\sin {\mathrm{\θ}}_2& \cos (\mathrm{\β}-{\mathrm{\α}}_2)\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\\ 1& 0& 0\\ \cos \left(\mathrm{\θ}\right)& \cos (\frac{\mathrm{\π}}{2}-\mathrm{\β})& 0\end{array}\right]$

其中，θ₁为第一地基干涉雷达入射角，θ₂为第二地基干涉雷达入射角，α₁为像元 斜距在观测墩A的第一坐标系中xy平面内的投影与y轴的夹角，α₂为像元斜距在观测墩B的第二坐标系中xy平面内的投影与y轴的夹角，可从雷达极坐标图像中读取。β为第一坐标系的x轴与第二坐标系的x轴之间的夹角，即为观测墩A和观测墩B之间的夹角。 

根据最小二乘法δ^TPδ＝min，权阵P为单位矩阵，通过下式计算得到三维形变： 

r＝(K^TPK)^-1K^TPR。 

本发明中雷达系统在非连续监测模式下工作，系统安装简便快捷；本发明通过对多幅图像进行时间滤波预处理，消弱了斑点噪声的影响；本发明基于幅度信息求解三维形变，不受大气扰动影响，且不需要进行相位解缠；本发明可获得监测区域三维形变，克服了单一视线向变形无法反映监测区域真实变形的缺点。 

Claims (6)

1.一种基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，包括以下步骤：1）、选择监测区域，并在所述监测区域内设置观测墩A和观测墩B，求解所述观测墩A和观测墩B之间的夹角；

2）、在所述监测区域内均匀设置若干角反射器，所述观测墩A上设置有第一地基干涉雷达，非连续监测模式下，在t₁和t₂时刻分别采集SAR复图像

和幅；

3）、对同一时刻采集的多幅图像进行时间滤波预处理，得到非相干平均图像A_incoh，其中，非相干平均图像A_incoh通过下式计算得到：

$A_{\mathrm{incoh}}(x,y)=\frac{\underset{k=0}{\overset{N_A}{\mathrm{\Σ}}}{A}_k(x,y)}{M_A}$

式中，(x,y)为像元坐标，A_k(x,y)为第k幅图像像元(x,y)幅度值，A_incoh(x,y)为处理后像元(x,y)幅度值，对观测墩A在t₁和t₂时刻采集的多幅图像进行预处理后，得到主图像M_A和辅图像S_A；

4）、在主图像M_A上选取Q个像元，在这些像元上进行配准，得到Q对偏移量，Q对偏移量中一部分位于稳定区域，其余位于变形区域，每一对偏移量包括距离向和方位向两个方向的位移(R_ro,R_ao)，主图像M_A和辅图像S_A的配准转换模型可用下式表示：

$\left\{\begin{array}{c}{x}_S=x_M+R_{\mathrm{ro}}\\ y_S=y_M+R_{\mathrm{ao}}\end{array}\right.$

式中，(x_M,y_M)为主图像上的像元坐标，(x_S,y_S)为辅图像上的像元坐标；

5）、利用稳定区域的偏移量及其像元坐标(x,y)构建如下双二次多项式模型：

$\left\{\begin{array}{c}{R}_{\mathrm{ro}}^0=a_0+a_{1}x+a_{2}y+a_3\mathrm{xy}+a_4{x}^2+a_5{y}^2\\ R_{\mathrm{ao}}^0=b_0+b_{1}x+b_{2}y+b_3\mathrm{xy}+b_4{x}^2+b_5{y}^2\end{array}\right.$

式中，

为距离向偏移量，为方位向偏移量，a_i(i＝0,1,…,5)和b_i(i＝0,1,…,5)为多项式系数，采用最小二乘准则进行求解系数；利用双二次多项式模型改正变形区域像元偏移量,得到完全由变形引起的距离向和方位向位移量

6）、以所述角反射器为公共点，通过所述公共点，将幅度图像投影到外部DEM上，结合观测墩A的高程，求得选取的像元与观测墩A之间的高程差，利用反正弦函数可求得第一地基干涉雷达入射角θ₁；

7）、所述观测墩B上设置有第二地基干涉雷达，非连续监测模式下，在t₁和t₂时刻,观测墩B上分别采集SAR复图像

和幅，采用与步骤3）-6）相同的方法，得到完全由变形引起的距离向和方位向位移量

和第二地基干涉雷达入射角θ₂；以所述角反射器为公共点建立观测墩A和观测墩B之间的图像转换关系,把观测墩A和观测墩B中的像元一一对应，每个像元有四个不同方向的位移量：距离向和

方位向

和 $R_{\mathrm{ao}}^B;$

8）、设(r_x,r_y,r_z)为像元的三维形变量，步骤7）中四个不同方向的位移量：距离向

和

方位向

和根据雷达扫描成像几何关系及观测墩A和观测墩B之间的几何关系，构造下面的误差方程：

δ＝K·r-R

式中， $R={[R_{\mathrm{ro}}^A,R_{\mathrm{ro}}^B,R_{\mathrm{ao}}^A,R_{\mathrm{ao}}^B]}^T,$ r＝[r_x,r_y,r_z]^T， $\mathrm{\δ}={[{\mathrm{\δ}}_{\mathrm{ro}}^A,{\mathrm{\δ}}_{\mathrm{ro}}^B,{\mathrm{\δ}}_{\mathrm{ao}}^A,{\mathrm{\δ}}_{\mathrm{ao}}^B]}^T$ 为相应的观测误差，K为投影矩阵：

$K=\left[\begin{array}{ccc}\cos (\frac{\mathrm{\π}}{2}-{\mathrm{\α}}_1)& \cos {\mathrm{\α}}_1\sin {\mathrm{\θ}}_1& \cos {\mathrm{\θ}}_1\\ -\cos (\frac{\mathrm{\π}}{2}+{\mathrm{\α}}_2-\mathrm{\β})\sin {\mathrm{\θ}}_2& \cos (\mathrm{\β}-{\mathrm{\α}}_2)\sin {\mathrm{\θ}}_2& \cos {\mathrm{\θ}}_2\\ 1& 0& 0\\ \cos \left(\mathrm{\θ}\right)& \cos (\frac{\mathrm{\π}}{2}-\mathrm{\β})& 0\end{array}\right]$

其中，θ₁为第一地基干涉雷达入射角，θ₂为第二地基干涉雷达入射角，α₁为像元斜距在观测墩A的第一坐标系中xy平面内的投影与y轴的夹角，α₂为像元斜距在观测墩B的第二坐标系中xy平面内的投影与y轴的夹角，β为第一坐标系的x轴与第二坐标系的x轴之间的夹角，β与观测墩A和观测墩B之间的夹角相同，

根据最小二乘法δ^TPδ＝min，权阵P为单位矩阵，通过下式计算得到三维形变：

r＝(K^TPK)^-1K^TPR。

2.如权利要求1所述的基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，步骤1）中，所述监测区域为长期缓慢变形区域，所述观测墩A和观测墩B的上表面均位于同一水平面上，所述观测墩A和观测墩B之间的夹角大于90度。

3.如权利要求1所述的基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，步骤4）中所述Q个像元包括角反射器、裸露的岩石或人工构筑物。

4.如权利要求1所述的基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，步骤8）中观测墩A的第一坐标系的x轴为所述第一地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面；观测墩B的第二坐标系的x轴为所述第二地基干涉雷达滑轨方向，y轴与x轴垂直且位于水平面内，z轴垂直于水平面。

5.如权利要求1所述的基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，步骤2）和7）中每次采集的SAR复图像为15-20幅。

6.如权利要求1所述的基于幅度信息的地基干涉雷达三维形变提取方法，其特征在于，步骤2）中每次观测结束后均拆下所述第一地基干涉雷达；步骤7）中每次观测结束后均拆下所述第二地基干涉雷达。

Images