顧及地形復雜度因子權重的多波束點云抽稀算法

2

(1.山東科技大學 測繪與空間信息學院，山東 青島 266590； 2.自然資源部 海洋測繪重點實驗室，山東 青島 266590)

海底地形測量是海洋測繪的基礎性任務，其獲取的海底地形地貌信息是支撐海洋開發活動的基本資料[1]。多波束測深系統以寬覆蓋、高精度的優勢，大幅提高了海底地形的測量效率和質量，已成為水深測量任務中首選的探測設備[2]。由于多波束點云數據具有海量性的特點，大量冗余的水深點在參與地形分析運算時會造成時間和硬件資源的浪費，為海底地形數據的管理、存儲與成圖工作帶來運算壓力。因此，進行多波束點云數據的簡化研究具有重要的現實意義。

目前常用的點云抽稀算法大致可以歸納為不顧及地形特征和顧及地形特征兩類。不顧及地形特征的算法主要包括基于規則格網的抽稀和基于系統規則的抽稀[3]，基本原理是依據一定的隨機采樣規則進行數據簡化。此類算法側重于提升運算效率而忽略了地形特征，地形精度損失較大，容易導致局部細節的丟失。許多學者提出了多種顧及地形特征的算法，如基于不規則三角網(triangulated irregular network，TIN)的抽稀算法[4]，通過計算局部所有三角形的法向量，根據最大與最小法向量的夾角大小進行數據點取舍。也有學者將地形特征因子與TIN結合，如基于TIN與坡度的算法[5]通過比較某點關聯三角面的最大坡度差與設定閾值的大小進行抽稀，其閾值判斷標準還可以是高差[6]、曲率[7]等地形特征因子?；谄露褥氐某橄8]在此基礎上引入熵理論，將局部區域的坡度熵作為表征地形豐富和變化程度的指標。但是，僅采用單一特征參數進行地形變化的分析較為片面，閾值選定對于抽稀質量的影響較大，并且在一些特殊地形點會判斷錯誤[9]?；诙嘣匦翁卣鞯某橄10]利用主成分分析與多因子理論，提取幾類貢獻度較大的特征因子來描述地形復雜度，抽稀質量相對于采用單一地形因子算法有所改善，但復雜度經驗公式存在優化空間，地形精度仍有待提高。

因此，提出一種顧及地形復雜度因子權重的多波束點云數據抽稀算法，本算法基于點云構建的不規則三角網結構提取局部地形起伏度、坡度和粗糙度因子，引入改進CRITIC(criteria importance though intercrieria correlation)法確定三者權重，從而構建一種地形復雜度指標作為抽稀判定標準，并將局部最深最淺點、邊界點視為特征點予以保留。通過對比實驗表明，本算法的抽稀結果在精度上優于傳統簡化算法，提升了地形表達的真實性和完整性。

1 地形復雜度指標構建

為了減小數據壓縮過程中的地形精度損失，點云取舍標準的確定至關重要。本研究將地形復雜度指標作為抽稀依據，首先考慮提取多個地形特征因子，以克服單一參數量化的局限性；然后遵循獨立、有效、易于計算[11]的原則篩選所需因子，保證復雜度指標的構建可行性與表達準確性；最后融合地形因子，關鍵在于權重系數的分配，本研究引入CRITIC多因子綜合評價法，以獲得最優定權結果。

1.1 地形因子篩選

地形復雜度是反映整體或局部地球表面起伏和褶皺程度的一類指標，可用地形特征因子表征，常見因子包括坡度、坡向、地表曲率、高程標準差、地形起伏度、地形粗糙度、地形分形維數[11]等。在篩選時，可以遵循相對獨立性、因子有效性、易于局部計算三個原則。其中，局部地形的起伏程度受極值影響較大，地形起伏度代表局部最大高程與最小高程之差，無均值參與運算，能夠有效反映地形變化；基于地形坡面特征計算出的坡度、坡向存在較大相關性，坡向多用于氣象、生態領域研究，而坡度用來表征地表曲面的傾斜程度，被認為是最重要且有效的地形因子；地表曲率、地形粗糙度可對地形曲面特征進行有效描述，而TIN結構中曲率一類的參數難以計算，因此不宜選擇。綜上考慮，本研究選用地形起伏度、地形坡度、地形粗糙度作為構建地形復雜度指標的三個評價因子。

