余 洋 靳国旺 熊 新 童成功 张国平

(信息工程大学地理空间信息学院 郑州 450001)

1 引言

合成孔径雷达(Synthetic Aperture Radar,SAR)具备全天时、全天候工作的能力，在地形测绘、灾害监测、目标侦察等领域具有突出优势[1]。中国首颗高分辨率C波段多成像模式SAR卫星高分三号(GF-3)，受卫星位置测量误差、传感器授时精度和测距精度影响，影像系统级定位误差约为40 m[2]。利用控制点对几何模型参数进行精化可达到提高定位精度的目的，包括少量高质量地面控制点[3,4]、已有数字正射影像图(Digital Orthophoto Map,DOM)和数字高程模型(Digital Elevation Model,DEM)等高精度地理信息数据中提取的内业控制点[5]。地面控制点布设耗时费力，且不能应用于境外等难以布设区域，内业控制点依赖于地理信息数据自身的精度。

近年来，将星载激光测高数据作为控制数据对光学遥感影像几何模型参数进行精化，可有效提高定位精度，星载激光测高数据全球控制点提取、联合遥感影像定位逐步成为研究热点[6]。冰、云和陆地高程卫星-2 (The Ice,Cloud and land Elevation Satellite-2,ICESat-2)是世界上首颗光子计数激光测高卫星，搭载高级地形激光测高系统(Advanced Topographic Laser Altimeter System,ATLAS)，以10 kHz的重复频率发射激光，可获得相邻光子沿轨距离为0.7 m的全球观测数据[7,8]。其中，名为全球地理定位光子数据的ATL03级产品提供了沿轨光子经纬度、高程、沿轨距离、置信度标识[9]等信息，并以沿轨距离约20 m对接收到的光子进行分组，是其他产品的基础数据[10]。由于ATLAS传感器灵敏度高，且发射波长为532 nm的绿光进行探测，数据中包含较多因大气散射与太阳光噪声引起的噪点，常用不同形状的滤波核进行去噪处理[11]。ATLAS数据光斑大小约为17 m，高程测量精度受光斑内地形坡度、地表植被覆盖影响，在地形平坦、无云且地表植被稀疏等理想情况下高程精度达到0.1 m，与地面控制点高程精度相当，是目前高程精度最高的星载激光测高数据[12,13]。

Rosiek等人[14]、何钰等人[15]、耿迅等人[16]针对行星表面缺乏控制数据，遥感影像测绘产品精度低问题，利用激光测高数据与光学遥感影像进行光束法区域网平差，行星表面地形测绘成果精度得到有效提升，初步验证了星载激光测高数据用于提升光学遥感影像定位精度的可行性。王晋等人[17]、张鑫磊等人[18]、Zhang等人[19]将星载激光测高数据与高分辨率光学立体影像联合处理，通过提取激光点为光学影像做高程控制，有效地消除了几何模型中的系统误差，定位精度得到了大幅提升。谭建伟等人[20]、王密等人[21]、郑迎辉等人[22]分别利用ICESat-1,ICESat-2测高数据，通过属性参数约束等策略筛选，提取出可满足中比例立体测图高程控制要求的激光高程点[23]。曹宁等人[24]、唐新明等人[25]利用中国搭载了激光测高系统的光学立体测绘卫星影像，进行了激光测高数据与遥感影像联合平差实验，为中国1:10000立体测图需求提供了技术支持，综上，星载激光测高数据与光学遥感影像联合应用具有重要价值，但ATLAS数据与SAR影像联合几何处理相关研究较少。

为精化SAR影像几何参数并提高立体定位精度，本文借鉴星载激光测高数据光学遥感影像高程控制点提取思路和秩自相似(Rank-based Self-Similarity,RSS)描述子匹配算法[26]，设计了一种从ATLAS数据中提取SAR高程控制点方法。利用中国登封市和日本横须贺市两个区域的ATLAS数据进行了高分三号SAR高程控制点提取实验。

2 多策略ATLAS数据优选与影像匹配相结合的SAR高程控制点提取

