快速LZD匹配異源DEM空洞修復(fù)應(yīng)用研究

李小虎，蘆 穎，田 壯

中國電子科技集團公司第十五研究所，北京 100080

地形數(shù)據(jù)信息處理技術(shù)是數(shù)字地形分析（digital terrain analyse，DTA）的基礎(chǔ)，一直以來作為測繪學(xué)、地貌學(xué)和計算機圖形學(xué)的研究熱點，被廣泛應(yīng)用于道路規(guī)劃、精準農(nóng)業(yè)、數(shù)字戰(zhàn)場、實景三維等眾多領(lǐng)域[1-2]。隨著全球信息化、數(shù)字化的快速發(fā)展，作為一種重要的地理信息資源，數(shù)字高程模型（digital elevation model，DEM）產(chǎn)品的數(shù)據(jù)采集方式更加多元，可根據(jù)用戶的個性化需求定制不同精度、形態(tài)復(fù)雜度、參考框架和覆蓋范圍的數(shù)據(jù)產(chǎn)品[3-4]。其中，機載激光雷達（airborne light detection and ranging，ALiDAR）憑借高密度、高精度的三維數(shù)據(jù)采集優(yōu)勢已成為獲取高質(zhì)量DEM的主要數(shù)據(jù)來源[5-7]。然而，受云霧遮擋、高程地形（陡峭的斜坡和深谷）和信號干擾等因素對立體成像的影響[8]，導(dǎo)致構(gòu)建的ALiDAR航測點云在某些區(qū)域散布較多的數(shù)據(jù)空洞-像素高程值為零，從而嚴重影響了地表數(shù)字化模擬產(chǎn)品的可用性[9-10]。

如何有效填補上述DEM 數(shù)據(jù)空缺，是學(xué)者們進行技術(shù)性研究任務(wù)的關(guān)鍵難點。目前主流方法有兩種[11-13]：一是利用空缺區(qū)域鄰域信息的直接插值法；二是結(jié)合三維數(shù)據(jù)融合與空間插值的綜合法。改進的插值方法有：基于高斯權(quán)值-向量基（Gauss weight-vector base，GW-VB）的插值方法，在明確稀疏與稠密點云的劃定界限后，能夠有效提升離散點間的平滑處理效果[14]；等高線間隔重差值法（contour interval redifference，CIR），該方法能夠較好地捕捉到總體的地形趨勢，但是由于缺乏詳細的描述特征，常應(yīng)用于水文區(qū)域的空缺填充[15]；高精度曲面模型（high accuracy surface modelling，HASM）插值法：只需借助空缺周圍數(shù)據(jù)進行填補，能夠快速填補位于簡單地形表面的小尺度的數(shù)據(jù)空洞[16]。改進后的自適應(yīng)Kriging 插值法：能夠克服離散點分布不均勻條件下所引入的建模誤差，通過最優(yōu)無偏估計提高插值精度[17]。

插值法在空洞范圍較小且平坦的數(shù)據(jù)修復(fù)中表現(xiàn)較為理想，對于山脊、山谷等大面積區(qū)域缺失的情況下，常采用三維數(shù)據(jù)融合與空間插值綜合法，以滿足高程產(chǎn)品的質(zhì)量要求。其中，岳國棟等[18]提出了干系數(shù)與垂直基線高度加權(quán)融合的方法，利用一組主輔影像去除側(cè)視異常區(qū)域；李世明[19]提出了多源點云漸進配準和拉普拉斯融合算法，在面向城市三維空間精細建模中完善了局部復(fù)雜的結(jié)構(gòu)表面；基于泊松融合[20]的無縫修復(fù)法，通過統(tǒng)計分析極大地增加了填充可用候選集，且具備復(fù)制紋理的效果；陳雪榮[21]探索出了優(yōu)化后的徑向基神經(jīng)網(wǎng)絡(luò)（radial basis function，RBF）方法，可有效去除大尺度噪聲，并用網(wǎng)絡(luò)預(yù)測值代替實測值有效解決非線性空洞修補；余婷婷等[22]提出的將隨機森林回歸算法（random forest regression，RFR）應(yīng)用在DEM 空洞修補精度校正上，適用于處理各種地形異質(zhì)程度的高程數(shù)據(jù)，泛化能力強。上述方法雖在去除大尺度數(shù)據(jù)空洞方面做了大量工作，但當數(shù)據(jù)分辨率較高且空洞覆蓋面積較大時，海量的訓(xùn)練樣本導(dǎo)致模型不易收斂，從而影響了修復(fù)效率。

