InSAR应用中角反射器研究

InSAR中人工角反射器研究摘要：传统的差分干涉测量（DInSAR）受到时间和空间失相干以及大气效应影响，在应用上受到很大限制，人工角反射器（CR）在SAR图像上能够显示稳定的、清楚的、较高的振幅信息，特别适合于低相干区域的形变监测，近年来得到了广泛的应用。

本文基于西安市布设的三面角反射器，研究了城市中角反射器的识别问题，同时应用CR对地理编码后的影像通过仿射模型进行了校正，最后求解了CR点间的形变量。

关键词：InSAR 人工角反射器地理编码 LAMBED 地表形变1 引言差分合成孔径雷达干涉测量技术（DInSAR)在最近10来年有了很大的发展。

但是，由于时间和空间的失相干以及大气效应的影响，传统的差分干涉测量的应用受到很大的限制。

基于此，一些学者提出了利用离散的、相位稳定的目标点作为研究对象的新技术[1]。

人工角反射器（CR）由于可人为地控制其几何形状、尺寸、结构和安放位置，因此在SAR图像上显示出稳定的、清楚的、较高的振幅信息，能够实现在低相干区域进行InSAR技术监测地表微量形变的潜力，近年来得到了广泛的应用和发展。

CR被安装在研究区域，雷达入射光线照射到CR相互垂直的两个或三个表面，经过几次反射，入射光线将沿原路径的逆方向反射回去，在图像上形成十字丝形状的亮点，亮点尺寸和亮度宽度均为一个分辨单元[2]。

目前，许多InSAR机构已经开始研究利用CR来探测城市地表微量形变、滑坡变形监测等，也相继布设了一系列的CR点。

然而很多时候在SAR图像上找不到相应的CR点，特别是在城市布设的CR点，由于周围地物的强反射特性，很难准确识别出真正的CR点。

为此，需要专门对城市布设的CR点做探测研究。

同时，CR本身被当做地面控制点，还可以对地理编码后的结果进行校正，以满足SAR图像平面上的精度要求。

基于CR原理求解形变量的理论目前还处于研究阶段，其中的难点就是CR点的相位解缠问题，LAMBED法可以快速准确求解出GPS整周模糊度,在CR 解算中同样可以采用这种方法来解决相位缠绕问题[3]。

2 CR的探测研究表明，即使三面角反射器物理大小大于SAR的空间分辨率，它在SAR图像上仍然表现为一个单独的点目标[4]。

为了准确地找出这个点目标，可以根据CR安装前后的SAR影像的强度信息进行对比研究，进而确定CR点[5]。

为此，需要对获得的包含研究区域的6景Envisat ASAR影像进行处理以探测CR点，具体影像以及CR安装情况如下表所示：表1 六幅影像CR安装情况CR点名20061125 20070310 20071006 20071110 20090103 20090314XJ07 没安装没安装安装安装安装安装XJ08 没安装没安装没安装没安装安装安装从表中可以看到，影像20090103、20090314上全部安装了角反射器，20061125、20070310影像两个角反射器都没有安装，20071006和20071110影像只安装了XJ07。

下图1为安装在西安污水处理厂的XJ07和XJ08两个角反射器在Google earth上的大体位置，图2为具体的CR示意图。

CR的初始坐标由GPS精确测定，单视复数影像图（SLC）经过裁剪、配准、重采样、坐标转换、提取强度偏差后[6],得到了XJ07和XJ08两个角反射器在上述六景SAR图像上的位置(图3)。

图 1 Google earth 上CR 的位置 图 2 安装的CR20061125 20070310 2007100620071110 20090103 20090314图3：CR 点在强度图上的影像3 利用CR 做地理编码校正经过成像和干涉处理得到的InSAR 数据在方位-距离成像几何坐标系中可以表示为(i, j,H,σ),其中: i 表示方位向的序号; j 表示距离向的序号;H 表示该像素对应的高程;σ表示雷达的后向散射量,即图像的灰度。

地学编码处理就是将成像几何坐标系中的数据变成(x,y,H,σ),x,y 表示通用的参考坐标系(我国一般都采用高斯平面坐标)[7]。

目前地理编码采用的方法很多，多项式、共线方程和雷达成像几何模型是经常采用的方法。

多项式模型和共线方程模型都必须依赖地面控制点来反演出模型参数,XJ07XJ07北北XJ07XJ08北北北北XJ07XJ08然后才能进行图像的定位纠正,是一种“地-空-地”的地理编码处理方法,不能直接进行地学编码处理,必须依赖一定的地面控制点。

在InSAR 地理编码中由于很难找到一定数量的地面控制点，因此经常采用距离-多普勒的地理编码方法[8]。

这种方法主要由椭球方程、斜距方程、多普勒三个方程建立地面点和SAR 图像之间的关系。

