以Doris进行雷达干涉技术产生数值地形模型之  文件类型：PDF/Adobe Acrobat   文件大小：字节

1
以Doris进行雷达干涉技术产生数值地形模型之
研究
The Study of DEM derived from SAR Interferometry with Doris
谢嘉声1,史天元2
Chia-Sheng Hsieh,Tian-Yuan Shih
摘要
合成孔径雷达干涉量测技术(SAR Interferometry, INSAR) 系利用雷达影像的
相位(phase)值来推求地表面的空间位置,是一项可以获得高精度数值高程模型及
地表位移的方法.该项技术虽具有应用的潜力,但因处理的过程复杂且条件要求
严格,目前除应用高价的商业软体外,一般研究人员不易投入该项技术的研发及
应用.
本研究以荷兰Delft大学所发展的公用软体Doris来进行雷达影像产生数值
高程模型的研究.实验影像采用欧洲太空总署ERS卫星在台湾中部地区所拍摄的
雷达影像,目前实作至去除地球曲率的干涉图.该软体的功能相当完整,但因影
像格式的不同,所以须针对差异处进行修改,才能完成所有步骤.
关键词:合成孔径雷达,数值高程模型,雷达干涉量测技术
Abstract
Key Words : SAR, DEM, INSAR
-----------------------------------------------------------------------------------------------------
1国立交通大学土木工程系博士生
2国立交通大学土木工程系教授
2
一,前言
雷达干涉技术(Radar Interferometry),利用不同时间或不同位置的雷达天线所
获得的两幅或多幅复数影像,藉由相位值的差异量来获取地表的三维资讯.目前
主要的应用包含数值地形模型(DEM,Digital Elevation Model)的产生[Zebker et.al.,
1994],地表变形量的研究等[Massonnet et.al., 1993].因为主动微波系统能全天候
的搜集地表资讯,配合相位值高精度的量测技术,因此可用以进行快速的制图及
有效的侦测地表的变形[Madsen and Feigl, 1998][Rosen et al., 2000].
本研究以荷兰Delft大学所发展的公用软体Doris来进行雷达影像产生数值高
程模型的研究,因为该软体主要就是以研究获取地表的三维资讯及地表变形而设
计,且该其为一完全免费的软体,开放软体的程式码供研究者使用[Kampes and
Usai, 1999].以该软体来进行干涉的处理,除可了解该项技术的实用性,更可配
合各种研究项目进行功能的修改,适合对雷达影像及干涉处理有兴趣者使用.
因为Doris只能接受欧洲太空总署ERS-1/2卫星所获取的雷达影像,本研究
则以中大太遥中心所接收的雷达影像来进行实验,实验区为台湾中部地区,目前
处理的进度为地球曲率修正的干涉图.
二,研究方法及流程
以雷达干涉量测技术来产生地形资料的理论如图一所示[Zebker and
Goldstein,1994].
B A2 α
A1
θ γ+δ
γ
h
Z(y)
y
图一,雷达干涉几何示意图
由两个天线所接收的同一地面点讯号的相位差值为 ,因相位差为距离差值δ
所产生的影响;且在此发射及接收的模式中,其讯号行走之路径差为2δ,故其相
关式为δ=λ /4π,由图一的几何关系配合三角公式的推导可得(1)
()
()
Zyh
B
B
=-
-
·--
·
λ
π
θα
λ
π
θ
4
2
4
2
2
sin
cos (1)
其中:
Z(y) : 地形高度 h : 卫星高度
λ : 雷达波波长 : 相位差
B : 基线长 θ : 视角
α : 倾角
从(1)式中可知,当像对间几何关系中的相关参数值已知,便可计算出其对应
的高程值.
3
因为雷达的成像原理为斜距投影,与传统的光学影像不同,其方位(azimuth)
方向的位置与雷达拍摄的时间相关,斜距(range)方向的像元位置与雷达到地面间
的距离有关[Rantakokko and Rosenholm,1999].Doris在进行处理时,将每一个步
骤当成一独立的计算模组,每一模组处理完后会有输出档纪录处理的结果及数
据,下一个模组在执行前会先读取所需的参数数据,所以每一个模组都可以单独
进行测试,以了解其运作的情形及选择最合适的参数值.其整体的作业流程如图
二.
Master SAR image
(SLC)
Slave SAR image
(SLC)
1.读取档头资料
(M_READFILES)
2.精确轨道计算
(M_PORBITS)
3.选取影像范围
(M_CORP)
基线计算
(COARSEORB)
1.初步套合(COARSECORR)
2.方位滤波处理
(M_FILTAZI,S_FILTAZI)
3.精密套合(FINE)
4计算转换参数(COREGPM)
5.重新取样(RESAMPLE)
斜距滤波处理
(FILTRANGE)
计算干涉图
(INTERFERO)
消除地球曲率
(COMPREFPHA)
数值地形模型
DEM
地理座标转换
(GEOCODE)
高程计算
(SLANT2H)
全相位回复
(UNWRAP*)
差分干涉量测
(DINSAR)
滤除相位杂讯
(FILTPHASE)
产生相干图
(COHERENCE)
干涉图
(INTERFERO)
以DEM消除地球曲率
(COMPREFDEM)
精确轨道
资料
1.读取档头资料
(S_READFILES)
2.精确轨道计算
(S_PORBITS)
3.选取影像范围
(S_CORP)
图二,Doris作业流程图,图中括弧内名称为Doris的模组名称.
目前Doris的模组以可从一开始的SLC影像处理至 DEM的结果,但中间的
全相位回复(Phase Unwrapping)模组仍未完成,所以必须配合其他的软体来完成.
在整个处理模组中,其必须按照流程的顺序来处理,但几个滤波的处理,则可以
选择是否进行,这些处理可以针对不同的杂讯於以滤除,提升处理的结果.
2-1影像格式
Doris所能处理的SAR影像为ERS 1/2 SLC(Single Look Complex)的影像,储
存的格式以符合SAR CEOS格式的ESA SAR.SLC 规格[Ola Grabak,1998].
2-2基线计算
基线的计算为计算影像对间的各种几何参数值,如图一所示,地面的高程值
就必须由相关的几何参数来计算,所以各参数的计算结果将影直接影影想成果的
品质.
4
基线的计算方法:因为为重复轨道的方式,所以两卫星轨道间的几何位置为
接近於平行的两曲线,且基线为拍摄同一地面点时两卫星的空间位置,依SAR
的成像原理可知,其斜距方向与轨道方向为接近垂直之情形,其位置为两轨道间
之最短距离.如以master的轨道为例,选择像中心所对应的地面点为共同点,则
以此影像中心的获取时间,配合精密轨道可以计算出卫星的瞬间位置.而在
Slave的影像上,其共轭的卫星的空间位置为在 Slave的轨道中,求出与前卫星
轨道位置距离最短的位置,即为相对应的卫星空间座标.以此两座标就可以计算
出各种基线参数.
1.实际的计算共轭点座标:先以 参考影像(master image)的中心点为主,求出中心
点的卫星位置,再以此卫星位置,影像的座标(Line,pixel)值依 Doppler方程式来
求所对应的地面点座标.有此地面点座标后,再依Doppler方程式来求出所应的
套合影像(Slave image)的方位时间,有此方位时间即可依精密轨道求出其卫星的
空间座标.有此两座标,就可以计算像对间的各种基线参数.
2.经确定其两点的座标后,其相对应的各种基线值如图三所示.其计算式为:
)(22
VHBBB+==)(2
//
2BB+⊥,
2
)(tan)(tan
//
11π
θα-+==⊥--
B
B
B
B
H
V
)cos(α·=BBH,)sin(α·=BBV,)cos(αθ-·=⊥BB,)sin(//αθ-·=BB (2)
图三,基线各项参数的几何关系
2-3 套合处理
套合处理的目的在确定两影像对间的空间(spatial)及光谱(spectral)的重叠区
域,一般而言套合处理在计算两影像间的转换参数及重新取样.其主要的处理步
骤:1.计算干涉像对的参考影像(master image)及套合影像(slave image)位移量.2.
计算预估轨道及地球的几何关系.3.计算转换参数.4.将套合影像转换至参考影
像坐标.
为确保套合处理能自动化且非常准确,所以将处理的部份分为三个部份,第
一先以精密的轨道资料来计算两影像间的大约位移量.计算步骤为:
(1) 由参考的中心像元(line,pixel)计算点P在椭球上的座标(X,Y,Z)
(2) 以 Doppler 方程式来计算点P在套合影像中的(line,pixel).
(3) 两影像间的line及pixel差即为其位移量(offset).
因为Doppler 方程式为非线性的方程式,所以必须给予适当的起始值,该值
5
可以由影像档头内影像中心的经纬度值来计算.如果最后的结果没有收敛,可以
将两影像以人工量取的方式来计算其大约的位移量,其精度只要在30个像元内
即可.
第二步骤为计算两影像的相关系数来决定其匹配的位置.该方法以不同的整
数位移值来求取两影像间的相关系数,并以最大的相关系数值所在的位移量为偏
移量.该方法可以将套合的位置提升至1个像元的准确度.其计算方式可以选择
在空间域的方式来计算,或将影像转换至频率域来计算.
雷达影像因为随机性的斑驳,使参考影像的特徵在套合影像上未必是特徵;
且因为地物背向散射系数的敏感性,使地物在两影像上可能有不同的轮廓,无法
以一般的特徵方法来进行套合.而较常使用的方法为将两影像中重叠区域区分为
若干个网格点,每个网格点选取邻近的像元成为一区块(patch),而后以两区块的
灰阶亮度(power)计算最大相关值来求区块位移量.
第三步骤为精密的套合,该方法可将两影像精密的套合至次像元的精度.该
方法采用较大的视窗,仍然以相关系数来计算,并求出次像元的精度.在以该法
来实施套合的处理前,必须确定其初始位置在几个像元之内.
在计算出点位次像元的套合精度后,便可多项式来代表两者间的转换关系.
计算前将相关系数值较低的点位去除,再以最小二乘法求出该转换参数.求出转
换参数后须将套合影像依前述的多项式,予以重新取样,经过取样后再与参考影
像进行精密套合的处理.
2-4 全相位回复
经地形效应纠正后之相位图,其所对应的高程变化仍然超过2π,表现在相位
图上形成一圈延著高程的等相位曲线.当我们要计算每一点的高程时,就必须加入
每一相位量测中的整数波数,这个要解决实际加回多少波长(2π)的问题,就称为全相
位回复.要在相位图上将波数辨认出来,这与干涉条纹的品质有关,如果有杂
讯,斑驳及地物变迁等破坏干涉的情形存在,则会有部份区域无法计算.
在doris的软体中,目前仍没有全相位回复处理模组.建议在此项的处理方
法可由其他的软体来处理,Two-dimensional Phase Unwrapping: Theory, Algorithms,
and Software书中所附的软体来进行处理.
三,研究成果探讨
3.1 输入影像资料
因为输入的SAR影像须为 CEOS的影像格式,且Doris会读取影像的档头
及影像的资料,因此必须具备Image,Leader及Volume等三种以上的影像资料.
本研究所使用的影像为向中央大学太空及遥测中心所购买的1996年3月6日
ERS-1及1996年3月7日ERS-2的SLC(Single Look Complex)影像,影像的像元
解析度为斜距方向(range)为7.904890 m,方位方向(azimuth)为3.947885 m .所以
全幅影像的范围为斜距方向(range)38544.243935 m,方位方向(azimuth)为
105759.885855 m .
因为中央大学太遥中心所接收的SAR SLC影像资料,其储存的格式为CCRS
的格式,虽然仍符合CEOS的标准,但其详细的格式内容与ESA的格式有些许
的不同.本研究为避免程式读取资料错误,将档头资料予以部份更改,使其格式
与ESA的格式相符.其中在VOL档的部份,更改的项目为将第三个纪录(records)
6
拿掉,只存4个纪录档.
在Lea档部份,因为ERS的SLC Lea档中包含有SAR影像所获取的时间及
处理的时间,根据此时间可计算其影像的空间位置,以供套合计算的起始值.实
验的影像资料lea档并无此项纪录,所以必须将此一时间加入.但因影像档头中
无此一资料,所以仅能由其他的参数中推估,以计算的方法将此时间值加入Lea
的档中.
经读取档头资料后,便须将影像的主体由资料档中取出,可以选取全部的范
围,亦可选取局部的区域.经选取全幅影像,并以Doris软体的工具程式来展现
影像的灰值图,如图四所示.从图中可以看出两幅影像的范围相当接近,影像所
在的区域为台中,彰化一带.
(a) (b)
图四,(a)1996年3月6日ERS-1影像.(b) 1996年3月7日ERS-2影像.
3.2 影像的基线
以Doris来计算之基线值,因所采用的方法为先以影像的中心来计算地面座
标,而地面座标的近似值则以影像档的中心经纬度值来计算,但因影像档头的格
式有些不同,所以计算出来的近似值有误,以致计算式无法收敛而产生大误差.
为解决此一问题,则由方法二来计算,所计算得之各项基线值如表一.
表一,影像对的基线资料
项目名称 数值
Bperp [m]: -58.50
Bpar [m]: -47.12
Bh [m]: -70.13
Bv [m]: 26.92
B [m]: 77
alpha [deg]: 159
theta [deg]: 17.68
Height_amb [m]: -116.4
Btemp: [days]: 1
7
3.3 雷达影像套合处理
套合处理的步骤可以分成三个部份,首先大约位移量粗估的部份经计算在方
位方向为-105,在斜距方向为-9个像元.经第二步骤的影像相关计算后,修正为
(-94,-6).以此为起始值,进行精密套合的计算,得到的点位及相关系数如表二,
其显示图形如图五所示.
表二,套合点位的相关系数及位移量
点号 line pixel 相关系
数
位移量
_L
位移量
__P
0 121 33 1 3.69 3.29
5 1788 3405 0.93 -1.63 -1.39
19 10121 2249 0.85 -1.01 -1.41
22 11787 2345 0.61 -1.39 -1.19
25 13454 2442 0.75 -1.27 -0.97
28 15121 2538 0.7 -1.02 -0.81
32 16787 4272 0.55 -0.47 0.12
35 18454 4369 0.63 -0.35 0.54
38 20121 4465 0.64 0.41 0.26
40 21787 2923 0.64 -0.35 -0.28
43 23454 3019 0.71 -0.41 -0.18
44 23454 4657 0.62 5.59 1.83
45 25120 1478 0.51 -0.47 -0.53
46 25120 3116 0.54 -0.41 0.04
图五,影像套合处理的点位分布,相关
系数(圆圈)及位移量(箭头)
3.4 产生干涉图
所产生的干涉图如图六.从图中可以明显的看出,未经地球曲率改正的干涉
图其图形的干涉条纹并不清楚,而在经地球曲率改正后即可看出较明确类似等高
线的干涉条纹产生.
(a) (b)
图六,比较是否经地球曲率改正产生干涉图之差异.(a)未经改
正之干涉图 (b)经地球曲率改正之干涉图
8
四,结论与建议
从研究中可知,Doris确实可以有效的进行各个模组的处理,且其说明文件
亦叙述的相当清楚,非常适合雷达干涉影像的处理.但因仅读取ESA SLC SAR
的影像格式,与我们较易获得的中大影像CCRS格式有些不同,造成处理上的困
难.但因所有的程式码都公开,所以只要能确实了每一个步骤处理的细节,便可
以更改参数的方式来解决此一问题.
五,志谢
本研究所使用的处理软体Doris及其工具程式系由Delft大学(Delft University
of Technology)DEOS(Delft Institute for Earth-Oriented Space Research)研究中心所提
供之免费公用软体.并感谢该中心所提供之精确卫星轨道.
参考文献
Kampes, B., and Usai, S. (1999). Doris: the Delft Object-oriented Radar Interferometric
Software, Proceedings ITC 2nd ORS symposium, Netherlands.
Leberl , F. W. (1990). Radargrammetric image processing, Artech House, Norwood.
Massonnet, D., Rossi, M, . Carmona, C,. Adragna, F,. Peltzer, G.,. Feigi, K. and
Rabaute, T. (1993) The displacement field of the Landers earthquake mapped by
radar interferometry, Nature, Vol. 364, pp. 138-142.
Massonet, D.and Feigl, K.L. (1998). Radar Interferometry and its application to changes
in the earth's surface. Reviews of Geophysics. 36(4). pp.441-500.
Ola Grabak. (1998). ERS SAR.SLC CCT and EXABYTE, ESA, ER-IS-EPO-GS-5902.3
Rantakokko, H, and Rosenholm, D. (1999). Rectification of Slant Range Imagery
Through a Direct Image to Ground Relationship, Photogrammetric Record, Vol. 16,
No. 94, pp. 685-694.
Rosen, P.A., Hensley, S., Joughin, I. R., Li, F.K., Madsen, S.N., Rodriguez, E.,
Goldstein, R.M., 2000. Synthetic Aperture Radar Interferometry. Proceeding of
the IEEE. 88(3). pp333-382.
Zebker, H.A., Farr, T.G., Salazar, R.P., and Dixon, T.H. (1994). Mapping the world's
topography using radar interferometry: the TOPSAT mission, IEEE Proceedings,
Vol. 82, No. 12, pp 1774-1786.