地形起伏度是定量描述地表形態的重要指標，等于區域內的最大高差，具體算式為：

R=Hmax-Hmin。

(1)

式中，R代表地形起伏度，Hmax、Hmin分別為TIN結構中某局部中心點與其關聯鄰近點中的最大、最小高程。

地形坡度等于坡面垂直高度和水平方向距離的比。在TIN結構下，局部區域的坡度可以由中心點關聯的每個三角面坡度的均值表示，三角面坡度的定義為：

(2)

式中，nx、ny、nz為標準矢量n三個方向的法向量，標準矢量n是方向向外并垂直于平面的直線。

地形粗糙度又稱局部褶皺度，一般由區域實際面積與相應投影面積之比計算得到[12]。由于地表實際面積難以計算，因此在TIN結構中，常用中心點關聯的所有三角面面積代表區域實際面積，則地形粗糙度Kr定義為：

(3)

式中，Ss為局部實際面積，S為局部投影面積，cosa由坡度經三角函數轉換求得，N為中心點關聯的三角面個數。

1.2 地形因子權重分配

確定參與評價的三個地形因子后，各因子的權重分配將極大影響地形復雜度模型的可靠性。賦權方法主要分為主觀和客觀賦權法，主觀賦權法脫離實測數據，而客觀賦權法通過數學方法從原始數據中提取指標間的信息，避免主觀因素干擾，具有嚴密的數學理論基礎，適用于地形復雜度指標的評價。其中，大多客觀賦權法憑借單一標準進行權數確定，如主成分分析法、變異系數法、熵權法等，而CRITIC法同時考慮了評價指標間的對比強度與沖突性，其評價結果相對其他方法更為全面客觀[13]。

CRITIC法采用Pearson積矩相關系數描述指標間的沖突性，相關系數越小，沖突性越高，所占權重越大；對比強度反映樣本在同一指標下的離散程度，采用標準差σ表示，σ越大，對應權重越大。雖然傳統CRITIC法的賦權方式已經較為完善，但一方面使用Pearson系數只限于評價指標之間的線性關系，并不能準確衡量實際相關關系[14]，且樣本需滿足正態分布，數據適用性低；另一方面，三種地形因子的量綱、單位、均值不同，常規采用標準差衡量對比強度失去意義，不再適用。因此，本研究對傳統CRITIC法進行兩方面改進。

1) 采用距離相關系數(distance correlation)代替Pearson相關系數進行衡量沖突性的改進。距離相關系數是一種新的度量指標相關性的方法，其優勢在于能夠實際描述指標間的任意相關關系[15]，且無需滿足任何模型假設和參數要求。

距離相關系數基于距離協方差dCov(u,w)和距離方差dCov(u,u)dCov(w,w)計算，則指標u與指標w之間的距離相關系數定義為：

(4)

式中，dCov2(u,w)=S1+S2-2S3，S1、S2、S3的計算公式分別為：

(5)

(6)

(7)

同理可計算dCov(u,u)和dCov(w,w)。

同時，為提高針對海量多波束點云的計算效率，可采用快速算法[16]計算該系數。

(8)

則基于改進CRITIC法構建地形復雜度指標的步驟如下。

i. 根據式(4)與式(8)計算出各地形因子的距離相關系數dCorr和變異系數v；

ii. 根據CRITIC賦權法的基本思想，構造融合對比強度與沖突性的綜合系數

(9)

式中，k為地形因子個數，dCorrij為第i個與第j個地形因子間的距離相關系數，vj為第j個地形因子的變異系數，Cj越大表示第j個地形因子對地形復雜度的貢獻越大，應賦予更大的權重。

iii. 求取各地形因子的歸一化權重

(10)

iv. 構建地形復雜度指標時，由于變異系數消除了各地形因子的量綱影響，為了表示地形因子原始值與復雜度指標的關系，需要進行反無量綱化處理，從而得到地形復雜度指標T與各地形因子原始值的關系式：

(11)

2 點云抽稀

點云抽稀的主要目的是降低數據冗余，但實現這一前提需要盡可能保證簡化后的數據質量，最大程度保持原始地形特征與細節，從而便于快速構建高精度的海底數字高程模型(digital elevation model)，提高海底地形可視化效果。另外，針對多波束測深數據的特性，在抽稀過程中需要考慮保留部分特征點以兼顧航行安全與海圖生產業務的需求。因此，本研究首先根據地形復雜度進行初步抽稀，然后對局部最深與最淺點[3]、邊界點[17]予以保留，實現多波束數據的高保真抽稀，使結果數據滿足上述實際應用要求。

