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

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

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

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

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

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

具体公式如下：
（多普勒方程）（斜距方程）（椭球方程）)()(2
)
()(1)(2
22
12
2121s r s r D s r s r ij p
e V V R R R
f 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点作为已知点，采用最小二乘法求解参数。

模型公式如下：
2
21321b 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.94326
108.94106
108.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）表示为最终的需要求解的参数值。

}{,,,,,,φ
φ
φ
φ
φ
φ
φ
k
n o i s e
x k
o r b i t
x k
a t m o s
x k obj
x k defo
x k topo
x k x W +++++=
φk
topo x ,=h h R B k x k inc
x k x K X
∆=∆⊥βθλπ
,sin 4
φ
k defo
x ,=
k x r ∆λ
π
4 （3）
m k m k k obj x x x ,,,ηξφφφ
+==x m x m
k x k
x m x m
k x k
ηϑλ
π
ϑλ
π
ξϑλ
π
ϑλ
π
⋅-
+⋅-
)sin 4sin 4(
)sin 4sin 4(
φ
k
o r b i t
x ,=ηξ⋅+⋅+c b a
ijnoi i
j ij ij t v k φλ
ππφ++-=42 （4） 相位解缠的目的就是要解得方程中的整周未知数ij k 。

LAMBED 法可以快速准确搜索出GPS 整周模糊度，这里同样可以用来求解CR 的整周未知数ij k 。

为此，建立相位观测方程的随机模型和函数模型，用最小二乘法解得整周未知数ij k 和沉降速率j v 的浮点解，代入LAMBED 程序确定ij k 的整数解并最终确定沉降速率j v 。

利用西安的三景影像，采用LAMBED 方法分别求取了XJ02-XJ03和XJ07-XJ08间的整周未知数和年沉降速率。

干涉对 点号 整周未知数
090103-090314 Xj02-xj03 1
090103-090314 Xj07-xj08 0 090103-090627 Xj02-xj03 -1 090103-090627
Xj07-xj08 0
表4将解算的CR 点年沉降速率与GPS 测得结果做了比较，发现CR 点间相对沉降量总体趋势与GPS 测得结果是一致的。

要证明CR 结果的可靠性，需要比较同影像获取时间一致的水准测量结果。

5 总结
人工角反射器技术（CR ）作为InSAR 学科的一个分支，可用于监测城市地面沉降、山体滑坡等地质灾害，已经被越来越多的研究机构所重视。

城市中布设的CR 点，由于周围存在太多的强散射体，以致在SAR 图像中不容易分辨出来，本文提出了利用CR 安装前后影像强度偏差的思想准确识别出了安装在西安的CR 点。

把CR 点作为地面控制点，采用仿射变换模型，对地理编码后的SAR 影像做校正，提取了三个GPS 点做检验，得出校正后的结果明显优于校正前。

用CR 技术做形变监
点号 07年12月-06年12月 08年6月-07年6月 08年12月-07年12月 CR 点年沉降率
XJ07-XJ08 0.84(cm) 0.6(cm) -1.07(cm) 0.095cm XJ02-XJ03
-1.93(cm)
-4.52(cm)
-3.56(cm)
-2.215cm
表4 CR 点间形变量与GPS 比较
表3 LAMBED 方法求解整周未知数
测是其最终目的，针对目前存在的结算模型及方法，用LAMBED解缠思想，得到两CR点的相对形变量，和GPS比较具有较一致性。

参考文献
[1] 向荣. CRInSAR大气校正及监测地表形变研究[D].湖南长沙：中南大学，2007.
[2] 程滔.CR、PS干涉形变测量联合解算算法研究与应用[D].北京：中国地震局地质研究所，2007.
[3] 高成发,赵毅,万德均. 用LAMBDA改进算法固定GPS整周模糊度[J]. 武汉大学学报. 信息科学版, 2006,31 (8) : 744～747.
[4] Jos.Groot,Matern.Otten,”SAR Imaging of Corner Reflectors Larger Than the Spatial Resolution,”IEEE TRANSACTIONS ON GEOIENCE AND REMOTE SENSING, 1994,3(20):721-723. [5] TaKako Sakurai_Amano. Deteetion of L一band Comer Refleetors in JaPanese Residential Envirotunent from JERs一1 Sar Imagery.
[6] 单世铎,赵拥军. 干涉合成孔径雷达(INSAR)复图像配准方法[J].测绘学院学报，2005，2(22):131-133.
[7] 尤红建,向茂生. 机载干涉雷达数据的直接地学编码处理[J].现代雷达，2006，1(28):38-41.
[8] 陈尔学，李增元.星载SAR影像地理编码算法研究[J].高技术通讯，2000，2(10):53-58;；
[9] 刘秀芳,刘佳音,洪文.SAR图像地理编码处理流程的研究[J].遥控遥测，2005:12-17.
[10]Y. Xia;Hermann,Kaufmann;xiaofang,guo.Differential SAR Interferometry Using Corner Refletors[J]. IGARSS 2002,Toronto,24-28 June 2002,pp.1243-1246.
[11]Y. Xia, H. Kaufmann, X. F. Ge, and X. D. Yang.．CR-Based TerraSAR-X Interferometry for Landslide Monitoring[J]; IGARSS 2008.
[12]X.L. Ding.Acquisition and Analysis of SAR Images for Studying Slope Movement and Ground Deformation in Hong Kong[R].July 2007;
[13]BERT M.KAMPES.RADAR INTERFEROMETRY Persistent Scatterer Technique[M].German Aerospace Center (DLR):2006;。