CN101706577B - InSAR监测高速公路路面沉降方法 - Google Patents

Abstract

本发明涉及基于遥感影像的大地测量领域，是一种InSAR监测高速公路路面沉降的方法。该方法包括：首先，对SAR数据进行预处理、配准和干涉，得到InSAR干涉相位和幅度影像；其次，对InSAR干涉相位进行平地效应、地形效应和轨道残余趋势相位消除，得到只包含有地表形变信息的相位值；再次，通过高速公路在InSAR幅度影像中的条带特征，识别高速公路在SAR影像中的条带位置和坐标。然后利用该坐标提取干涉相位图中相应位置的相位值，并采用条带特征目标的滤波和解缠算法恢复其真实相位值；最后，对其进行地理编码和形变值转换，得到高速公路路面沉降值。本方法具有实现简单、费用低、监测精度高、监测范围大、自动化程度高等优点。

Description

InSAR监测高速公路路面沉降方法

技术领域

本发明属于基于遥感影像的大地测量领域，涉及一种利用InSAR(合成孔径雷达干涉测量)技术监测高速公路路面沉降方法。

技术背景

作为一种全天候，大范围，高精度的大地测量技术，合成孔径雷达干涉测量(InSAR-Interferometric Synthetic Aperture radar)一直是近年来国内外关注的热门之一。随着SAR数据的极大丰富，该技术获得了前所未有的发展和应用，如监测由于地下水抽取、煤矿开采、石油开采、地下铁路修建、填海等引起的地表变形。

InSAR技术最基本的原理是借助于覆盖同一个地区的两幅或两幅以上SAR图像，利用包含在SAR图像中的相位信息提取出雷达天线到地表之间的距离，进行相位干涉处理，结合雷达的姿态参数重建地表的数字高程模型(DEM)，或者是测量地表的形变(可以达到厘米到毫米级精度)等。

图1为二轨法进行差分干涉测量的几何示意图。重复轨道干涉测量的相位Δφ由以下几部分组成：

Δφ＝φ_flt+φ_topo+φ_def+φ_atm+φ_noise+k·2π         (1)

其中φ_flt为平地相位；φ_topo为地形引起的相位；φ_def为视线方向(LOS)的形变引起的相位；φ_atm为两次成像期间大气延迟不同引起的相位；φ_noise为噪声相位，k为整周模糊度。

对(1)式，平地相位φ_flt可根据基线进行估计，地形引起的相位φ_top可根据已有的DEM估计，若忽略大气相位和噪声相位的影响，则可估计出地表形变引起的相位。对形变相位进行解缠，即可得：

$D=\frac{\mathrm{\λ}}{4\mathrm{\π}}{\mathrm{\φ}}_{\mathrm{def}}---\left(2\right)$

其中λ为雷达波长。对(2)微分，可得：

$\mathrm{\ΔD}=\frac{\mathrm{\λ}}{4\mathrm{\π}}\mathrm{\Δ}{\mathrm{\φ}}_{\mathrm{def}}---\left(3\right)$

当Δφ_def＝2π时，可得形变模糊度：

$\mathrm{\ΔD}=\frac{\mathrm{\λ}}{2}---\left(4\right)$

由(4)可知：2π相位变化所对应的形变量为波长的一半。而对于工作于C波段的ERS和ENVISAT卫星，形变模糊度为2.8cm，波长相对较长的L波段的ALOS卫星PALSAR传感器也仅为11.8cm。因此，InSAR监测地表形变的灵敏度比较高，理论上可达到cm到mm级的监测精度。

InSAR技术虽然已经被研究用于监测大范围地表沉降，但目前针对高速公路路面这一条带目标沉降监测的研究很少，因为InSAR技术本身是一种基于面的处理方法，用来监测线状条带窄目标，有很多核心算法需要突破。目前，高速公路沉降监测方法主要是水准测量和GPS测量，但是它们都存在着一些缺陷：①水准点(标石点)的稳定性问题；②测量的是沉降点、线，构成沉降面还必须经过数据的内插，如果测量点过于稀疏，就无法给出整个区域的变化趋势，而如果要得到整个区域的变化趋势，就需要观测很多观测点，这又需要耗费大量的人力和物力；③需要事先预计到沉降的大致方位和范围，才能布置下一步的测量工作。