具体公式如下：（多普勒方程）（斜距方程）（椭球方程）)()(2)()(1)(222122121s r s r D s r s r ij pe V V R R Rf R R R R R R z h R y x -⋅--=-⋅-==+++λ （1）只要给定SAR 图像上每个像素的斜距ij R 和多普勒中心频率D f ,联立求解椭球方程、多普勒方程、斜距方程,即可求得每个像素对应的位置地球惯性坐标系中的坐标)(1,1,11z y x R 。

其中，e R 表示地球半径，))(/11(h R f R e p +-=，h 为该点的地面高程值,f 为平坦度因子，S S V R 、为卫星飞行时间的函数，可以通过卫星提供的轨道数据得到[9]。

然而，通过这样的地理编码计算，得到的WGS84坐标系下的地理坐标(B,L)与实际的地理坐标相差较大，这主要是由不精确的轨道数据产生的。

在西安布设的四个角反射器，这些角反射器都是安装在GPS 观测墩上的，因此很容易通过GPS 精确测定这些角反射器的WGS84坐标，同时可以得到经地理编码后的CR 点WGS84坐标，对比后的结果如下图4所示。

从图中可以看到，纬度偏移量约为0.72秒，经度偏移量约为：-7.4秒，经度上的偏移量相对来说较大。

图4：CR 点上地理编码误差角反射器点最初的应用目的在于作为地面控制点，进而对影像图进行校正。

这里可以借用CR 点作为控制点，对地理编码后的SAR 影像图做校正。

为此，引用仿射变换模型对纬度和经度进行拟合，利用XJ02 、XJ03、XJ07和XJ08点作为已知点，采用最小二乘法求解参数。

模型公式如下：221321b L b B b L a L a B a B SAR SAR GPS SAR SAR GPS +⨯+⨯=+⨯+⨯= （2）其中GPS GPS L B 、表示由GPS 测得的CR 点经纬度，认为是真值，SAR SAR L B 、表示地理编码后SAR 影像上的CR 点经纬度，321321b b b a a a 、、、、、表示相应的参数，代入相应的值可解得参数值为：⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡5759780.00607545 8360361.00045853 5806430.00157584- 1437010.00984951 6607100.00016232- 8821651.00022322 321321b b b a a a 为检验校正后的效果，根据GPS 坐标求取精确的CR 点坐标在SLC 图像上的统计结果，提取了三个GPS 点在SLC 图像上的斜距多普勒坐标，经地理编码，再利用仿射变换模型进行校正，最后的到得结果如表2所示：表2：地理编码校正结果点号GPS 坐标（度）地理编码坐标（度）仿射变化后坐标（度） 地理编码差值（秒） 仿射校正后差值（秒） 1纬度 34.232010 34.232213 34.2320202 0.7305999 0.036555297 经度108.94251 108.94035 108.942439 -7.7637999 -0.26104641 2纬度 34.232015 34.232206 34.23201313 0.68749999 -0.006822 经度108.94293 108.94082 108.9429053 -7.6142999 -0.110737 3纬度 34.232086 34.232275 34.23208210 0.68029999 -0.0141086 经度108.94326108.94106108.94314732-7.9269-0.4233294从表中看出经过校正后点位的地理编码结果得到明显改善，能够满足SAR 影像平面精度的要求。

4 CR 点间提取形变量角反射器点在InSAR 中的主要应用还是提取CR 上的形变量。

国内外许多学者对CR 的模型以及算法进行了研究，其中应用较为典型的是德国波茨坦地学中心xiaye 博士、希腊NPA 团队以及香港理工大学丁晓利教授提出的模型及其算法[10]。

xiaye 博士利用安装在三峡地区的二十多个三面角反射器监测了新滩、树平、卡子湾等地区的滑坡形变，得出树平月最大沉降量可达到50mm ，卡子湾月最大沉降量可达46—53mm 的结论[11]。

丁晓利教授基于荷兰Delf 大学KAMPES 博士关于PS 的解算思想提出了一套新的关于CR 的解算模型和算法，并且成功应用在香港地区的地面沉降监测中[12]。

公式（3）是通用的差分干涉相位方程。

在进行CR 处理的过程中，我们只把地形相位误差和形变相位作为模型参数参与解算，其他相位都作为随机相位的一部分参与处理[13]。

在实际的差分干涉相位方程中，大气延迟相位是一个很难确定的量，但是通过二次差分相邻的CR 点，可以有效地减弱甚至是消除大气延迟相位。

CR 点的高程可以由GPS 精确测定，因此地形相位贡献认为是确定的。

轨道相位和目标点散射相位通过一定的措施可以减弱。

因此方程中只剩下形变相位贡献和噪声相位，噪声相位认为是均值为0的高斯噪声。

公式（4）表示为最终的需要求解的参数值。

}{,,,,,,φφφφφφφkn o i s ex ko r b i tx ka t m o sx k objx k defox k topox k x W +++++=φktopo x ,=h h R B k x k incx k x K X∆=∆⊥βθλπ,sin 4φk defox ,=k x r ∆λπ4 （3）m k m k k obj x x x ,,,ηξφφφ+==x m x mk x kx m x mk x kηϑλπϑλπξϑλπϑλπ⋅-+⋅-)sin 4sin 4()sin 4sin 4(φko r b i tx ,=ηξ⋅+⋅+c b aijnoi ij ij ij t v k φλππφ++-=42 （4） 相位解缠的目的就是要解得方程中的整周未知数ij k 。