本文提出了一種基于異源輔助數(shù)據(jù)的快速最小高差三維匹配算法，配合統(tǒng)一計算架構(gòu)（compute unified device architecture，CUDA）實現(xiàn)GPU 多核策略并行計算，提高LZD 空間匹配模型求解和“臺階”融合算法對空洞邊緣離散點平滑處理的效率和效果。修復(fù)后得到的DEM數(shù)據(jù)完整度高、過程可靠，利用該方法形成的軟件成果已在工程應(yīng)用中取得了較好的實效，具有一定的推廣價值。

1 方法原理

1.1 LZD表面匹配原理

為了進行快速、高精度的空洞修復(fù)，首先需要統(tǒng)一輔助數(shù)據(jù)與待修補數(shù)據(jù)間的坐標基準，將兩個結(jié)構(gòu)、位姿差異的三維空間進行表面匹配。目前廣泛使用的方式有最近點迭代算法（iterative closest point，ICP）和LZD 算法等[23]。鑒于待匹配兩數(shù)據(jù)規(guī)則格網(wǎng)的地形表述特點，選用LZD 算法進行待轉(zhuǎn)換空間映射模型的計算，如圖1所示。

圖1 空間坐標關(guān)系示意圖Fig.1 Spatial coordinate relationship diagram

依據(jù)兩組DEM 數(shù)據(jù)組織形式，設(shè)待匹配模型為Zc=f(Xc,Yc)，其中Xc、Yc是待修復(fù)表面投影至二維平面的坐標值；同理，標準模型為Zs=f(Xs,Ys)。待修復(fù)表面上一點Pc=[Xc,Yc,Zc]T與基準數(shù)據(jù)上對應(yīng)點Ps=[Xs,Ys,Zs]T間的幾何關(guān)系可進行如下表示：

其中，T=[tx ty tz]T為Pc相對Ps的空間三個維度平移矩陣，s為比例縮放系數(shù)，R是含有參數(shù)θ、?、ψ的表示旋轉(zhuǎn)偏移的單位正交矩陣。由于DEM格網(wǎng)地理坐標作為一種相對坐標，不存在偏扭、彎曲等復(fù)雜情況，根據(jù)布爾莎[24]仿射變換模型得到其簡化表達式如下：

式中，θ、?、ψ分別為繞O-Zs、O-Ys、O-Xs軸的滾轉(zhuǎn)角、偏航角和俯仰角。將式（2）帶入式（1）中可建立待匹配坐標系O-XcYcZc向基準坐標系O-XsYsZs進行空間轉(zhuǎn)換的矩陣模型：

其中，樣本個數(shù)i=m×n，m、n分別為待匹配柵格的橫、縱坐標像素量。依據(jù)最小二乘（ordinary least squares，OLS）迭代估計方法，取Z方向高差平方和最小值建立目標方程。則最小高差極值的計算公式為：

存在平面同名點時，格網(wǎng)值的權(quán)重因子wsi=1，否則wsi=0。當基準與待匹配數(shù)據(jù)分辨率不一致時，利用雙線性插值法（bilinear interpolation，BI）進行重采樣，構(gòu)造j個離散的虛擬像素pfj(xj,yj,zj)，其中zj=0。

將局部虛擬空間近似為連續(xù)二次線性曲面：

通過4 個臨近參考點確定橫縱兩個維度的曲線參數(shù)a0、a1、b0、b1，將pfj帶入式（5）計算等距插值zj=hfj(xj,yj)，將數(shù)據(jù)修改為相同尺寸。

