圆迹SAR子孔径图像序列联合相关DEM提取方法

张金强，索志勇，李真芳，保 铮

(西安电子科技大学雷达信号处理国家重点实验室，陕西西安710071)

0 引 言

合成孔径雷达(synthetic aperture radar, SAR)获取观测场景地形高程模型(digital elevation model, DEM)主要通过两种技术途径:雷达立体摄影技术[1-2]和雷达干涉测高技术[3-5]。基于不同雷达下视角下获取的两幅SAR图像,雷达立体摄影技术利用视差提取目标高程,而雷达干涉测高技术利用相位差提取目标高程。由于单一照射方向获取的SAR图像存在叠掩或阴影区域,导致利用上述技术途径无法获取观测场景全方位DEM[6-7]。

圆迹SAR(circular SAR, CSAR)随雷达平台以观测场景为中心作360°圆周运动,波束始终照射同一地面场景,实现对其全方位观测[8-11]。利用CSAR可以提取观测场景全方位DEM[12-14]。文献[12-14]提出了基于CSAR子孔径图像序列的DEM提取方法。首先,将360°圆环划分为多段圆弧,每段圆弧划分为多个子孔径,对各子孔径回波进行聚焦成像处理,获取子孔径图像序列。然后,针对每段圆弧,利用圆弧内部分子孔径图像序列间的相关性估计观测场景DEM。文献[12]利用圆弧内各子孔径图像与中心子孔径图像间相关系数的和值作为测度函数,每段圆弧提取一幅DEM。文献[13-14]利用圆弧内各子孔径图像与中心子孔径图像间的相关系数分别作为测度函数,获取多幅DEM后进行平均,最终获取一幅DEM。相比于文献[13],文献[14]在对圆弧内各子孔径图像与中心子孔径图像进行配准时,采用了金字塔分层搜索的方法,提高了配准的稳健性。由于叠掩或阴影的影响,不同段圆弧提取的DEM中存在不同的无效高程值区域。最后,融合不同段圆弧DEM获取全方位DEM。

针对每段圆弧,文献[12-14]只利用了各子孔径图像与中心子孔径图像之间的相关性信息,非中心子孔径图像之间的相关性信息没有得到利用。为了解决上述问题,本文提出了联合相关法。该方法将圆弧内所有子孔径图像之间的联合相关系数作为测度函数。该测度函数利用了圆弧内全部子孔径图像之间的相关性信息,能够提高测度函数灵敏度和DEM提取精度。同时,将子孔径图像序列向三维空间投影以校正其几何形变,提高子孔径图像序列之间的相关性。本文首先对CSAR子孔径图像序列进行分析,然后给出联合相关法处理流程,最后利用CSAR实测数据验证所提算法的有效性和精确性。

1 CSAR子孔径图像序列分析

CSAR成像几何模型如图1所示。根据分辨率要求确定子孔径方位角宽度,然后将360°圆环划分为多个子孔径,利用后向投影(back projection, BP)算法[11]对各子孔径回波进行聚焦成像处理,获取CSAR子孔径图像序列。不失一般性,假设各子孔径图像成像于地面坐标系。

图1 CSAR成像几何模型Fig.1 CSAR imaging geometry

成像平面高程与目标实际高程不一致,将导致目标在子孔径图像上的平面位置与实际不一致,引起子孔径图像发生几何形变[12]。如图1所示,假设子孔径图像#i对应的雷达轨迹中心为Ci,雷达下视角为θi,方位角为φi,观测场景中一目标点P实际高程与聚焦平面高程之差为Δh,则其在子孔径图像#i上的平面位置相对于实际位置偏移量为

(1)

根据式(1)可知,目标实际高程与成像平面高程差值越大,子孔径图像几何形变越大。不同子孔径图像几何形变不同,与其雷达下视角和方位角有关。如图1所示,假设子孔径图像#j对应的雷达轨迹中心为Cj,雷达下视角为θj,方位角为φj,目标点P在子孔径图像#j上的平面位置相对于其在子孔径图像#i上的位置偏移量为

(2)

