多策略ATLAS數據優選與影像匹配相結合的SAR高程控制點提取

余 洋 靳國旺 熊 新 童成功 張國平

(信息工程大學地理空間信息學院 鄭州 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提供，在此表示感謝。