2

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

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

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

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

1 地形復(fù)雜度指標構(gòu)建

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

1.1 地形因子篩選

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

地形起伏度是定量描述地表形態(tài)的重要指標，等于區(qū)域內(nèi)的最大高差，具體算式為：

R=Hmax-Hmin。

(1)

式中，R代表地形起伏度，Hmax、Hmin分別為TIN結(jié)構(gòu)中某局部中心點與其關(guān)聯(lián)鄰近點中的最大、最小高程。

地形坡度等于坡面垂直高度和水平方向距離的比。在TIN結(jié)構(gòu)下，局部區(qū)域的坡度可以由中心點關(guān)聯(lián)的每個三角面坡度的均值表示，三角面坡度的定義為：

(2)

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

地形粗糙度又稱局部褶皺度，一般由區(qū)域?qū)嶋H面積與相應(yīng)投影面積之比計算得到[12]。由于地表實際面積難以計算，因此在TIN結(jié)構(gòu)中，常用中心點關(guān)聯(lián)的所有三角面面積代表區(qū)域?qū)嶋H面積，則地形粗糙度Kr定義為：

(3)

式中，Ss為局部實際面積，S為局部投影面積，cosa由坡度經(jīng)三角函數(shù)轉(zhuǎn)換求得，N為中心點關(guān)聯(lián)的三角面?zhèn)€數(shù)。

1.2 地形因子權(quán)重分配

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

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

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

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

(4)

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

(5)

(6)

(7)

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

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

(8)

則基于改進CRITIC法構(gòu)建地形復(fù)雜度指標的步驟如下。

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

ii. 根據(jù)CRITIC賦權(quán)法的基本思想，構(gòu)造融合對比強度與沖突性的綜合系數(shù)

(9)

式中，k為地形因子個數(shù)，dCorrij為第i個與第j個地形因子間的距離相關(guān)系數(shù)，vj為第j個地形因子的變異系數(shù)，Cj越大表示第j個地形因子對地形復(fù)雜度的貢獻越大，應(yīng)賦予更大的權(quán)重。

iii. 求取各地形因子的歸一化權(quán)重

(10)

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

(11)

2 點云抽稀

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

2.1 基于地形復(fù)雜度抽稀

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

2.2 特征點保留

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

1) 設(shè)計一種用來保留局部范圍內(nèi)最深點與最淺點的自適應(yīng)格網(wǎng)分區(qū)算法，具體步驟為：

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

(12)

式中，a為計算保留點個數(shù)的經(jīng)驗系數(shù)，本研究取值0.02。

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

(13)

iii. 遍歷所有點云的坐標，按照其位置索引存入所在網(wǎng)格單元內(nèi)；

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

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

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

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

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

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

2.3 算法流程

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

1) 利用離散點云數(shù)據(jù)構(gòu)建TIN結(jié)構(gòu)，獲得水深點與鄰近關(guān)聯(lián)點的空間關(guān)系；

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

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

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

5) 通過自適應(yīng)格網(wǎng)分區(qū)，保留局部范圍內(nèi)最深點與最淺點；

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

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

3 實驗與分析

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

3.1 實驗區(qū)域概述

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

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

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

3.2 實驗結(jié)果分析

利用本算法構(gòu)建兩個實驗區(qū)域的地形復(fù)雜度模型，得到兩個區(qū)域的地形復(fù)雜度指標T與各地形因子原始值之間的關(guān)系分別為：

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)

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

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

表2 區(qū)域1抽稀結(jié)果誤差統(tǒng)計Tab. 2 Error statistics of thinning results in area 1

表3 區(qū)域2抽稀結(jié)果誤差統(tǒng)計Tab. 3 Error statistics of thinning results in area 2

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

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

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

圖5 區(qū)域1實驗結(jié)果暈渲圖Fig. 5 Hill shading of area 1 experimental results

圖6 區(qū)域2實驗結(jié)果暈渲圖Fig. 6 Hill shading of area 2 experimental results

4 結(jié)論

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