发明内容

本发明要解决的技术问题是：克服当前InSAR技术的缺陷和不足，提供一种InSAR技术针对高速公路路面的沉降监测方法，从而降低高速公路的沉降监测成本，提高监测的精度和准确性。

本发明的技术解决方案如下：

一种InSAR监测高速公路路面沉降方法，其特征在于，包括以下步骤：

1)InSAR初处理：该步骤包括主辅SAR影像【主辅SAR影像分别对应于图1中卫星在位置A₁和A₂时获取的SAR影像】的选取、预处理、配准和干涉，从而得到监测地区的干涉相位影像和幅度影像(主辅影像应该选取空间基线尽量小而时间跨度涵盖地表形变发生的过程的影像对，主辅影像的选取为现有技术中的常规技术)；干涉过程为：通过对配准后的主辅SAR影像进行共轭相乘操作，生成干涉图；

2)InSAR干涉影像的处理：该步骤包括平地相位、地形相位和由卫星轨道不精确引起的趋势相位的去除，从而获取了代表地表形变的含有噪声和整周模糊度的相位图；

平地相位利用以下公式计算：

${\mathrm{\φ}}_{\mathrm{flt}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}B\sin ({\mathrm{\θ}}^0-\mathrm{\α});$

λ为雷达波长、B为空间基线、θ⁰为雷达入射角、α为空间基线与水平方向的夹角；

地形相位利用下式计算：

${\mathrm{\φ}}_{\mathrm{topo}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}^0}{R_1\sin {\mathrm{\θ}}^0}h;$

B_⊥ ⁰为空间基线的垂直分量、R₁为主影像获取时卫星的斜距、h为地面点高程；

趋势相位φ_ramp的提取方法为：先找出只包含有趋势相位的点的区域，提取该区域的点的相位值，利用一个双二次多项式模型来进行拟合：

φ_ramp＝c₁+c₂i+c₃j+c₄ij+c₅i²+c₆j²

其中c_m(m＝1，2，...，6)为多项式系数，(i，j)为提取相位值点的影像坐标，得到多项式模型后，完成整幅SAR影像的趋势相位的拟合；即根据该多项式模型计算整幅SAR影像的趋势相位；

具体过程如下：在去除由于卫星轨道不精确引起的趋势相位时，首先通过先验信息(即通过水准或GPS测量获取的研究区域的形变信息)，将整个干涉图分为形变区和稳定区。再在整幅干涉图上均匀取若干个点(如32×32)，剔除位于形变区(即由水准或GPS监测的地表产生位移的区域)和相干性较低(＜0.7)的点，提取剩余点的相位值，此时这些点的相位值就只包含有趋势相位φ_ramp，利用上述的双二次多项式模型来进行拟合。得到多项式模型后，完成整幅SAR影像的趋势相位的拟合，并从重复轨道干涉测量的相位

中去除。

3)InSAR幅度影像的处理：该步骤利用高速公路在InSAR幅度影像的条带(线性)目标特征，自动识别高速公路在影像中所对应的像素坐标，从而在干涉相位影像中提取相应的代表高速公路形变的相位值；

假设x₀为幅度影像A中的任意一个像素，利用x₀所在的中心区(如图3中的区域1，具体来说，是指x₀所在的一个像素区，如大小为3*9像素的区域)和比邻区(如图3中的区域2和区域3，具体来说，是指上述中心区左右两边的2个小区域，如中心区左右两边各有一个3*9个像素的区域)的幅度比来计算两个边缘检测响应，其中较小的那个就是线特征检测比值r。采用最小线检测比值的阈值r_min，当r＞r_min时，则认为x₀为高速公路上的点，否则就不是高速公路上的点，阈值r_min通过虚警率确定。然后利用Hough转换将已检测出来的公路点构建成连续的高速公路【自动识别高速公路在影像中所对应的像素坐标步骤为本领域的现有技术】。

