李海森，陸丹，魏玉闊，么彬

(哈爾濱工程大學(xué) 水聲技術(shù)國家級重點實驗室，黑龍江 哈爾濱 150001)

多波束測深聲吶作為當前海洋勘測的重要工具之一，被廣泛應(yīng)用于海洋資源調(diào)查、海底地形測量、水下打撈等領(lǐng)域.該聲吶獲得的厘米級分辨率測深數(shù)據(jù)經(jīng)濾波處理后，可進行海底地形建模(digital terrain modeling，DTM)，實現(xiàn)海底地形可視化［1-3］.其中等深線圖能夠有效反映海底地形的高低起伏特征，是一種兼顧海底地形定性與定量信息的高層次表達手段.從海底地形中生成等深線圖是多波束測深聲吶海底成圖軟件的最重要功能之一.多波束測深聲吶具有寬覆蓋、高分辨、大面積掃海測量特點，隨著多波束測深數(shù)據(jù)量的不斷增大，如何提高大數(shù)據(jù)量海底地形的等深線生成速度成為必須面對的問題，也是國產(chǎn)多波束測深聲吶必須解決的緊迫問題之一.

海底地形等深線作為一種深度屬性等值線，其生成方法與常規(guī)等值線生成方法類似.多年來，眾多學(xué)者對等值線快速生成方法進行了研究［4-8］.在格網(wǎng)或三角網(wǎng)拓撲關(guān)系已知條件下，等值線生成過程通常可分為2部分:1)遍歷格網(wǎng)或三角網(wǎng)確定等值線搜索起點;2)從起點開始利用拓撲關(guān)系進行等值線跟蹤直至生成完整的閉合或非閉合等值線.這兩部分共同決定了等值線生成效率.研究人員通常采用各種索引結(jié)構(gòu)來提高第一部分中搜索起點確定的效率［4-7］，采用顧及方向的策略簡化網(wǎng)格邊上等值點存在的判斷條件來提高第二部分中等深線跟蹤的效率［7］.然而，索引結(jié)構(gòu)雖然能夠提高搜索起點的確定效率但其構(gòu)建和查詢比較復(fù)雜并且占用空間較大;顧及方向的策略雖然在一定程度上簡化了等值點判斷條件，但是這種策略只適用于規(guī)則格網(wǎng)DTM，對不規(guī)則三角網(wǎng)(triangulated irregular network，TIN)結(jié)構(gòu)的DTM并不具適用性.為同時提高這兩部分的效率，實現(xiàn)等深線快速生成，本文提出一種基于等深值索引序列的等深線快速生成算法.該算法在提高等深線生成速度的同時更加適合計算機編程實現(xiàn)，并且能同時適用于規(guī)則格網(wǎng)和不規(guī)則三角網(wǎng)結(jié)構(gòu)DTM的等深線快速生成.目前，該算法已在國產(chǎn)首臺便攜式多波束測深聲吶的海底成圖軟件中成功實現(xiàn)，并在吉林松花湖以及湖北清江平洛湖水下地形等深線圖顯示中得到應(yīng)用.

1 常規(guī)等深線生成

常規(guī)等深線可以從規(guī)則格網(wǎng)和TIN2種結(jié)構(gòu)的DTM中生成.由于多波束測深聲吶獲得的測深數(shù)據(jù)呈不規(guī)則離散分布，且TIN結(jié)構(gòu)的海底地形模型直接利用原始深度數(shù)據(jù)，避免了內(nèi)插的精度損失，更適合建立TIN海底地形模型，因此本文重點研究從TIN海底地形模型生成等深線.文獻［3］給出了使用野值剔除后的離散多波束測深數(shù)據(jù)進行TIN海底地形建模的方法.以此為基礎(chǔ)，常規(guī)的等深線生成方法按三角形順序進行搜索，基本過程如下［9］:

1)對給定的等深值h，與所有數(shù)據(jù)點zi(i=1，2，…，n)進行比較，若zi=h，則將zi作降一微量(如10－4)處理［10-11］.

2)設(shè)立三角形標志數(shù)組，設(shè)置其初始值為0，每一個元素與一個三角形對應(yīng)，凡處理過的三角形將標志置為1，以后不再處理，直至等深值改變.

