一种基于矢量地图的SAR图像几何精校正方法

吴 戌，叶 伟，劳国超

(1.航天工程大学 研究生管理大队，北京 101416；2.航天工程大学 信息装备系，北京 101416)

由于合成孔径雷达(Synthetic Aperture Radar，SAR)系统的工作模式为侧视成像，并且受到传感器姿态误差、平台星历误差、地形起伏等因素的影响，SAR图像会存在不同程度的几何形变[1]。通常利用基于R-D模型[2-4]的校正法来实现系统级的几何校正。基于R-D模型的校正方法受到目标距离误差、平台星历误差以及目标高度误差等[5]影响，精度较低，不能满足高精度图像应用需求。通过实测地面控制点(Ground Control Point，GCP)修正R-D模型参数可提高几何校正精度，但这种方法需要平台星历数据、成像参数以及R-D方程解算的具体细节，较复杂，并且人工选取控制点费时费力，该方法应用受到限制[6]。所以，寻找精度较高且容易获取的控制信息，提高校正精度同时减少工作量是几何精校正中亟待解决的一个问题。

矢量地图精度较高且来源广泛，可以作为SAR图像几何精校正的控制信息。提取SAR图像与矢量地图中的道路线特征，利用基于特征的图像匹配方法对二者进行匹配，可实现SAR图像的几何精校正。本文首先提取SAR图像的道路线特征，并对其进行参数化描述；然后建立线特征匹配的最优化模型[8]并设计模型的优化求解算法。

1 SAR图像道路线特征提取

矢量地图是对地物形状及位置的记录，不具有地物的灰度或色彩等信息，主要包括点实体、线实体、面实体[9]等，具体如道路线特征，区域边缘线特征等。本文利用道路线特征实现二者的匹配。图1与图2为待精校正高分辨率SAR图像与对应地区1∶5 000比例尺成图的矢量地图。

图1 某地待精校正高分辨率SAR图像

图2 对应地区1∶5 000成图矢量地图

1.1 SAR图像预处理

在高分辨率SAR图像中，建筑区域主要表现为亮线和亮点的集合[10]，由于其灰度值明显高于其他周围地区，此区域的边缘强度值较高。本文在进行边缘检测前，首先利用该图割方法[11-12]对建筑区域进行分割与限幅处理；然后采用直方图均衡化技术[13]，增强边缘信息。预处理结果如图3所示。原始图像中的建筑区域得到明显抑制，且道路边缘得到增强。

图3 经预处理后的SAR图像

1.2 SAR图像道路边缘检测

SAR图像与光学图像不同，由于其固有的相干斑噪声，使得道路呈现多边缘形态[14]。为克服这一影响，本文利用ROEWA(Ratio of Exponentially Weighted Averages，指数加权均值比率)算子[15]联合OTSU最大类间方差[16]自适应阈值分割对SAR图像进行边缘检测，边缘检测结果如图4所示；然后针对高分辨SAR图像边缘检测结果中虚警较多的问题，根据道路边缘与非道路边缘几何特征的差异，将边缘的面积[17]与最小外接矩形的改进长宽比[18]作为过滤因子对非道路边缘进行过滤。对边缘检测结果进行过滤并经数学形态学[19]腐蚀运算后的图像如图5所示，与图4相比道路边缘得以保留而其他地物的边缘被有效过滤。

图4 ROEWA算子联合OTSU边缘检测

1.3 SAR图像边缘线特征提取

为完成后续的匹配工作，需要将SAR图像中的边缘特征抽象成直线段。由于高分辨率SAR图像中其他地物的干扰，边缘存在着断裂、不规则等情况，本文利用Hough变换对边缘检测结果进行线特征提取，并有效解决了边缘断裂等问题[20-21]。经Hough变换后得到的线特征描述：

line((x1,y1),(x2,y2),l,ρ,θ).

其参数分别为线段的两个端点坐标，线段长度，原点到线段的距离与线段的法线与x轴正向的夹角。Hough变换线特征提取结果如图6所示。

图5 经过滤后的边缘检测

图6 Hough变换线特征提取结果