4)高速公路形变相位的滤波处理：该步骤根据高速公路的特征，对提取的高速公路形变相位值进行条带特征目标滤波，从而消除或减弱高速公路形变相位中的噪声的影响；

计算公式为：

$\hat{m}=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{w}_n{m}_n$

其中，为融合估计值，即滤波后的相位值；m_n和w_n分别为第n个线性窗口的相位平均值和权值；N为所选线性窗口个数，取值范围为1-8，具体由干涉图的相干性决定。权重w_n的计算公式为：

$w_n=\frac{1/{{\mathrm{\σ}}_n}^2}{\underset{K=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{\σ}}_k^2};n=\mathrm{1,2},...,N$

其中，σ_n ²为第n个线性窗口的方差；

5)高速公路形变相位的解缠处理：该步骤根据高速公路的特征，对滤波后的高速公路形变相位值进行条带特征目标解缠，从而恢复高速公路形变相位的真实值：

先计算高速公路上的残数点，并将高速公路的边界(如图3中区域1和区域2、3的分界线)设置为解缠的限制边界，从高速公路起始的地方需找第一个残数点作为起始点，以该点为中心3×3像素的窗口搜索下一个残差点，若找到残数点则将其与当前残数点相连接；若没找到残数点，则将窗口按高速公路方向从左至右，从上至下依次寻找，直至遍历待检测的整条高速公路；连接所有残数点，并建立枝切线；随后，沿高速公路避开枝切线进行逐点积分，则任意一点p的解缠后的相位为：

其中

为相位梯度，C为避开枝切线的解缠路径，

为解缠起始点p₀的相位。直至搜索完全部残数点，完成相位解缠，得到只含有形变信息的相位值Δφ_road。

6)形变转换和地理编码处理：利用 $D_{\mathrm{road}}=\frac{\mathrm{\λ}}{4\mathrm{\π}}\mathrm{\Δ}{\mathrm{\φ}}_{\mathrm{road}}$ 将相位值Δφ_road转换为高速公路路面沉降量D_road。然后对其进行地理编码。首先确定每个像元的笛卡尔坐标，再将时间参数和主影像的轨道信息作为地理地位的输入数据，计算每个像元的WGS 84坐标，将影像坐标系中的高速公路路面沉降量投影到国际标准地理参考系中。最后利用若干个控制点(如GPS或水准资料)将相对形变值转换为绝对形变值，即得出了WGS 84坐标系下的高速公路在主辅SAR影像获取时间之内的路面沉降值。

有益效果：

根据这一思路，可以得到如图2所示的InSAR监测高速公路路面沉降方法的实现框图。从图中可以看出，整个方法充分融合了InSAR和高速公路的特点，即将传统InSAR面目标监测拓展到了符合高速公路这一类的条带目标监测，又克服了传统高速公路沉降监测方法费时费力且无法得到整个高速公路的形变趋势的缺点。整个流程结构清晰，具有实现简单、费用低、监测精度高、监测范围大、自动化程度高等优点。

附图说明

图1是InSAR监测地表形变的基本原理图。

图2是InSAR监测高速公路路面沉降方法流程图。

图3是InSAR幅度影像中条带结构的检测模型。

具体实施方式

以下将结合图和具体实施过程对本发明做进一步详细说明：

实施例1：

本发明的具体实施方案包含以下几个步骤：

(1)数据的选择和预处理。根据监测区域坐标范围定制SAR raw影像，挑选SAR影像时要考虑如下几个因素：传感器的特性(包括波长、带宽、SNR、轨道和重复周期等)、数据的可用性、时间和空间基线、地面特征以及大气等等。随后对已选好的raw影像进行预处理，将其转成单视复数影像(SLC)。