3)按順序判斷每一個三角形的三邊中的2條邊是否有等深線穿過.若三角形一邊的2端點為P1(x1，y1，z1)，P2(x2，y2，z2)，則采用如下公式判斷該邊上是否存在等深點:

若上式成立，則說明該條邊上存在等深點.獲得搜索起點，通過線性內(nèi)插計算等深點在該邊上的平面坐標.

4)搜索該等深線在該三角形的出口邊，也就是相鄰三角形的入口邊，并內(nèi)插該等深點的平面坐標.

5)進入相鄰三角形，重復(fù)第4)步，進行等深線跟蹤.若等深點回到起點，則形成閉合等深線;若等深點到達三角網(wǎng)邊界，則從搜索起點開始反向跟蹤，直至到達另一邊界，形成非閉合等深線.

6)當一條等深線全部跟蹤完后，將其輸出.然后繼續(xù)三角形的搜索，直至全部三角形處理完成.改變等深值，重復(fù)以上過程，直到全部等深線繪制完成為止.

2 基于等深值索引序列的等深線快速生成算法

在常規(guī)等深線生成算法中，要初始化三角形標志數(shù)組，首先要判斷三角形中包含哪些等深值.這就需要查找哪些等深值在三角形頂點的最小和最大深度值范圍內(nèi).由于三角網(wǎng)中一個點總是被多個三角形共用，顯然，常規(guī)算法在獲取三角形中所含等深值時，同一數(shù)據(jù)點的深度值會在不同的三角形中被重復(fù)計算.另外，由于實測的深度值為浮點型數(shù)據(jù)，當DTM中數(shù)據(jù)量很大時，必然需要大量浮點計算，以致等深線生成速度無法滿足多波束海底成圖系統(tǒng)的實時顯示需求.因此，本文提出一種基于等深值索引序列的等深線快速生成算法.算法先將浮點型深度屬性三角網(wǎng)映射為整型索引屬性三角網(wǎng)，根據(jù)索引序列和深度范圍的映射關(guān)系避免了標志數(shù)組初始化中同一數(shù)據(jù)點在不同三角形中的重復(fù)計算過程，然后提出一種基于0/1異或的等深線出口邊判斷方法，將等深線跟蹤過程建立在整型運算的基礎(chǔ)上，通過降低計算機運算量來實現(xiàn)等深線快速生成.由于算法運算量的降低與DTM的具體結(jié)構(gòu)無關(guān)，該算法同時適用于規(guī)則格網(wǎng)和不規(guī)則三角網(wǎng)結(jié)構(gòu)DTM的等深線快速生成.

2.1 深度與索引屬性的映射

DTM數(shù)據(jù)深度與索引屬性的映射是將絕大部分運算從浮點運算轉(zhuǎn)化為整型運算的關(guān)鍵，也是等值線快速生成算法的基礎(chǔ)，主要包括設(shè)置整型索引屬性和對每個數(shù)據(jù)點增加索引屬性值2部分，處理過程如下:

1)根據(jù)等深值設(shè)置深度范圍索引.假設(shè)Zmin和Zmax為DTM的最小和最大深度值，V1，V2，…，Vn(以升序排列)為待顯示的n個等深值.顯然，Zmin≤V1＜V2＜…≤Zmax，可得到個遞增的深度范圍，用一個從0到n的整數(shù)序列表示:C={0，1，2…，n}.

2)建立等深值與深度范圍索引的直接映射關(guān)系.如圖1所示，建立索引序列C={0，1，2，…，n}與深度范圍的直接映射關(guān)系:C0=0，映射［Zmin，V1］深度范圍;Ci=i(i=1，2，…，n－1)，映射［Vi，Vi+1］深度范圍;Cn=n，映射(Vn，Zmax］深度范圍.

圖1 索引序列與深度范圍的映射關(guān)系Fig.1 Relation of the index array and depth ranges

3)給DTM中每個數(shù)據(jù)點添加索引屬性，并將浮點型深度屬性三角網(wǎng)轉(zhuǎn)化為整型索引屬性三角網(wǎng).例如，圖2(a)中最大深度Zmin=22.3，最小深度Zmax= 27.2，待顯示等深值為 V1=23.0，V2=25.0，V3= 27.0，則可得索引序列C={0，1，2，3}，轉(zhuǎn)化為整型索引屬性三角網(wǎng)后如圖2(b)所示.由于索引序列與等深值存在直接映射關(guān)系，因此無需任何計算過程即可從索引屬性三角網(wǎng)中直接得到任意三角形邊上包含的等深值.三角形A的edge1上包含3個等深值分別為V1，V2和V3;edge2上包含2個等深值V2和V3; edge3上包含一個等深值V1.而三角形B的3個頂點索引值相等，因此該三角形中不包含任何等深值.顯然，將深度屬性三角網(wǎng)轉(zhuǎn)化為整型索引屬性三角網(wǎng)后，無需繁瑣地對每個三角形頂點的深度值進行比較以判斷其中是否包含等深值或包含了哪些等深值，大大減少了計算機運算量和程序復(fù)雜度.