由于受到道路附近的一些干扰，道路边缘极不规则，并且在高分辨率SAR图像中，道路表现为双边缘特征(如图5的右上方边缘)，因此，在一条道路附近Hough变换将检测出多条线特征，需要将这些线特征拟合为一条线段作为道路的控制线。本文对这些线特征根据参数ρ,θ进行分类拟合。

图7 拟合后的道路特征

图8 SAR图像道路特征流程

2 SAR图像几何精校正问题的最优化模型及求解算法

本文利用SAR图像同矢量地图配准的方式实现SAR图像的几何精校正，因此图像变换模型参数的求解是其关键所在。通常，可以利用最小二乘法求解变换模型参数[22]，但该方法需要在两幅图像上人工选取大量分布均匀的同名点才能得到较为精确的校正结果，工作量大。本文根据SAR图像与矢量地图中提取的线特征基于空间关系的相似性度量，建立一个最优化模型，通过求解最优化模型的极值，得到变换模型参数的最优解，使SAR图像几何校正达到全局最优。

2.1 最优化模型

2.1.1 基于空间关系的相似性度量

在地形平坦，不考虑高程信息的情况下，本文采用图像配准中常用的仿射变换(Affine Transform)模型来实现SAR图像的几何精校正[23]。仿射变换

变换模型参数向量为p=(a,b,c,d,e,f)，(x,y)为SAR图像中某一点的坐标，(x′,y′)为仿射变换后对应点的坐标。

当SAR图像与矢量地图中两个线段相匹配时，SAR图像中的线段端点经过仿射变换后应落在矢量地图中对应线段所在的直线上，由此利用线特征的参数ρ,θ形式来定义二者距离。

距离d就是线特征基于空间关系的相似性度量。求解SAR图像所有线段与矢量地图所有线段的之间距离，构成相似性度量矩阵D。D(i,j)表示SAR图像中的第i个线段与矢量地图的第j个线段的相似性度量。D满足关系式：

2.1.2 最优化模型的目标函数

为保证参与计算的特征均是正确匹配的特征，提高计算精度，引入一个特征关系矩阵M，控制可参与计算的特征对。

特征关系矩阵M的定义如下，M(i,j)表示SAR图像的第i个特征与矢量地图的第j个特征之间的关系，M(i,j)=1表示特征匹配，0为不匹配。由于SAR图像的一个特征至多只与矢量地图中的一个特征相匹配，并且矢量地图中的一个特征至多只与SAR图像的一个特征相匹配。因此M矩阵满足如下约束条件。

因此，最优化模型的目标函数[8]，通过引入高斯函数，使得相似性度量矩阵与特征关系矩阵起到了相互约束的作用：当M矩阵存在错误的匹配特征对时，高斯函数使得对应的相似性度量趋近于0；而当特征对(i,j)不是匹配的特征对时，M矩阵使之不参与计算，提高求解精度。

2.2 优化求解算法

最优化模型的目标函数中有两组未知的变量M与P。可通过交替固定其中一组，优化另一组的双步迭代策略求解，直至目标函数收敛。

首先估计变换模型参数作为初始值开始迭代求解。6参数仿射模型至少需要3对同名点来求解变换模型参数。根据图2与图7，选取SAR图像与矢量地图中相对应3组线段交叉点坐标来估计仿射变换模型参数的初值P0。利用P0求相似性度量矩阵D。

当相似性度量矩阵D已知，方程E变化为

E的未知参数为N1×N2个0-1变量，其最优化问题为0-1整数规划问题，可利用解决0-1整数规划问题经典的原始对偶内点法(Primal Dual Interior Point Method)[24]进行求解。由于M缺少等式约束，直接进行求解将始终得到解的边界值。因此本文定义一个等式约束条件，并给出了该约束条件n的求解方法。

由于0-1规划问题解决的是求最小值问题，因此求E的最大值转化为求-E的最小值。

计算-E(1)，由于当n=1时，只有一个M值为1，则其对应的D′应满足D′→-1，若不满足则证明初值选取误差过大需要重新选取初值；

循环计算ΔE=-E(i)-(-E(i-1)),i=2,3,4，…,n；