(2)主辅SAR影像配准。SAR影像配准的精度直接影响后续的处理和最终结果的质量，要生成高质量的干涉条纹图，配准精度要优于0.2个像素。一般先利用卫星轨道信息进行粗配准，使得配准精度达到几个像素之间。然后在此基础上计算主辐SAR影像的幅度影像的相关性，按照一定分布在整幅影像上找出若干个同名点(如64×64)。再利用这些同名点确定几何变换多项式模型，该多项式模型一般取为双二次模型。最后，利用确定的几何变换模型对辅影像进行重采样。重采样的插值算法包括最近邻法、双线形插值法、三次样条插值等。在InSAR数据处理时，要特别注意需对SLC影像的实部和虚部分别进行重采样。至此就完成了主辅SAR影像的配准工作。配准步骤为现有技术。

(3)生成干涉图。SAR SLC影像u为一个m×n的矩阵，其中每一个像素点上有一个复数值，包含实部u_real和虚部u_image，则该点的幅度值 $\left|u\right|=\sqrt{{u_{\mathrm{real}}}^2+{u_{\mathrm{image}}}^2},$ 相位值通过对配准后的主辅SAR SLC影像u₁、u₂的共轭相乘，生成干涉图：

其中A为干涉图的幅度，为干涉图的相位。但是在实际操作中，两幅SLC影像在空间域共轭相乘对应于它们各自频谱的卷积，因此截止频率不一定能满足奈奎斯特条件，从而出现频谱混迭的现象，这时需要在干涉之前对SLC影像进行过采样。

(4)形变相位提取。如式(1)所示，

包括平地相位、地形相位、形变相位、大气相位、噪声以及整周模糊度，为了从中提取出形变相位，必须首先计算出剩余的相位贡献。其中平地相位利用以下公式计算：

${\mathrm{\φ}}_{\mathrm{flt}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}B\sin ({\mathrm{\θ}}^0-\mathrm{\α})$

这里λ为雷达波长、B为空间基线、θ⁰为雷达入射角、α为基线与水平方向的夹角。而地形则利用下式计算：

${\mathrm{\φ}}_{\mathrm{topo}}=-\frac{4\mathrm{\π}}{\mathrm{\λ}}\frac{B_{\⊥}^0}{R_1\sin {\mathrm{\θ}}^0}h$

这里B_⊥ ⁰为空间基线的垂直分量、R₁为主影像获取时卫星的斜距、h为地面点高程。其中地面点高程h需通过外部数据来获取，例如SRTM数据或TanDEM数据。

从上面两式可以看出，空间基线B估计的好坏直接影响到平地相位和地形相位去除精度。但是在实际的数据处理时，基线B的精确估计往往是一个难点，因为基线B的估计依赖于对卫星轨道状态的计算。而卫星轨道状态的计算需要有足够多的地面控制点，这在实际应用中是很难满足的，因此在干涉图中，往往会因为基线估计的不精确而出现趋势相位，特别是目前分辨率较高的PALSAR数据，这种趋势相位会更加明显，往往会掩盖真实的地表形变，因此需要进行剔除。本发明提出一种去除趋势相位的方法。通过先验信息(即通过水准或GPS测量获取的研究区域的形变信息)，将整个干涉图分为形变区和稳定区。再在整幅干涉图上均匀取若干个点(如32×32)，剔除位于形变区(即由水准或GPS监测的地表产生位移的区域)和相干性较低(＜0.7)的点，提取剩余点的相位值，此时这些点的相位值就只包含有趋势相位φ_ramp，利用一个双二次多项式模型来进行拟合：

φ_ramp＝c₁+c₂i+c₃j+c₄ij+c₅i²+c₆j²

这里c_i(i＝1，2，...，6)为多项式系数，(i，j)为提取相位值点的影像坐标。得到多项式模型后，完成整幅SAR影像的趋势相位的拟合，并从重复轨道干涉测量的相位

中去除。