圖2 深度屬性三角網(wǎng)轉(zhuǎn)換為索引屬性三角網(wǎng)Fig.2 Turn depth attribute triangulation net to index attribute triangulation net

2.2 基于0/1異或的等深線走向判斷方法

在深度屬性映射為索引屬性的基礎(chǔ)上，進一步研究如何快速獲取等深線在三角形中走向的問題.

首先考慮最簡單的情況，即三角形頂點索引值只為0或1的情況:假設(shè)三角形頂點的索引值分別為Cdot1，Cdot2Cdot3，且最小索引值為0，最大索引值為1.顯然，只要三角形邊端點索引值的異或結(jié)果為1，該邊上必然存在等深點.因此，在已知等深線入口邊條件下，只需判斷三角形第三個頂點的索引值與入口邊任一端點索引值的異或結(jié)果是否為1即可判斷等深線出口邊.如圖3所示，假設(shè)等深線入口邊為三角形的edge1，若Cdot1∧Cdot3=1，則判斷等深線出口邊為三角形的edge3;否則，等深線出口邊為三角形的edge2.

圖3 基于0/1異或的等深線走向判斷Fig.3 Contour trend decision based on 0/1 exclusive-OR

然后延伸到一個三角形內(nèi)通過多條等深線的情況:假設(shè)三角形頂點的最小索引值為Cmin=i，最大索引值為Cmax=j，通過索引值與深度范圍的直接映射關(guān)系可知該三角形內(nèi)包含Vi+1到Vj的(j－i)個等深值.要獲得Vk(i+1≤k＜j)等深線在三角形內(nèi)的走向，需進行如下處理將頂點索引值轉(zhuǎn)化為0或1:

圖4 3條等深線在三角形內(nèi)的走向判定Fig.4 Trend decision of three contours in a triangle

2.3 等深線快速生成

本文提出的等深線生成算法主要由2部分組成:一是確定三角形中的等深線起點并初始化三角形標志數(shù)組;二是等深線跟蹤.為提高等深線生成速度，算法先將深度屬性映射為索引屬性，在索引屬性三角網(wǎng)上直接獲取等深線起點并初始化標志數(shù)組，然后從起點開始使用0/1異或方法快速判斷等深線走向進行等深線跟蹤.算法基本流程如圖5所示.

圖5 等深線快速生成算法流程圖Fig.5 Flow chart of fast algorithm for contours generation

在整型索引屬性三角網(wǎng)上，遍歷所有三角形，當三角形頂點的最大最小索引值之差d=Cmax－Cmin≠0時，表示該三角形內(nèi)含有d個等深值.將當前三角形號、所含等深值數(shù)量d、d個等深值索引及其相應(yīng)的標志按順序存入標志數(shù)組中.然后，以標志數(shù)組中對應(yīng)的三角形為起始三角形進行等深線跟蹤.依次遍歷當前三角形內(nèi)的d個等深值，若當前等深值標志為1，表示該等深值已被搜索過，跟蹤該三角形的下個等深值;若當前等深值標志為0，表示該等深值未被搜索過，進行等深線跟蹤:先在起始三角形內(nèi)獲得起始等深點，使用0/1異或方法判斷等深線在三角形內(nèi)走向，獲得出口點，以該出口點作為出口邊相鄰三角形的新入口點搜索新出口點，如此反復(fù)，直至等深線到達邊界或返回起點.若等深線到達邊界，則從搜索起點進行反向跟蹤，直至到達另一邊界，形成一條非閉合等深線;若等深線返回起點，則形成一條閉合等深線.三角形標志數(shù)組中的所有三角形均處理完畢后，該DTM的等深線生成完成.

