SAR Elevation Control Point Extraction Combining Multistrategy ATLAS Data Preference and Image Matching
-
摘要: 为了精化星载SAR影像几何参数并提高立体定位精度,借鉴星载激光测高数据光学遥感影像高程控制点提取思路,设计了一种多策略高级地形激光测高系统(ATLAS)数据优选与影像匹配相结合的SAR高程控制点提取方法。该方法采用非夜间观测光子滤除、高置信度光子选取、SRTM DEM辅助的粗差剔除、大偏心率椭圆滤波核平坦区域光子筛选等多种策略,从ATLAS数据ATL03级产品中提取高质量、平坦区域的激光高程点,再依据SRTM DEM对斜距SAR影像进行地理编码,按激光高程点的平面坐标选取局部谷歌地球影像作为足印影像,采用秩自相似描述子进行足印影像与SAR地理编码影像的匹配,得到与激光高程点对应的SAR影像像点坐标,从而提取SAR高程控制点。采用中国登封市、日本横须贺市两个区域的ATLAS数据进行了高分三号SAR高程控制点提取实验,利用提取的高程控制点进行SAR影像几何参数精化,大幅提升了立体定位精度,验证了该文高程控制点提取方法的可行性和有效性。
-
关键词:
- 合成孔径雷达 /
- 高分三号 /
- 冰、云和陆地高程卫星-2 /
- 有理函数模型 /
- 定位
Abstract: To refine the geometric parameters of satellite-based SAR images and improve the stereo positioning accuracy, a method of SAR elevation control point extraction combining multistrategy ATLAS data preference and image matching has been developed. This method is based on the concept of optical remote sensing image elevation control point extraction from satellite-based laser altimetry data. The method employs various strategies, such as non-night observation photon filtering, high confidence photon selection, SRTM DEM-assisted coarse difference rejection, and large eccentricity elliptical filtering kernel flat area photon screening. To extract laser elevation points with high quality and flat area from ATLAS data for ATL03 level products. Then the geocoding of the slant range SAR images is performed using the SRTM DEM. The local Google Earth images are selected as the footprint images according to the plane coordinates of the laser elevation points. The rank self-similarity descriptor is used to match the footprint images with the SAR geocoded images. The coordinates of the SAR images corresponding to the laser elevation points are obtained. Thus, SAR elevation control points are extracted. The extraction of GF-3 SAR elevation control points was performed using ATLAS data from two regions: Dengfeng, China, and Yokosuka, Japan. The geometric parameter refinement of SAR images using extracted elevation control points significantly improved the accuracy of stereo positioning and verified that the method for extracting elevation control points described in this paper is feasible and effective. -
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:10 000立体测图需求提供了技术支持,综上,星载激光测高数据与光学遥感影像联合应用具有重要价值,但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影像高程控制点。
2.1 多策略ATLAS数据优选
多策略ATLAS数据优选流程如图2所示,分为夜间高质量光子提取与大偏心率椭圆平坦区域光子提取两个步骤。夜间高质量光子提取中,根据光子三维坐标计算太阳高度角,从而滤除非夜间观测光子数据;由置信度信息提取位于地面的信号光子;通过比较光子高程值与SRTM DEM高程值剔除粗差。最后,利用大偏心率椭圆滤波核提取位于平坦区域的光子作为激光高程点。算法需要读取的参数如表1所示。
表 1 算法使用到的ATL03参数Table 1. ATL03 parameters used by the algorithm参数 含义 本文符号 lat_ph 每个接收光子纬度(WGS84) B lon_ph 每个接收光子经度(WGS84) L h_ph 每个接收光子高程(WGS84椭球高) H signal_conf_ph 每个接收光子置信度 C data_start_utc 数据开始时间 t0 dist_ph_along 每个接收光子在该组中的沿轨距离 li segment_length 每组沿轨距离(19.8~20.2 m) Li segment_ph_cnt 每组内光子数 wi (1) 夜间高质量光子提取
为减小太阳光噪声对高程控制点提取的影响,依据光子三维坐标计算的太阳高度角
Hθ 提取夜间观测数据。Hθ 为太阳光线与地平面之间的夹角,该值小于零时可认为数据获取时刻在当地为夜间,计算式为Hθ=arcsin(sinB⋅sinδ+cosB⋅cosδ⋅cosΩ) (1) 式中,
δ 为太阳赤纬,Ω 为时角,δ=0.006918−0.399912cos(2π(N−1)365)+0.070257sin(2π(N−1)365)−0.006758cos(4π(N−1)365)+0.000907sin(4π(N−1)365)−0.002697cos(6π(N−1)365)+0.00148sin(6π(N−1)365) (2) N为观测时间与该年第1天之间的天数差,
Ω=15⋅(t′0−12) (3) t′0 为由t0 计算的真太阳时[27]。ATL03产品为每个光子提供不同地形下基于泊松算法的信号置信度标识,分别为0~4,在被标记为陆地的数据中提取C大于0的光子。
根据拉依达准则,比较光子高程值与相同平面坐标下SRTM DEM高程值,保留高程差值小于3倍SRTM DEM标准差的光子,即
ΔH=|Hice−(HSRTM+ξ)|<3σD (4) 式中,
Hice 为光子高程值,HSRTM 为相同平面坐标下SRTM DEM高程值,ξ 为egm96模型水准面差,σD 为SRTM DEM标准差。(2) 大偏心率椭圆平坦区域光子提取
地面有较大的坡度时,物像不一致问题会造成高程控制偏差问题。如图3所示,S为成像时刻SAR卫星位置,地物A对应的同名像点为a,A匹配到的像点为
a′ ,a′ 对应同名地物点为A′ ,a与a′ 在距离向、方位向上相差Δy 、Δx 个像素,当地面存在大小为α 的坡度时,会分别引起高程为ΔZ1 与ΔZ2 的差值:{ΔZ1≈Δy⋅Ry⋅sinαΔZ2=Δx⋅Rx⋅tanα (5) 式中,
Ry ,Rx 为SAR影像距离向、方位向分辨率。则当Δx=1 ,Δy=1 ,α=1.5∘ ,Ry=0.6m ,Rx=0.4m 时,由式(5)计算可得,将引起高程为ΔZ1=0.0157m ,ΔZ2=0.0105m 的误差。因此,筛选出位于平坦区域的光子即可保证激光点高程精度,也可减弱物像不一致问题造成的高程控制偏差问题。星载激光测高数据常用空间聚类算法进行地面光子的提取,椭圆形滤波核更符合地面信号光子的分布特征。如图4所示,对于某个光子p,定义以点p为中心的长轴沿轨椭圆形窗口,点q是否在该椭圆内可用式(6)判断:
d(p,q)=(lp−lq)2a2+(hp−hq)2b2 (6) 式中,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°,能够满足高程控制的要求。为了提取平坦区域裸地的光子,以当前光子为中心设置宽度为
WS 、长度应能覆盖所有光子的矩形窗口,计算窗口内提取的夜间高质量光子高程方差,剔除方差大于阈值δW 的光子数据,得到优选的激光高程数据。2.2 谷歌影像、SAR影像秩自相似描述子匹配
为使激光高程点与SAR影像匹配更具鲁棒性与准确性,本文通过影像匹配的方法实现二者的匹配。如图6所示,将提取出的激光高程点按平面坐标投影至GE中,并与SAR地理编码影像进行匹配,建立SAR影像像点与激光点高程值之间的对应关系,作为高程控制条件引入几何模型参数精化中,从而提升SAR影像定位精度。
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)定义为
RSq(x,y)=R(SADq(x,y)) (7) 式中,
SAD 为差绝对值和运算,R为由排序计算得到的秩值。将模板窗口与匹配窗口划分为若干块,通过对数极坐标网格(分区网格)将块内的像素划分为12个区间(径向区间数量为2,方向区间数量为4, 8),每个区间的秩值连接构成该块的描述子,窗口内所有块描述子连接构成RSS描述子向量。利用NCC计算描述子间的相似性值:
RSSNCC=l∑k=1(V1(k)−ˆV1)(V2(k)−ˆV2)√l∑k=1(V1(k)−ˆV1)2l∑k=1(V2(k)−ˆV2)2 (8) 式中,
V1 ,V2 是模板窗口与匹配窗口上的RSS描述符,ˆV1 ,ˆV2 为V1 ,V2 的平均值。在搜索窗口内,逐像素计算匹配窗口与模板窗口间的相似性值,相似性最大的像素(xm,ym) 即为匹配点。(3) 亚像素坐标计算
如图8所示,为了得到亚像素级的匹配结果,利用相似性值最大的匹配点位
(xm,ym) 与其邻域8个点的相似性度量结果,拟合亚像素匹配模型参数值,模型函数偏导为0的位置为亚像素坐标。亚像素匹配模型为
RSSNCC(x,y)=k0+k1x+k2y+k3xy+k4x2+k5y2 (9) 式中,
k0 ,k1 ,k2 ,k3 ,k4 ,k5 为模型参数,由最大的匹配点位(xm,ym) 及邻域8个点的相似性值,通过最小二乘平差方法求解。亚像素坐标(xa,ya) 由式(10)计算得到:{k1+2k3xa+k5ya=0k2+k5xa+2k4ya=0 (10) (4) 坐标反算与高程控制点获取
在得到SAR地理编码影像中的匹配像点坐标
(xa,ya) 后,由RPC模型反解斜距SAR影像像点坐标(x,y) 。像点(xa,ya) 对应平面坐标(X,Y) 为{X=X0+xa⋅SY=Y0+ya⋅S (11) 式中,
(X0,Y0) 为SAR地理编码影像原点地面坐标,由(X,Y) 插值SRTM DEM相同坐标下的高程值Z。将地面坐标与影像坐标正则化到–1与1之间,避免在计算过程中引入舍入误差,正则化公式为{Xn=X−XoXsYn=Y−YoYsZn=Z−ZoZs (12) (Xn,Yn,Zn) 为正则化的地面坐标,Xo ,Xs ,Yo ,Ys ,Zo ,Zs 为地面坐标的正则化参数, RPC模型多项式为{xn=3∑u1=13∑u2=13∑u3=1pu1u2u31Xu1nYu2nZu3n3∑u1=13∑u2=13∑u3=1pu1u2u32Xu1nYu2nZu3nyn=3∑u1=13∑u2=13∑u3=1pu1u2u33Xu1nYu2nZu3n3∑u1=13∑u2=13∑u3=1pu1u2u34Xu1nYu2nZu3n (13) (xn,yn) 为正则化的像点坐标,pu1u2u3v(v=1,2,3,4) 为RPC模型多项式系数,像点坐标(x,y) 为{x=xn⋅xs+xoy=yn⋅ys+yo (14) 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所示。
实验收集了获取时间为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影像数据基本参数Table 2. Basic parameters of SAR image data区域 序号 成像日期 轨道类型 入射角 尺寸(像素) 中心经纬度 区域1 影像1 20191021 升轨 23.31°~24.58° 32966×12576 (113.114°E, 34.514°N) 影像2 20191206 降轨 42.52°~43.18° 30944×15648 (113.120°E, 34.512°N) 影像3 20190829 降轨 46.24°~46.82° 28224×16056 (113.130°E, 34.511°N) 区域2 影像4 20190116 升轨 38.65°~39.38° 30558×14624 (139.645°E, 35.290°N) 影像5 20190116 降轨 25.16°~26.32° 31920×12576 (139.645°E, 35.296°N) 影像6 20190121 降轨 35.46°~36.25° 30638×13600 (139.642°E, 35.296°N) 3.2 SAR高程控制点提取与分析
(1) 激光高程点优选
由表3可知利用夜间高质量光子提取算法筛选光子数据,其中,区域1、区域2分别提取出夜间观测光子208 935个,309 832个,相比于原始数据剔除率为97.11%, 94.52%,表明原始数据中存在较多噪点;通过置信度提取得到光子171 242个,276 938个,相较于夜间观测数据剔除率为18.04%, 10.62%,表明夜间观测数据整体质量较高;通过SRTM DEM进行粗差剔除,提取到光子157 635个,276 890个,剔除率7.95%, 0.02%,表明被标记为较高置信度的光子中仍有少量粗差点;最终通过大偏心率椭圆滤波核(阈值
T=20 )并利用矩形窗口(WS=30 ,δW=0.2 )剔除误提取点后,得到平坦区域光子6 609个,25 641个,即激光高程点,剔除率95.81%, 90.74%,表明平坦区域光子数量相对原始数据较少。表 3 不同筛选条件下光子数量Table 3. The number of photons under different screening conditions数据类型 区域1 区域2 光子数量(个) 数据剔除率(%) 光子数量(个) 数据剔除率(%) 原始数据 7 221 634 0 5 651 808 0 夜间观测数据 208 935 97.11 309 832 94.52 高置信度数据 171 242 18.04 276 938 10.62 粗差剔除数据 157 635 7.95 276 890 0.02 平坦区域数据 6 609 95.81 25 641 90.74 为验证所提激光高程点的准确性,通过由机载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数据优选算法得到了数量充足的光子数据作为激光高程点。
(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所示。
(3) SAR高程控制点提取结果
激光高程点高程信息与匹配后的像点坐标构成SAR高程控制点,区域1、区域2 SAR高程控制点分布范围如图12所示,高程信息如表4、表5所示。
表 4 区域1 SAR高程控制点高程信息Table 4. Elevation information of SAR elevation control points in area 1序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 1 457.09 11 402.28 21 426.60 31 359.66 41 368.92 2 430.53 12 397.58 22 373.99 32 355.80 42 375.77 3 435.54 13 404.51 23 351.45 33 352.67 43 377.07 4 435.32 14 403.89 24 348.68 34 343.63 44 373.36 5 435.43 15 385.88 25 333.22 35 336.50 45 372.15 6 435.10 16 384.23 26 412.32 36 336.41 46 371.26 7 439.42 17 382.38 27 390.60 37 322.30 47 370.16 8 408.47 18 407.05 28 389.62 38 334.37 48 368.57 9 406.28 19 403.99 29 363.33 39 334.76 10 404.04 20 407.59 30 360.98 40 351.15 表 5 区域2 SAR高程控制点高程信息Table 5. Elevation information of SAR elevation control points in area 2序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 1 95.58 8 95.69 15 39.42 22 40.07 29 39.47 2 38.71 9 204.95 16 39.32 23 40.05 30 69.25 3 107.14 10 205.64 17 40.42 24 39.73 31 104.99 4 107.19 11 71.39 18 39.08 25 39.95 32 71.12 5 107.19 12 71.53 19 40.33 26 38.90 33 71.23 6 57.84 13 222.95 20 40.44 27 78.16 34 60.62 7 83.25 14 38.98 21 39.63 28 39.44 35 60.64 (4) SAR高程控制点有效性验证
为验证提取出的SAR高程控制点的有效性,由在航空摄影获取的区域1分辨率为0.50 m DOM和1.00 m DEM上采集11个检查点作为定位结果验证数据。如图13所示,检查点主要位于道路交叉口、水库监测站等明显、易于判读的位置,且均位于3景SAR影像重叠区域内;其在3景SAR影像上的像点坐标均由人工量测获取,检查点信息如表6所示,经纬度坐标基准为WGS84,数值高位以“***”代替,高程基准为WGS84椭球高,分布情况如图14所示。
表 6 检查点信息Table 6. Information of checkpoints序号 经度(°) 纬度(°) 高程(m) 1 *.**3996 *.**7348 459.38 2 *.**0267 *.**2819 487.50 3 *.**6372 *.**4978 438.96 4 *.**2031 *.**1293 399.94 5 *.**9026 *.**1702 384.47 6 *.**7216 *.**9004 405.98 7 *.**4967 *.**1815 430.68 8 *.**3398 *.**2909 374.35 9 *.**8198 *.**0609 323.37 10 *.**9957 *.**5764 383.76 11 *.**2113 *.**4711 343.14 分别利用GF-3影像原始RPC模型参数(方法1)、提取出的10个SAR高程控制点精化模型参数(方法2)、提取出的48个SAR高程控制点精化模型参数(方法3)进行立体定位,精度统计如表7所示,其中平面误差为高斯3°投影带下的平面坐标差,高程误差为WGS84椭球高坐标差,各检查点高程残差如图15所示。
表 7 定位结果Table 7. Results of positioning定位
方法SAR高程控制
点数(个)检查
点数(个)平面中
误差(m)高程中
误差(m)方法1 0 11 10.42 30.42 方法2 10 11 7.35 3.77 方法3 48 11 7.34 3.69 由图15可知,利用方法1的高程误差具有明显系统性,方法2与方法3高程误差较小,且无明显系统性。由表7可知,采用提取出的SAR高程控制点进行几何模型参数精化后,立体定位精度明显提升,且选择少量分布合理的点也可实现较高精度的立体定位。通过不同数量SAR高程控制点的立体定位方法,验证了本文SAR高程控制点提取方法的有效性。
4. 结语
本文设计了一种由ATLAS数据提取SAR高程控制点方法,通过多策略ATLAS数据优选,提取出了数量充足的高质量、平坦区域激光高程点,通过谷歌影像、SAR地理编码影像秩自相似描述子匹配,实现了激光点高程点高程值与斜距SAR影像像点坐标之间的对应,有效提取了SAR高程控制点。采用中国登封市和日本横须贺市两个区域的ATLAS数据进行高分三号SAR高程控制点提取,验证了方法的有效性。
目前仅利用了覆盖范围较小的聚束模型SAR影像,后续将进行大区域ATLAS数据SAR高程控制点提取与区域网平差实验,进一步验证本文方法的可靠性与普适性。
致谢 实验中的高分三号SAR影像数据由中国资源卫星应用中心提供,谷歌地图影像由GGGIS提供,在此表示感谢。
-
表 1 算法使用到的ATL03参数
Table 1. ATL03 parameters used by the algorithm
参数 含义 本文符号 lat_ph 每个接收光子纬度(WGS84) B lon_ph 每个接收光子经度(WGS84) L h_ph 每个接收光子高程(WGS84椭球高) H signal_conf_ph 每个接收光子置信度 C data_start_utc 数据开始时间 t0 dist_ph_along 每个接收光子在该组中的沿轨距离 li segment_length 每组沿轨距离(19.8~20.2 m) Li segment_ph_cnt 每组内光子数 wi 表 2 SAR影像数据基本参数
Table 2. Basic parameters of SAR image data
区域 序号 成像日期 轨道类型 入射角 尺寸(像素) 中心经纬度 区域1 影像1 20191021 升轨 23.31°~24.58° 32966×12576 (113.114°E, 34.514°N) 影像2 20191206 降轨 42.52°~43.18° 30944×15648 (113.120°E, 34.512°N) 影像3 20190829 降轨 46.24°~46.82° 28224×16056 (113.130°E, 34.511°N) 区域2 影像4 20190116 升轨 38.65°~39.38° 30558×14624 (139.645°E, 35.290°N) 影像5 20190116 降轨 25.16°~26.32° 31920×12576 (139.645°E, 35.296°N) 影像6 20190121 降轨 35.46°~36.25° 30638×13600 (139.642°E, 35.296°N) 表 3 不同筛选条件下光子数量
Table 3. The number of photons under different screening conditions
数据类型 区域1 区域2 光子数量(个) 数据剔除率(%) 光子数量(个) 数据剔除率(%) 原始数据 7 221 634 0 5 651 808 0 夜间观测数据 208 935 97.11 309 832 94.52 高置信度数据 171 242 18.04 276 938 10.62 粗差剔除数据 157 635 7.95 276 890 0.02 平坦区域数据 6 609 95.81 25 641 90.74 表 4 区域1 SAR高程控制点高程信息
Table 4. Elevation information of SAR elevation control points in area 1
序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 1 457.09 11 402.28 21 426.60 31 359.66 41 368.92 2 430.53 12 397.58 22 373.99 32 355.80 42 375.77 3 435.54 13 404.51 23 351.45 33 352.67 43 377.07 4 435.32 14 403.89 24 348.68 34 343.63 44 373.36 5 435.43 15 385.88 25 333.22 35 336.50 45 372.15 6 435.10 16 384.23 26 412.32 36 336.41 46 371.26 7 439.42 17 382.38 27 390.60 37 322.30 47 370.16 8 408.47 18 407.05 28 389.62 38 334.37 48 368.57 9 406.28 19 403.99 29 363.33 39 334.76 10 404.04 20 407.59 30 360.98 40 351.15 表 5 区域2 SAR高程控制点高程信息
Table 5. Elevation information of SAR elevation control points in area 2
序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 序号 高程(m) 1 95.58 8 95.69 15 39.42 22 40.07 29 39.47 2 38.71 9 204.95 16 39.32 23 40.05 30 69.25 3 107.14 10 205.64 17 40.42 24 39.73 31 104.99 4 107.19 11 71.39 18 39.08 25 39.95 32 71.12 5 107.19 12 71.53 19 40.33 26 38.90 33 71.23 6 57.84 13 222.95 20 40.44 27 78.16 34 60.62 7 83.25 14 38.98 21 39.63 28 39.44 35 60.64 表 6 检查点信息
Table 6. Information of checkpoints
序号 经度(°) 纬度(°) 高程(m) 1 *.**3996 *.**7348 459.38 2 *.**0267 *.**2819 487.50 3 *.**6372 *.**4978 438.96 4 *.**2031 *.**1293 399.94 5 *.**9026 *.**1702 384.47 6 *.**7216 *.**9004 405.98 7 *.**4967 *.**1815 430.68 8 *.**3398 *.**2909 374.35 9 *.**8198 *.**0609 323.37 10 *.**9957 *.**5764 383.76 11 *.**2113 *.**4711 343.14 表 7 定位结果
Table 7. Results of positioning
定位
方法SAR高程控制
点数(个)检查
点数(个)平面中
误差(m)高程中
误差(m)方法1 0 11 10.42 30.42 方法2 10 11 7.35 3.77 方法3 48 11 7.34 3.69 -
[1] WANG Taoyang, LI Xin, ZHANG Guo, et al. Large-scale orthorectification of GF-3 SAR images without ground control points for China’s land area[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1–17. doi: 10.1109/TGRS.2022.3142372 [2] 润一. 高分三号卫星在轨几何定标及与高分二号光学卫星影像联合定位[D]. [硕士论文], 武汉大学, 2017.RUN Yi. On-orbit geometric calibration of GF-3 satellite and joint-positioning of GF-3 and GF-2 satellite images[D]. [Master dissertation], Wuhan University, 2017. [3] 张红敏, 靳国旺, 徐青, 等. 利用单个地面控制点的SAR图像高精度立体定位[J]. 雷达学报, 2014, 3(1): 85–91. doi: 10.3724/SP.J.1300.2014.13138ZHANG Hongmin, JIN Guowang, XU Qing, et al. Accurate positioning with stereo SAR images and one ground control point[J]. Journal of Radars, 2014, 3(1): 85–91. doi: 10.3724/SP.J.1300.2014.13138 [4] 魏钜杰, 张继贤, 赵争, 等. 稀少控制下TerraSAR-X影像高精度直接定位方法[J]. 测绘科学, 2011, 36(1): 58–60, 50. doi: 10.16251/j.cnki.1009-2307.2011.01.006WEI Jujie, ZHANG Jixian, ZHAO Zheng, et al. High-precisely direct geo-location method for TerraSAR-X image with sparse GCPs[J]. Science of Surveying and mapping, 2011, 36(1): 58–60, 50. doi: 10.16251/j.cnki.1009-2307.2011.01.006 [5] 张祖勋, 陶鹏杰. 谈大数据时代的“云控制”摄影测量[J]. 测绘学报, 2017, 46(10): 1238–1248. doi: 10.11947/j.AGCS.2017.20170337ZHANG Zuxun and TAO Pengjie. An overview on “cloud control” photogrammetry in big data era[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1238–1248. doi: 10.11947/j.AGCS.2017.20170337 [6] 方勇, 龚辉, 张丽, 等. 从全球激光点云到三维数字地球空间框架: 全球精确测绘进阶之路[J]. 激光与光电子学进展, 2022, 59(12): 1200002. doi: 10.3788/LOP202259.1200002FANG Yong, GONG Hui, ZHANG Li, et al. From global laser point cloud acquisition to 3D digital geospatial framework: The advanced road of global accurate mapping[J]. Laser &Optoelectronics Progress, 2022, 59(12): 1200002. doi: 10.3788/LOP202259.1200002 [7] NEUMANN T A, MARTINO A J, MARKUS T, et al. The ice, cloud, and land elevation satellite-2 mission: A global geolocated photon product derived from the advanced topographic laser altimeter system[J]. Remote Sensing of Environment, 2019, 233: 111325. doi: 10.1016/j.rse.2019.111325 [8] XIE Huan, XU Qi, YE Dan, et al. A comparison and review of surface detection methods using MBL, MABEL, and ICESat-2 photon-counting laser altimetry data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14: 7604–7623. doi: 10.1109/JSTARS.2021.3094195 [9] ICE, CLOUD, and land elevation satellite-2 (ICESat-2) algorithm theoretical basis document (ATBD) for global geolocated photons (ATL03)[EB/OL]. https://nsidc.org/sites/default/files/icesat2_atl03_atbd_r005_0.pdf, 2021. [10] LIN Xiaojuan, XU Min, CAO Chunxiang, et al. Estimates of forest canopy height using a combination of ICESat-2/ATLAS data and stereo-photogrammetry[J]. Remote Sensing, 2020, 12(21): 3649. doi: 10.3390/rs12213649 [11] 张帅台, 李国元, 周晓青, 等. 基于多特征自适应的单光子点云去噪算法[J]. 红外与激光工程, 2022, 51(6): 20210949. doi: 10.3788/IRLA20210949ZHANG Shuaitai, LI Guoyuan, ZHOU Xiaoqing, et al. Single photon point cloud denoising algorithm based on multi-features adaptive[J]. Infrared and Laser Engineering, 2022, 51(6): 20210949. doi: 10.3788/IRLA20210949 [12] ZHU Xiaoxiao, NIE Sheng, WANG Cheng, et al. A ground elevation and vegetation height retrieval algorithm using micro-pulse photon-counting Lidar data[J]. Remote Sensing, 2018, 10(12): 1962. doi: 10.3390/rs10121962 [13] MARKUS T, NEUMANN T, MARTINO A, et al. The ice, cloud, and land elevation satellite-2 (ICESat-2): Science requirements, concept, and implementation[J]. Remote Sensing of Environment, 2017, 190: 260–273. doi: 10.1016/j.rse.2016.12.029 [14] ROSIEK M R, KIRK R L, ARCHINAL B A, et al. Utility of Viking orbiter images and products for Mars mapping[J]. Photogrammetric Engineering & Remote Sensing, 2005, 71(10): 1187–1195. doi: 10.14358/PERS.71.10.1187 [15] 何钰, 吴绍民, 邢帅. 基于RFM的嫦娥一号CCD影像区域网平差研究[J]. 测绘科学, 2013, 38(6): 5–6, 15. doi: 10.16251/j.cnki.1009-2307.2013.06.034HE Yu, WU Shaomin, and XING Shuai. Block adjustment of Chang’e-1 CCD images based on RFM[J]. Science of Surveying and Mapping, 2013, 38(6): 5–6, 15. doi: 10.16251/j.cnki.1009-2307.2013.06.034 [16] 耿迅. 火星形貌摄影测量技术研究[D]. [博士论文], 解放军信息工程大学, 2014.GENG Xun. Research on photogrammetric processing for Mars topographic mapping[D]. [Ph. D. dissertation], Information Engineering University, 2014. [17] 王晋, 张勇, 张祖勋, 等. ICESat激光高程点辅助的天绘一号卫星影像立体区域网平差[J]. 测绘学报, 2018, 47(3): 359–369. doi: 10.11947/j.AGCS.2018.20170425WANG Jin, ZHANG Yong, ZHNAG Zuxun, et al. ICESat laser points assisted block adjustment for mapping Satellite-1 stereo imagery[J]. Acta Geodaetica et Cartographica Sinica, 2018, 47(3): 359–369. doi: 10.11947/j.AGCS.2018.20170425 [18] 张鑫磊, 邢帅, 徐青, 等. ATLAS数据与资源三号02星影像联合区域网平差[J]. 红外与激光工程, 2020, 49(S2): 20200194. doi: 10.3788/IRLA20200194ZHANG Xinlei, XING Shuai, XU Qing, et al. Joint block adjustment for ATLAS data and ZY3-02 stereo imagery[J]. Infrared and Laser Engineering, 2020, 49(S2): 20200194. doi: 10.3788/IRLA20200194 [19] ZHANG Guo, JIANG Boyang, WANG Taoyang, et al. Combined block adjustment for optical satellite stereo imagery assisted by spaceborne SAR and laser altimetry data[J]. Remote Sensing, 2021, 13(16): 3062. doi: 10.3390/rs13163062 [20] 谭建伟, 程春泉. 建筑影像高程控制点的激光测高全波形分解提取[J]. 测绘科学, 2021, 46(8): 1–7, 13. doi: 10.16251/j.cnki.1009-2307.2021.08.001TAN Jianwei and CHENG Chunquan. Extracting building image elevation control points by decomposing full waveform of laser altimetry[J]. Science of Surveying and Mapping, 2021, 46(8): 1–7, 13. doi: 10.16251/j.cnki.1009-2307.2021.08.001 [21] 王密, 韦钰, 杨博, 等. ICESat-2/ATLAS全球高程控制点提取与分析[J]. 武汉大学学报:信息科学版, 2021, 46(2): 184–192. doi: 10.13203/j.whugis20200531WANG Mi, WEI Yu, YANG Bo, et al. Extraction and analysis of global elevation control points from ICESat-2/ATLAS data[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 184–192. doi: 10.13203/j.whugis20200531 [22] 郑迎辉, 张艳, 王涛, 等. 基于ICESat-2数据的高程控制点提取和精度验证[J]. 地球信息科学学报, 2022, 24(7): 1234–1244. doi: 10.12082/dqxxkx.2022.210667ZHENG Yinghui, ZHANG Yan, WANG Tao, et al. Elevation control points extraction and accuracy validation based on ICESat-2 data[J]. Journal of Geo-Information Science, 2022, 24(7): 1234–1244. doi: 10.12082/dqxxkx.2022.210667 [23] 中华人民共和国国家质量监督检验检疫总局, 中国国家标准化管理委员会. GB/T 12341-2008 1: 25000 1: 50000 1: 100000 地形图航空摄影测量外业规范[S]. 北京: 中国标准出版社, 2008.General Administration of Quality Supervision, Inspection and Quarantine of the People's Republic of China and Standardization Administration of China. GB/T 12341-2008 Specifications for aerophotogrammetric field work of 1∶25000 1∶50000 1∶100000 topographic maps[S]. Beijing: Standards Press of China, 2008. [24] 曹宁, 周平, 王霞, 等. 激光测高数据辅助卫星成像几何模型精化处理[J]. 遥感学报, 2018, 22(4): 599–610. doi: 10.11834/jrs.20187252CAO Ning, ZHOU Ping, WANG Xia, et al. Refined processing of laser altimeter data-aided satellite geometry model[J]. Journal of Remote Sensing, 2018, 22(4): 599–610. doi: 10.11834/jrs.20187252 [25] 唐新明, 刘昌儒, 张恒, 等. 高分七号卫星立体影像与激光测高数据联合区域网平差[J]. 武汉大学学报:信息科学版, 2021, 46(10): 1423–1430. doi: 10.13203/j.whugis20210417TANG Xinming, LIU Changru, ZHANG Heng, et al. GF-7 satellite stereo images block adjustment assisted with laser altimetry data[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1423–1430. doi: 10.13203/j.whugis20210417 [26] XIONG Xin, XU Qing, JIN Guowang, et al. Rank-based local self-similarity descriptor for optical-to-SAR image matching[J]. IEEE Geoscience and Remote Sensing Letters, 2020, 17(10): 1742–1746. doi: 10.1109/LGRS.2019.2955153 [27] 王国安, 米鸿涛, 邓天宏, 等. 太阳高度角和日出日落时刻太阳方位角一年变化范围的计算[J]. 气象与环境科学, 2007, 30(S1): 161–164. doi: 10.16765/j.cnki.1673-7148.2007.s1.031WANG Guoan, MI Hongtao, DENG Tianhong, et al. Calculation of the change range of the sun high angle and the azimuth of sunrise and sunset in one year[J]. Meteorological and Environmental Sciences, 2007, 30(S1): 161–164. doi: 10.16765/j.cnki.1673-7148.2007.s1.031 期刊类型引用(1)
1. 徐源,罗运华,赵耀. 基于参数重建的机载SAR图像目标高精度定位. 电子设计工程. 2024(13): 176-180 . 百度学术
其他类型引用(1)
-