多策略ATLAS数据优选与影像匹配相结合的SAR高程控制点提取流程如图1所示。对ATLAS数据采用多种策略组合优选激光高程点，将激光高程点平面坐标投影至谷歌地球(Google Earth,GE)影像中，选取局部区域作为激光足印影像(Laser Footprint Image,LFI)；利用SRTM DEM对斜距SAR影像进行地理编码，由RSS描述子实现SAR地理编码影像与LFI的匹配，从而建立激光高程点高程值与SAR地理编码影像坐标之间的对应关系，再通过坐标反算得到斜距SAR影像高程控制点。

图1 多策略ATLAS数据优选与影像匹配相结合的SAR高程控制点提取流程Fig.1 Process of SAR elevation control point extraction combining multistrategy ATLAS data preference and image matching

2.1 多策略ATLAS数据优选

多策略ATLAS数据优选流程如图2所示，分为夜间高质量光子提取与大偏心率椭圆平坦区域光子提取两个步骤。夜间高质量光子提取中，根据光子三维坐标计算太阳高度角，从而滤除非夜间观测光子数据；由置信度信息提取位于地面的信号光子；通过比较光子高程值与SRTM DEM高程值剔除粗差。最后，利用大偏心率椭圆滤波核提取位于平坦区域的光子作为激光高程点。算法需要读取的参数如表1所示。

表1 算法使用到的ATL03参数Tab.1 ATL03 parameters used by the algorithm

图2 多策略ATLAS数据优选流程Fig.2 Process of multistrategy ATLAS data preference

(1) 夜间高质量光子提取

为减小太阳光噪声对高程控制点提取的影响，依据光子三维坐标计算的太阳高度角Hθ提取夜间观测数据。Hθ为太阳光线与地平面之间的夹角，该值小于零时可认为数据获取时刻在当地为夜间，计算式为

式中，δ为太阳赤纬，Ω为时角，

N为观测时间与该年第1天之间的天数差，

为由t0计算的真太阳时[27]。

ATL03产品为每个光子提供不同地形下基于泊松算法的信号置信度标识，分别为0～4，在被标记为陆地的数据中提取C大于0的光子。

根据拉依达准则，比较光子高程值与相同平面坐标下SRTM DEM高程值，保留高程差值小于3倍SRTM DEM标准差的光子，即

式中，Hice为光子高程值，HSRTM为相同平面坐标下SRTM DEM高程值，ξ为egm96模型水准面差，σD为SRTM DEM标准差。

(2) 大偏心率椭圆平坦区域光子提取

地面有较大的坡度时，物像不一致问题会造成高程控制偏差问题。如图3所示，S为成像时刻SAR卫星位置，地物A对应的同名像点为a，A匹配到的像点为a′，a′对应同名地物点为A′，a与a′在距离向、方位向上相差 Δy、Δx个像素，当地面存在大小为α的坡度时，会分别引起高程为ΔZ1与ΔZ2的差值：

图3 物像不一致造成的高程差示意图Fig.3 Schematic diagram of height difference caused by difference between object and image

式中，Ry,Rx为SAR影像距离向、方位向分辨率。则 当 Δx=1,Δy=1,α=1.5°,Ry=0.6 m,Rx=0.4 m 时，由式(5)计算可得，将引起高程为ΔZ1=0.0157 m，ΔZ2=0.0105 m的误差。因此，筛选出位于平坦区域的光子即可保证激光点高程精度，也可减弱物像不一致问题造成的高程控制偏差问题。

星载激光测高数据常用空间聚类算法进行地面光子的提取，椭圆形滤波核更符合地面信号光子的分布特征。如图4所示，对于某个光子p，定义以点p为中心的长轴沿轨椭圆形窗口，点q是否在该椭圆内可用式(6)判断：

图4 椭圆滤波核提取分析示意图Fig.4 Schematic diagram of ellipse filter kernel extraction

式中，a,b为椭圆长、短半径，lp,lq为p,q的沿轨距离，由光子对应li与前序所有组Li的和相加得到，hp,hq为p,q的高程，θ为椭圆旋转角度，当d(p,q)≤1时q在椭圆内；即对于任意光子pi，根据式(6)遍历位于以pi为中心的椭圆内光子数W，若W大于某一阈值T，则保留pi。