利用最小二乘迭代求解（4）中E接近限差時的7個約束方程參數(shù)，從而構(gòu)建三維匹配的空間坐標轉(zhuǎn)換模型。

1.2 臺階融合

受不同時相、傳感幾何輻射差異等因素的影響，空洞銜接邊緣會產(chǎn)生不可避免的高程值突變[25]，針對存在的“臺階”式拼接縫，利用圖像融合技術(shù)對其進行平滑處理。

經(jīng)臺階融合算法處理得到的高程坐標Zfusion，是根據(jù)其自身每個格網(wǎng)處的高程值Zdem與同名點高程值Zref按如下規(guī)則加權(quán)平均得到：

式中，pdem表示待處理DEM 數(shù)據(jù)點的權(quán)重，與其對應(yīng)的同名參考點權(quán)值之和為1。結(jié)合離散圖像相鄰區(qū)域像素值變化小于區(qū)域內(nèi)部的特點，可利用柵格距離衰減規(guī)律對每個格網(wǎng)處的pdem進行計算。

若取窗口大小為n×n，中心點權(quán)重p(xc,yc)=0.5，其余8個方向n2-1 個點的權(quán)重和之為0.5，則待修復(fù)數(shù)據(jù)鄰域內(nèi)各點的高程梯度倒數(shù)絕對值為：

式中，像素偏移量i、j∈[0,n-1]，且i、j不同時為0，若Zdem(xc+i,yc+j)=Zdem(xc,yc)，則gdem(xc,yc,i,j)=0。

以待修復(fù)點作為掩模中心進行加權(quán)鄰域求梯度平均，可得到待處理高程權(quán)值系數(shù)：

將pdem帶入式（6），經(jīng)臺階融合輸出空洞填充邊緣的高程坐標信息，使格網(wǎng)像素能夠自然、平滑地銜接在一起。

實現(xiàn)過程如下：

（1）設(shè)定每個格網(wǎng)點初始化的權(quán)值p0如下：

i、j分別為DEM數(shù)據(jù)的行號和列號，有效初值為215，為確定待修復(fù)點的高程權(quán)值分布，需要對圖像順序、逆序兩次進行計算。

（2）定義一個偏移量offset：?。ǜ鞣较颍?quán)值重新計算的格子數(shù)默認為20。

（3）以柵格平面的左上角點作為原點，行號、列號依次遞增，遍歷尋找空洞邊緣點，取該像元鄰域內(nèi)左、左上、上、右上方向上各offset∈[1,(n-1)/2]個點并根據(jù)梯度計算依次設(shè)置權(quán)重。

（4）更改原點為像平面右下角，行號、列號依次遞減，在左下、下、右下、右方向上逆序重復(fù)上述操作，更新上一步重新計算的像元的權(quán)重pdem=min{p(xc,yc),pdem} 。

（5）取滿足pdem(x,y)＜1 的所有點，利用加權(quán)平均的方法重新計算高程值，以消除待修復(fù)數(shù)據(jù)存在的臺階問題。

2 快速DEM空洞修復(fù)方案

考慮到大尺度數(shù)據(jù)空洞修復(fù)時存在的實時性問題，首先，基于LZD 原則將表面匹配中三個維度的殘差分析簡化為取Z 方向極小值，進而利用GPU 并行處理機制，將同名點搜索和格網(wǎng)點像元值的填補過程分解至多個線程塊同步進行，以縮短處理過程的總體運算周期。

2.1 匹配模型參數(shù)獲取

為了得到最優(yōu)的表面匹配結(jié)果，通過最小二乘迭代逐步優(yōu)化模型數(shù)值，確定各參數(shù)近限差的增量矩陣dX。對變換陣Trans的7 參選取初始值X0=[R0,t0,s0]T，其中令Rψ0、R?0、Rθ0、tx0、ty0、tz0=0，s0=1 表示空間兩坐標系重合，無相對旋轉(zhuǎn)、平移，默認等比例縮放。

遍歷待匹配樣本中的每個高程點Pci=[Xci,Yci,Zci]T，帶入式（3）進行Trans變換后，在標準DEM中尋找格網(wǎng)坐標滿足xci=xsi、yci=ysi的同名點Psi。