根据式(2)可知,子孔径图像间方位夹角增大,导致子孔径图像间几何形变变化增大,引起子孔径图像间相关性降低。同时,目标散射特性随方位角的变化,也将导致子孔径图像间相关性降低[13]。因此,子孔径图像间相关性随着方位夹角的增大而降低。下面利用美国空军研究实验室(air force research laboratory, AFRL)公开发布的CSAR实测数据(gotcha volumetric SAR data set, Version 1.0)验证该结论。

子孔径图像序列生成过程与本文第3节相同。在子孔径图像序列中根据方位夹角随机选择两幅图像,逐像素估计其相关系数,相关系数估计窗口大小为15像素×15像素。每一方位夹角均随机选择5对子孔径图像,计算所有像素相关系数的平均值,其随方位夹角的变化如图2所示,方位夹角变化范围为3°～42°。根据图2可知,子孔径图像对间相关系数随着方位夹角的增大而减小。

图2 CSAR 子孔径图像对间相关系数随方位夹角的变化Fig.2 Cross-correlation coefficients between subaperture imagepairs of CSAR versus azimuth angle intervals

CSAR子孔径图像序列的相关性将影响DEM提取精度,因此对CSAR子孔径图像序列相关性的分析结果对于下文联合相关DEM提取方法中圆环的划分具有重要参考价值。

2 联合相关DEM提取方法

CSAR具有360°全方位观测能力,但是对于不同的方位观测角,子孔径图像中存在不同的叠掩或阴影区域。为了获取观测场景全方位DEM,将整个圆环划分为多段圆弧,利用每段圆弧内的子孔径图像序列提取观测场景DEM。不同段圆弧提取的DEM中由于叠掩或阴影导致的无效高程值区域不同,融合不同段圆弧提取的DEM即可获取全方位DEM。联合相关DEM提取方法流程图如图3所示,下面对主要步骤进行详细介绍。

图3 联合相关DEM提取方法流程图Fig.3 Flowchart of joint cross-correlation DEM extraction method

联合相关DEM提取方法以CSAR子孔径图像序列作为输入。

为了消除CSAR子孔径图像间几何形变变化对相关性的影响,将子孔径图像向三维空间投影,如图4所示。首先,根据观测场景的先验平面位置和目标高程范围,在地面坐标系下建立三维格网;然后,针对每一格网点根据式(1)进行SAR反向定位,获取格网点在子孔径图像中的像素坐标;最后,通过插值提取该像素的幅值赋予格网点。当格网点对应的高程与目标高程一致时,子孔径图像几何形变得到校正,子孔径图像间相关性最高。

根据第1节分析可知,子孔径图像间方位夹角越大相关性越低,相关性降低将导致DEM提取精度降低。根据式(2)可知,子孔径图像间方位夹角较小时,目标平面位置偏移量对目标高程不敏感,DEM提取精度较低。因此,在进行子孔径图像序列划分时需要根据其相关性选择圆弧方位角宽度。

图4 CSAR子孔径图像向三维空间投影示意图Fig.4 Project subaperture images of CSAR on 3D space

联合相关系数用公式表示为

(3)

式中,M为圆弧内子孔径图像数量;sm(l,k)为所选窗口内图像#m的像素幅度值;μm为相应的幅度平均值;(2L+1)×(2K+1)为所选窗口大小。所选窗口为以待估计像素为中心的矩形窗,并且窗口内其他像素与中心待估计像素的高程值相同。对于每一平面坐标(x,y),沿高程方向(即图4中Z轴方向)计算JC随h的变化。然后选择使JC取最大值时的h作为估计高程,即

(4)

文献[13-14]利用圆弧内各子孔径图像与中心子孔径图像间的相关系数作为测度函数,获取多幅DEM后进行平均,最终获取一幅DEM。两幅子孔径图像间的相关系数用公式表示为

ρ1,m=

(5)

式中,下标1表示中心子孔径图像编号;下标m∈[2,3,…,M]表示其他各子孔径图像编号。文献[12]利用圆弧内各子孔径图像与中心子孔径图像间相关系数的和值作为测度函数,估计得到一幅DEM。

ρ1,m

(6)

对比式(3)、式(5)和式(6)可知,文献[12-14]仅利用了圆弧内部分子孔径图像间的相关性信息,而本文方法利用了圆弧内全部子孔径图像间的相关性信息。

将所有圆弧提取到的DEM进行融合[12],获取全方位DEM并输出。

3 实测数据处理及分析