相对于整体信号光子提取时将椭圆滤波核长短轴之比设为4:1的做法，本文将滤波核长轴a设为25 m，短轴b设为0.5 m，并将滤波核转角θ设为0，以直接提取平坦区域处的地面光子。该滤波核的提取结果如图5所示，此时椭圆滤波核的偏心率为0.9998，所提光子在沿轨方向上的坡度不超过1.15°，能够满足高程控制的要求。

图5 平坦区域光子提取结果示意图Fig.5 Schematic diagram of laser points extraction result in flat area

为了提取平坦区域裸地的光子，以当前光子为中心设置宽度为WS、长度应能覆盖所有光子的矩形窗口，计算窗口内提取的夜间高质量光子高程方差，剔除方差大于阈值δW的光子数据，得到优选的激光高程数据。

2.2 谷歌影像、SAR影像秩自相似描述子匹配

为使激光高程点与SAR影像匹配更具鲁棒性与准确性，本文通过影像匹配的方法实现二者的匹配。如图6所示，将提取出的激光高程点按平面坐标投影至GE中，并与SAR地理编码影像进行匹配，建立SAR影像像点与激光点高程值之间的对应关系，作为高程控制条件引入几何模型参数精化中，从而提升SAR影像定位精度。

图6 GE影像与SAR地理编码影像匹配示意图Fig.6 Schematic diagram of matching between GE image and SAR geocoding image

GE影像与SAR影像的匹配属于异源遥感影像匹配工作，匹配前利用GF-3影像导航信息(几何模型参数)与外部DEM对SAR影像进行地理编码，即将SAR影像纠正到与GE影像同一坐标基准下，从而消除GE影像与斜距SAR影像间的旋转与尺度差异，且可获得含有较小误差的几何对应关系。为了得到精确的匹配结果，利用模板匹配算法进行，但此类算法效率低，可以根据几何对应关系预估激光高程点在SAR地理编码影像上的位置，然后在局部窗口内逐像素计算匹配算子相似性度量结果。

(1) SAR影像地理编码

将图像坐标系的SAR影像纠正到地理坐标系下，即由SAR地理编码影像像点坐标 (x′,y′)，根据分辨率S计算地面点平面坐标(X,Y)，由SRTM DEM内插得到相应高程值Z，利用SAR影像有理多项式系数(Rational Polynomial Coefficients,RPC)模型计算斜距影像像点坐标 (x,y)，根据斜距影像灰度分布内插 (x,y)处的灰度值赋给SAR地理编码影像中的相应像点，逐像点赋值后得到SAR地理编码影像。

(2) 秩自相似描述子匹配

由于平坦区域激光点往往对应影像特征稀疏区域，基于特征的影像匹配结果误差较大，这里使用RSS描述子进行匹配工作，并针对效率和精细匹配需求，进行了减少分区网格数量和基于二次多项式的亚像素坐标计算的改进。如图7所示，将激光高程点投影至GE影像中某个像元，以该像元为中心，大小为(2R1+1)×(2R1+1)的窗口提取LFI(模板窗口)；将激光高程点投影至SAR地理编码影像中，选取大小为(2R1+1)×(2R1+1)和(2R2+1)×(2R2+1)的区域为匹配窗口与搜索窗口。影像秩表面RSq(x,y)定义为

图7 匹配窗口示意图Fig.7 Matching window diagram

式中，SAD为差绝对值和运算，R为由排序计算得到的秩值。

将模板窗口与匹配窗口划分为若干块，通过对数极坐标网格(分区网格)将块内的像素划分为12个区间(径向区间数量为2，方向区间数量为4,8)，每个区间的秩值连接构成该块的描述子，窗口内所有块描述子连接构成RSS描述子向量。利用NCC计算描述子间的相似性值：

式中，V1,V2是模板窗口与匹配窗口上的RSS描述符，,为V1,V2的平均值。在搜索窗口内，逐像素计算匹配窗口与模板窗口间的相似性值，相似性最大的像素(xm,ym)即为匹配点。

(3) 亚像素坐标计算

如图8所示，为了得到亚像素级的匹配结果，利用相似性值最大的匹配点位(xm,ym)与其邻域8个点的相似性度量结果，拟合亚像素匹配模型参数值，模型函数偏导为0的位置为亚像素坐标。