當同名點不存在時，本算法中使用了雙線性插值構(gòu)造虛擬同名點的高程坐標，數(shù)據(jù)量小時也可考慮反距離加權(quán)平方（inverse distance weight，IDW）算法。

給定一個足夠小的步長delta，當前點在X、Y方向的一階導(dǎo)數(shù)計算公式如下：

式中，待匹配所有點滿足Zc=f(Xc,Yc)。對誤差函數(shù)e(x,y)進行一階泰勒展開，得到線性化誤差函數(shù)：

每次迭代各點誤差函數(shù)的過程就是求最優(yōu)增量參數(shù)矩陣dX*的過程：

經(jīng)系數(shù)反推，可建立增量方程組B=AdX，設(shè)計觀測系數(shù)矩陣Ai=[dzxdzy-1]T，得到Bi：

將Ai*Bi的結(jié)果存儲到A矩陣的第i行，Zci-Zsi的值存儲到B矩陣的第i行。

遍歷結(jié)束后，使用奇異值分解（singular value decomposition，SVD）方法解超定方程以計算dX增量：

并將dX*累加到變換矩陣中X=X+dX*，如果dX小于迭代限差Vx或達到最大迭代次數(shù)Va時，跳出迭代過程完成參數(shù)計算。

在匹配過程中，迭代終止條件由各參數(shù)相鄰迭代值之差來判定，即七個轉(zhuǎn)換參數(shù)的增量，七個轉(zhuǎn)換參數(shù)的限差值和設(shè)定的最大迭代次數(shù)由各個觀測值的精度決定。

匹配模型參數(shù)獲取的具體流程如圖2所示。

圖2 算法流程圖Fig.2 Algorithm flow chart

2.2 空洞修復(fù)任務(wù)的GPU并行方案設(shè)計

引入最小高差算法雖然避免了三個維度最小距離計算造成的冗余問題，但在進行同名點搜索、鄰域加權(quán)時，窮舉式的距離查找過程和大規(guī)模的DEM 數(shù)據(jù)量依然影響了空間匹配的算法效率。

因此，加入GPU并行處理技術(shù)，通過PCIe總線協(xié)助CPU 進行運算控制操作，充分發(fā)揮CPU 的復(fù)雜邏輯控制能力和GPU在處理圖像高密集計算和超長流水線任務(wù)的結(jié)構(gòu)優(yōu)勢，形成的CPU-GPU異構(gòu)模型如圖3所示。

圖3 CPU-GPU異構(gòu)模型Fig.3 CPU-GPU heterogeneous model

NVIDA 推出的標準運算平臺CUDA，代碼層面更易編寫，在性能、成本維護方面表現(xiàn)突出。本文利用其并行計算特性，提出一種多核多線程并發(fā)的DEM 匹配修復(fù)方案。其中，并行化編程最需考慮的線程同步和數(shù)據(jù)歸并問題，利用CUDA平臺線程調(diào)度和內(nèi)存的訪問管理原則，將空洞修復(fù)計算過程以核-格網(wǎng)-線程快-線程（KERNAL-GRID-BLOCK-THREAD，KGBT）的三級組織架構(gòu)進行調(diào)度配置。

構(gòu)造順序的處理邏輯步驟如下：

（1）獲取所有待匹配格網(wǎng)像元同名點坐標值；

（2）迭代計算殘差極值并求解空間映射模型的特征向量；

（3）利用加權(quán)融合算法重新計算空洞邊緣高程坐標。

獲取同名點坐標和加權(quán)融合過程需要進行大量重復(fù)計算，因此分別映射至knl_search、knl_restore 兩個核函數(shù)中，并采用單指令多數(shù)據(jù)流（single instruction multiple data，SIMD）以異步的形式執(zhí)行相同的算法指令：

首先，以程序初始化和影像數(shù)據(jù)GeoTiff（tiff），Erdas Imagine Images（img），ASCII DEM（dem）等格式的解析作為起始模塊。利用GDAL（geospatial data abstraction library）庫讀取元數(shù)據(jù)和高程信息數(shù)組取至主機端內(nèi)存。

