星载全方位SAR联合图像配准和DEM提取方法

张金强, 秦 强, 宫 俊, 刘 伟, 苏皎阳

(1.上海无线电设备研究所,上海201109;2.上海目标识别与环境感知工程技术研究中心,上海201109)

0 引言

星载合成孔径雷达(Synthetic Aperture Radar,SAR)因其全天时、全天候的全球观测能力,在军事侦察、国民经济建设和科学研究中得到了广泛的应用[1]。为了实现目标的精细化探测、高精度定位和全方位描述,近年来人们提出了星载全方位SAR系统。SAR卫星按照四组不同轨道根数的轨迹绕同一地面场景飞行,通过控制雷达波束指向实现对目标360°全方位观测[2-4]。每一次航过的雷达波束斜视角范围在-45°～+45°之间,合成孔径时间可达百秒量级。将全孔径回波数据进行子孔径划分,对各子孔径回波数据进行聚焦成像处理,获取多角度SAR图像序列。为了利用单颗SAR卫星单航过录取数据实现地形数字高程模型(Digital Elevation Model,DEM)的高精度提取,本文基于星载全方位SAR获取的多角度图像序列,提出了多角度SAR图像序列配准和DEM提取联合处理方法,并进行了仿真验证。

1 联合配准和DEM提取方法

1.1 星载全方位SAR几何构型

星载全方位SAR几何构型如图1所示。图中所有参数均定义在地心地固坐标系中。雷达天线相位中心(Antenna Phase Center,APC)位置矢量和速度矢量分别用P和V表示。P在地球表面的投影点称为其星下点,用P′表示,P至P′的连线是地球表面的法线。雷达APC与目标最接近的时刻称为零多普勒时刻。零多普勒时刻雷达APC所在位置称为最近位置,用P0表示,P0的星下点用P′0表示。T表示地面目标位置矢量,P到T的矢量R称为斜距矢量,P0到T的矢量R0称为最短斜距矢量。垂直于雷达APC速度矢量并且包含雷达APC的平面称为零多普勒面。R与零多普勒面之间的夹角称为斜视角,用φ表示。当R与V的夹角小于90°时,φ的值为正;反之,φ的值为负。

为了描述不同航过雷达APC位置间的关系,定义方位角和方位夹角这两个参数。在地球表面建立如图1(b)所示的二维直角坐标系,其中坐标原点o为观测场景中心目标所在位置,x轴和y轴相互垂直。o至P′的矢量与x轴的夹角称为方位角,用φ表示。o到任意两个雷达APC位置的星下点P′i和P′j的矢量之间的夹角称为方位夹角,用Δφij表示。

图1 星载全方位SAR几何构型

1.2 多角度SAR图像序列联合配准方法

利用多角度SAR图像序列提取观测场景DEM时,图像对间方位夹角越大,配准误差导致的高程提取误差越小。但是,图像对间方位夹角越大,其相关性越低,配准偏移量估计精度越低[5]。针对该问题,基于星载全方位SAR获取的多角度SAR图像序列方位角连续这一特性,提出了多角度SAR图像序列联合配准方法。该方法利用小方位夹角图像对之间的配准偏移量,估计大方位夹角图像对之间的配准偏移量。

不失一般性,假设所有图像均聚焦成像于斜平面,聚焦多普勒频率等于子孔径回波数据的多普勒中心频率,并且聚焦多普勒频率均匀分布。设聚焦多普勒频率间隔为ΔfD,0,图像数量为2N+1,图像编号用向量表示为[0,1,2,…,2N],则多角度图像序列的聚焦多普勒频率可以用向量表示为[-N,-(N-1),…,0,…,(N-1),N]×ΔfD,0。根据图像对间聚焦多普勒频率差值的不同(即图像对间方位夹角的不同),将其分为小方位夹角图像对和大方位夹角图像对。对于聚焦多普勒频率为fD,i的图像#i和聚焦多普勒频率为fD,j的图像#j,如果两者的聚焦多普勒频率差与判别门限CΔfD之间的关系满足

则称之为小方位夹角图像对,反之则称之为大方位夹角图像对。CΔfD需要根据具体的图像分辨率、信噪比以及观测场景内地形起伏和地物类型等参数确定。小方位夹角图像对数量用K表示,并且K≥2N。

图2给出了N=3,CΔfD=2ΔfD,0和CΔfD=3ΔfD,0时,小方位夹角图像对连接示意图及其数量。随着CΔfD的增大,所选择的小方位夹角图像对的数量增加。小方位夹角图像对间的配准偏移量可以利用相关函数法等传统SAR图像对配准方法进行估计[6]。

图2 小方位夹角图像对连接示意图