2.1 基于地形復雜度抽稀

將1.2節構建的地形復雜度指標T作為點云取舍的主要判斷準則，依次求取所有水深點的地形復雜度T，由于T值代表該點所在局部區域的地形起伏和褶皺程度，因此按T值大小對所有水深點進行排列，抽稀時根據實際需要設置合適的閾值，優先剔除T值較小的冗余點，從而保留并輸出對地形變化貢獻較大的水深點，實現考慮地形特征的水深點初步簡化。

2.2 特征點保留

在水深點簡化過程中，一些對地形變化貢獻度不大的點會被剔除掉，而此類數據對構建真實海底地形模型、滿足海圖制圖條件以及后續數據分析都具有一定作用。結合多波束水深數據的特點，將局部范圍內的最深與最淺點、測區邊界點視為兩類特征點，并進行提取和保留。

1) 設計一種用來保留局部范圍內最深點與最淺點的自適應格網分區算法，具體步驟為：

i. 根據式(12)，由總點云數量p動態計算二維網格行列數q，然后遍歷所有點云坐標數據求出X、Y方向最大與最小坐標值(xmax、xmin、ymax、ymin)。

(12)

式中，a為計算保留點個數的經驗系數，本研究取值0.02。

ii. 依據式(13)求出網格單元分別在X、Y方向的坐標間隔，并建立q×q大小的矩形區域；

(13)

iii. 遍歷所有點云的坐標，按照其位置索引存入所在網格單元內；

iv. 遍歷所有網格單元，當網格單元至少存在兩組數據時，保留網格內的最深與最淺點，若網格單元只包含一組水深點數據則直接保留，其余水深點刪除。

圖1 提取邊界點示意圖Fig. 1 Diagram of extracting boundary points

圖2 抽稀算法具體流程Fig. 2 Specific flow of thinning algorithm

該算法視測區數據量與點云分布情況不同，對原始測深點云進行自適應格網劃分，從而將需要保留的最深、最淺點數量控制在合理范圍內，有助于補充完善整體區域的細節特征，同時避免了大面積數據缺失的情況，兼顧了多波束抽稀的地形完善性準則[3]。

2) 采用Alpha-shapes算法提取整體區域的邊界點(圖1)。在離散水深點集內，設置一個判別半徑r，過任意兩點繪制半徑為r的滾動圓，判斷此圓內部是否存在其他點云，若不存在，則認為這兩個點是邊界點。

統計提取的兩類地形特征點，將特征點集與2.1節計算的初步簡化點集合并，并剔除重復點，輸出最終抽稀結果。

2.3 算法流程

圖2展示了改進算法的流程，可概括為以下七個主要步驟：

1) 利用離散點云數據構建TIN結構，獲得水深點與鄰近關聯點的空間關系；

2) 遍歷所有水深點，提取各點的三種地形因子信息；

3) 基于改進CRITIC法構造地形復雜度指標作為點云取舍的判定準則；

4) 計算各點地形復雜度，保留復雜度大于特定簡化率下的閾值的點，實現數據初步簡化；

5) 通過自適應格網分區，保留局部范圍內最深點與最淺點；

6) 提取原始點云數據的邊界點，與步驟5)結果合并為特征點集；

7) 將步驟4)、6)的結果數據合并，輸出所有保留的水深點。

3 實驗與分析

選擇由多波束系統獲取的真實海底點云作為實驗數據，將本算法和其他傳統算法的抽稀結果對比分析，驗證提出算法的可行性和適用性。

3.1 實驗區域概述

實驗選取的海底地形數據采用R2Sonic-2024型號的多波束測深系統采集，于2017年4月完成外業測量，測區位于中國舟山群島東南部附近海域，實驗選取地形變化較為明顯的兩個局部區域作為研究對象(圖3)。其中，截取的區域1地形起伏相對較大，區域2地形較為粗糙，二者的基本信息如表1所示。實驗所用數據均為經過原始數據解析、預處理的點云結果數據。

圖3 實驗區域示意圖Fig. 3 Schematic illustration of the experimental area

表1 實驗點云數據基本信息表Tab. 1 Basic information of experimental point cloud data

3.2 實驗結果分析

