面向高分三号全球正射影像生成的无控定位精度提升方法

李 欣,蒋博洋,汪韬阳,张 过,崔 浩,程 前

1.武汉大学测绘与遥感信息工程国家重点实验室,武汉 430079; 2.武汉大学遥感信息工程学院,武汉 430079

合成孔径雷达卫星能进行全天时全天候对地成像,且具备一定地表穿透能力,广泛应用于海洋、水利、气象、测绘和军事等方面。在全球范围内进行空间环境研究和地表覆盖制图的最有用方法之一是将SAR卫星影像与光学影像相结合,国内很多研究团队已经研制了各类光学/SAR数字正射影像图(digital orthophoto map,DOM)[1-4],但目前仍缺少基于国产SAR影像进行全球尺度、快速高精度正射影像生成的技术方法和产品研制。GF-3是我国首颗C频段多极化高分辨率SAR卫星[5-6],于2016年8月10日在太原搭载长征四号丙运载火箭成功发射,该卫星系统具备高分辨率、大成像幅宽、多成像模式和高质量成像特点。其中精细条带二(fine strip map model Ⅱ,FSⅡ)模式影像数据具有兼顾大幅宽和高分辨率的特点,可满足全球尺度的测绘制图需求[7-8]。

在星载SAR影像的大区域制图方面,SAR影像匹配和区域网平差方法是消除区域影像几何定位误差的关键核心技术。通常情况下,引入高精度地面控制点是提升区域网平差相对和绝对几何精度的前提条件。文献[9—15]针对多源卫星影像的高精度几何定位开展研究,试验结果表明高精度地面控制信息在区域制图中仍然不可或缺。但是,在境外、无人区或纹理匮乏区域,地面控制点难以获取,因此,在缺少地面控制点条件下的大区域SAR影像无控定位精度提升是亟须解决的技术难题。面向全球尺度的星载SAR正射影像生成,如何摆脱对地面控制点依赖实现高精度绝对定位是一大难点;其次,面向大区域SAR影像数据匹配,如何保障快速高精度可靠的影像连接点提取;然后,针对大区域SAR影像区域网平差,需要通过快速稳健大规模平差模型解算,实现影像间高精度相对定位;最后,海量数据的规模化全流程处理依赖于高效率的多节点并行计算策略。

本文提出的方法可为大区域高精度和高效率的SAR影像制图提供参考,与以往的研究相比,有以下创新点:①采用几何再定标来解决SAR影像初始定位精度低的问题,通过标定SAR影像成像时在距离向和方位向的时间项误差,提升SAR影像的初始无控定位精度;②针对大区域SAR影像提出了一种基于级联匹配的SAR影像快速匹配策略(SAR-SGFM),以及无控大规模自适应区域网平差方法,保障大区域SAR影像间的高精度相对定位;③全流程通过采用多节点并行和GPU/CPU协同提高SAR影像匹配、平差、正射纠正处理效率。

1 星载SAR影像对地定位误差分析

星载SAR通过计算雷达观测信号的传播时间能够获取天线相位中心与地面目标点之间的距离,地面目标在SAR影像方位向上的位置与平台飞行记录的时间有关,在距离向上的位置与地面目标回波到达接收天线的时间先后有关。图1为SAR系统工作示意,结合距离多普勒模型和SAR影像的成像过程分析可知SAR影像的定位误差来源可以分为4类:SAR载荷、卫星平台、观测环境和地面处理[16]。根据误差随时间和空间的变化特性,将误差分为动态误差和静态误差两类,动态误差主要包括卫星平台、观测环境和地面处理引入的误差。静态误差主要包括SAR载荷引入的误差。

图1 侧视成像雷达系统