下面根据小方位夹角图像对间配准偏移量与大方位夹角图像对间配准偏移量之间的关系构建线性方程组。不失一般性,选择聚焦多普勒频率等于0Hz的SAR图像作为参考图像,参考图像编号为#0。在参考图像中选择任意像素,设其坐标为(aref,rref),该像素对应的目标在其他2N幅图像上的像素坐标可以用向量表示为

式中:a为方位向坐标向量;r为距离向坐标向量;T表示矢量或矩阵的转置操作。在对小方位夹角图像对进行配准处理时,需要指定主辅图像,其编号可以用向量形式表示为

式中:M为主图像编号向量;S为辅图像编号向量;k=1,…,K,表示小方位夹角图像对编号。

K个小方位夹角图像对间方位向和距离向的配准偏移量可以用向量表示为

式中:δa为方位向配准偏移量;δr为距离向配准偏移量。每一小方位夹角图像对之间的方位向和距离向配准偏移量可以表示为

式(9)和式(10)定义了两个方程组,每一个方程组包含K个等式、2N个未知数,可以用矩阵表示为

式中:Aa和Ar均为K×2N维的矩阵;aref和rref均为K×1维的向量。对于∀k=1,…,K和∀n=1,…,2N,Aa、Ar、aref和rref对应的元素可以表示为

式中:n表示矩阵Aa和Ar的列索引。

矩阵Aa和Ar均为稀疏矩阵,其元素值取决于所选择的小方位夹角图像对。由稀疏矩阵的特性以及K≥2N可知,Aa和Ar均为列满秩矩阵。因此,其他各图像相对于参考图像的配准偏移量可以利用最小二乘算法[7]估计得到

利用该方法可以同时估计得到所有图像相对于参考图像的配准偏移量。

1.3 联合图像配准和DEM提取方法

所有图像相对于参考图像的配准偏移量,还可以根据SAR成像几何,利用几何法计算得到。几何法流程包括两步:首先,通过正向定位获取参考图像像素对应的目标三维位置矢量;然后,通过反向定位获取上述目标在图像#n上的坐标。针对参考图像上某一像素(aref,rref),将成像参数和轨道数据代入距离-多普勒模型[8],求解目标三维位置矢量,即正向定位。距离-多普勒模型表达式为

利用目标三维位置矢量T、图像#n成像参数和轨道数据,首先根据式(20)计算该目标对应的图像#n的聚焦多普勒时刻tn,将不同时刻的Pref和Vref代入式(20)计算公式左右两边之间的差值,最小差值所对应的时刻即为tn;然后根据式(19)计算时刻tn目标到雷达的斜距Rn;最后计算该目标在图像#n上的二维坐标(an,rn)。上述过程称为反向定位。其中an和rn的计算式为

式中:tn0表示图像#n方位向第一个像素对应的聚焦多普勒时刻;fp表示脉冲重复频率;Rn0表示图像#n距离向第一个像素对应的雷达斜距;fs表示脉冲采样频率;c表示光速。

利用几何法计算图像对间配准偏移量时,需要输入目标高程h。若输入先验目标高程存在误差,将导致几何法计算得到的配准偏移量存在误差。下面分析几何法计算得到的配准偏移量误差随输入先验目标高程误差的变化。

根据正向定位模型,将式(19)～式(21)中目标三维位置矢量对目标高程求一阶偏导数,得到

根据反向定位模型,将式(22)和式(23)中的聚焦多普勒时刻tn和雷达斜距Rn分别对目标三维位置矢量T求一阶偏导数,得到

其中

式中:Pn、Vn和An分别表示T对应的图像#n的聚焦多普勒时刻的雷达APC位置矢量、速度矢量和加速度矢量。

根据式(22)～式(28),得到目标在图像#n上的二维坐标相对于目标高程的一阶偏导数为

由于星载SAR雷达斜距较大,在局部范围内可以认为,目标在图像对间的二维偏移量随其高程线性变化。

将其他各图像相对于参考图像的几何配准偏移量用向量形式表示为

将利用式(17)和式(18)估计得到的配准偏移量与几何配准偏移量对比,可以得到

其中

因此,几何配准时所用先验目标高程误差可以利用最小二乘算法估计得到

将其与先验目标高程相加,即为目标估计高程。先验目标高程可以设置为0m,通过迭代处理最终实现高精度DEM提取。

观察式(11)、式(12)和式(33),可以利用一个线性方程组将其联合

式中:Ea和Er均为2N阶单位矩阵。

因此,可以利用最小二乘算法同时估计得到所有图像相对于参考图像的配准偏移量和先验目标高程误差

2 仿真数据处理及结果分析

2.1 多角度SAR图像序列仿真