图8 匹配点及邻域点示意图Fig.8 Schematic diagram of matching points and neighborhood points

亚像素匹配模型为

式中，k0,k1,k2,k3,k4,k5为模型参数，由最大的匹配点位(xm,ym)及邻域8个点的相似性值，通过最小二乘平差方法求解。亚像素坐标 (xa,ya)由式(10)计算得到：

(4) 坐标反算与高程控制点获取

在得到SAR地理编码影像中的匹配像点坐标(xa,ya)后，由RPC模型反解斜距SAR影像像点坐标(x,y)。像点(xa,ya)对 应平面坐标(X,Y)为

式中，(X0,Y0)为SAR地理编码影像原点地面坐标，由 (X,Y)插值SRTM DEM相同坐标下的高程值Z。将地面坐标与影像坐标正则化到-1与1之间，避免在计算过程中引入舍入误差，正则化公式为

(Xn,Yn,Zn)为 正则化的地面坐标，Xo,Xs,Yo,Ys,Zo,Zs为地面坐标的正则化参数，RPC模型多项式为

xo,xs,yo,ys为影像坐标正则化参数，计算出与激光高程点高程值H对应的斜距SAR影像像点坐标(x,y) 后，得到相应的SAR高程控制点(x,y,H)。

3 实验与分析

3.1 数据概况

本文选取中国登封市与日本横须贺市两个区域进行实验，区域1位于嵩山山脉中部，地形复杂，坡度多在6°以上，高程范围为229～1454 m。区域2位于东京湾西南岸的三浦半岛，地形起伏较小，高程范围为21～281 m，各区域地形如图9所示。

图9 各区域地理位置及影像范围图Fig.9 Geographical location and topographic map of each area

实验收集了获取时间为2018年11月至2021年10月覆盖区域1、区域2的ATL03级数据16,18轨，各区域分别选取3景GF-3聚束模式SAR影像，每景覆盖范围约为10 km×10 km，影像数据参数如表2所示。进行光子粗差剔除与地理编码时使用的DEM为30 m分辨率SRTM DEM (SRTM 1″)。

表2 SAR影像数据基本参数Tab.2 Basic parameters of SAR image data

3.2 SAR高程控制点提取与分析

(1) 激光高程点优选

由表3可知利用夜间高质量光子提取算法筛选光子数据，其中，区域1、区域2分别提取出夜间观测光子208935个，309832个，相比于原始数据剔除率为97.11%,94.52%，表明原始数据中存在较多噪点；通过置信度提取得到光子171242个，276938个，相较于夜间观测数据剔除率为18.04%,10.62%，表明夜间观测数据整体质量较高；通过SRTM DEM进行粗差剔除，提取到光子157635个，276890个，剔除率7.95%,0.02%，表明被标记为较高置信度的光子中仍有少量粗差点；最终通过大偏心率椭圆滤波核(阈值T=20)并利用矩形窗口(WS=30,δW=0.2)剔除误提取点后，得到平坦区域光子6609个，25641个，即激光高程点，剔除率95.81%,90.74%，表明平坦区域光子数量相对原始数据较少。

表3 不同筛选条件下光子数量Tab.3 The number of photons under different screening conditions

为验证所提激光高程点的准确性，通过由机载LiDAR获取的区域1分辨率为1.00 m的高精度DEM作为高程验证数据。将提取出的激光高程点按平面坐标确定1.00 m DEM中对应像素的高程值，直接计算两类数据高程残值，所有光子高程残差直方图如图10所示，其中，高程残差为2～4 m的光子数据主要位于水体。统计高程平均绝对误差(Mean Absolute Error,MAE)为1.43 m、均方根误差(Root Mean Square Error,RMSE)为2.65 m。利用多策略ATLAS数据优选算法得到了数量充足的光子数据作为激光高程点。

图10 高程残差直方图Fig.10 Elevation difference histogram

(2) 影像匹配