(5)高速公路像素的提取。在遥感影像中，高速公路是一类比较特殊的地物，具有很大空间跨度的条带特征，本身就具有重要的意义。目前对于从遥感影像中提取线形特征的研究有很多，包括差分算子、Canny算子、零交叉算子等，但是这些传统的算法对光学影像可以取得很好的效果，对InSAR幅度图而言则难以奏效，这是因为幅度图中斑点(Speckle)噪声的影响。因此需要引入一种关于被检测条带结构的大范围先验知识的全局方法，以融合上述算法提取的条带特征片断，从而组成具有较大空间跨度的条带。

如图3所示，x₀为幅度影像A中的任意一个像素，μ₁为x₀所在的中心区域的幅度均值，μ₂和μ₃分别为两个比邻区域的幅度均值，则有：

${\mathrm{\μ}}_i=\frac{1}{n_i}{\mathrm{\Σ}}_{s\∈i}{A}_s$

这里A_s代表像素S的幅度值，i＝1，2，3代表图3中不同区域，n_i代表区域i的像素总数。则区域i和j之间的边缘检测响应可以定义为r_ij：

$r_{\mathrm{ij}}=1-\min (\frac{{\mathrm{\μ}}_i}{{\mathrm{\μ}}_j},\frac{{\mathrm{\μ}}_j}{{\mathrm{\μ}}_i})$

对于图3中的三个区域，则可得到两个边缘检测响应r₁₂和r₁₃，因此线特征检测比值就是这两个边缘检测响应较小的那个：

r＝min(r₁₂，r₁₃)

设定一个最小线检测比值的阈值r_min，当r＞r_min时，则认为x₀为高速公路上的点，否则就不是。按上述步骤检测不同方向窗口下的r，响应最大的窗口中所示方向即为高速公路的方向。

而r_min可通过虚警率来确定。假设斑点噪声充分发展(fully developed)，得的L视InSAR幅度影像的强度概率密度函数(pdf)：

$f_A\left(t\right|<I>)=\frac{2}{\mathrm{\&Gamma;}\left(L\right)}{\left(\frac{L}{<I>}\right)}^L{t}^{2L-1}{e}^{-({\mathrm{Lt}}^2/<I>)}$

其中t为任意像素的强度值，L为视数，I为整个灰度影像的强度值，<·>代表求均值运算，г(·)代表Gamma运算。

可推导出线检测比值的pdf：

$f_r\left(t\right|c_{12},c_{13})=\frac{4\left[(n_1+n_2)L\right]\mathrm{\&Gamma;}\left[(n_1+n_3)L\right]}{\mathrm{\&Gamma;}{\left(n_{1}L\right)}^2\mathrm{\&Gamma;}\left(n_{2}L\right)\mathrm{\&Gamma;}\left(n_{3}L\right)}\&times;{n_1}^{2n_{1}L}{n_2}^{n_{2}L}{n_3}^{n_{3}L}$

$\&times;[g\left(t\right|c_{12}){\&Integral;}_t^{1}g\left(x\right|c_{13})\mathrm{dx}+g\left(t\right|c_{13}){\&Integral;}_t^{1}g\left(x\right|c_{12})\mathrm{dx}]$

其中n₁，n₂和n₃分别对应于区域1，2和3的像素个数；幅度比 $c_{12}=\sqrt{\frac{<I_1>}{<I_2>}},$ $c_{13}=\sqrt{\frac{<I_1>}{<I_3>}};$ $g\left(x\right|c_{1i})=\frac{c_{1i}^{2n_{i}L}{(1-x)}^{2n_{1}L-1}}{{[{(1-x)}^2{n}_1+n_i{c}_{1i}^2]}^{L(n_1+n_i)}}+\frac{{\left(\frac{1}{c_{\mathrm{li}}^2}\right)}^{n_{1}L}{(1-x)}^{2n_{i}L-1}}{{[{(1-x)}^2{n}_i+\frac{n_1}{c_{1i}^2}]}^{L(n_1+n_i)}}.$ 则对应于某个阈值r_min和幅度比c₁₂、c₁₃，线检测的概率Pdf为：

$P_d(r_{\min },c_{12},c_{13})={\&Integral;}_{r_{\min }}^1{f}_r\left(t\right|c_{12},c_{13})\mathrm{dt}$