3 算法運算量分析與測試結(jié)果

3.1 算法運算量分析

假設(shè)海底DTM中包含N個數(shù)據(jù)點，則最多能構(gòu)成2N－5個三角形［12］，當N較大且數(shù)據(jù)點均勻分布時，取三角形個數(shù)T≈2N，等深值數(shù)量為k，k值可根據(jù)顯示效果需求來選取.

在三角形標志數(shù)組初始化階段，常規(guī)算法可分為2個步驟:1)查找數(shù)據(jù)點深度值zi(i=1，2，…，N)，是否等于等深值hj(j=1，2，…，k)，若相等，對zi作降微量處理.設(shè)降微量處理操作的運算量為c1，則該步驟的運算量為Nc1lbk.2)遍歷T個三角形，根據(jù)每個三角形的最小最大深度值zmin和zmax查找升序等深值序列 {hi}(i=1，2，…)，k中滿足zmin＜hj＜zmax條件的hj，其中h=a，a+1，…，b且1≤a≤b≤k，然后標記該三角形所含等深值并初始化.設(shè)在k個等深值中查找hj的運算量為c2lbk，初始化該三角形等深值標志的運算量為c3，則該步驟的運算量共為T(c2lbk+c3).取T≈2N，常規(guī)算法進行三角形標志數(shù)組初始化的總運算量近似為N(c1lbk+2c2lbk+ 2c3).

本文算法的三角形標志數(shù)組初始化也可分為2個步驟:1)將N個數(shù)據(jù)點的深度屬性轉(zhuǎn)化為索引屬性，即查找升序等深值序列 {hi}(i=1，2，…，k)中滿足z≤hj(j=1)，hj＜z≤hj+1(j=1，2，…，k－1)或z＞ hj(j=k)條件的等深值索引j，其中z為當前數(shù)據(jù)點深度值.設(shè)在k個等深值中查找z位置的運算量為c4lbk，則該步驟的運算量為Nc4lbk.2)遍歷T個三角形，根據(jù)三角形頂點的最大最小索引值直接獲取等深值索引并初始化標志數(shù)組.設(shè)初始化一個三角形等深值標志的運算量為c5，顯然，c5=c3，因此該步驟的運算量減少為Tc3.取T≈2N，本文算法的三角形標志數(shù)組初始化總運算量近似為N(c4lbk+2c3).

由于c1、c2、c3和c4都大于0，且在升序等深值序列{hi}中查找滿足zmin＜hj＜zmax(j=a，a+1，…，b且1≤a≤b≤k)的 hj所需運算量必然大于在{hi}(i=1，2，…，k)中查找滿足z≤hj(j=1)，hj＜z≤hj+1(j=1，2，…，k－1，)或(z＞hj(j=k))的j所需運算量，即c2lbk＞c4lbk，因此，N(c1lbk+2c2lbk+ 2c3)大于N(c4lbk+2c3).由此可知，三角形標志數(shù)組初始化運算量隨數(shù)據(jù)量的增加呈線性增長，但本文算法運算量的增長斜率小于常規(guī)算法，這意味著本文算法的低運算量優(yōu)勢隨著數(shù)據(jù)量的增加而更加突出.

在等深線跟蹤階段，常規(guī)算法根據(jù)式(1)需要2次浮點數(shù)減法，1次浮點數(shù)乘法，1次關(guān)系運算來判斷等深線在三角形內(nèi)的出口邊.而本文算法利用數(shù)據(jù)點的整型索引屬性僅使用2次整型條件表達式操作，1次異或操作即可判斷出口邊.由于每次出口邊判斷都能節(jié)省部分運算量且計算機處理整型數(shù)據(jù)的速度高于處理浮點型數(shù)據(jù)的速度，因此在等深線跟蹤階段，本文算法所需運算量必然小于常規(guī)算法，且所節(jié)省運算量隨數(shù)據(jù)量的增加而增大.

為驗證本文算法正確性與低運算量優(yōu)勢，根據(jù)運算量與運行時間的正比關(guān)系，使用實際測深數(shù)據(jù)對兩種算法的運行時間進行測試并給出測試結(jié)果.