卫星平台引入的误差一般通过GNSS定位系统来减小,经过精密定轨技术[17]实现厘米级的定轨精度,对SAR影像几何定位精度影响非常小,可忽略不计。观测环境引起的误差主要为雷达信号在传播路径中由于大气环境(图2)所导致的传播延迟[18],这种影响随传播路径的增加而变大,且随大气环境的变化而变化。SAR信号传播路径的改变影响着距离向测量精度,导致几何定位精度受到严重影响,因此需要通过引入大气传播延迟模型补偿该误差项[19]。地面处理引入的误差指在进行SAR影像对地定位时地形高程误差引起的平面投影误差,可通过引入高精度的数字高程模型(digital elevation map,DEM)进行消除。

图2 大气延迟效应

SAR载荷引起的误差包括系统时延和方位向时间误差,方位向时间误差主要是指星载SAR影像辅助文件中各扫描行起始时间与实际系统记录时间的同步误差,是较为稳定的系统量。因此可通过几何再定标对相关参数进行估计并补偿。

综上,影响星载SAR正射影像无控定位精度的主要因素为观测环境和SAR载荷引起的几何定位误差,在几何再定标过程中,可以通过引入大气传播延迟模型[18]和几何再定标模型来提升“无控”绝对定位精度。通过高精度SAR影像匹配算法和区域网平差算法,提升区域SAR影像间“无控”相对定位精度。最后,在正射影像生成时需要引入一个外部高精度全球DEM数据来降低地形起伏处高程误差引起的平面投影误差。因此本文提出的“无控”定位精度提升方法,主要是在星载SAR影像几何处理模型参数精化阶段,实现不依赖地面控制点的定位精度提升。

2 超大规模SAR影像无控定位精度提升方法

2.1 整体架构

本文的总体技术流程如图3所示。首先通过几何再定标有效消除绝大部分系统定位误差,提高星载SAR影像的绝对几何定位精度;其次,基于改进的SAR影像配准算法SAR-SGFM逐级匹配提取可靠的连接点;然后,利用自适应区域网平差方法实现大规模平差解算,从而实现超大规模SAR影像“无控”定位精度提升。最后,通过目前较为成熟的正射纠正、强度一致性处理和镶嵌技术实现全球SAR正射影像生成。

图3 无控制点支持的SAR正射影像制作流程

2.2 GF-3几何再定标

根据初步精度验证结果[7],GF-3 FSⅡ模式的影像(10 m分辨率)的初始定位精度约为30 m,难以直接满足全球大区域图像处理的精度要求。文献[20]提出了星载SAR定标模型和4个定标参数:距离时间、方位角时间、脉冲重复频率和采样频率。文献[21—22]进行了多模混合几何再定标试验,考虑了大气传播延迟和基于提出的4个几何再定标参数的几何交叉校准方法,试验结果[23]表明经过几何再定标后,SAR影像的定位精度可以达到1个像素。文献[19,21]验证了成像时间跨度为5个月的多时相SAR数据,证明了GF-3 SAR卫星的几何再定位性能相对稳定。根据第2节中星载SAR影像对地定位误差分析可知,影响SAR影像几何定位的系统误差项可以通过几何再定标来消除。系统误差包含两个方面:①SAR载荷引入的误差;②观测环境引入的误差。其中SAR载荷引入的误差包括系统时延和方位向时间误差,观测环境引入的误差主要为大气传播延迟引起的误差。本文通过采用顾及大气传播延迟改正的几何再定标方案提升几何定位精度[24]。

几何再定标模型的本质是建立顾及大气时延的SAR系统时延误差和方位向时间误差的补偿模型,具体形式如下

(1)

式中,tf、ts分别为距离向的快时间和方位向的慢时间;tf0、ts0分别为距离向起始时间的测量值和方位向起始时间测量值;tdelay为大气传播延迟时间;Δtf、Δts为系统时延误差;fp为脉冲重复频率;fs为距离向采样频率;x、y为像素坐标,width、height为SAR影像的宽和高,x∈[0,width-1],y∈[0,height-1]。