那么当确定高速公路方向后，发生错误检测的情况有两种：①c₁₂＝c₁₃＝1；②c₁₂＝1或c₁₃＝1。上述两种情况下，都可以给出如下虚警率：

$P_{\mathrm{\&phi;}}(r_{\min },c)=P_d(r_{\min },1,c)=P_d(r_{\min },c,1)={\&Integral;}_{r_{\min }}^1{f}_r\left(t\right|1,c)\mathrm{dt}$

利用Kolmogorov-Smirnov检验，可得虚警率为0.01。因此，通过上式和虚警率就可以反推确定最终的最小线形特征检测比值的阈值r_min【反推求阈值r_min部分为现有技术】。记录下所检测的高速公路点的像素坐标和高速公路方向，但是这些点都是高速公路基元，在大多数的情况下，这些基元是不连续的，需要将它们构建成连续的高速公路。

首先，对所检测到的像素中独立的点进行压缩，因为它们属于高速公路的可能性较小。利用这些独立点的高速公路方向d_k，k∈{0，...，N_d}，如果其它像素的高速公路方向接近于d_k(如d_k-1，d_k或d_k+1)，则保留这个独立的点，否则就将其删除。其次，对小块的线形或条带结构进行压缩。在邻域中确定一个正确的高速公路方向，利用一个20×20像素大小的窗口对这些小块结构进行局部Hough转换，保留符合正确高速公路方向的像素。最后是像素连接，将高速公路相近且距离不超过4个像素的点连接在一起，最终形成具有很大空间跨度的高速公路条带，记录其所在的像素坐标。

(6)高速公路的形变相位滤波。利用幅度影像上高速公路所在像素坐标，提取相位图中对应位置的相位值，即为高速公路的形变相位。值得注意的是，此时的相位值仅是真实相位的主值部分，要得到真实的相位值，需要进行相位解缠。但是在相位解缠之前，必须先对相位进行滤波以去除噪声的影响。目前，针对InSAR干涉图存在的相位噪声，目前已有多种滤波方法，包括空域滤波和频域滤波，这些方法均能在一定程度上减少干涉图的噪声，提高干涉图的质量。对于高速公路这一特殊地物而言，由于路面的地形起伏一般是非常缓变的，路面上的相位相对而言也是比较稳定的，且具有很强的相关性，而噪声则是统计独立。在这种情况下，本发明提出一种条带特征目标的相位滤波方法。针对高速公路延伸方向和对应的条纹都表现出的条带(线性)特征，选取线性形状的同态象元(相位接近)求均值，然后按照最优化融合的方法将线性窗口的均值进行融合达到滤波的效果。最优化融合是一种多传感器数据融合技术，该技术将多元数据按照方差倒数的权重进行加权平均，实现数据的融合。最优化融合方法使本技术得到的滤波最终结果是最优的(方差最小)，因此经过融合后的结果具有最好的相位平滑度。

考虑到噪声分布不均匀，根据相干性的大小来选择线性窗口的数量，在相干性高的区域其噪声较弱，此时采用较少的线性窗口进行融合滤波，减少信息损失；在相干性低的区域其噪声较强，选择较多线性窗口进行融合滤波，增强去噪能力。根据相干值大小确定融合窗口的数量，象元方差大小选择线性窗口，然后对所选线性滤波窗口的均值以方差倒数为权重，求线性窗口的加权平均值来代替像元的值。其计算公式为：

$\hat{m}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_i{m}_i$

其中，

为融合估计值，即滤波后的相位值；m_i和w_i分别为第i个线性窗口的相位平均值和权值；N为所选线性窗口个数，取值范围为1-8，具体由干涉图的相干性决定。权重w_i计算公式为：

$w_i=\frac{1/{{\mathrm{\&sigma;}}_i}^2}{\underset{K=1}{\overset{N}{\mathrm{\&Sigma;}}}1/{\mathrm{\&sigma;}}_k^2},i=1,2,...,N$

