第30卷 第1期测 绘 学 报
V o l
.30,N o .1 2001年2月
A CTA GEODA ET I CA et CA R TO GRA PH I CA S I N I CA
Feb .,2001
文章编号:100121595(2001)0120060207中图分类号:P 237 文献标识码:A
刘国祥，丁晓利，李志林，陈永奇，章国宝
(香港理工大学土地测量与地理资讯学系，香港)
Co -Reg istration of Satellite SAR Com plex I mages
L I U Guo 2x iang ,D I N G X iao 2li ,L I Zh i 2lin ,CH EN Yong 2qi ,ZHAN G Guo 2bao
(D ep t .of L and S u rvey ing and Geo 2Inf or m atics ,T he H ong K ong P oly technic U niversity ,H ong K ong ,Ch ina )
Abstract :Spacebo rne Synthetic A perture R adar Interferom etry (InSA R )and D ifferential InSA R (D InSA R )
have been demonstrated to be the very po tential too ls fo r generating large 2area even o r globalD E M as w ell as m easuring and monito ring Earth surface defo r m ati on ,respectively .Fo r the purpo se of computing the accu 2rate interferogram ,the co 2registrati on of satellite SA R comp lex i m ages is one of the mo st i m po rtant p rocess 2ing p rocedures invo lved in satellite SA R interferom etry .Based on the discussi on about transfo r m ati ons a 2mong several spatial coo rdinate system s associated w ith interferom etric p rocesses ,th is paper p ropo ses the ge 2om etric co 2registrati on schem e of SA R i m ages depending on satellite o rbit state vecto rs and SA R i m aging ge 2om etry .O n the o ther hand ,the co rrelative co 2registrati on algo rithm that m akes full use of amp litude info r 2m ati on of i m age ch i p s is also designed on the basis of initial results of geom etric co 2registrati on .F inally ,the co 2registrati on experi m ents are carried out using th ree types of comp lex SA R i m age pairs that w ere co llected by ER S 21
option602,JER S 21,and RADA R SA T ,respectively .Som e significant comparisons and analyses are m ade .Key words :SA R comp lex i m age ;geom etric co 2registrati on ;co rrelative co 2registrati on ;radar interferom e 2
try
摘 要：星载干涉雷达测量已被证实用于生成大范围(乃至全球)的D E M 或监测地表位移具
有广阔的应用前景，其中SA R 复数图像序列的精确配准是保证生成有效干涉相位图的重要环节之一。本文分析了卫星InSA R 系统中所涉及的几种坐标系统变换问题，并提出了基于此的几何配准方案，即基于卫星轨道状态矢量、SA R 成像几何进行配准的算法；同时，还设计了基于几何配准结果和能量影像相关技术的配准算法。最后使用ER S 21 2、JER S 21和RADA R SA T 卫星获取的3种SA R 复数图像对分别做几何配准和相关配准试验和比较，并对 卫星SA R 图像质量和轨道数据质量作了定性分析。
收稿日期:2000201204；修回日期:2000204229
基金项目：香港理工大学研究设备基金(G 29025)；香港理工大学研究基金(G 2V 747)
作者简介：刘国祥(19682)，男，湖南临澧人，西南交通大学测量工程系讲师，在读博士生，主要从事InSA R 研究。
关键词:SA R 复数图像；几何配准；相关配准；雷达干涉测量
1 引 言
近十多年来，国际上掀起了对合成孔径雷达
干涉技术(Syn thetic A pertu re R adar In terferom 2etry ，简称InSA R )理论和应用研究的热潮。In 2SA R 具有从覆盖同一地区的SA R 复数图像对提
取干涉相位图，并借助于雷达成像时的姿态数据重建地表3维模型(即数字高程模型)的巨大潜力[1～4]，特别地，基于多幅SA R 图像处理的差分干涉技术(D ifferen tial InSA R ，简称D 2InSA R )可以用于监测地表形变，精度可达厘米级甚至更高，其监测空间分辨率是前所未有的[5～7]。与机载SA R 系统相比，星载SA R 系统具有可操作性强、方便、SA R 图像相对廉价的优势，因而星载SA R 图像是InSA R 应用的主要数据源，故卫星SA R 干涉处理也是目前研究的焦点。星载SA R 系统一般使用单天线采集信号，故对某个局部地区来说一次卫星通过只能获得一幅影像，卫星以一定的时间间隔和轻微的轨道偏离重复对该地区成像，相邻轨道间隔构成合成干涉基线B (可分解为平行于侧视方向的分量B 平和垂直于侧视方向的分量B 垂，如图1所示)，一般为几十米至一公里左右。一般情况下，卫星轨道数据和传感器参数附在
SA R 图像上一并提供给用户[2,3]
。
图1 InSA R 几何
F ig .1 InSA R geom etric configurati on
一般来说,InSA R 的后处理必须面临一个不可逾越而且是影响干涉相位图合成质量的至关重要的问题，即重复轨道复数图像的空间配准(co 2
registrati on )和重取样[2,3]
。已有一些学者建议完全采用基于影像匹配的方法来进行SA R 图像的配准，只是他们各自所采用的相关测度不同而已[1,8,9]，但配准效率和可靠性方面还有待提高。而本文首先提出了完全通过卫星轨道的状态矢量(瞬时位置和速度)和成像几何作SA R 图像空间配准的算法(后面称几何配准)；随后设计了基于几何配准结果和能量影像信息相关技术的配准算法，即利用Fou rier 正反变换求互相关函数并以
此寻求相关峰值点，从而确定配准像素对(后面称相关配准)；最后使用ER S 21 2、JER S 21和
RADA R SA T 卫星[10]
获取的3种SA R 复数图像对分别做几何配准和相关配准的试验及比较，对轨道数据质量作了定性分析。另外，通过相关配准点对的影像坐标差统计特性讨论卫星SA R 图像的一般特征。值得指出的是，相关配准结果仍是最可靠、也是最终的配准结果。
2 干涉处理中坐标系的转换
这里首先定义3种坐标系:①地球固定坐标系(Earth Cen tered Ro tating ,ECR )，如图1中的O 2X Y Z ，这是卫星轨道状态矢量(瞬时位置s i 和速度v )的标准框架;②卫星SA R 平台坐标系s -xy z ，以天线中心为坐标原点，坐标轴指向如图1中矢量表达式所定义，注意，这些矢量的参考系均是ECR 系;③SA R 图像坐标系，斜距向(即影像行向)垂直于卫星飞行方向，相邻像素具有一定的斜距间隔，如ER S 21SA R 的斜距取样间隔约为7.9m ，而方位向(即影像列向)与卫星飞行方向平行，实际上对应于卫星的一个历元，如ER S 21SA R 的方位向取样时间间隔约为0.000595s ，对应于实地约3.9m 。对于一幅SA R 图像第一行对应于一个确定的卫星历元(时间T 0)，第一列对应于一个确定的从卫星到地面目标的斜距R 0。实际处理中，除了大地坐标和空间直角坐标的相互转换外，还有SA R 影像坐标和大地坐标的相互转换，下面将以SA R 平台坐标系为过渡说明SA R 图像坐标和大地坐标的相互转换。
2.1 侧视方向(给定单位矢量r )与地面的
交点坐标计算
如图1所示，已知传感器中心S 在ECR 系下的位置矢量为s =(s x ,s y ,s z )，若已知雷达侧视方向(如图1中的R 1或R 2)在ECR 系下的单位矢量为r =(r x ,r y ,r z )，则沿该方向与地面有且仅有一个交点P (设地面大地高为h ),P 的ECR 坐标(X ,Y ,Z )近似满足如下关系
X a +h
2
+
Y a +h
2
+
Z c +h
2
=1(1)
其中,a ,c 分别为椭球的长短半轴(如可以选用W GS 84椭球)。另一方面，由矢量关系有
P =(X ,Y ,Z )=S +l r =
(S x +lr x ,S y +lr y ,S z +lr z )(2)其中,l 为传感器中心S 至地面点P 的斜距。联合
1
6第1期 刘国祥等：星载SA R 复数图像的配准
式(1)、式(2)可容易解得P的ECR坐标(X,Y, Z)。
2.2 SAR图像坐标到大地坐标的正反算
由于反算处理与正算类似，在此只讨论正算问题。为计算方便，定义一个以卫星SA R传感器中心为原点的直角坐标系s-xy z，使卫星速度矢量v和其位置矢量s(如图1所示)位于x z平面内。根据SA R成像的基本原理,SA R图像的某一像素的Dopp ler频率F满足下列等式[11]
F=2 V co s8 Κ(3)其中,Κ为雷达波长,8为卫星瞬时速度v方向与雷达侧视方向R的夹角(如图2所示)，对于零多普勒投影有F=0，对于非零多普勒投影F应根据Dopp ler中心系数求解[3,7,11]。设从传感器中心至地面点的空间矢量为R，根据图2中的矢量分解和矢量运算法则，有
R=(x R x,y R y,z R z)=
R ·(x co s8y sin8co s;z sin8sin;)(4)
x=V V y=V×S V×S z=x×y(5
)
其中,；为R在y z面上的投影与y轴的夹角，一般称为俯角，平台系三个坐标轴指向的单位矢量分别为x,y和z，这里(×)表示矢量的交叉积。
图2 平台坐标系及雷达侧视矢量分解
F ig.2 P latfo r m coo rdinate system and decompo siti on of
radar side2look ing vecto r
给定像素点坐标(设沿法线方向n的大地高为h，相应斜距 时间为(R,T),Dopp ler频率为F)，从式(3)不难求出卫星瞬时速度v方向与雷达侧视方向R的夹角8。对于俯角；需迭代求解，将SA R图像坐标转化为大地坐标(Β,L)的处理可按如下步骤完成:
1.由卫星轨道状态矢量传播器(satellite o r2
b it vecto r p rop agato r)计算在ECR系下卫星T 时刻的瞬时位置和速度矢量(s,v)，按式(5)可计算卫星平台局部坐标系3轴的单位矢量x,y和z，并根据SA R传感器参数给定初始俯角;; 2.根据式(4)计算对应的侧视方向单位矢量r;
3.由式(1)、式(2)求沿r从传感器中心到地面点的斜距l;
4.按一定的规律改变俯角;，重复2、3两步直到l与给定斜距R之差足够小为止;
5.按式(2)求对应地面点的ECR坐标(X, Y,Z)，并可继而转化为大地坐标(Β,L)。
3 基于卫星轨道数据和成像几何的SAR复数图像配准
如图1所示，在S1轨道上获取的图像称为主图像(m aster)，在S2轨道上获取的图像称为从图像(slave)。选定主图像为目标空间，将从图像配准到主图像空间的几何配准方案如下所述。
首先确定主从图像地理空间覆盖范围，并计算出覆盖区域内主图像的斜距 时间坐标(R,T)范围{[R m in m,T m in m];[R m ax m,T m ax m]}，并在此范围内设置一个M×N的格网(如6×6)，这些格网点被称为主点(m aster po in ts)，设主点在主图像中的斜距 时间为(R m,T m)，对应的像素坐标为(r1,t1)，已知斜距向取样间隔为∆R，方位向时间间隔为∆3，则有
R m=R m in m+r1∆R T m=T m in m+t1∆T(6)实际处理时，设主点的大地高为0。然后，根据上述算法将(R m,T m)转化为大地坐标(Β,L)，最后由大地坐标反求在从图像中对应的斜距 时间(R s,T s)，继而求解出相应的像素坐标(r2,t2)(一般为非整数值)，配准点对的斜距向和方位向像素坐标偏差分别为
∃r1=r2-r1
∃t1=t2-t1(7)当所有主点(r1,t1)对应的∃r1、∃t1、全部求解出来后，采用两个多项式(如双三次)分别拟合斜距向和方位向像素坐标偏差，并可利用最小二乘算法求解多项式系数。
由于联系点高程数据一般难以得到，上述处理中假设主点大地高为0。实际上，这一假定只会对斜距向的几何配准产生影响，而对方位向配准不会产生影响，对斜距像素偏差影响大小应是基线倾角a、基线长度B、地面高程h、雷达侧视角和斜距的函数。如图3所示，当地面点高程为0时，雷达瞬时位置S1和S2到该点的斜距差∆R可近似为
26测 绘 学 报 第30卷
∆R =R 2-R 1≈-B sin (Η
0-a )(8)
当地面点高程变为h 时，引起侧视角增量∆Η1(近
似等于h (R 1sin Η0))，斜距差∆R ′可表示为∆R ′=R 2′-R 1≈-B sin (Η0+∆Ηt -a )≈-B sin (Η0-a )-B co s (Η0-a )
h R 1sin Η0
(9)
式(8)和式(9)相差可估计出因忽略高程造成的几何配准斜距偏差的误差+∃r 1
+∃r 1
=[B co s (Η0-a )
h R 1sin Η0
] r c (10)
这里,r c 表示雷达斜距取样间隔。以ER S 21SA R
系统为例，设基线长度B 和倾角a 分别为200m ,25°，在最远、最近斜距处(雷达侧视角和斜距分别取24.5°,874km 、21.5°,852km ),+∃r 1随高程h 变化的趋势如图4所示，尽管+∃r 1与h 有线性关系，但斜率很小。随斜距增大,h 对+∃r 1的影响变小；比如，高程为1km 时，误差接近0.08像素；高程达到10km 时，误差接近0.8像素。因此，对于大地高较小的地区，几何配准时假定联系点高程为0，对斜距向配准的影响是可以忽略的
。
图3 斜距差与高程的关系
F ig .3 R elati on betw een range difference and heigh
t
图4 高程变化对斜距向配准误差的影响
F ig .4 Effects of heigh t variati on on co 2registered erro r
of range di m ensi on
4 基于能量影像的相关配准
相关配准单纯从影像能量互相关的角度来考察。在主从图像地理空间覆盖范围内的主图像空间设置一个K ×L 的格网(一般比几何配准时所
采用的格网要密，如8×8)，各格网点的影像坐标
为(r 1,t 1)，利用几何配准时所估计的相应拟合多项式计算出在从图像中的影像坐标(r 2,t 2)，对应
的初始坐标偏差为∃r 1、∃t 1
，这对点称为初始联系点(tie po in t )；以此为基础采用下面的影像相关
算法作配准，使配准精度达到子像素级。
以联系点(r 1,t 1)、
(r 2,t 2)为中心分别从主从SA R 图像提取一定窗口大小(如256×256)的影像块(ch i p s )，设其能量描述函数分别为f m (x ,y )、f s (x ,y )，对应的频谱函数分别为F m (U ,V )、F s (U ,V )，根据信号相关定理(或称维纳2辛钦定理)[12]，两个影像块的互相关函数和互能谱是一个Fou rier 变换对，即
Ν[R f m f s (x s ,y s )]=F m F 3
s
(11)R f m f s (x s ,y s )=
∑+M
-M
∑+N
-N
f
m
(x ,y )f s (x -x s ,y -y s )
(12)中心架
式中Ν表示Fou rier 变换,(3)表示共轭复数,(x s ,
y s )表示信号的空间位移(实际处理时有一个限定范围)，采用快速Fou rier 变换将极大地提高求解相关函数的效率。基于信号相关定理配准的具体算法设计如下:
1.分别在主从图像空间以初始联系点为中
心提取两个复数影像块，分别计算各自对应的能量影像；对每一块作这样的处理：统计影像块的的平均能量值和其标准偏差，并将能量值大于3倍标准偏差的像素点以平均能量值代替(即平滑处理)，这样做会减少噪声点对相关结果的影响；根据下式计算相关系数[11],
Coh =
2<R f m f s >
R f m f m R f s f s
-1(13)
如相关系数低于阈值(th resho ld ，如0.2)，则此对
人工鱼礁联系点无效，移至下一对处理。相关系数低的原因可能是影像块落入阴影区(shadow ing )，雷达迭掩区(layover )或水域区等;
2.当联系点对有效时，分别对经上述处理后
的主从能量图像块做快速Fou rier 变换，并做共轭相乘，然后做快速Fou rier 反变换，即获得相关函数;
3.寻求极值点(p eak po in t )，即获得近似最大相关配准点，并作可靠性检验，采用双二次多项式拟合法内插出子像素级精度的准确峰值点；值得指出的是，如果影像信号中包含高频窄带随机
3
6第1期 刘国祥等：星载SA R 复数图像的配准
噪声或包含特别强的高频信号，都会使互相关函数出现多个极值点[12]，因此，对求得的峰值作可靠性检验可防止错相关，这里采用峰宽和相对峰值两个标准控制。如果相关结果可靠，则主图像像素(r1,t1)在从图像中的精确配准点为(r2+∃r2,t2 +∃t2)，故精确配准后联系点对的斜距向和方位向像素坐标总偏差分别为:
∃r=∃r1+∃r2
∃t=∃t1+∃t2(14)
4.重复1～3步，基于所有其他初始联系点对作相同的相关配准运算。
5 配准实例及相应结果分析
为了验证上述几何配准、相关配准算法的可靠性，考察这两种配准结果的差异以及比较不同卫星SA R图像配准的情况，这里选取3种卫星SA R系统获取的图像对作为研究对象，其成像时间、覆盖区域、名义干涉基线分量和主图像中选择的用于配准的范围如表1所示。图5显示了ER S12 2联合飞行方式获取的降轨(descending o rb it)图像(参数见表1)，为节约篇幅，显示影像经过了一次减少数据量的重取样，图像块覆盖香港元朗区，主从图像右上角和左下角是海区，信噪比很低且没有相关性。由于篇幅所限，其他的图像这里未给出。再次指出，相关配准结果是最终的配准结果，也是这里的评判标准。
按上述算法对3类图像对进行几何配准时均选择6×6的格网，即均匀分布在图像对重叠区域内的36个主点，而进行相关配准时选择8×8的格网，即64个联系点，两次配准统计结果如表2所示。表2中统计了36个主点经几何配准后，斜距向和方位向像素坐标差均值(∃r1,∃t1)、像素坐标差标准偏差(∆∃r1,∆∃t1)；同时也统计了相关配准后的上述各项指标。值得注意的是，相关配准中依赖于联系点对的有效性和可靠性检验已剔除出不合格的联系点，有效的配准点对会小于64。
表1 3种SAR图像对试验数据
Tab.1 Test dat a:three types of SAR i m age pa irs
SA R图像对类型覆盖区域
成像时间名义基线分量 m主图像中被选择的范围 (m·s-1)
主图像从图像B平B垂R m in
m,T m in
m R m ax
m,T m ax
m
ER S21 2Yuen L ong
(香港)
(ER S21)
1996203218
(ER S22)
1996203219
80.0999.05
val-031
860300.443m,
10500.91799s
864804.138,
10503.43556s
JER S21M ount Fuji
(日本)
19932072071993208220644.16510.48
717933.260m,
5033.00027s
729046.905m,
5034.15594s
RADA R2
SA T
Bathurst岛
(加拿大西北部)
19962032041996203228-909.04-1061.67
1096267.306m,
46960.34965s
1107994.184m,
46962.14110s 表2 3种SAR图像几何配准和相关配准结果比较
Tab.2 Co m par ison of co-reg istration results am ong three types of SAR i m ages P ixels
SA R 图像对
几何配准联系点统计相关配准后联系点统计
∃r1
∃t1
∆∃r1
∆∃t1∆∃1
∃r
∃t
∆∃r
∆∃t
有效
点对
∆∃
ER S21 2-5.881
-65.124
延时电路
±0.045
±0.021
±0.026
-5.978
-66.230
±0.014
±0.017
54±0.056
JER S21-3.503
-30.123
±0.355
±0.176
±0.10
-5.156
22.380
±0.279
±0.280
56±0.191
RADA R SA T -3.940
-7.494
±0.772
±0.773
±0.095
0.088
-11.862
±0.546
wcdl±0.575
49±0.134
对于主图像中的某一待配点来说，分别采用几何配准和相关配准算法在从图像中会得到不同的配准点。从斜距向和方位向像素坐标差均值之差(∃r-∃r1,∃t-∃t1)可以看出，对于ER S21 2的SA R图像对，斜距向差值为0.097像素，方位向差值为1.106像素；对于JER S21的SA R图像对分别为1.653像素、-52.503像素；对于RADA R SA T的SA R图像对分别为-4.028像素、-19.356像素，显然ER S21 2的差异最小。从几何配准点对像素坐标差标准偏差(∆∃r1,∆∃t1)
46测 绘 学 报 第30卷
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议