通过距离多普勒定位模型的间接定位算法,结合几何再定标模型就可以建立起地面点坐标与方位向时间和星地距离的对应关系,而方位向时间和星地距离又是与SAR影像方位向行号和距离向列号一一对应的,在定标景SAR影像内提取角反射器位置作为控制点,通过最小二乘平差就可以精确计算方位向和距离向起始时间改正值,基于解算得到的几何再定标参数(方位向开始时间误差、距离向起始时间延迟误差、脉冲重复频率、距离向采样频率),提升星载SAR影像数据无控条件下的绝对几何定位精度。

有理多项式模型(rational function model,RFM)由于其计算效率较高且能够高精度拟合距离多普勒方程模型[23,25-26],越来越多地应用于星载SAR影像几何处理。基于文献[23,25—26]求解方法来获得SAR影像的(rational polynomial coefficient,RPC)模型参数,其拟合精度优于0.01像素,说明RFM模型可代替距离多普勒方程模型进行后处理。

2.3 星载SAR影像匹配

本文在SAR-SIFT基础上改进,通过采用类梯度位置方向直方图(gradient location and orientation histogram,GLOH)的描述子构建方法建立SAR影像特征描述符,通过快速样本一致和几何模型约束的级联匹配方法,实现区域SAR影像高精度配准,本文将其命名为SAR-SGFM。星载SAR影像匹配技术路线如图4所示。

图4 SAR影像匹配技术路线

首先基于RFM模型构建影像拓扑结构关系及重叠区域估计,基于几何约束下进行重叠区域估计,然后通过改进的SAR-SIFT算子,实现相邻影像特征提取,通过邻近匹配方法获取初始匹配对后利用快速样本一致性算法实现特征点的粗匹配,最后建立基于几何模型约束的粗差剔除方法进行匹配精化。在获取每个分块区域的配准后,逐块处理直至所有匹配任务结束。

2.3.1 SAR影像特征提取

在SIFT[27]的基础上,文献[28]通过改进像素点提取计算方式和改进最终描述子生成方式,提出了性能更好的SAR-SIFT算法。SAR-SIFT构建SAR-Harris尺度空间实现极值检测,通过方向直方图的梯度大小和方向进行特征描述。

SAR-Harris矩阵和响应值为

(2)

RSAR(x,y,α)=Det(HSAR(x,y,α))-

d·Trace(HSAR(x,y,α))2

(3)

为进一步保持局部特征点的稳定性,使用Laplace空间(图5)进行特征点提取,进一步加强了局部特征的仿射不变性,Laplace响应值为

图5 Laplace空间极值点检测示意

图6 原形邻域和对数极坐标网格

F(x,y,α)=α2|Lx(x,y,α)+Ly(x,y,α)|

(4)

将高于SAR-Harris响应阈值的点在Laplace的3×3×3空间邻域内进行峰值检测,进一步加强了局部特征的仿射不变性。对于特征描述符的构建,本文在SAR-SIFT的圆形描述符基础上使用基于GLOH的半径为12圆形邻域和有17个子区域的对数极坐标网格生成特征描述子[29]。

2.3.2 SAR影像特征匹配

(1) 基于快速样本一致算法的点匹配。最近邻匹配法通过在特征空间内两个影像对应特征点的相似度实现匹配。首先需要求取两图对应特征点特征向量的欧氏距离,欧氏距离值最小的两点可定义为匹配点。但由于特征点较多特征向量较为复杂,单一欧氏距离计算可能出现特征点相似特征干扰,为降低干扰影响,除比较欧氏距离最邻近点外,还需要分析第二邻近点的欧氏距离大小,当两者数值比值小于给定阈值dratio时,就判定该特征点对正确,二维空间特征点对欧氏距离计算公式如下

(5)

本文从参考图像和待配准图像中分别提取两个点特征集{p1,p2,…,pm}和{q1,q2,…,qn},并从两组点中找出最佳候选点构建两景影像的对应关系,将这种对应关系定义为ti={pi,qi},pi的坐标(x1,y1)表示参考图像中的一个特征点,qi的坐标(x2,y2)表示待配准图像中的一个特征点。T={t1,t2,…,ti,…,tn}表示临时对应点集。近邻匹配法处理后仍会出现部分误匹配点对,通过随机采样一致性算法(RANSAC)进行误匹配点对剔除。

