CN104360332A - 基于地基合成孔径雷达干涉的大气相位屏提取方法 - Google Patents

Abstract

本发明公开了一种基于地基合成孔径雷达干涉的大气相位屏提取方法，依据地基雷达差分干涉相位模型，基于最小二乘法拟合线性形变和时空滤波方法提取大气相位，该干涉相位模型充分考虑地基雷达成像机理以及差分干涉原理，从而提取大气相位屏，本发明提出振幅离差指数、相干指数以及相位稳定性三级PS点选取方法，最后经Kriging插值获得每景雷达数据的大气相位屏，准确地提取受大气影响产生的相位对雷达干涉测量精度的提高具有重要作用。通过大坝实测数据验证本方法的正确性，具有一定的实际应用价值和通用性。

Description

基于地基合成孔径雷达干涉的大气相位屏提取方法

技术领域

本发明涉及测量学领域技术领域，具体涉及一种基于地基合成孔径雷达干涉的大气相位屏提取方法。

背景技术

地基合成孔径雷达干涉(Ground Based雷达,GB雷达，简称地基雷达)技术是基于微波探测主动成像方式获取监测区域二维影像,通过合成孔径技术和步进频率技术实现雷达影像方位向和距离向的高空间分辨率,克服了星载SAR影像受时空失相干严重和时空分辨率低的缺点,通过干涉技术可实现优于毫米级微变形监测，采用GB雷达技术能精确测定被测物表面沿雷达视线向(LOS)的微量变形信息,其基本原理是：通过合成孔径雷达技术获取监测区域的二维影像,利用SF-CW技术提高雷达的距离向分辨率，通过比较影像中目标点的电磁波相位信息，采用干涉技术求取监测区域的变形量。

地基雷达进行工程监测具有以下优势：(1)利用干涉技术，地基雷达提取的形变精度高；(2)不受重访周期限制，地基雷达提取的形变具有连续性；(3)地基雷达可根据大型建筑的变形特点，选择合适的时间和位置进行观测，具有灵活性；(4)地基雷达的设站和周边环境易于控制和测量，测量误差具有可控性。但是同时还具有以下缺点，(1)地基雷达数据处理过程相对复杂；(2)地面干涉测量精度虽高，但易受到大气扰动或系统噪声影响，而其中大气扰动造成的相位影响最大，导致测量精度低，有必要对其进行提取处理。

发明内容

本发明所解决的技术问题是克服现有的地基雷达数据处理方法中存在数据处理过程相对复杂，易受到大气扰动或系统噪声影响的问题。本发明的基于地基合成孔径雷达干涉的大气相位屏提取方法，依据地基雷达差分干涉相位模型，通过最小二乘法拟合线性形变和时空滤波方法提取大气相位，该干涉相位模型充分考虑地基雷达成像机理以及差分干涉原理，从而提取大气相位屏，准确地提取受大气影响产生的相位对雷达干涉测量精度的提高具有重要作用，具有良好的应用前景。

为了解决上述目的，本发明所采用的技术方案是：

一种基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：包括以下步骤，

步骤(1)，基于地基雷达差分干涉原理，建立差分相位模型，则第n景影像的第k个像素点的差分干涉相位如公式(1)所示，

其中，为第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位、为第n景影像的第k个像素点在两次观测下大气延迟差的贡献值、为第n景影像的第k个像素点在两次观测像点的反射特性变化产生的相位贡献、为第n景影像的第k个像素点在干涉测量中，由于时间和空间去相干影响而产生的噪声以及SAR影像本身所产生的热噪声；

步骤(2)，对雷达数据进行定标后，进行幅度信息比较，将地基雷达数据进行归一化处理，得到归一化像元值r_ij，如公式(2)所示，

$r_{\mathrm{ij}}=\frac{u_{\mathrm{ij}}-\min \left\{u_{\mathrm{ij}}\right\}}{\max \left\{u_{\mathrm{ij}}\right\}-\min \left\{u_{\mathrm{ij}}\right\}}---\left(2\right)$

其中，u_ij为原始影像第i行、第j列像素的幅值，并根据公式(2)，得到归一化矩阵R，如公式(3)所示，