本文选择AFRL公开发布的CSAR实测数据中航过1、HH极化方式的数据验证所提算法的有效性和精确性。为了对比分析所提算法的性能,同时利用文献[12]中的算法(称为和相关法)和文献[14]中的算法(称为改进和相关法)进行DEM提取。上述实测数据由X波段、640 MHz的雷达系统录取。美国地理学会(US geography society, USGS)获取的观测场景的光学正射影像如图5(a)所示。将360°圆环划分为120个子孔径,每个子孔径对应的方位角宽度为3°。利用BP算法对每个子孔径的回波数据在地面坐标系下进行聚焦成像,成像平面高程为0 m,格网间距为0.2 m,图像大小为501像素×501像素。对所有子孔径图像进行非相干叠加[15]获取的SAR图像如图5(b)所示。

图5 AFRL CSAR系统的观测场景Fig.5 Imaged area of AFRL CSAR system

3.1 测度函数灵敏度分析

利用测度函数主瓣3 dB宽度(称之为主瓣宽度)和第一旁瓣峰值与主瓣峰值的比值(称之为峰值旁瓣比)评价其灵敏度[12]。主瓣宽度和峰值旁瓣比越小,表明测度函数灵敏度越高,DEM提取精度越高。

根据方位角宽度在子孔径图像序列中随机选择若干幅方位角连续的子孔径图像组成一段圆弧。针对每一平面坐标(x,y)计算和相关系数和联合相关系数随高程h的变化,然后分别估计其主瓣宽度和峰值旁瓣比。每一方位角宽度均随机选择5组子孔径图像,窗口大小选为5像素×5像素和11像素×11像素。图6给出了和相关函数和联合相关函数的主瓣宽度和峰值旁瓣比均值随圆弧方位角宽度的变化,圆弧方位角宽度的变化范围为21°～87°。由图6可知,与和相关函数法相比,本文方法测度函数的主瓣宽度和峰值旁瓣比更小,灵敏度更高。

图6 测度函数灵敏度随圆弧方位角宽度的变化Fig.6 Sensitivities of measure function versus azimuthangle intervals of arcs

3.2 DEM提取结果及分析

根据图1和图2的分析结果,选择圆弧方位角宽度为60°,将360°圆环划分为相互重叠的24个圆弧[12]。选择大小为5像素×5像素和11像素×11像素的两种窗口计算相关系数随目标高程的变化,将得到的两组相关系数在对应高程处相乘,然后基于此估计目标高程。上述处理流程能够提高高程估计性能[12]。观测场景内目标高程范围为-2～2 m,格网间隔选择0.2 m。和相关法提取DEM如图7所示,改进和相关法提取DEM如图8所示,联合相关法提取DEM如图9所示。对比图7～图9的结果可知,本文方法提取的DEM中汽车轮廓更清晰。由于观测场景内大部分地物后向散射系数较小、子孔径图像序列间相关性较低,导致DEM提取精度较低。

图7 和相关法提取DEMFig.7 DEM extracted by accumulation cross-correlation

图8 改进和相关法提取DEMFig.8 DEM extracted by improved accumulation cross-correlation

图9 联合相关法提取DEMFig.9 DEM extracted by joint cross-correlation

为了定量评估算法性能,选择图5(b)中所标注的7辆汽车评估所提取的DEM精度,图中箭头指向汽车的车头方向。7辆汽车的图片如图10所示,7辆汽车的实际长、宽和高如表1所示。汽车#F的DEM提取结果如图11所示,选择图11(a)中内框所示的区域评估DEM提取精度。统计所提取目标高程的均值和均方根误差[16],其中均方根误差的计算公式为

(5)

图10 停车场中汽车的图片Fig.10 Photos of vehicles in parking lot

图11 汽车#F的DEM提取结果Fig.11 Extracted DEM of car #F

m

表2 汽车的高程估计结果

4 结 论

本文通过对CSAR子孔径图像序列相关性的分析发现,子孔径图像间的相关性随着方位夹角的增大而降低。基于上述分析,本文提出利用子孔径图像序列间的联合相关系数作为测度函数提取观测场景DEM。该测度函数利用了保持相关性的一段圆弧内所有子孔径图像之间的相关性信息,提高了测度函数灵敏度和DEM提取精度。实测数据处理结果和分析验证了本文算法的有效性和精确性。