本文基于RANSAC算法提出一种快速样本一致性算法以提高算法可靠性和效率。从集合T中提取正确匹配率高的子集Th进行快速样本一致性处理,然后在集合T中找到最大的共识集。提取子集Th的过程如下。本文设置了两个dratio参数dmax和dmin,dmin是一个相对较小的值,Th中的对应关系通过dmin的比例进行匹配;dmax是一个相对较大的值,T中的对应关系通过dmax的比例进行匹配。Th表示本文取样的样本集,T表示本文找到最大共识集的共识集。最大共识集中的对应关系通过其转换误差来判断,ci={pi,qi}的转换误差表示为

式中,θ是由Th的样本决定的转换模型参数。本文方法对参考影像和待配准影像进行精匹配,相较于RANSAC算法,本文方法的计算代价和稳定性均较优。

(2) 基于几何模型约束的精匹配。建立星载SAR影像的RPC模型[23],形式如下

(6)

式中,(P,L,H)为标准化的地面坐标;(X,Y)为标准化的影像坐标。

基于RPC模型实现坐标正反变换,选取基准影像Fbase和待匹配影像Fmatch,分别将基准影像和待匹配影像上的同名点计算其地面点坐标,当差值小于某一阈值即为同名点,否则就视为粗差点

Fbase(Latitude,Longitude,Height)=

T(line,sample)

(7)

Fmatch(Latitude,Longitude,Height)=

T(line,sample)

(8)

分别获取其对应地面点坐标后,二者差值Δ,若满足阈值则为同名点,否则为误差点

(9)

式中,1为同名点;0为误差点。

在山区,由于SAR影像透视收缩、叠掩、阴影等问题引起的影像局部纹理信息丢失导致匹配困难,匹配时通过叠掩自动定位算法[30]规避叠掩区域。对于缺少纹理特征的沙漠区域,直接对几何再定标后的SAR影像进行正射纠正,镶嵌时可通过强度一致性处理和边缘羽化来进行辐射平滑。

2.4 区域网平差

2.4.1 平差模型构建

构建大区域自适应区域网平差模型,结合基于像方仿射变换的系统误差模型,提升区域SAR影像间的相对定位精度,形式如下

(10)

式中,(sample,line)为正则化的像方坐标;(e0,e1,e2,…,f0,f1,f2,…)表示系统误差补偿参数,即为平差所需求解的未知数。

对未知数全微分构建误差方程如下

(11)

式中,Fx表示方程在列方向的误差;Fy表示方程在行方向的误差;Fx0、Fy0表示泰勒级数0次项;e0、e1、e2和f0、f1、f2为仿射变换参数;line、sample为影像行列坐标;x、y为精化后的影像行列坐标;∂Fx/∂lat、∂Fx/∂lon、∂Fx/∂h分别表示纬度、经度和高程方向上的一阶偏导数。

误差方程的矩阵形式为

V=Bt+AX-lP

(12)

式中,B为仿射变换改正数系数矩阵;t为仿射变换改正数;A为地面点坐标改正数系数矩阵;X为地面点改正数;l为常数项;V为观测值残差向量。

2.4.2 大规模平差解算

误差方程法化得到法方程

(13)

式(13)可以表示为

(14)

式中,N矩阵为外方位元素改正数t的系数矩阵;M为地面点坐标改正数X的系数矩阵;K为协方差矩阵;L1和L2分别为外方位元素和地面点坐标改正数的常数向量。当进行大规模区域网平差时,整个区域内影像连接点数量常常到千万甚至上亿量级,法方程矩阵规模巨大,如果直接求解不论是内存开销还是运行效率都无法满足快速求解的需求。本文通过采用消元改化法方程的策略来进行平差解算。