3.2 算法測試結(jié)果

算法使用VC++和OpenGL實現(xiàn)，實驗的硬件環(huán)境為Intel(R)Core(TM)i3，2.93 GHz，3.46 GB內(nèi)存.實驗數(shù)據(jù)為吉林松花湖區(qū)域和湖北清江平洛湖區(qū)域的實際水深數(shù)據(jù).

圖6和圖7給出了本文算法對上述區(qū)域?qū)崪y水深數(shù)據(jù)進行等深線生成得到的等深線圖.圖6中松花湖區(qū)域DTM數(shù)據(jù)量為12×104，最小與最大深度值為7.13 m和62.67 m，設(shè)置等深值范圍為8～62 m，等深間隔1 m，等深值數(shù)量k=55，等深線生成耗時約112.8 ms.圖7中平洛湖區(qū)域DTM數(shù)據(jù)量為14.5×104，最小與最大深度值為 10.3 m 和110.9 m，設(shè)置等深值范圍為15～110 m，等深間隔5 m，等深值數(shù)量 k=20，等深線生成耗時約51.3 ms.由于圖6中設(shè)置的k值大于圖7中的k值，也就是說圖6中等深線分布比圖7更加密集，因此雖然松花湖區(qū)域DTM的數(shù)據(jù)量小于平洛湖區(qū)域，但是其等深線生成耗時大于平洛湖區(qū)域.

圖6 松花湖部分區(qū)域等深線圖Fig.6 Contour chart of part area of Songhua Lake

圖7 平洛湖部分區(qū)域等深線圖Fig.7 Contour chart of part area of Pingluo Lake

圖8給出了6組不同數(shù)據(jù)量實驗數(shù)據(jù)的本文算法和常規(guī)算法運行時間測試結(jié)果.測試數(shù)據(jù)的深度值在6.64～66.23 m之間，設(shè)置最小與最大等深值范圍為8 m與66 m，等深間隔為2 m，即等深值數(shù)量k=30.

圖8 本文算法和常規(guī)算法運行時間比較Fig.8 Runtime comparison of the algorithm in this paper and a general algorithm

測試結(jié)果表明本文算法在三角形標志數(shù)組初始化階段和等深線跟蹤階段所需運行時間均小于常規(guī)算法.從圖8(a)可以看出，三角形標志數(shù)組初始化運行時間呈線性增長，本文算法運行時間的增長斜率小于常規(guī)算法的增長斜率，符合運算量線性增長和斜率c4lbk+2c3＜c1lbk+2c2lbk+2c3的規(guī)律.從圖8(b)可以看出，等深線跟蹤運行時間也隨數(shù)據(jù)量呈線性增長，本文算法的運行時間總是低于常規(guī)算法.圖8(c)給出的2種算法的等深線生成總運行時間表明:1)本文算法的等深線生成總運行時間低于常規(guī)算法;2)等深線生成運行時間隨數(shù)據(jù)量呈線性增長;3)本文算法的運行時間增長斜率小于常規(guī)算法.顯然，與常規(guī)算法相比，本文算法的優(yōu)勢隨著DTM數(shù)據(jù)量的增加而不斷增大，尤其適合大數(shù)據(jù)量海底地形的等深線快速生成.

4 結(jié)論

多波束掃海測量數(shù)據(jù)量已達到百萬、千萬乃至海量量級，本文針對大數(shù)據(jù)量海底地形等深線生成算法進行了深入研究，提出一種基于等深值索引序列的等深線快速生成算法.通過實測海底地形數(shù)據(jù)處理分析，得出以下結(jié)論:

1)通過深度屬性到索引屬性的映射，將等深線生成過程建立于整型索引屬性三角網(wǎng)的基礎(chǔ)上，降低了算法的復(fù)雜度并為基于0/1異或的等深線走向快速判定提供了條件，從而同時降低了等深線搜索起點確認和等深線跟蹤的運算量.

2)算法的低運算量優(yōu)勢隨著數(shù)據(jù)量的增大更加突出，因此該算法尤其適用于數(shù)據(jù)量達到百萬級、千萬級乃至海量的多波束測深地形等深線快速生成.