其中，σ_i ²为第i个线性窗口的方差，N为所选线性窗口个数。

在对线状或条带地物进行滤波时，同时考虑线性窗口和相干性使该滤波方法具有很好的保真性和自适应性。能够在抑制噪声的同时极大程度的减少信息的丢失，较好的保持干涉图的边缘和细节信息。

(7)高速公路的形变相位解缠。正如前面所说，目前的相位值只是真实相位的主值部分，其取值范围在(-π，+π)之间，而真实的相位是这个值与2π的整数倍之和。从理论上来说，提取相位的偏导数，利用路径积分就可以反演出真实相位值，但在一般的InSAR数据处理中，噪声等因素往往使得相位解缠变得困难，尤其对于高速公路等条带式地物，因为可供选择的积分路径比大范围干涉图少。

考虑到本次解缠的对象是在高速公路条带上的相位，地表起伏缓慢，而且已经进行了符合高速公路条带特征的滤波处理，本发明提出一种适用于条带或线性特征目标的相位解缠方法。首先计算高速公路上的残数点。由于在沿着二维的某一区域的闭合路径对缠绕的相位差进行积分时，沿不同的路径的积分结果不一致的现象，而这些现象只发生于一些孤立的点或小区域，这些点就是残数点，在解缠的时候需要避开这些残数点。然后将高速公路的边界(如图3中区域1和区域2、3的分界线)设置为解缠的限制边界。从高速公路起始的地方需找第一个残数点作为起始点，以该点为中心3×3像素的窗口搜索下一个残数点。若找到残数点则将其与当前残数点相连接。若没找到残数点，则将窗口按高速公路方向从左至右，从上至下依次寻找，直至到达高速公路边界。连接所有残数点，并建立枝切线。随后，沿高速公路避开枝切线进行逐点积分，则任意一点p的解缠后的相位为：

其中

为相位梯度，C为避开枝切线的解缠路径，

为解缠起始点p₀的相位。直至搜索完全部残数点，完成相位解缠。值得注意的是，在高速公路形变微小的特殊情况下(即形变量小于波长的一半)，相位值不会超出(-π，+π)的范围，可不进行解缠。

(8)形变转换及地理编码。完成上述工作后，得到了只含有形变信息的相位值Δφ_road，则可以利用下面公式将其转换为高速公路路面沉降量D_road：

$D_{\mathrm{road}}=\frac{\mathrm{\&lambda;}}{4\mathrm{\&pi;}}\mathrm{\&Delta;}{\mathrm{\&phi;}}_{\mathrm{road}}$

此时的形变量是在影像坐标系下的相对形变值，为了得到最终用户可以使用的产品，必须对其进行地理编码，即将影像坐标系中的高速公路路面沉降量投影到国际标准地理参考系中(如WGS 84坐标系)，使其能够和其它在同一参考系中的地理图件进行比较。

首先确定每个像元的笛卡尔坐标，再将时间参数和主影像的轨道信息作为地理地位的输入数据，计算每个像元的WGS 84坐标。此外，为了提高地理编码的精度，必须利用外部的SRTM数据或GPS数据提供的坐标进行校正。由于所得的WGS 84坐标一般是不规则分布的，最后将其插值到一个规则网格中。最后利用若干个控制点(如GPS或水准资料)将相对形变值转换为绝对形变值。

Claims (3)

1.一种InSAR监测高速公路路面沉降方法，其特征在于，包括以下步骤：

1)InSAR初处理：该步骤包括主辅SAR影像的选取、预处理、配准和干涉，从而得到监测地区的干涉相位影像和幅度影像；干涉过程为：通过对配准后的主辅SAR影像进行共轭相乘操作，生成干涉图；

2)InSAR干涉相位影像的处理：该步骤包括平地相位、地形相位和由卫星轨道不精确引起的趋势相位的去除，从而获取了代表地表形变的含有噪声和整周模糊度的相位图；

平地相位利用以下公式计算：

λ为雷达波长、B为空间基线、θ°为雷达入射角、α为空间基线与水平方向的夹角；