$R={\left(r_{\mathrm{ij}}\right)}_{m\×n}=\left(\begin{array}{ccc}{r}_{11}& \·\·\·& r_{1n}\\ \·& & \·\\ \·& & \·\\ \·& \·\·\·& \·\\ r_{m1}& \·\·\·& r_{\mathrm{mn}}\end{array}\right)---\left(3\right)$

其中，m为矩阵R的行数，n为矩阵R的列数；

步骤(3)，利用最小二乘法对连续观测的地基雷达数据进行配准，如公式(4)所示，

$\left\{\begin{array}{c}v=M_{x,y}-S_{x_2,y_2}\\ x_2=a_0-a_1{x}_1+a_2{y}_1\\ y_2=b_0-b_1{x}_1+b_2{y}_1\end{array}\right.---\left(4\right)$

其中，M_x,y为主影像像元(x,y)强度值、为从影像像元(x₂,y₂)强度值、a₀，a₁，a₂，b₀，b₁，b₂为几何畸变参数，利用强度差的平方和最小原理进行配准，即Σvv＝min，其中，v为主从影像强度值的较差；

步骤(4)，将配准后的地基雷达进行差分干涉处理，如公式(5)所示，

其中，为主影像对应像点的相位值、为从影像对应像点的相位值，式中由SAR复数影像获取相位公式为Im(u)和Re(u)分别为雷达观测数据u的虚部和实部；

步骤(5)，通过振幅离差指数法，提取永久散射体候选点D_A，振幅离差的数学表达式，如公式(6)所示，

$D_A=\frac{{\mathrm{\σ}}_A}{{\mathrm{\μ}}_A}---\left(6\right)$

其中，μ_A、σ_A分别为振幅的均值和标准偏差，若D_A的值不大于0.25，作为提取初始的永久散射体点目标的阈值；

步骤(6)，将第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位根据公式(7)，分为线性视线向形变相位和非线性视线向形变相位的之和

其中，l为雷达波长，t为雷达干涉时间间隔,结合公式(1)，则第n景影像的第k个像素点的差分干涉相位表示为公式(8)，

其中，为残余相位；

步骤(7)，通过时间域低通滤波去除噪声，去除估算线性形变后的残余相位为各差分干涉残余相位均值为作为主大气相位值，由于噪声在时间序列上是一个高频信号，因此，将为主大气相位从残余相位值中去除减掉在时间序列上通过低通滤波的剥离噪声；

步骤(8)，将在空间上进行平滑结果，作为时间T上获取的从影像的大气相位APS_slave，再将APS_slave加上主大气相位值得到对应每个差分干涉图对应的大气相位；

步骤(9)，当从差分干涉相位中去除大气相位屏后，根据公式(9)对影像上的每个像素点进行时序分析，

来重新计算时序影像上每个点的整体相干性，通过不大于阈值0.8以选择最终的永久散射体PS点；

步骤(10)，利用每一个最终的永久散射体PS点的大气相位，通过Kriging插值法估计出对应于每一景雷达图像中每一个像元的大气相位，插值公式如公式(10)所示，

$\mathrm{APS}\left(s_0\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{APS}\left(s_i\right),\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i=1---\left(10\right)$

其中，APS(s_i)为第i个永久散射体PS点的大气相位、s₀为未知像点位置、s_i为已知像点位置、M为估算大气相位值所需的永久散射体PS的点个数、λ_i为距离反比加权函数，表达式为d_i为待估像点与永久散射体PS点的距离值。

前述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(2)将地基雷达数据进行归一化处理，其中地基雷达数据的影像取景数量不少于30景。

前述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(7)，通过时间域低通滤波去除噪声，采用均值滤波器。

前述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(5)，通过振幅离差指数法，D_A等于0.25，作为提取永久散射体候选点目标的阈值。

本发明的有益效果是：本发明的基于地基合成孔径雷达干涉的大气相位屏提取方法，依据地基雷达差分干涉相位模型，基于最小二乘法拟合线性形变和时空滤波方法提取大气相位，该干涉相位模型充分考虑地基雷达成像机理以及差分干涉原理，从而提取大气相位屏，本发明提出振幅离差指数、相干指数以及相位稳定性三级PS点选取方法，最后经Kriging插值获得每景雷达数据的大气相位屏，准确地提取受大气影响产生的相位对雷达干涉测量精度的提高具有重要作用，通过大坝实测数据验证本方法的正确性，具有一定的实际应用价值和通用性。

附图说明

图1是本发明的基于地基合成孔径雷达干涉的大气相位屏提取方法的流程图。

图2是本发明的相位离差的数学模拟结果的示意图。

图3是本发明的一实施例坝顶点P1的大气改正前后的位移时间序列图。

图4是本发明的一实施例坝底点P2的大气改正前后的位移时间序列图。

具体实施方式

下面将结合说明书附图，对本发明作进一步的说明。

本发明的基于地基合成孔径雷达干涉的大气相位屏提取方法，依据地基雷达差分干涉相位模型，基于最小二乘法拟合线性形变和时空滤波方法提取大气相位，该干涉相位模型充分考虑地基雷达成像机理以及差分干涉原理，从而提取大气相位屏，本发明提出振幅离差指数、相干指数以及相位稳定性三级PS点选取方法，最后经Kriging插值获得每景雷达数据的大气相位屏，准确地提取受大气影响产生的相位对雷达干涉测量精度的提高具有重要作用，如图1所示，具体包括以下步骤，

步骤(1)，基于地基雷达差分干涉原理，建立差分相位模型，则第n景影像的第k个像素点的差分干涉相位如公式(1)所示，

其中，为第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位、为第n景影像的第k个像素点在两次观测下大气延迟差的贡献值、为第n景影像的第k个像素点在两次观测像点的反射特性变化产生的相位贡献、为第n景影像的第k个像素点在干涉测量中，由于时间和空间去相干影响而产生的噪声以及SAR影像本身所产生的热噪声，实验中地基雷达的参数设置如下表1所示。

表1 地基雷达的参数设置表

步骤(2)，对雷达数据进行定标后，进行幅度信息比较，将地基雷达数据进行归一化处理，这里的地基雷达数据，需要连续取景不少于30景，得到归一化像元值r_ij，如公式(2)所示，

$r_{\mathrm{ij}}=\frac{u_{\mathrm{ij}}-\min \left\{u_{\mathrm{ij}}\right\}}{\max \left\{u_{\mathrm{ij}}\right\}-\min \left\{u_{\mathrm{ij}}\right\}}---\left(2\right)$

其中，u_ij为原始影像第i行、第j列像素的幅值，并根据公式(2)，得到归一化矩阵R，如公式(3)所示，

$R={\left(r_{\mathrm{ij}}\right)}_{m\×n}=\left(\begin{array}{ccc}{r}_{11}& \·\·\·& r_{1n}\\ \·& & \·\\ \·& & \·\\ \·& \·\·\·& \·\\ r_{m1}& \·\·\·& r_{\mathrm{mn}}\end{array}\right)---\left(3\right)$

其中，m为矩阵R的行数，n为矩阵R的列数；

步骤(3)，利用最小二乘法对连续观测的地基雷达数据进行配准，如公式(4)所示，

$\left\{\begin{array}{c}v=M_{x,y}-S_{x_2,y_2}\\ x_2=a_0-a_1{x}_1+a_2{y}_1\\ y_2=b_0-b_1{x}_1+b_2{y}_1\end{array}\right.---\left(4\right)$

其中，M_x,y为主影像像元(x,y)强度值、为从影像像元(x₂,y₂)强度值、a₀，a₁，a₂，b₀，b₁，b₂为几何畸变参数，利用强度差的平方和最小原理进行配准，即Σvv＝min，这里的v为主从影像强度值的较差；

步骤(4)，将配准后的地基雷达进行差分干涉处理，如公式(5)所示，

其中，为主影像对应像点的相位值、为从影像对应像点的相位值，式中由SAR复数影像获取相位公式为Im(u)和Re(u)分别为雷达观测数据u的虚部和实部；

步骤(5)，通过振幅离差指数法，提取永久散射体候选点D_A，振幅离差的数学表达式，如公式(6)所示，

$D_A=\frac{{\mathrm{\σ}}_A}{{\mathrm{\μ}}_A}---\left(6\right)$

其中，μ_A、σ_A分别为振幅的均值和标准偏差，若D_A的值不大于0.25，作为提取初始的永久散射体点目标的阈值，振幅离差方法是PS技术中的提取相干点目标的方法中的经典方法，它的主要思想是利用振幅离差和相位标准偏差统计特性上的关系，当即在高信噪比g/σ_n>4的条件下，相位标准偏差近似等于振幅离差D_A，如图2所示，D_A等于0.25，作为提取永久散射体候选点目标的阈值，

步骤(6)，将第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位根据公式(7)，分为线性视线向形变相位和非线性视线向形变相位的之和

其中，l为雷达波长，t为雷达干涉时间间隔,结合公式(1)，则第n景影像的第k个像素点的差分干涉相位表示为公式(8)，

其中，为残余相位；

步骤(7)，通过时间域低通滤波去除噪声，去除估算线性形变后的残余相位为各差分干涉残余相位均值为作为主大气相位值，由于噪声在时间序列上是一个高频信号，因此，将(减号，是将主大气相位从残余相位值中去除减掉)在时间序列上通过低通滤波剥离噪声；

步骤(8)，将在空间上进行平滑结果作为时间T上获取的从影像的大气相位APS_slave，然后再将APS_slave加上主大气相位值(这个就是残余相位均值，也被估计为主大气相位值)，得到对应每个差分干涉图对应的大气相位；

步骤(9)，当从差分干涉相位中去除大气相位屏后，根据公式(9)对影像上的每个像素点进行时序分析，

来重新计算时序影像上每个点的整体相干性，通过不大于阈值0.8以选择最终的永久散射体PS点；

步骤(10)，利用每一个最终的永久散射体PS点的大气相位，通过Kriging插值法估计出对应于每一景雷达图像中每一个像元的大气相位，插值公式如公式(10)所示，

$\mathrm{APS}\left(s_0\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{APS}\left(s_i\right)\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i=1---\left(10\right)$

其中，APS(s_i)为第i个永久散射体PS点的大气相位、s₀为未知像点位置、s_i为已知像点位置、M为估算大气相位值所需的永久散射体PS点个数、λ_i为距离反比加权函数，表达式为d_i为待估像点与永久散射体PS点的距离值。

本发明的利用matlab程序处理，经过以上步骤对地基雷达进行大气相位提取后，可以对雷达视线向距离进行改正，进而得到更为准确的观测值，下面以具体实施例，经过永久散射体的选取后，取坝顶点P1及坝底点P2两点进行位移时间序列分析如下，如图3所示，坝顶点P1的大气改正前后的位移时间序列图，如图4所示，坝底点P2大气改正前后的位移时间序列图，大气的精确相位无法用实测数据直接获取，但由于测量误差的主值为大气扰动引起的相位值，因此，可以通过大坝实测的垂线数据来验证大气相位提取的精度，地基雷达数据经过大气改正后得到的视线向位移同大坝垂线数据的较差，明显地要小于未经改正的结果，从而也验证本方法所提取大气相位是准确的，具有一定的实用性。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (4)

1.基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：包括以下步骤，

步骤(1)，基于地基雷达差分干涉原理，建立差分相位模型，则第n景影像的第k个像素点的差分干涉相位如公式(1)所示，

其中，为第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位、为第n景影像的第k个像素点在两次观测下大气延迟差的贡献值、为第n景影像的第k个像素点在两次观测像点的反射特性变化产生的相位贡献、为第n景影像的第k个像素点在干涉测量中，由于时间和空间去相干影响而产生的噪声以及SAR影像本身所产生的热噪声；

步骤(2)，对雷达数据进行定标后，进行幅度信息比较，将地基雷达数据进行归一化处理，得到归一化像元值r_ij，如公式(2)所示，

$r_{\mathrm{ij}}=\frac{u_{\mathrm{ij}}-\min \left\{u_{\mathrm{ij}}\right\}}{\max \left\{u_{\mathrm{ij}}\right\}-\min \left\{u_{\mathrm{ij}}\right\}}---\left(2\right)$

其中，u_ij为原始影像第i行、第j列像素的幅值，并根据公式(2)，得到归一化矩阵R，如公式(3)所示，

$R={\left(r_{\mathrm{ij}}\right)}_{m\×n}=\left(\begin{array}{ccc}{r}_{11}& ...& r_{1n}\\ .& & .\\ .& & .\\ .& ...& .\\ r_{m1}& ...& r_{\mathrm{mn}}\end{array}\right)---\left(3\right)$

其中，m为矩阵R的行数，n为矩阵R的列数；

步骤(3)，利用最小二乘法对连续观测的地基雷达数据进行配准，如公式(4)所示，

$\left\{\begin{array}{cc}v=M_{x,y}-S_{x_2,y_2}& \\ x_2=a_0-a_1{x}_1+a_2{y}_1& \\ y_2=b_0-b_1{x}_1+b_2{y}_1& \end{array}\right.---\left(4\right)$

其中，M_x,y为主影像像元(x,y)强度值、为从影像像元(x₂,y₂)强度值、a₀，a₁，a₂，b₀，b₁，b₂为几何畸变参数，利用强度差的平方和最小原理进行配准，即Σvv＝min，其中，v为主从影像强度值的较差；

步骤(4)，将配准后的地基雷达进行差分干涉处理，如公式(5)所示，

其中，为主影像对应像点的相位值、为从影像对应像点的相位值，式中由SAR复数影像获取相位公式为Im(u)和Re(u)分别为雷达观测数据u的虚部和实部；

步骤(5)，通过振幅离差指数法，提取永久散射体候选点D_A，振幅离差的数学表达式，如公式(6)所示，

$D_A=\frac{{\mathrm{\σ}}_A}{{\mathrm{\μ}}_A}---\left(6\right)$

其中，μ_A、σ_A分别为振幅的均值和标准偏差，若D_A的值不大于0.25，作为提取初始的永久散射体点目标的阈值；

步骤(6)，将第n景影像的第k个像素点两次雷达观测中像素在视线向形变相位根据公式(7)，分为线性视线向形变相位和非线性视线向形变相位的之和

其中，l为雷达波长，t为雷达干涉时间间隔,结合公式(1)，则第n景影像的第k个像素点的差分干涉相位表示为公式(8)，

其中，为残余相位；

步骤(7)，通过时间域低通滤波去除噪声，去除估算线性形变后的残余相位为各差分干涉残余相位均值为作为主大气相位值，由于噪声在时间序列上是一个高频信号，因此，将为主大气相位从残余相位值中去除减掉在时间序列上通过低通滤波的剥离噪声；

步骤(8)，将在空间上进行平滑结果，作为时间T上获取的从影像的大气相位APS_slave，再将APS_slave加上主大气相位值得到对应每个差分干涉图对应的大气相位；

步骤(9)，当从差分干涉相位中去除大气相位屏后，根据公式(9)对影像上的每个像素点进行时序分析，

来重新计算时序影像上每个点的整体相干性，通过不大于阈值0.8以选择最终的永久散射体PS点；

步骤(10)，利用每一个最终的永久散射体PS点的大气相位，通过Kriging插值法估计出对应于每一景雷达图像中每一个像元的大气相位，插值公式如公式(10)所示，

$\mathrm{APS}\left(s_0\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{APS}\left(s_i\right),\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\λ}}_i=1---\left(10\right)$

其中，APS(s_i)为第i个永久散射体PS点的大气相位、s₀为未知像点位置、s_i为已知像点位置、M为估算大气相位值所需的永久散射体PS的点个数、λ_i为距离反比加权函数，表达式为d_i为待估像点与永久散射体PS点的距离值。

2.根据权利要求1所述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(2)将地基雷达数据进行归一化处理，其中地基雷达数据的影像取景数量不少于30景。

3.根据权利要求1所述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(7)，通过时间域低通滤波去除噪声，采用均值滤波器。

4.根据权利要求1所述的基于地基合成孔径雷达干涉的大气相位屏提取方法，其特征在于：步骤(5)，通过振幅离差指数法，D_A等于0.25，作为提取永久散射体候选点目标的阈值。