針對基準模型大小遠大于待匹配區(qū)域的現(xiàn)象，利用OpenCv庫對存儲DEM的數(shù)據(jù)模型文件進行裁剪處理，調(diào)用RasterIO 將大圖分割為512×512 大小的塊后，將柵格數(shù)組與元數(shù)據(jù)一同拷貝至全局顯存（gloabl memory，GM），完成程序的初始化。

設(shè)置臨近點計算的輻射距離為0.2 m，采用CUDA編程模型并調(diào)用核函數(shù)knl_search，以線程id 作為檢索項，將遍歷像點的任務(wù)映射成若干并發(fā)執(zhí)行的子任務(wù)。在某一線程網(wǎng)格（Grid）的多個線程塊（Block）中對區(qū)域分割后的子圖進行同名點搜索，利用異步流迭代方法將同名點對應(yīng)的坐標集合存儲至當前Block（32×32）的共享存儲器中，供塊內(nèi)線程共用，減少訪問延遲以提高讀寫效率。在knl_search計算模型中，對于每個線程，分別搜索當前待匹配格網(wǎng)距離相近的所有點，最近距離小于等于閾值時確定為同名點對。同名點不存在時，設(shè)置橫縱縮放比例和幾何中心利用BI 插值函數(shù)構(gòu)造虛擬點，resize調(diào)整數(shù)據(jù)至相同尺寸。待所有樣本完成查詢操作后，輸出參數(shù)控制所有點對信息進行坐標反算。

迭代求解映射參數(shù)向量的過程，是采用奇異值分解（singular value decomposition，SVD）對稠密矩陣進行降維，其計算代價大，但由于運算先后關(guān)聯(lián)性較強，內(nèi)部邏輯耦合度高，不適合進行并行化改造，因此該矩陣運算仍在CPU串行域中進行。

待模型數(shù)據(jù)解析完備，針對圖中多個待修復(fù)區(qū)域進行的8鄰域權(quán)值計算，并將多個空洞分配給多個線程塊block同步進行“臺階”融合。正、逆遍歷所有待融合點，計算中心點p(xc,yc,zc)與鄰域各點高程權(quán)重關(guān)系。對拼接邊緣像素進行平滑處理，以函數(shù)返回值形式合并修復(fù)后數(shù)據(jù)至CPU內(nèi)存，將結(jié)果模型寫入輸出的DEM文件中，完成空洞修復(fù)運算。

任務(wù)基本設(shè)置：

（1）分配2個順序進行的核函數(shù)（kernel），進行同名點插值計算和臺階像素融合。將knl_search函數(shù)設(shè)置為逐點搜索，n個點對應(yīng)n個線程。等待圖像所有點遍歷結(jié)束，通知主線程迭代求解單應(yīng)性矩陣；對于knl_restore函數(shù)，若相鄰Z坐標缺失個數(shù)＞1，即判定為待修復(fù)空洞，采用流水線作業(yè)實行多個空洞邊緣并發(fā)修復(fù)。

（2）knl_search 的m個子圖分別對應(yīng)m個格網(wǎng)，將格網(wǎng)劃分為多個pixel 的線程塊（block），每個線程塊分別對應(yīng)256個CUDA線程。

（3）子圖中第k個空洞映射至knl_restore 的第k個線程塊中。

采用程序并行化關(guān)鍵解決的幾點問題：

（1）在大樣本點集查詢?nèi)蝿?wù)中，采用K近鄰法（Knearest neighbor，KNN）查詢單元格間的拓撲關(guān)系，通過兩步法確定法向最小距離點。KNN算法的搜索過程可有以下三種常用方式：窮舉法、kd樹（k-dimensional-tree）和球樹（ball-tree），根據(jù)預(yù)測平均時間隨樣本數(shù)量變化曲線知[26]，窮舉的方式更適用于容量較大的DEM數(shù)據(jù)。