地形相位利用下式计算：

为空间基线的垂直分量、R₁为主影像获取时卫星的斜距、h为地面点高程；

趋势相位利用下式计算：先找出只包含有趋势相位的点的区域，提取该区域的点的相位值，利用一个双二次多项式模型来进行拟合：

φ_ramp＝c₁+c₂i+c₃j+c₄ij+c₅i²+c₆j²

其中(i，j)为提取相位值点的干涉相位影像坐标，c_m为多项式系数，m＝1，2，...，6，其数值通过上述的拟合过程来确定，得到多项式模型后，完成整幅SAR干涉相位影像的趋势相位的拟合；即根据该多项式模型计算整幅SAR干涉相位影像的趋势相位；

3)InSAR幅度影像的处理：该步骤利用高速公路在InSAR幅度影像的条带目标特征，自动识别高速公路在干涉相位影像中所对应的像素坐标，从而在干涉相位影像中提取相应的代表高速公路形变的相位值；

定义x₀为幅度影像A中的任意一个像素，利用x₀所在的中心区和比邻区的幅度比来计算两个边缘检测响应，其中较小的那个就是线特征检测比值r，采用最小线特征检测比值的阈值r_min，当r＞r_min时，则认为x₀为高速公路上的点，否则就不是高速公路上的点，其中阈值r_min通过虚警率确定；然后利用Hough转换将已检测出来的公路上的点构建成连续的高速公路；

4)高速公路形变相位的滤波处理：该步骤根据高速公路的特征，对提取的高速公路形变相位值进行条带特征目标滤波，从而消除或减弱高速公路形变相位中的噪声的影响；

5)高速公路形变相位的解缠处理：该步骤根据高速公路的特征，对滤波后的高速公路形变相位值进行条带特征目标解缠，从而恢复高速公路形变相位的真实值；

6)形变转换及地理编码处理：利用公式

将公路形变相位值Δφ_road转换为高速公路路面沉降量D_road，然后对其进行地理编码：首先确定干涉相位影像中高速公路涵盖的每个像元的笛卡尔坐标，再将时间参数和主影像的轨道信息作为地理地位的输入数据，计算干涉相位影像中高速公路涵盖的每个像元的WGS 84坐标，将影像坐标系中的高速公路路面沉降量投影到国际标准地理参考系中。

2.根据权利要求1所述的InSAR监测高速公路路面沉降方法，其特征在于：在进行高速公路形变相位滤波时，采用了以下适用于条带或线性特征目标的相位滤波方法：

计算公式为：

$\hat{m}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_n{m}_n$

其中，为融合估计值，即滤波后的相位值；m_n和w_n分别为第n个线性窗口的相位平均值和权值；N为所选线性窗口个数，取值范围为1-8，权值w_n的计算公式为：

$w_n\text{=}\frac{1/{{\mathrm{\&sigma;}}_n}^2}{\underset{K=1}{\overset{N}{\mathrm{\&Sigma;}}}1/{\mathrm{\&sigma;}}_k^2};n=\mathrm{1,2},...,N$

其中，

为第n个线性窗口的方差。

3.根据权利要求1所述的InSAR监测高速公路路面沉降方法，其特征在于：在进行高速公路形变相位解缠时，采用了一种适用于条带或线性特征目标的相位解缠方法：首先计算高速公路上的残数点，并将高速公路的边界设置为解缠的限制边界，从高速公路起始的地方寻找第一个残数点作为起始点，以该起始点为中心3×3像素的窗口搜索下一个残数点；若找到残数点则将其与当前残数点相连接；若没找到残数点，则将窗口按高速公路方向从左至右，从上至下依次寻找，直至遍历待检测的整条高速公路；连接所有残数点，并建立枝切线，随后，沿高速公路避开枝切线进行逐点积分，则任意一点p的解缠后的相位为：

其中为相位梯度，C为避开枝切线的解缠路径，为解缠起始点p₀的相位，直至搜索完全部残数点，完成相位解缠。

Images