如图7所示,法方程结构分为3个部分:左上角是外方位元素的系数矩阵,右下方占据绝大多数空间的是地面点坐标项,两边对称分布是协方差矩阵。由于其形状为带状,也被称为带状矩阵,可利用法方程特有的稀疏结构进行优化求解。由于大规模区域网平差过程中,地面点的数量远远大于影像的数量,可以考虑消去地面点坐标X,构建包含外方位元素改正数t的方程

图7 星载SAR影像大规模区域网平差法方程的稀疏结构

[N-KM-1KT]t=L1-KM-1L2

(15)

本文引入Schur消元法进行了法方程矩阵改化[31],相比于直接矩阵求逆,其优势在于:

(1) 在消元过程中,M矩阵由对角的多个小方阵构成,因此对于M-1的结果可以由对内部每个小方阵求逆后组合而成,这个过程可并行化计算,大大节约内存且显著提升计算效率;

(2) 在求解了t后,地面点坐标的改正数X可由式(16)快速求解

X=M-1(L2-KTt)

(16)

3 试验结果与分析

3.1 试验设计

为了验证本文方法的有效性,通过设计系列试验开展验证,具体如下。

(1) 几何再定标试验:通过随机选择多个试验区域不同成像时间的SAR影像数据进行几何再定标处理,通过高精度控制数据验证几何再定标前后的几何定位精度提升效果。

(2) SAR影像匹配试验:利用本文提出的SAR-SGFM方法在不同地形条件(平地、丘陵、山地、高山地)下进行SAR影像匹配试验,验证匹配方法的配准效率和精度。匹配精度通过中误差(root-mean-square error,RMSE)、匹配点数、正确匹配点数来统计,配准效率通过时间进行评估。

(3) 区域网平差试验:根据本文提出的大区域影像平差方法实现区域影像的相对几何定位精度提升。区域网平差的精度评价通过连接点中误差进行评估。

(4) 正射影像精度评估:通过高精度地面控制作为检查点进行正射影像的几何精度评价,统计检查点物方坐标残差的中误差作为精度指标,计算方法如下

(17)

3.2 试验数据

本文利用GF-3 FSII模式影像数据进行了全球正射影像生成,系统地验证了本文方法的效果。整个测区跨越了全球南北纬70°之间的陆地和主要岛屿,包括了平原、丘陵、山地、高山地等多种地形。影像分布图如图8所示,数据相关参数见表1。

表1 测区基本参数

图8 高分三号影像全球覆盖

3.3 试验结果

3.3.1 几何再定标精度验证

全球范围内选取了8景影像进行精度验证,中国区域的验证数据为高精度光学正射影像,国外区域的验证数据源为Google Earth底图。数据分布情况如图9所示,表2给出了测试数据的详细参数。

表2 测试数据基本参数

图9 几何再定标精度验证测试数据分布情况

本文统计了每个区域的几何定位精度中误差,图10为各个区域的几何再定标前后定位精度对比图。结合表3和图10分析可知,几何再定标前后GF-3 FSⅡ模式影像几何定位精度从30 m左右提升至10 m以内。

表3 几何再定标精度验证试验结果

图10 几何再定标前后影像几何定位精度对比

根据国内外区域几何再定标试验可知,8个测试区域单景影像几何定位精度均得到明显提升,达到了无控制点条件下优于10 m的几何定位精度。以上结果验证了本文提出的几何再定标方法的适用性和稳定性,可以适用于全球SAR正射影像的几何处理。

3.3.2 SAR匹配算法验证

为了评估本文提出的SAR影像匹配SAR-SGFM方法的准确性和效率,分别在湖北和西藏两个区域开展算法验证。测试环境为一台6核i7-10875 CPU、16 GB内存、GPU为NVIDIA GeForce RTX 2070的便携式笔记本。试验区域相关参数见表4。

表4 星载SAR影像匹配试验数据情况

利用湖北省5景SAR影像,可在4 min内提取大约21 984个可靠的连接点(图11);利用西藏4景SAR影像,可在2 min内提取大约14 590个可靠的连接点。