（2）設(shè)Q為大小為u的n維參考點集，R為大小為v的查詢點集，通過線性掃描對k個最近鄰域點進行表決。表面匹配算法只關(guān)心距離最小點的坐標值，在此基礎(chǔ)上，對v個長度為k點集的分別進行插入排序。這種暴力實現(xiàn)的方式非常耗費運算資源，而結(jié)果集合的搜索與有序表的生成在匹配點之間不存在任何依賴關(guān)系，因此，可充分利用GPU 編程將遍歷格網(wǎng)數(shù)據(jù)進行同名點查尋詢的任務(wù)映射成若干并發(fā)執(zhí)行的子任務(wù)，提高算法效率。

（3）融合算法是根據(jù)所有點鄰域加權(quán)取均值，進行Z坐標平滑。這一過程在程序中，遍歷待修復(fù)DEM 所有點獲取權(quán)重是獨立且高重復(fù)性的，因此將任務(wù)分配至GPU運算單元重疊執(zhí)行，在并行層面分別對所有待處理點高程進行重算，將循環(huán)遍歷的過程原子化，顯著降低了內(nèi)存占用和運算時間。

CPU并行架構(gòu)下，設(shè)3維參考點集Q的樣本數(shù)為待匹配柵格大小i×j，查詢點集為v，由于核心較少導(dǎo)致能夠同時開啟的線程數(shù)目有限，可得鄰域搜索的時間復(fù)雜度為O(ijuv)。結(jié)合圖形處理單元加速運行該算法，每個樣本單獨開啟一個核函數(shù)同步處理，在相同數(shù)據(jù)量的情況下，時間復(fù)雜度減少至O(uv)。然后，對于搜索到的臨近點集，進行v個長度為k的插入排序，該操作的時間復(fù)雜度為O(vk2)，使用GPU 并開啟v個線程并行執(zhí)行可將時間復(fù)雜度降低為O(k2)。在邊緣平滑算法的優(yōu)化中，逐點高程值重算被簡化為同步進行的窗口(n×n)卷積，計算時間縮短樣本個數(shù)倍，時間復(fù)雜度為O(n2)。從算法理論角度分析該策略具備明顯的優(yōu)化性能，而系統(tǒng)實際運行的總時長，不僅需要對算法執(zhí)行指令的條數(shù)進行評估，還涉及到了線程工作量后的步驟復(fù)雜度。

圖4 為GPU 多核架構(gòu)與多線程作業(yè)下的具體實現(xiàn)流程。

圖4 GPU多核多線程空洞修復(fù)程序流程Fig.4 GPU multi-core multi-thread hole repair program

3 實驗內(nèi)容

3.1 實驗數(shù)據(jù)及平臺

本文選取的基準數(shù)據(jù)來源于奮進號航天飛機的雷達地形測繪（shuttle radar topography mission，SRTM）數(shù)據(jù)，全圖數(shù)據(jù)包含1 200×1 200個采樣點高程值，水平和垂直精度分別為±20 m和±16 m，格網(wǎng)間距90 m（3″），投影基準采用WGS84 坐標系統(tǒng)，經(jīng)緯范圍（102.4 E～104 E，29 N～30.3 N）。待修復(fù)數(shù)據(jù)均為中國境域采集的SRTM-DEM 30數(shù)據(jù)，格網(wǎng)間距30 m（1″）。

分別對待修復(fù)數(shù)據(jù)中海拔較高的青藏高原（區(qū)域1）和地勢相對平緩的長江三角洲沖積平原（區(qū)域2）的山脊數(shù)據(jù)缺失和條狀云層遮擋進行算法驗證，其中區(qū)域1空洞占比4.5%，連通空洞個數(shù)為8；區(qū)域二占比2.2%，空洞個數(shù)為4。

實驗環(huán)境：軟件平臺采用Visual Studio2010和C/C++語言進行軟件設(shè)計與編譯，通過觀測模型的收斂速度、修復(fù)總耗時、加速比等，在搭建的三個異構(gòu)混合平臺上依次進行對照組實驗，以驗證本文方法可用性和實現(xiàn)效率。工作站所用處理器（CPU）及顯卡（GPU）型號說明和預(yù)先設(shè)置匹配模型收斂條件分別如表1、表2所示。