3)算法同時適用于規(guī)則格網(wǎng)數(shù)據(jù)結(jié)構(gòu)和不規(guī)則三角網(wǎng)數(shù)據(jù)結(jié)構(gòu)DTM的等深線深快速生成，對DTM數(shù)據(jù)結(jié)構(gòu)具有更廣泛的適用性.

此外，由于海底地形等深線是一種深度屬性的等值線，其生成方法與大部分等值線生成方法類似，因此本文算法不僅適用于海底地形的等深線快速生成，同樣能夠應(yīng)用于其他屬性等值線，如等高線、等壓線、等溫線等的快速生成，具有廣泛的適用性.

［1］李海森，陳寶偉，么彬，等.多子陣高分辨海底地形探測算法及其FPGA和DSP陣列實現(xiàn)［J］.儀器儀表學(xué)報，2010，31(2):281-286.

LI Haisen，CHEN Baowei，YAO Bin，et al.Implementation of high resolution sea bottom terrain detection method based on FPGA and DSP array［J］.Chinese Journal of Scientific Instrument，2010，31(2):281-286.

［2］么彬，李海森，周天，等.多子陣超寬覆蓋海底地形探測方法試驗研究［J］.哈爾濱工程大學(xué)學(xué)報，2008，29(10): 1076-1081.

YAO Bin，LI Haisen，ZHOU Tian，et al.Test study of a multiple sub-array super-wide-coverage seabed terrain detection method［J］.Journal of Harbin Engineering University，2008，29(10):1076-1081.

［3］LU D，LI H S，WEI Y K，et al.An improved merging algorithm for Delaunay meshing on 3D visualization multibeam bathymetric data［C］//IEEE ICIA 2010.Harbin，China，2010:1170-1176.

［4］KREVELD M V.Efficient methods for isoline extraction from a TIN［J］.International Journal of Geographical Information Science，1996，10(5):523-540.

［5］王濤，劉紀平，毋河海.基于排序預(yù)處理的等高線生成算法［J］.測繪學(xué)報，2006，35(11):390-394.

WANG Tao，LIU Jiping，WU Hehai.The extraction of contour lines from grid DEM based on sorting［J］.Acta Geodaetica et Cartographica Sinica，2006，35(11):390-394.

［6］顧濱兵，楊兆海，高宇，等.基于網(wǎng)格和端點量化的等值線生成算法［J］.吉林大學(xué)學(xué)報:信息科學(xué)版，2010，28 (1):89-94.

GU Binbing，YANG Zhaohai，GAO Yu，et al.Algorithm of contour lines drawing based on grid using technique of grid scan and point digital［J］.Journal of Jilin University:Information Science Edition，2010，28(1):89-94.

［7］王濤，毋河海，劉紀平.基于區(qū)間樹索引的等高線生成算法［J］.武漢大學(xué)學(xué)報:信息科學(xué)版，2007，32(2):131-134.

WANG Tao，WU Hehai，LIU Jiping.An algorithm for extracting contour lines based on interval tree from grid DEM［J］.Geomatics and Information Science of Wuhan University，2007，32(2):131-134.

［8］JONES N L，KENNARD M J，ZUNDEL A K.Fast algorithm for generating sorted contour strings［J］.Computers＆Geosciences，2000，26:831-837.

［9］鄔倫，劉瑜，張晶，等.地理信息系統(tǒng)——原理、方法和應(yīng)用［M］.北京:科學(xué)出版社，2005:215-216.

［10］遲寶明，李治軍，葉勇，等.基于GIS的地下水水位等值線圖自動生成算法研究［J］.吉林大學(xué)學(xué)報:地球科學(xué)版，2007，37(2):261-265.

CHI Baoming，LI Zhijun，YE Yong，et al.Study on the arithmetic of automatically drawing isoline of groundwater level based on GIS［J］.Journal of Jilin University:Earth Science Edition，2007，37(2):261-265.

［11］張立華.港口工程水深測量中的不規(guī)則區(qū)域、任意方向、不等尺度網(wǎng)格法追蹤等深線研究［J］.水運工程，2005(4):9-13.

ZHANG Lihua.A study on drawing depth contours using grids based on irregular area，arbitrary direction and different scale in bathymetric survey for port engineering［J］.Port＆Water Engineering，2005(4):9-13.

［12］周培德.計算幾何——算法分析與設(shè)計［M］.3版.北京:清華大學(xué)出版社，2008:89-91.