GE影像使用18级数据，SAR地理编码影像与GE影像地面分辨率S约为1.07 m。选取位于3景SAR影像重叠范围内分布良好的激光高程点投影至GE影像中与SAR地理编码影像匹配，即按每个激光高程点平面坐标为中心，选取GE影像与SAR地理编码影像的局部区域为模板窗口与匹配窗口，窗口大小为201像素×201像素，搜索窗口大小为281像素×281像素，分区网格的半径为3像素，相邻两个分区网格块之间在行与列方向上的重叠大小为3像素，则模板窗口与匹配窗口的RSS描述子维度为 (4+8)×48×48。在搜索窗口内逐像素计算模板窗口与匹配窗口描述子相似性值，得到相似性值最大的像素即为匹配像素，由该像素与邻域8个像素的行列号、描述子相似性值，通过亚像素匹配模型计算最终匹配结果。

通过人工目视判读的方法检查匹配结果，每个区域3景影像全部匹配正确作为SAR高程控制点，否则予以剔除，区域1、区域2分别得到48个与35个，各区域随机选取3个激光高程点对应SAR影像像点匹配结果如图11所示。

图11 匹配结果示意图Fig.11 Schematic diagram of matching results

(3) SAR高程控制点提取结果

激光高程点高程信息与匹配后的像点坐标构成SAR高程控制点，区域1、区域2 SAR高程控制点分布范围如图12所示，高程信息如表4、表5所示。

表4 区域1 SAR高程控制点高程信息Tab.4 Elevation information of SAR elevation control points in area 1

表5 区域2 SAR高程控制点高程信息Tab.5 Elevation information of SAR elevation control points in area 2

图12 激光高程点分布示意图Fig.12 Schematic diagram of laser elevation point distribution

(4) SAR高程控制点有效性验证

为验证提取出的SAR高程控制点的有效性，由在航空摄影获取的区域1分辨率为0.50 m DOM和1.00 m DEM上采集11个检查点作为定位结果验证数据。如图13所示，检查点主要位于道路交叉口、水库监测站等明显、易于判读的位置，且均位于3景SAR影像重叠区域内；其在3景SAR影像上的像点坐标均由人工量测获取，检查点信息如表6所示，经纬度坐标基准为WGS84，数值高位以“***”代替，高程基准为WGS84椭球高，分布情况如图14所示。

表6 检查点信息Tab.6 Information of checkpoints

图13 区域1部分检查点示意图Fig.13 Schematic diagram of some checkpoints in the area 1

图14 检查点分布情况示意图Fig.14 Schematic diagram of the distribution of checkpoints

分别利用GF-3影像原始RPC模型参数(方法1)、提取出的10个SAR高程控制点精化模型参数(方法2)、提取出的48个SAR高程控制点精化模型参数(方法3)进行立体定位，精度统计如表7所示，其中平面误差为高斯3°投影带下的平面坐标差，高程误差为WGS84椭球高坐标差，各检查点高程残差如图15所示。

由图15可知，利用方法1的高程误差具有明显系统性，方法2与方法3高程误差较小，且无明显系统性。由表7可知，采用提取出的SAR高程控制点进行几何模型参数精化后，立体定位精度明显提升，且选择少量分布合理的点也可实现较高精度的立体定位。通过不同数量SAR高程控制点的立体定位方法，验证了本文SAR高程控制点提取方法的有效性。

图15 检查点误差情况Fig.15 Checkpoint error condition

表7 定位结果Tab.7 Results of positioning

4 结语

本文设计了一种由ATLAS数据提取SAR高程控制点方法，通过多策略ATLAS数据优选，提取出了数量充足的高质量、平坦区域激光高程点，通过谷歌影像、SAR地理编码影像秩自相似描述子匹配，实现了激光点高程点高程值与斜距SAR影像像点坐标之间的对应，有效提取了SAR高程控制点。采用中国登封市和日本横须贺市两个区域的ATLAS数据进行高分三号SAR高程控制点提取，验证了方法的有效性。

目前仅利用了覆盖范围较小的聚束模型SAR影像，后续将进行大区域ATLAS数据SAR高程控制点提取与区域网平差实验，进一步验证本文方法的可靠性与普适性。

致谢实验中的高分三号SAR影像数据由中国资源卫星应用中心提供，谷歌地图影像由GGGIS提供，在此表示感谢。