表1 CPU/GPU參數(shù)Table 1 CPU/GPU parameters

表2 空間映射模型7參增量限差Table 2 Seven parameters for spatial mapping model

根據(jù)兩個DEM（或DSM 等）的元數(shù)據(jù)信息，設(shè)定Geotiff格式數(shù)據(jù)點平面的坐標。以標準DEM的左上角點作為原點Os(0,0)，Xs（行）方向每行間隔為經(jīng)度方向分辨率×111 000 m，Ys（列）方向間隔為緯度方向的分辨率×111 000 m。待匹配DEM 數(shù)據(jù)設(shè)定為Os(0,0)-XcYc原理同標準數(shù)據(jù)。

在修復(fù)前進行數(shù)據(jù)準備：

（1）將待匹配DEM數(shù)據(jù)用GDAL庫讀入內(nèi)存，包括元數(shù)據(jù)信息和高程數(shù)組；

（2）將標準DEM中對應(yīng)待匹配DEM區(qū)域（略大）提取為算法中的標準DEM；

（3）根據(jù)精度需求進行異源數(shù)據(jù)等間隔采樣；

（4）切分地形區(qū)域為若干（m）子圖；

3.2 實驗結(jié)果與分析

LZD 算法的模型解算精度取決于各參數(shù)限差和迭代次數(shù)，鑒于兩組地區(qū)地理面積相近，此處均采用最大迭代次數(shù)300，采樣大小500×500 pixel。當所有參數(shù)增量滿足迭代限差時，算法終止，并確定模型參數(shù)向量。

分別取待匹配DEM 中青藏高原山區(qū)、長江三角洲沖積平原兩組典型的數(shù)據(jù)模型作為對照組進行格網(wǎng)數(shù)據(jù)空洞修復(fù)，以驗證該方法在不考慮實現(xiàn)策略時，不同地勢條件下算法的可行性。其中，青藏高原地區(qū)地勢變化較大，受傳感范圍限制，在數(shù)據(jù)的采集和重構(gòu)過程中易丟失高海拔山脊處的三維信息；長江三角洲平原位于長江入?？?，該地區(qū)平均坡度較小，但由于厚云遮擋導(dǎo)致內(nèi)含兩條狹長的數(shù)據(jù)空缺。兩組待修復(fù)數(shù)據(jù)大小均為26 MB，全國基準DEM 大小為5 GB，經(jīng)模型求解和高程坐標估值計算后，得到修復(fù)后的結(jié)果如圖5所示。

圖5 漏洞修補對比圖Fig.5 Comparative diagram of vulnerability patching

由圖5可知，對于兩地區(qū)分布的不規(guī)則空洞的修復(fù)效果表現(xiàn)為：均能夠較好地恢復(fù)局部紋理特征、幾何性質(zhì)穩(wěn)定，在源數(shù)據(jù)的基礎(chǔ)上保持現(xiàn)有分辨率，并擴大了數(shù)據(jù)的空間覆蓋，實現(xiàn)了影像質(zhì)量的有效提升。

在盡可能保證精度的前提下，取區(qū)域1的數(shù)據(jù)進行計算性能測試。分別采用CPU多線程和CUDA多線程調(diào)度策略，設(shè)置大小相同的數(shù)據(jù)區(qū)域（區(qū)域1）和采集密度（500×500 pixel），依次記錄模型的匹配耗時與空洞修復(fù)總耗時。配備的三組硬件實驗環(huán)境為，平臺1：CPU為Xeon CPU E3-1230V6，顯卡采用NVIDA GeFouce GTX 980；平臺2：CPU 為Xeon CPU E3-1230V6，顯卡采用NVIDA GeFouce GTX 1080 Ti；平臺3：CPU 為Intel Core i7-7700，顯卡采用NVIDA GeFouce GTX 1080 Ti。在多種主機設(shè)備組合配置下，對照并行算法的執(zhí)行效果。

實驗結(jié)果表明，部分任務(wù)調(diào)用CPU-GPU 異構(gòu)的CUDA 編程模型進行實現(xiàn)后，修復(fù)總時長＜0.86 s，具體運行效率如表3所示。