多角度SAR图像序列仿真参数如表1所示。轨道数据利用STK软件生成,四组轨道的星下点轨迹如图3所示,图中箭头指向卫星的飞行方向。观测场景地形在航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission,SRTM)DEM库[9]中提取,高程变化范围为(1225～1537)m。利用分形插值方法[10]对其进行插值处理,地面格网间距设为0.5m,地物类型设为土壤或岩石,目标后向散射系数根据Ulaby模型[11]计算得到。插值后DEM如图4所示。

表1 多角度SAR图像仿真参数

航过一仿真得到57幅图像,聚焦多普勒频率在(-270～+290)kHz范围内均匀分布;航过二仿真得到50幅图像,聚焦多普勒频率在(-220～+270)kHz范围内均匀分布;航过三仿真得到50幅图像,聚焦多普勒频率在(-230～+260)kHz范围内均匀分布;航过四仿真得到47幅图像,聚焦多普勒频率在(-230～+230)kHz范围内均匀分布。四个航过图像序列聚焦多普勒频率间隔均为10kHz,图像信噪比均为20dB。观测场景中心目标对应的四个航过第一幅、最后一幅和斜视角为0°的图像的聚焦多普勒时刻雷达APC位置的星下点,分别如图3中黑色实点所示。在该仿真中,图1(b)中所定义的坐标系的x轴,由观测场景中心目标位置指向其航过一零多普勒时刻雷达APC位置的星下点,y轴方向与航过一雷达APC速度方向相反。

图3 卫星星下点轨迹及观测场景位置

图4 仿真观测场景分形插值后SRTM DEM

图5给出了方位角φ分别为0°,95.98°,176.97°和266.11°的仿真图像,图像大小均为19 000像素×10 500像素。

2.2 仿真数据处理及结果分析

图5 多角度SAR图像

利用本文方法对仿真得到的四个航过SAR图像序列分别进行处理,分别选择方位角为0°,176.97°,95.98°和266.11°的图像作为参考图像。选择小方位夹角图像对时,令式(2)中的α等于3。利用相关函数法逐像素估计小方位夹角图像对间的精配准偏移量。逐像素估计大方位夹角图像对间精配准偏移量时,剔除相关系数低于0.5的小方位夹角图像对间的配准偏移量。如果剩余小方位夹角图像对数量少于原数量的60%或导致式(37)中的观测矩阵为非列满秩矩阵,则将参考图像该像素对应的目标高程值置为无效值。利用几何法计算图像对间配准偏移量时,在仿真DEM上添加均值为0 m、标准差为17 m的随机高程误差,将添加高程误差后的DEM作为先验目标高程输入值。

利用四个航过的多角度SAR图像序列估计得到四组观测场景目标高程,图6(a)～图6(d)分别给出了四个航过目标高程估计误差的分布。高程估计误差为所估计DEM与仿真DEM的差值。统计得到四个航过目标高程估计误差范围分别为(-17.697 8～+16.020 3)m,(-12.329 7～+12.337 3)m,(-15.730 8～+15.337 4)m,(-12.320 8～+14.504 4)m;均方根值分别为1.891 1,1.604 0,1.387 4,1.327 5 m;无效值像素数量分别为1,0,0,0。由上述结果可知,基于星载全方位SAR多角度图像序列,利用本文方法可以高精度提取观测场景DEM。

由于雷达波束照射方向不同,不同角度SAR图像中叠掩区域或低信噪比区域不同,导致利用不同航过多角度SAR图像序列所提取的目标高程的无效值像素以及估计误差大小的分布不同。为了提高DEM提取精度、降低无效值像素的数量,基于相关系数的高低将四个航过获取的DEM进行融合处理。图6(e)给出了融合后目标高程估计误差的分布。融合后目标高程估计误差范围为(-4.199 8～+3.699 7)m,均方根值为0.443 4 m,无效值像素数为0。由上述结果可知,通过融合基于不同角度SAR图像序列获取的DEM,可以提高观测场景DEM提取精度,同时降低无效值像素的数量。

图6 高程误差分布

3 结论

针对星载全方位SAR获取的多角度图像序列,提出了多角度SAR图像序列配准和DEM提取联合处理方法。首先,给出了星载全方位SAR几何构型;然后,提出了多角度SAR图像序列联合配准方法,利用小方位夹角图像对间配准偏移量估计大方位夹角图像对间配准偏移量,解决了大方位夹角图像对由于相关性低而难以精确配准的问题;最后,通过构建图像对间配准偏移量与目标高程的关系,提出了多角度SAR图像序列配准和DEM提取联合处理方法,实现配准偏移量和目标高程的同时估计。仿真数据处理结果验证了本文方法的有效性和精确性。