依据匹配到的连接点通过区域网平差方法,验证匹配连接点精度见表5。湖北区域的21 984个连接点的中误差为0.658像素,西藏区域的14 590个连接点的中误差为0.871像素。两个区域的连接点中误差均优于1个像素,连接点的匹配效率约为30 s/景。试验结果表明本文提出的SAR-SGFM影像匹配方法可以适用于大规模SAR影像匹配。

此外,为了进一步验证SAR-SGFM匹配算法的稳健性,本文在测试区域内,选择了2000×2000大小的平地、丘陵和山区等3种不同地形条件的局部影像对进行匹配试验。

由表6结果可知,在不同地形条件下,本文的匹配方法均可以实现亚像素级的配准精度。

表6 不同地形下的SAR影像匹配结果

3.3.3 区域网平差试验验证

对全球覆盖的GF-3 FSII模式影像进行大区域试验验证,按照亚欧大陆、非洲、大洋洲、南美洲和北美洲5个不同区域的SAR影像区域网平差结果见表7。

表7 全球分区域GF-3区域网平差试验结果

由表7可知,不同区域平差后连接点中误差均优于1个像素,保证了影像间良好的接边精度。亚欧大陆区域8581景影像的连接点数量为40.7万,连接点中误差为0.763像素,平差解算时间为84.60 s;非洲区域3911景影像的连接点数量为18.6万,连接点中误差为0.738像素,平差解算时间为37.25 s;大洋洲区域1293景影像的连接点数量为6.2万,连接点中误差为0.823像素,平差解算时间为12.81 s;南美洲区域2318景影像的连接点数量为11.0万,连接点中误差为0.598像素,平差解算时间为21.20 s;北美区域3734景影像的连接点数量为17.7万,连接点中误差为0.451像素,平差解算时间为33.40 s。根据统计结果,平差解算效率约为3 s/万连接点,平差精度均优于1个像素,可满足后续区域影像无缝镶嵌的几何相对定位精度要求。

本文通过自主研发的SAR正射影像生产软件,并引入30 m格网的AW3D DEM数据[32]完成全球正射影像生产,正射影像生成效率优于50 MB/s(较CPU方法效率提升10倍),然后采用SAR影像强度一致性处理算法[33]进行匀色,最终研制的全球SAR正射影像产品如图13所示。

图13 全球SAR正射影像产品

针对SAR正射影像国内区域,采集273个地面控制点进行精度验证,统计检查点总体中误差为8.01 m。针对国外区域,通过日本部分地区DOM产品(平面精度±2 m)均匀选择100个检查点进行精度验证,统计检查点总体中误差为9.739 m。针对全球其他地区的精度验证,通过Google Earth影像采集检查点进行精度验证,其绝对定位精度均优于10 m。

4 结 论

本文面向GF-3 SAR影像全球正射影像生成过程所遇到的技术难题,系统提出了一套全球星载SAR影像无控几何精度提升方法,完成了全球主要陆地和岛屿10 m分辨率SAR正射影像产品研制。

(1) 通过几何再定标实现了单景影像几何定位精度提升,通过不同区域的数据验证结果可知,实现了GF-3 FSⅡ模式数据初始定位从30 m提升至10 m,可以为全球正射影像生产奠定较好的几何基础。

(2) 通过本文提出的星载SAR影像匹配和平差方法,实现各区域平差后连接点中误差均优于1个像素,保障了相邻SAR影像几何无缝的拼接效果。通过真实数据验证,全球SAR正射影像产品绝对定位精度优于10 m(1个像素)。

(3) 基于多节点并行计算集群系统,影像匹配平均时间为30 s/景,区域影像平差平均解算效率为3 s/万连接点,正射纠正效率优于50兆/秒,显著提升了SAR正射影像生成效率。

本文提出的无控定位精度提升方法有效支撑高分三号全球SAR正射影像的规模化生产,可为多源卫星载荷的大区域制图提供方法参考。下一步工作将重点研究基于迭代平差的全球数据更新,解决现有正射影像快速更新问题,以及面向全球的SAR正射影像的叠掩补偿方法解决现有正射影像信息缺失问题。