表3 不同策略下區(qū)域1修復(fù)耗時與加速比Table 3 Time-consuming and speed-up ratio of Region 1 repair under different strategies

采用GPU 技術(shù)協(xié)助CPU 處理重復(fù)計算，縮短了點對最優(yōu)距離計算和排序時間，然而隨像素尺寸增加數(shù)據(jù)在CPU-GPU 間拷貝的開銷增大，擴容后的多線程CPU策略（平臺2）較平臺1 總耗時明顯縮減，且匹配耗時占比由20.4%增加至26.1%，串行優(yōu)勢在充足內(nèi)存容量條件下更為顯著，同時，并行計算效率也提升了3.6倍。平臺1與平臺2采用相同CPU和相異的圖形處理器，利用并行異構(gòu)策略，使總執(zhí)行加速比分別提升至4.0 和4.5，優(yōu)化后的匹配、融合耗時較CPU 多線程最多可減少71.082%和76.006%。在平臺3上運行異構(gòu)程序，算法總體修復(fù)的執(zhí)行效率相較于單線程CPU 能夠達到9.7 倍以上的加速比，具有較好的修復(fù)性能。

為了檢驗本文算法在不同地形場景下的匹配精度，分別將待修復(fù)數(shù)據(jù)中的區(qū)域1和區(qū)域2與全國基準數(shù)據(jù)進行LZD 三維匹配。其中，青藏山區(qū)和沖積平原數(shù)據(jù)大小為各取50 個點作為測試樣本，與基準數(shù)據(jù)對照計算Z坐標估值殘差，得到的最小參數(shù)增量、誤差矩陣和Z方向精度均方誤差（單位：m）如表4所示。

表4 迭代最終參數(shù)增量與誤差Table 4 Iterate over final parameter increment and error

以區(qū)域1的山脊為例，可看到模型計算的高程值在空洞邊緣存在“臺階”形斷層，為完善修復(fù)效果，對該空洞區(qū)域進行了加權(quán)融合。取空洞邊緣點Z坐標值得加權(quán)平均，進行DEM高程值重算，再次賦值后達到奇異點去除和邊緣平滑的效果，空洞平滑修補前后對比如圖6所示。

圖6 空洞平滑修補前后對比圖Fig.6 Comparative diagram of cavity before and after smooth repair

如圖6（a）所示，融合前的高海拔數(shù)據(jù)缺失面，邊沿分布著非連續(xù)的像素條帶，在數(shù)據(jù)準確度和視覺效果層面都存在明顯缺陷。經(jīng)加權(quán)融合糾正后的該區(qū)域見圖6（b），保留了主要地形特征，未見明顯修復(fù)痕跡，達到趨近于基準數(shù)據(jù)的完善要求。

4 結(jié)束語

機載雷達（airborne radar，AR）在高海拔、崎嶇地勢和極端氣候時出現(xiàn)的測繪盲區(qū)，導(dǎo)致衍生出的數(shù)字高程數(shù)據(jù)（DEM）存在數(shù)據(jù)空洞等的問題，現(xiàn)有的空間插值和表面匹配算法為該空缺修復(fù)研究提供了有力支撐。為了進一步滿足準確、高效的應(yīng)用需求，本文提出了一種改進的快速LZD三維匹配算法結(jié)合“臺階”融合的修復(fù)方法。借助異源輔助數(shù)據(jù)，在搭建的不同平臺下，分別對2 組差異明顯的非結(jié)構(gòu)化地形樣本進行了方法驗證。實驗結(jié)果表明本文采用的方法，修復(fù)后的圖像信息豐富、空洞邊緣平滑，誤差控制在±37.725 m，符合地形數(shù)據(jù)的處理需求。在模型精度滿足固定限差的條件下匹配效率成倍提升，較CPU 多線程架構(gòu)運行耗時可減少76%，加速比最優(yōu)可達9.7倍，在大尺度數(shù)字高程模型修復(fù)領(lǐng)域具有明顯優(yōu)勢。