若|ΔE|<ε，停止循环，n=i-1。本文选取ε=0.5。

该算法的思想是首先选取匹配效果最好的一对特征计算其目标函数值作为初值，每次循环中，满足不等式约束的条件下，令n=n+1，即特征对的个数加1，若增加的特征对在当前仿射变换模型参数下匹配效果较好，则|ΔE|应无限接近于1或等于1，若增加的特征对匹配效果不好或是错误匹配特征对，则|ΔE|将趋于0。因此，通过|ΔE|的大小可以选出在当前仿射变换模型参数下最优的n个匹配特征对得到M矩阵。

若已知特征关系矩阵M，P表示对SAR图像线段特征的端点进行6参数仿射变换。

此时，E的未知参数为仿射变换模型的六参数，其最优化问题为一个无约束的多参数函数最优化的问题。本文利用求解无约束优化问题经典的L-BFGS算法(Limited-Memory BFGS)进行迭代求解。

利用原始对偶内点法与L-BFGS算法对目标函数进行迭代计算，设定迭代求解的终止条件为‖p1-p0‖<ε，p1和p0代表后一次迭代和前一次迭代求得的变换模型参数。通常ε取0.01～0.05。当变换模型参数收敛时，目标函数取得极值，几何校正效果达到最优。

几何校正模型参数求解过程如图9所示。

图9 优化求解流程图

最终求得的仿射变换模型参数如表1所示。

表1 仿射变换六参数

3 几何校正精度分析

本文从几何精校正后SAR图像与矢量地图的图像融合效果与在SAR图像与矢量地图中选取检验点计算相对误差两方面分析提出算法的精度。

在求得变换模型的参数后，对SAR图像进行仿射变换并对变换后的图像进行线性插值，将其绘制在矢量地图上如图10所示，融合效果较好。

图10 校正后的SAR图像与矢量地图融合结果

选取若干分布均匀的检验点计算SAR图像几何精校正的相对精度。检验点为SAR图像与矢量地图中相对应的12个道路交叉点，如图11所示中圆圈区域内的道路交叉中心，按从左至右，从上至下的顺序编为1-12号。同时，在矢量地图上选择对应的道路交叉点作为其同名点，这12个同名点的经纬度坐标利用ArcGIS软件读取。

计算SAR图像检验点经仿射变换后与矢量地图中的同名点在WGS84坐标系下的相对距离，其中φ1,λ1为SAR图像中某一检验点的经纬度，φ2,λ2为矢量地图中同名点的经纬度，Δ为两点之间的距离。表2第二栏为只进行系统校正的SAR图像中检验点的相对误差；表2第三栏为在系统校正的基础上采用本文方法对SAR图像进行精校正后检验点的相对误差。

Δ=111.199[(φ1-φ2)2+

图11 检验点分布

由表2可以得出，本文提出的精校正方法相对平均误差可达到10.49 m，与系统级几何校正相比，精度提高了约50 m左右。但本文的校正方法仍存在一定误差，并且误差较为分散，最大误差和最小误差相差15 m左右。这是由于SAR图像道路线特征提取得到的控制线不是完全意义上的道路中心线，在角度和位置上均存在一定偏差，导致后续的校正存在一定误差；并且本文提出的校正方法是一种全局最优方法，全局性的最优必然会存在某一位置上的误差相对较大的问题。从整体上来看，本文提出的精校正算法对比系统级校正在精度上较有了明显提高，并且与文献[3]利用控制点修正变换模型参数的校正精度相当。

表2 系统校正误差与本文精校正方法误差对比

4 结 论

本文研究基于矢量地图的SAR图像几何精校正方法，该方法利用来源广泛且精度较高的矢量地图，通过自动提取SAR图像中的道路线特征，利用线特征基于空间关系的相似性度量构建最优化模型，并设计优化算法计算模型参数的最优解，实现SAR图像的几何精校正，有效解决在SAR图像的几何精校正中人工选取控制点费时费力以及R-D模型参数解算复杂的问题。精校正后的平均相对误差在10 m左右，较系统级校正精度有明显提高，并且与利用控制点修正R-D模型的几何校正精度相当，满足后续的高精度目标定位等需求。