利用本算法構建兩個實驗區域的地形復雜度模型，得到兩個區域的地形復雜度指標T與各地形因子原始值之間的關系分別為：

T1=0.548 9×R+6.016 2×S+0.018 5×Kr，

(14)

T2=0.578 9×R+2.394 5×S+0.040 5×Kr。

(15)

為檢驗本算法的抽稀效果，實驗設置五檔簡化率(被剔除的點云數目與點云總數目的比值)代表不同的點云簡化程度，選取未改進算法、基于坡度[5]和基于高差[6]算法進行對比。其中，未改進算法與本算法的區別在于采用傳統CRITIC法構建地形復雜度指標。

首先，采用檢查點法的思想對抽稀結果進行精度評定，即按照隨機原則均勻抽取部分原始數據作為檢查點，將這些點位在TIN結構下的線性內插值與實際水深值逐一作差得到誤差，并統計所有點位的均方根誤差(root mean square error，RMSE)，RMSE值越小說明精度越高。兩個實驗區域均隨機抽取了原始數據的2%作為檢查點，實驗結果如表2、表3所示，圖4可以直觀地反映出各算法在不同簡化率條件下的精度變化情況。

表2 區域1抽稀結果誤差統計Tab. 2 Error statistics of thinning results in area 1

表3 區域2抽稀結果誤差統計Tab. 3 Error statistics of thinning results in area 2

圖4 不同算法精度對比Fig. 4 Accuracy comparison of different methods

比較表2、表3精度評定數據發現，針對兩個實驗區域，各算法的抽稀結果在整體趨勢上，隨著簡化率的增大，RMSE值逐漸增大，點云質量與簡化程度呈負相關關系；在所有簡化率條件下，本算法相比未改進算法的RMSE值均有所減小，點云抽稀質量更好；在簡化率最小(簡化率=19.0%)，即簡化程度最低時，本算法與坡度算法在兩個實驗區域下的RMSE數值大小非常接近，抽稀質量基本相當，而高差算法的誤差相對較大，在其他四檔簡化率條件下，本算法的RMSE值均為最小，特別是簡化率達到70%以上時，抽稀質量相比高差算法與坡度算法具有顯著優越性；進一步分析誤差曲線(圖4)，隨著簡化率增大，本算法的誤差值上升趨勢相對高差與坡度算法較為平緩。

其次，從主觀視覺角度觀察和分析本算法在地形真實性表達上的優劣。本研究以簡化率76.3%條件下的抽稀結果為例，基于Surfer軟件對框選區域i、ii、iii、iv(圖5(a)、圖6(a))的點云進行插值生成1 m×1 m格網數據，并渲染得到局部地形暈渲圖(圖5(b)、圖6(b))，可以發現，在局部區域ii、iii、iv內，基于高差抽稀的地形圖表現出較為嚴重的失真現象，坡度算法的抽稀結果在局部區域i、iii中也出現一定程度的異常。此外，對于其他局部區域，坡度與高差算法的地表較為光滑，部分地形特征的表達不明顯甚至消失，這是由于二者僅采用單一抽稀判定標準，部分對地形變化貢獻大的點被誤刪，且未能全面考慮特征點的影響，其抽稀結果存在局部點云缺失的情況，致使生成的格網模型與實際地形不符；未改進算法的整體表達較好，但對于部分特征的呈現相較本算法仍有所欠缺。綜上，本算法的圖像能夠更好地描述細節特征，在保持原始地表形態方面優于其他算法。

圖5 區域1實驗結果暈渲圖Fig. 5 Hill shading of area 1 experimental results

圖6 區域2實驗結果暈渲圖Fig. 6 Hill shading of area 2 experimental results

4 結論

點云抽稀是多波束測深數據處理的重要環節之一，為了提高點云簡化質量，以滿足后續地形表達與處理的需要，提出一種顧及地形復雜度因子權重的多波束點云抽稀算法。實驗分析表明，本研究通過綜合評價法對多個因子進行綜合分析，構建的地形復雜度指標能夠很好地匹配地形特征，相比采用單一地形因子，該算法的點云取舍標準更為合理。此外，算法考慮局部最深最淺點與邊界點的影響，確保了點云數據的整體性，并在一定程度上起到完善地形細節的作用。本算法在降低數據冗余的同時，相較傳統算法提高了地形精度，其結果在保持地形真實性上具有優越性，且算法在不同地形都表現出良好的適用性。