張志衡，彭認(rèn)燦，黃文騫,3，董 箭，劉國輝

1. 海軍大連艦艇學(xué)院海洋測繪系，遼寧 大連 116018； 2. 海軍大連艦艇學(xué)院海洋測繪工程軍隊重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116018； 3. 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054； 4. 海軍出版社，天津 300450

當(dāng)前，多波束測深系統(tǒng)已成為海底地形探測的主要設(shè)備，但在測量過程中由于儀器自噪聲、海況因素、設(shè)備參數(shù)設(shè)置不合理及海中生物的影響，導(dǎo)致多波束測深數(shù)據(jù)不可避免地存在較多粗差，嚴(yán)重影響了海底地形的真實(shí)表達(dá)[1-2]。為了提高多波束測深數(shù)據(jù)的精度，需要對其進(jìn)行濾波處理。依據(jù)多波束測深數(shù)據(jù)粗差探測手段的不同，多波束測深數(shù)據(jù)濾波方法可大致劃分為交互式濾波和自動濾波兩類[2-3]。交互式濾波是一個主觀識別和標(biāo)定錯誤的過程，即在三維可視的環(huán)境下，由作業(yè)人員人工標(biāo)定測深數(shù)據(jù)中的異常值，濾波效率相對較低，且嚴(yán)重依賴作業(yè)人員的主觀經(jīng)驗(yàn)；為克服交互式濾波的缺點(diǎn)，國內(nèi)外學(xué)者通過將數(shù)理統(tǒng)計、誤差處理、地形分析等理論和方法融入粗差探測技術(shù)中，提出了諸多自動濾波方法[4]。文獻(xiàn)[5]提出了一種基于加權(quán)平均移動模型的深度濾波方法，在假設(shè)海底地形連續(xù)的條件下識別多波束條帶內(nèi)標(biāo)準(zhǔn)偏差異常高于平均值的水深點(diǎn)，并運(yùn)用2σ或3σ準(zhǔn)則對多波束數(shù)據(jù)進(jìn)行濾波。文獻(xiàn)[6]將趨勢面濾波的思想引入到多波束測深數(shù)據(jù)異常值判定中，利用一定范圍內(nèi)的水深值z和平面坐標(biāo)(x,y)進(jìn)行多項式趨勢面擬合，并結(jié)合2σ或3σ準(zhǔn)則對異常值進(jìn)行判定。文獻(xiàn)[2]針對基于函數(shù)/統(tǒng)計推值的比較判別法在檢驗(yàn)異常數(shù)據(jù)時出現(xiàn)的異常值“遮蔽”現(xiàn)象，將抗差估計理論引入到檢驗(yàn)異常值的推值比較法中，提出了一種基于抗差估計的選權(quán)迭代插值比較檢驗(yàn)法，采用簡單加權(quán)平均數(shù)模型作為判別異常值的推值模型，在被檢測點(diǎn)鄰域內(nèi)選取正確觀測數(shù)據(jù)，拒絕不正確或者帶有粗差的數(shù)據(jù)，提高了檢測方法的抗差性。文獻(xiàn)[7—12]提出了一種基于不確定度理論的多波束測深數(shù)據(jù)濾波算法——CUBE(combined uncertainty and bathymetry estimator)算法，該算法主要利用測深數(shù)據(jù)的深度信息和不確定度信息，通過卡爾曼濾波方法剔除數(shù)據(jù)中可能存在的粗差，濾波效果較為理想。文獻(xiàn)[3]提出了一種基于中值濾波、局部方差檢測和小波分析理論的聯(lián)合濾波方法，該方法中的局部方差檢測模型避免了中值濾波對海底地形細(xì)節(jié)的破壞，小波分析理論增強(qiáng)了數(shù)據(jù)的抗差性和算法的可靠性。文獻(xiàn)[13]引入了觀測值識別變量的概念，構(gòu)建了識別變量均值漂移模型，提出了一種基于識別變量后驗(yàn)概率的粗差探測Bayes方法。在此基礎(chǔ)上，文獻(xiàn)[14]將該粗差探測Bayes方法應(yīng)用于多波束測深數(shù)據(jù)異常值探測中，通過構(gòu)造水深觀測識別變量判別模型對異常值的后驗(yàn)概率進(jìn)行分析，判斷和識別水深觀測數(shù)據(jù)中的異常值。文獻(xiàn)[15]針對最小二乘支持向量機(jī)(LS-SVM)在構(gòu)造海底趨勢面過程中存在的數(shù)值解稀疏性較差以及樣本數(shù)據(jù)偏差引起的擬合曲面形態(tài)失真等問題，提出了一種基于水深不確定度調(diào)控的趨勢面濾波方法，利用水深不確定度調(diào)控模型對訓(xùn)練樣本數(shù)據(jù)進(jìn)行形態(tài)優(yōu)化，構(gòu)造出反映海底實(shí)際變化的趨勢面，實(shí)現(xiàn)對測深異常值的有效剔除。

在諸多自動濾波方法中，由于趨勢面濾波法設(shè)計思想簡單、預(yù)設(shè)參數(shù)較少、易于并行計算等特點(diǎn)，受到了海洋測繪學(xué)者的廣泛關(guān)注[7,14-17]。然而由于海底地形的復(fù)雜多變以及局部曲面擬合范圍的不確定性，趨勢面濾波法存在其固有的技術(shù)缺陷，主要體現(xiàn)在：①由于海底地形復(fù)雜多變，局部曲面擬合范圍不易確定，單一的曲面擬合函數(shù)無法全面反映海底的實(shí)際地形，極易出現(xiàn)粗差濾除不徹底或誤將正常水深判定為粗差點(diǎn)進(jìn)而被濾除的情況；②當(dāng)局部曲面擬合范圍內(nèi)存在較大的突變粗差點(diǎn)時，將會導(dǎo)致正常水深點(diǎn)與所擬合的趨勢面之間產(chǎn)生較大偏差，進(jìn)而影響正常水深點(diǎn)的判定結(jié)論，當(dāng)局部曲面擬合范圍內(nèi)存在的粗差點(diǎn)相對集中且成簇出現(xiàn)時，趨勢面濾波法將無法檢測出所有粗差點(diǎn)，相反可能出現(xiàn)將粗差點(diǎn)誤判為正常水深點(diǎn)保留，而將正常水深點(diǎn)誤判為粗差點(diǎn)濾除的情況；③趨勢面濾波法的構(gòu)建前提是海底地形的連續(xù)性，當(dāng)海底出現(xiàn)斷裂點(diǎn)(如斷崖)或障礙物(如暗礁或沉船)時，趨勢面濾波法可能導(dǎo)致斷裂點(diǎn)或障礙物被不合理濾除的情況。為解決上述問題，本文提出了一種基于自然鄰點(diǎn)影響域的多波束測深數(shù)據(jù)趨勢面濾波改進(jìn)算法。該算法在引入散亂水深點(diǎn)局部最小范圍——自然鄰點(diǎn)影響域的概念的基礎(chǔ)上，通過構(gòu)造影響域內(nèi)特定局部坐標(biāo)系下的統(tǒng)一曲面擬合函數(shù)，進(jìn)行傳遞式迭代趨勢面濾波，逐步濾除影響正常水深點(diǎn)判定的粗差數(shù)據(jù)，并針對突變地形邊界點(diǎn)難以確定的問題，根據(jù)邊界點(diǎn)在其相鄰鄰域地形內(nèi)連續(xù)性不一致的特性，建立了面向突變地形邊界點(diǎn)的判斷準(zhǔn)則，最后通過實(shí)驗(yàn)驗(yàn)證了改進(jìn)算法的有效性。

1 傳統(tǒng)趨勢面濾波法及其存在問題

1.1 傳統(tǒng)趨勢面濾波法構(gòu)建原理

趨勢面濾波法是一種以曲面擬合函數(shù)為數(shù)值分析基礎(chǔ)的誤差處理算法[16-18]。其算法構(gòu)建原理為：依據(jù)波束腳印的深度值和平面位置坐標(biāo)，構(gòu)造出反映海底地形變化趨勢的多項式曲面函數(shù)，計算和統(tǒng)計實(shí)測水深值與所構(gòu)造趨勢面間的深度偏離量，結(jié)合誤差處理理論建立粗差數(shù)據(jù)的判定準(zhǔn)則，實(shí)現(xiàn)多波束測深數(shù)據(jù)的自動濾波[16-18]。趨勢面濾波法中曲面擬合函數(shù)的一般形式為[19]

(1)

式中，z表示波束腳印的深度值；(x,y)表示波束腳印的平面位置坐標(biāo)；alm表示各多項式系數(shù)；n表示多項式總階數(shù)；Q(xq,yq,zq)表示以待檢測點(diǎn)q(xq,yq,zq)為中心的曲面擬合函數(shù)f(x,y)的局部曲面擬合范圍。趨勢面濾波法的關(guān)鍵在于多項式總階數(shù)n與局部曲面擬合范圍Q的選擇[4]。相關(guān)研究表明：局部曲面擬合范圍Q選擇過小且多項式總階數(shù)n選擇過高，可能存在粗差點(diǎn)未被有效濾除的情況；范圍Q過大且總階數(shù)n過低，可能引起海底有效地形信息的丟失[4]。一般情況下，多項式總階數(shù)n與局部曲面擬合范圍Q由擬合范圍內(nèi)的海底地形復(fù)雜度f來確定[4]。

在假設(shè)局部海底地形復(fù)雜度適中且變化不大的前提下，局部曲面擬合范圍Q可以待檢測點(diǎn)q(xq,yq,zq)為中心且包含20～30個水深點(diǎn)的區(qū)域范圍代替，多項式總階數(shù)n控制在3階以內(nèi)(1≤n≤3)[15,20-21]。為更好地說明問題，本文以多項式總階數(shù)n=2為例，闡述趨勢面濾波法的濾波原理。

將n=2代入式(1)，并將其改寫為AX=L的誤差方程形式，其中

(2)

X=(ATA)-1ATL

(3)

根據(jù)局部曲面擬合范圍Q內(nèi)的平面坐標(biāo)(xi,yi)和水深值zi擬合出趨勢面后，趨勢面濾波的異常值判斷準(zhǔn)則為[4]

(4)

式中，k=2或3；σ表示局部曲面擬合范圍Q內(nèi)各實(shí)測水深值與所構(gòu)造趨勢面間的深度偏離量均方差。

1.2 傳統(tǒng)趨勢面濾波法存在的主要問題

由1.1的分析可知，趨勢面濾波法的關(guān)鍵在于多項式總階數(shù)n與局部曲面擬合范圍Q的匹配。傳統(tǒng)趨勢面濾波法中上述相關(guān)參數(shù)的經(jīng)驗(yàn)取值盡管在一定程度上簡化了算法的分析過程，提高了多波束測深數(shù)據(jù)的濾波效率，但不同經(jīng)驗(yàn)參數(shù)的取值將引起局部曲面擬合范圍Q內(nèi)水深點(diǎn)位分布與數(shù)量的較大差異，進(jìn)而導(dǎo)致曲面擬合函數(shù)與實(shí)際海底地形匹配的不確定性問題。因此，定量分析多項式總階數(shù)n與局部曲面擬合范圍Q之間的匹配關(guān)系，建立顧及海底地形復(fù)雜度的多項式曲面函數(shù)，將是該濾波法改進(jìn)方向之一。

其次，傳統(tǒng)趨勢面濾波法在構(gòu)造趨勢面過程中所采用的擬合數(shù)據(jù)可能含有粗差，該粗差數(shù)據(jù)點(diǎn)可能單個出現(xiàn)，也可能多個成簇出現(xiàn)[19,21]。如圖1和圖2所示，圖1為粗差點(diǎn)單個出現(xiàn)的情況，假設(shè)P1點(diǎn)為局部曲面擬合范圍Q內(nèi)的單個粗差點(diǎn)且其粗差值較大，則由范圍Q內(nèi)的水深點(diǎn)所構(gòu)造的趨勢面將在P1點(diǎn)處產(chǎn)生較大的形態(tài)突變，使P2和P3點(diǎn)與所構(gòu)造的趨勢面間產(chǎn)生較大的偏差ΔZ2和ΔZ3，從而誤將正常點(diǎn)P2和P3判定為粗差點(diǎn)；圖2為粗差點(diǎn)多個成簇出現(xiàn)的情況，假設(shè)P1、P2和P3點(diǎn)為局部曲面擬合范圍Q內(nèi)的多個成簇出現(xiàn)的粗差點(diǎn)群且P1點(diǎn)的粗差值遠(yuǎn)大于P2和P3點(diǎn)的粗差值，則所構(gòu)造的趨勢面將在P1點(diǎn)處產(chǎn)生較大的形態(tài)突變，可能造成P2和P3點(diǎn)與趨勢面間的偏差小于限差，被誤判為正常點(diǎn)，而正常點(diǎn)P4、P5和P6與所構(gòu)造的趨勢面間產(chǎn)生較大的偏差ΔZ4、ΔZ5和ΔZ6，從而誤將正常點(diǎn)P4、P5和P6判定為粗差點(diǎn)。這些問題產(chǎn)生的原因在于擬合海底地形表面的點(diǎn)中含有粗差，擬合出的趨勢面不能合理反映海底地形變化。因此，研究如何從所選范圍內(nèi)選取可信度較高的水深點(diǎn)，拒絕不正確或帶有粗差的可疑點(diǎn)，以此來擬合地形趨勢面，將是該濾波法改進(jìn)的另一方向。

圖1 單個粗差點(diǎn)影響濾波示意圖Fig.1 The situation of single gross error point

圖2 多個粗差點(diǎn)成簇影響濾波示意圖Fig.2 The situation of conglomerate multiple gross error points

最后，由于傳統(tǒng)趨勢面濾波法是在假設(shè)海底地形連續(xù)的前提下進(jìn)行的，當(dāng)海底出現(xiàn)局部突變地形(斷崖或障礙物)且其偏差大于限差時，可能出現(xiàn)突變地形整體或部分水深點(diǎn)被不合理剔除的情況。如圖3所示，P1和P2點(diǎn)為突變地形的邊界點(diǎn)，P3和P4點(diǎn)為突變地形上的局部連續(xù)點(diǎn)。當(dāng)局部曲面擬合范圍Q大于突變地形點(diǎn)群時，若突變地形點(diǎn)群占局部曲面擬合范圍Q內(nèi)總點(diǎn)數(shù)的比重較小，對所構(gòu)造的趨勢面影響較小，可能使得P1—P4點(diǎn)與所構(gòu)造的趨勢面間的偏差ΔZ1—ΔZ4均大于限差，則突變地形點(diǎn)群將被誤當(dāng)作成簇出現(xiàn)的粗差點(diǎn)群，導(dǎo)致突變地形被不合理剔除；若突變地形點(diǎn)群占擬合范圍Q內(nèi)總點(diǎn)數(shù)的比重較大，對所構(gòu)造的趨勢面影響較大，可能使得突變地形中部的P3和P4點(diǎn)與所構(gòu)造的趨勢面間的偏差ΔZ3和ΔZ4小于限差，而P1和P2點(diǎn)與所構(gòu)造的趨勢面間的偏差ΔZ1和ΔZ2大于限差，導(dǎo)致突變地形的邊界點(diǎn)P1和P2被不合理剔除，使海底地形不能精確地被表達(dá)。因此，有效區(qū)分和辨別成簇出現(xiàn)的粗差點(diǎn)群和突變地形點(diǎn)群及保留海底突變地形中的邊界點(diǎn)，進(jìn)一步提高海底地形表達(dá)精度，將是該濾波法改進(jìn)的又一方向。

圖3 突變地形點(diǎn)群影響濾波示意圖Fig.3 The situation of point group on the mutation terrain

2 趨勢面濾波改進(jìn)算法

2.1 局部曲面擬合范圍及擬合函數(shù)的確定

2.1.1 局部曲面擬合范圍的確定

參照文獻(xiàn)[23]，對于曲面S中任意一點(diǎn)p處的局部曲面規(guī)范形式定義為

(5)

式中，(x*，y*，z*)為局部曲面在局部坐標(biāo)系的三維坐標(biāo)；k1(p)和k2(p)分別為p點(diǎn)相應(yīng)的兩個主曲率；該局部坐標(biāo)系以p點(diǎn)為原點(diǎn)，以兩主曲率的方向?yàn)閤軸和y軸的方向，以局部曲面的法向量為z軸方向，o((x*)2+(y*)2)為高階無窮小量。

對式(5)分析可知，當(dāng)k1(p)和k2(p)同號且均不為零時，局部曲面近似為橢圓拋物面；當(dāng)k1(p)和k2(p)異號且均不為零時，局部曲面近似為雙曲拋物面；當(dāng)k1(p)和k2(p)有且只有一個為零時，局部曲面近似為拋物柱面。因此，任意曲面S在點(diǎn)p處的局部曲面為一個二次拋物曲面。

鑒于曲面S中任意一點(diǎn)p處的局部曲面規(guī)范形式相對簡單，且嚴(yán)格微分意義下的局部曲面地形復(fù)雜度近似為零，若采用任意一點(diǎn)周圍的若干點(diǎn)所組成的局部擬合范圍近似代表微分意義下的局部曲面，則該局部離散點(diǎn)擬合范圍應(yīng)具有最小的地形復(fù)雜度。根據(jù)文獻(xiàn)[19]中對地形復(fù)雜度f與間距為d的所有地形點(diǎn)的自相關(guān)系數(shù)R(d)間關(guān)系的描述可知，自相關(guān)系數(shù)R(d)越大，地形復(fù)雜度f越小。自相關(guān)系數(shù)R(d)的曲線可由高斯函數(shù)來描述，其公式為

(6)

式中，cov(d)是間距為d的所有地形點(diǎn)之間的協(xié)方差；V為由所有地形點(diǎn)計算得到的方差；c為cov(d)接近零時的相關(guān)間距。

由式(6)可知，自相關(guān)系數(shù)R(d)與各點(diǎn)間的間距d有關(guān)，d越小，R(d)越大。因此，為了降低擬合范圍內(nèi)的地形復(fù)雜度f，可在待檢測點(diǎn)q處選擇各方向距離最近的點(diǎn)組成局部曲面擬合范圍Q，該擬合范圍內(nèi)的地形復(fù)雜度f最小。而在Voronoi圖非結(jié)構(gòu)化網(wǎng)格中，q點(diǎn)各方向距離最近的點(diǎn)為其對應(yīng)的Voronoi單胞的鄰接離散點(diǎn)，被稱為自然鄰點(diǎn)[24-27]，由每一離散點(diǎn)的自然鄰點(diǎn)所構(gòu)成的局部范圍可稱為該點(diǎn)的自然鄰點(diǎn)影響域。根據(jù)多波束測深系統(tǒng)海底跟蹤控制技術(shù)的相關(guān)知識可知，一般認(rèn)為在任意方向至少有連續(xù)記錄的3～5個測深點(diǎn)才可以判斷海底是否存在凹(凸)地形[27]，因此自然鄰點(diǎn)影響域也可認(rèn)為是海底地形表達(dá)的最小局部范圍。因此，由自然鄰點(diǎn)影響域組成的擬合范圍可近似代表微分意義下的局部曲面。對于如何快速確定任意離散點(diǎn)的自然鄰點(diǎn)影響域，可參照文獻(xiàn)[28]或文獻(xiàn)[29]中描述的方法，根據(jù)Delaunay三角網(wǎng)構(gòu)建基本原則——最大最小角準(zhǔn)則或空外接圓特性來快速尋找局部范圍內(nèi)任意一點(diǎn)的自然鄰點(diǎn)。

由上述分析可知，若采用自然鄰點(diǎn)影響域作為趨勢面擬合范圍，可將復(fù)雜多變的海底地形表面分割成若干個由自然鄰點(diǎn)影響域構(gòu)成的局部簡單地形，且地形形態(tài)單一，均可采用算法相對簡單的低階固定多項式曲面函數(shù)進(jìn)行擬合。此外，自然鄰點(diǎn)影響域作為海底地形表達(dá)的最小局部范圍，使得海底地形易于精細(xì)化表達(dá)，更有利于海底微地貌或微小障礙物(如暗礁或沉船)的識別和判定。

2.1.2 曲面擬合函數(shù)的確定

通常自然鄰點(diǎn)影響域范圍內(nèi)的點(diǎn)數(shù)較少(4≤m<27)[21]，若仍采用完整的二次多項式來擬合自然鄰點(diǎn)影響域內(nèi)的二次拋物面，將因點(diǎn)數(shù)不夠而導(dǎo)致擬合過程終止。為解決此問題，需要減少曲面擬合函數(shù)中的多項式系數(shù)，對自然鄰點(diǎn)影響域范圍內(nèi)的曲面擬合函數(shù)進(jìn)行統(tǒng)一確定。

當(dāng)局部曲面所在的局部坐標(biāo)系的各坐標(biāo)軸方向按照式(6)定義后，曲面擬合函數(shù)可以被簡化成z*=a(x*)2+b(y*)2的形式，但在實(shí)際操作過程中，局部曲面的兩個主曲率方向很難直接確定，而其法向量方向可直接近似獲得，若在p點(diǎn)建立以p點(diǎn)為原點(diǎn)，以局部曲面的法向量為z軸方向，x軸和y軸方向任意設(shè)置的新局部坐標(biāo)系，假設(shè)新局部坐標(biāo)系下的三維坐標(biāo)為(x′,y′,z′)，則兩個局部坐標(biāo)系各點(diǎn)的關(guān)系為

(7)

式中，γ0為z軸旋轉(zhuǎn)角度。

不考慮式(5)中的高階無窮小量，將式(7)代入式(5)中，可得任意離散點(diǎn)的自然鄰點(diǎn)影響域內(nèi)的局部曲面在所構(gòu)建的新局部坐標(biāo)系下的二次拋物面多項式函數(shù)為

z′=a′(x′)2+b′x′y′+c′(y′)2

(8)

式中，a′、b′、c′表示整理后的各多項式系數(shù)，為待求參數(shù)。

由于式(8)所采用的新局部坐標(biāo)系與多波束測深數(shù)據(jù)所采用的全局坐標(biāo)系存在差異，因此，需要將多波束測深數(shù)據(jù)轉(zhuǎn)換至新局部坐標(biāo)系下的數(shù)據(jù)，才可采用式(8)進(jìn)行曲面擬合。

對于空間坐標(biāo)系轉(zhuǎn)換可采用變換矩陣來描述。假設(shè)全局坐標(biāo)系為A(x,y,z)，局部坐標(biāo)系為B(x*,y*,z*)，p點(diǎn)坐標(biāo)為(x0，y0，z0)，坐標(biāo)系轉(zhuǎn)換過程為將全局坐標(biāo)系A(chǔ)的原點(diǎn)平移至p點(diǎn)，使兩坐標(biāo)系原點(diǎn)重合，繞局部坐標(biāo)系原點(diǎn)p旋轉(zhuǎn)A的坐標(biāo)軸，使其z軸與B的z軸重合即可[30]。整個過程的數(shù)學(xué)描述為

B=TA+Δ

(9)

式中，Δ為坐標(biāo)系平移矩陣；T為坐標(biāo)系旋轉(zhuǎn)矩陣，T=Tx×Ty×Tz，各矩陣的計算公式詳見文獻(xiàn)[30]。

在旋轉(zhuǎn)過程中，z軸不進(jìn)行旋轉(zhuǎn)，而對于x軸和y軸的旋轉(zhuǎn)角度α和β，可采用局部曲面的法向量來求得。首先對自然鄰點(diǎn)影響域內(nèi)的數(shù)據(jù)用最小二乘法進(jìn)行平面擬合，用擬合平面的法向量N(nx,ny,nz)來近似代替p點(diǎn)的法向量，則可得旋轉(zhuǎn)角度α和β為

(10)

最后采用式(8)進(jìn)行曲面擬合，由于曲面擬合過程中只存在3個未知數(shù)，因此，自然鄰點(diǎn)影響域內(nèi)的點(diǎn)數(shù)(m≥4)完全滿足曲面擬合要求。

2.2 傳遞式迭代趨勢面濾波方法

針對擬合數(shù)據(jù)中含有粗差而導(dǎo)致所構(gòu)造的趨勢面與實(shí)際地形趨勢面之間存在較大偏差這一問題，借鑒抗差M估計濾波方法[2]，本文提出傳遞式迭代趨勢面濾波的方法。如圖4所示，由離散點(diǎn)集N0、N1和N2所組成的范圍分別是以P0、P1和P2點(diǎn)為中心點(diǎn)的連續(xù)3個自然鄰點(diǎn)影響域，假設(shè)N0點(diǎn)集的自然鄰點(diǎn)影響域已進(jìn)行趨勢面濾波處理，以N1點(diǎn)集的自然鄰點(diǎn)影響域?yàn)槔醾鬟f式迭代趨勢面濾波法的具體過程如下：

(1) 首先將N0點(diǎn)集中已經(jīng)確定含有粗差的水深點(diǎn)利用其擬合的趨勢面進(jìn)行改正，這樣可保證N0∩N1點(diǎn)集中的水深點(diǎn)不含粗差。

(2) 然后用N1點(diǎn)集的自然鄰點(diǎn)影響域內(nèi)所有的水深點(diǎn)來擬合地形表面，找出其中大于2倍(或3倍)中誤差的水深點(diǎn)，將其存入可疑點(diǎn)集Ω中。

(3) 用(N1-Ω)點(diǎn)集的水深點(diǎn)繼續(xù)擬合，重復(fù)步驟(2)，直到無水深點(diǎn)存入可疑點(diǎn)集Ω中，用(N1-Ω)點(diǎn)集所擬合的地形趨勢面達(dá)到穩(wěn)定狀態(tài)，得到最終的海底地形趨勢面。

(4) 最后計算可疑點(diǎn)集Ω中的水深點(diǎn)與最終海底地形趨勢面之間的最大偏差，以2倍(或3倍)中誤差為判斷標(biāo)準(zhǔn)，從可疑點(diǎn)集Ω中得出最終的粗差點(diǎn)集Ω′。

如圖4所示，通過步驟(1)可使得N0∩N1點(diǎn)集中的水深點(diǎn)不含粗差，這樣可以保證N1點(diǎn)集內(nèi)至少有4個不含粗差的水深點(diǎn)(m≥4)，因此在迭代過程中將不會出現(xiàn)由于點(diǎn)數(shù)不夠而導(dǎo)致擬合過程終止情況的發(fā)生。

2.3 突變地形邊界點(diǎn)判斷準(zhǔn)則

當(dāng)海底出現(xiàn)突變地形(斷崖或障礙物)且其偏差大于限差時，可能出現(xiàn)突變地形點(diǎn)群被誤認(rèn)為是成簇出現(xiàn)的粗差點(diǎn)群或者突變地形中的邊界點(diǎn)被不合理剔除等情況。對于突變地形點(diǎn)群被誤認(rèn)為是成簇出現(xiàn)的粗差點(diǎn)群這種情況，由于改進(jìn)算法采用自然鄰點(diǎn)影響域作為趨勢面擬合范圍，而自然鄰點(diǎn)影響域又是海底地形表達(dá)的最小局部范圍，其可將突出點(diǎn)群表面分割成若干個自然鄰點(diǎn)影響域范圍，通過判斷突出點(diǎn)群表面是否為局部連續(xù)的真實(shí)地形來區(qū)分成簇出現(xiàn)的粗差點(diǎn)群和突變地形點(diǎn)群，進(jìn)而完整保留突變地形點(diǎn)群上的局部連續(xù)點(diǎn)。但對于突變地形中的邊界點(diǎn)，由于其既具有突變地形點(diǎn)群局部連續(xù)的特性，又具有粗差點(diǎn)群單點(diǎn)突變的特性，采用上述改進(jìn)算法進(jìn)行濾波時，仍然將突變地形邊界點(diǎn)誤判為粗差點(diǎn)濾除，導(dǎo)致突變地形的完整性受損。

圖4 迭代趨勢面濾波過程示意圖Fig.4 The process of iterative trend surface filtering

針對這一問題，本文在上述趨勢面濾波改進(jìn)算法中加入突變地形邊界點(diǎn)判斷準(zhǔn)則，避免了突變地形邊界點(diǎn)被不合理剔除。如圖5所示，在上述迭代趨勢面濾波方法的濾波過程中，假設(shè)N0點(diǎn)集的自然鄰點(diǎn)影響域在正常的海底地形上，N1點(diǎn)集的自然鄰點(diǎn)影響域?yàn)檎５匦蔚酵蛔兊匦蔚倪^渡區(qū)域，P1點(diǎn)為突變地形邊界點(diǎn)，N2點(diǎn)集的自然鄰點(diǎn)影響域在突變海底地形上，P2點(diǎn)為突變地形上的局部連續(xù)點(diǎn)，所提突變地形邊界點(diǎn)的判斷準(zhǔn)則如下：

(2) 在N1點(diǎn)集的自然鄰點(diǎn)影響域內(nèi)進(jìn)行濾波時，將會形成一個曲率變化較大的趨勢面，以保證海底地形的連續(xù)性。

這樣在三次趨勢面濾波中對于改正前后的P1點(diǎn)是否為粗差點(diǎn)的判定便形成矛盾，而造成這一矛盾的原因是P1點(diǎn)為突變地形邊界點(diǎn)，因此根據(jù)邊界點(diǎn)在其相鄰鄰域地形內(nèi)連續(xù)性不一致這一特性，可對濾波過程中出現(xiàn)的矛盾點(diǎn)設(shè)定為突變地形邊界點(diǎn)，對其進(jìn)行保留操作。

圖5 突變地形邊界點(diǎn)判斷準(zhǔn)則示意圖Fig.5 The boundary point of mutation terrain detection

3 試驗(yàn)結(jié)果與分析

為驗(yàn)證上述改進(jìn)算法的有效性，本文通過VC++編程實(shí)現(xiàn)了基于自然鄰點(diǎn)影響域的趨勢面濾波改進(jìn)算法及傳統(tǒng)趨勢面濾波法，分別采用傳統(tǒng)趨勢面濾波法、基于自然鄰點(diǎn)影響域的趨勢面濾波改進(jìn)算法和當(dāng)前應(yīng)用較多的CUBE算法(HIPS & SIPS 6.1軟件自帶濾波算法)對包含粗差的實(shí)測數(shù)據(jù)和模擬數(shù)據(jù)進(jìn)行濾波試驗(yàn)。試驗(yàn)環(huán)境為Inter(R)core(TM)i3處理器，主頻為3.4 GHz，內(nèi)存為2 GB，最后利用Surfer8.0軟件對試驗(yàn)結(jié)果進(jìn)行了可視化顯示及對比分析。

3.1 實(shí)測算例

實(shí)測試驗(yàn)所采用的數(shù)據(jù)為我國東海某海區(qū)的多波束數(shù)據(jù)，由于邊緣波束的數(shù)據(jù)質(zhì)量較差，所以只選取靠近中央波束的數(shù)據(jù)進(jìn)行拼接，數(shù)據(jù)共包含716 516個水深點(diǎn)。在算法參數(shù)設(shè)置中，對于傳統(tǒng)趨勢面濾波法，采用完整的二次多項式，范圍半徑根據(jù)多波束測深數(shù)據(jù)密度保證所選范圍包含大約30個水深點(diǎn)，粗差判斷準(zhǔn)則中的k為2；對于其改進(jìn)算法，k同樣設(shè)定為2。兩種濾波方法耗時分別為2414 s和2643 s，CUBE算法耗時為78 s。最后利用Surfer8.0軟件中的反距離加權(quán)內(nèi)插法，對原始多波束測深數(shù)據(jù)、經(jīng)傳統(tǒng)趨勢面濾波法、改進(jìn)算法和CUBE算法濾波后的多波束測深數(shù)據(jù)進(jìn)行插值處理，形成229×200的格網(wǎng)，分別繪制各數(shù)據(jù)的海底地形表面及等深距為5 m的等深線圖，試驗(yàn)結(jié)果如圖6、圖7、圖8和圖9所示。此外，還對實(shí)測原始多波束數(shù)據(jù)(A)、經(jīng)傳統(tǒng)趨勢面濾波算法(B)、改進(jìn)算法(C)和CUBE算法(D)濾波后的實(shí)測多波束數(shù)據(jù)中的水深極大值、水深極小值、水深標(biāo)準(zhǔn)方差和數(shù)據(jù)點(diǎn)數(shù)進(jìn)行了統(tǒng)計分析，結(jié)果見表1。

表1 實(shí)測算例試驗(yàn)結(jié)果分析

由實(shí)測算例試驗(yàn)結(jié)果對比分析可知：①改進(jìn)算法和CUBE算法的濾波效果較好，而經(jīng)傳統(tǒng)趨勢面濾波法處理后的測深數(shù)據(jù)仍然包含部分粗差點(diǎn)，這些粗差點(diǎn)主要集中在復(fù)雜地區(qū)和粗差點(diǎn)成簇出現(xiàn)的區(qū)域；②經(jīng)傳統(tǒng)趨勢面濾波法、改進(jìn)算法和CUBE算法處理后的測深數(shù)據(jù)的水深值范圍均小于實(shí)測原始數(shù)據(jù)的水深值范圍，說明3種濾波方法均濾掉實(shí)測原始數(shù)據(jù)中的較大粗差點(diǎn)，而傳統(tǒng)趨勢面濾波法處理后的測深數(shù)據(jù)仍保留了一部分相對較小的粗差點(diǎn)，這主要是由于傳統(tǒng)趨勢面濾波法的抗差性較弱，較大粗差點(diǎn)嚴(yán)重影響所構(gòu)造的趨勢面，而產(chǎn)生的異常值“遮蔽”現(xiàn)象；③相比于改進(jìn)算法，傳統(tǒng)趨勢面濾波法和CUBE算法處理后的測深數(shù)據(jù)的海底地形表面相對更加平滑，對應(yīng)等深線圖中代表微小地貌的細(xì)部彎曲多數(shù)被濾除，相關(guān)統(tǒng)計數(shù)據(jù)諸如水深值范圍、標(biāo)準(zhǔn)方差和數(shù)據(jù)點(diǎn)數(shù)的數(shù)值表現(xiàn)均相對較小。

需要指出，在無法預(yù)先獲得真實(shí)海底地形的前提下，此類海底地形的相對平滑并不能代表各濾波算法的優(yōu)劣。主要原因在于平滑處理的海底地形信息低頻分量無法確定其是否為水深異常值。尤其在海底地形存在微小地貌的情況下，過度平滑可能造成海底地形表達(dá)精度的丟失。本文算法采用自然鄰點(diǎn)影響域作為趨勢面擬合范圍，且其又是海底地形表達(dá)的最小局部范圍，保證了微小地貌的準(zhǔn)確識別，傳遞式迭代濾波思想的應(yīng)用更是保證了算法的抗差性。盡管表1中的試驗(yàn)結(jié)論不足以驗(yàn)證改進(jìn)算法在保留微小地貌方面的優(yōu)勢，但數(shù)值統(tǒng)計規(guī)律表明改進(jìn)算法具有識別和保留海底地形信息低頻分量的特性。

圖6 實(shí)測原始數(shù)據(jù)的海底地形表面及等深線圖Fig.6 The submarine surface and contour map from the measured original data

圖7 傳統(tǒng)趨勢面濾波法濾波后的實(shí)測數(shù)據(jù)海底地形表面及等深線圖Fig.7 The submarine surface and contour map from the measured data after trend surface filtering

圖8 改進(jìn)算法濾波后的實(shí)測數(shù)據(jù)海底地形表面及等深線圖Fig.8 The submarine surface and contour map from the measured data after the improved filtering

圖9 CUBE算法濾波后的實(shí)測數(shù)據(jù)海底地形表面及等深線圖Fig.9 The submarine surface and contour map from the measured data after the method of CUBE

3.2 模擬算例

為進(jìn)一步驗(yàn)證改進(jìn)算法所具有的海底地形信息低頻分量保持特性，以及驗(yàn)證其在突變地形邊界點(diǎn)保留方面的優(yōu)勢，本文在上述海區(qū)人為設(shè)置粗差點(diǎn)及突變地形并進(jìn)行模擬濾波試驗(yàn)。首先，對上述海區(qū)的測深數(shù)據(jù)以2σ作為粗差判斷準(zhǔn)則進(jìn)行人工逐點(diǎn)檢查濾波處理，然后，選擇測區(qū)中的4塊小區(qū)域(長寬均為500 m)，對各小區(qū)域內(nèi)的400個水深點(diǎn)隨機(jī)加入大小為3σ～10σ不等的異常值各50個并進(jìn)行模擬濾波試驗(yàn)。如圖11所示，對4個等大的小區(qū)域(紅框區(qū)域)進(jìn)行編號，分別為I～I(xiàn)V，其中I和II區(qū)域?yàn)槠教沟匦?，III和IV區(qū)域?yàn)閺?fù)雜地形，在I和III區(qū)域內(nèi)均勻選取400個水深點(diǎn)并隨機(jī)加入異常值，以模擬平坦地形和復(fù)雜地形的擬合范圍內(nèi)單個粗差點(diǎn)出現(xiàn)的情況；在II和IV區(qū)域內(nèi)隨機(jī)選取25個水深點(diǎn)，在每個水深點(diǎn)的擬合范圍內(nèi)隨機(jī)加入大小為3σ～10σ不等的異常值各兩個，以模擬平坦地形和復(fù)雜地形的擬合范圍內(nèi)多個粗差點(diǎn)成簇出現(xiàn)的情況。此外，為了驗(yàn)證突變地形邊界點(diǎn)判斷準(zhǔn)則的可行性，在平坦區(qū)域設(shè)置了一個凸起的長方體(藍(lán)框區(qū)域，長為300 m，寬為200 m，水深值為-35 m)作為突變地形，編號為V。最后利用Surfer8.0軟件繪制人工處理后的原始數(shù)據(jù)(A)、加粗差后的模擬數(shù)據(jù)(B)、經(jīng)傳統(tǒng)趨勢面濾波法(C)、改進(jìn)算法(D)和CUBE算法(E)濾波后的模擬多波束測深數(shù)據(jù)的海底地形表面及等深距為5 m的等深線圖，試驗(yàn)結(jié)果如圖10—圖14所示。此外，還對各小區(qū)域內(nèi)的水深極大值、水深極小值、水深標(biāo)準(zhǔn)方差、數(shù)據(jù)點(diǎn)數(shù)和各方法粗差檢測正確率進(jìn)行了統(tǒng)計，結(jié)果見表2和表3所示。

由模擬算例試驗(yàn)結(jié)果對比分析可進(jìn)一步看出：

(1) 對于平坦海底處單個出現(xiàn)的粗差點(diǎn)(I區(qū)域)，由于海底地形較為簡單，各濾波方法所構(gòu)造的趨勢面與真實(shí)海底地形之間的偏差較小，3種濾波方法對不同程度的粗差點(diǎn)均具有較好的識別能力。

(2) 對于平坦海底處多個成簇出現(xiàn)的粗差點(diǎn)(II區(qū)域)，由于傳統(tǒng)趨勢面濾波法未考慮粗差數(shù)據(jù)對趨勢面構(gòu)造的影響，導(dǎo)致濾波后的多波束測深數(shù)據(jù)中仍有小部分粗差值較小的粗差點(diǎn)未被剔除，而改進(jìn)算法和CUBE算法在算法設(shè)計上增加了對測深數(shù)據(jù)的抗差性，減弱了粗差數(shù)據(jù)對趨勢面構(gòu)造的影響，濾波效果較好。

表2 模擬算例實(shí)驗(yàn)結(jié)果分析

表3 I～I(xiàn)V區(qū)域粗差檢測正確率

(3) 對于復(fù)雜海底處出現(xiàn)的單個和多個成簇出現(xiàn)的粗差點(diǎn)(III和IV區(qū)域)，由于傳統(tǒng)趨勢面濾波法在復(fù)雜地區(qū)局部范圍內(nèi)所構(gòu)造的趨勢面不能合理反映真實(shí)海底地形，導(dǎo)致多波束測深數(shù)據(jù)中的粗差點(diǎn)剔除不徹底，尤其在復(fù)雜海底處多個成簇出現(xiàn)的粗差點(diǎn)(IV區(qū)域)濾除情況最差；而

圖10 人工處理后的測深數(shù)據(jù)海底地形表面及等深線圖Fig.10 The submarine surface and contour map from the measured original data after manual handling

圖11 模擬原始數(shù)據(jù)的海底地形表面及等深線圖Fig.11 The submarine surface and contour map from the simulative original data

圖12 傳統(tǒng)趨勢面濾波法濾波后的模擬數(shù)據(jù)海底地形表面及等深線圖Fig.12 The submarine surface and contour map from the simulative data after trend surface filtering

圖13 改進(jìn)算法濾波后的模擬數(shù)據(jù)海底地形表面及等深線圖Fig.13 The submarine surface and contour map from the simulative data after the improved filtering

CUBE算法運(yùn)用水深信息局部相關(guān)性原則為每一個水深節(jié)點(diǎn)確定水深估值，所構(gòu)造的CUBE曲面為規(guī)則網(wǎng)格曲面，其與真實(shí)海底地形在微小地貌處仍存在一定的差異，特別在復(fù)雜海底處多個粗差點(diǎn)成簇出現(xiàn)的區(qū)域(IV區(qū)域)，可能出現(xiàn)少量粗差值較小的粗差點(diǎn)被保留；本文所提改進(jìn)算法采用自然鄰點(diǎn)影響域作為其擬合范圍，所構(gòu)造的趨勢面可以反映各種復(fù)雜程度的真實(shí)海底地形，特別對于復(fù)雜海底處多個成簇出現(xiàn)的粗差點(diǎn)，改進(jìn)算法濾波效果最好。

圖14 CUBE算法濾波后的模擬數(shù)據(jù)海底地形表面及等深線圖Fig.14 The submarine surface and contour map from the simulative data after the method of CUBE

(4) 在I～I(xiàn)V 4個區(qū)域中，經(jīng)改進(jìn)算法濾波后的測深數(shù)據(jù)的水深值范圍、標(biāo)準(zhǔn)方差和數(shù)據(jù)點(diǎn)數(shù)均大于傳統(tǒng)趨勢面濾波法和CUBE算法濾波后的測深數(shù)據(jù)，且與人工處理后的測深數(shù)據(jù)的各水深值統(tǒng)計數(shù)據(jù)較為接近，這說明經(jīng)改進(jìn)算法濾波后的測深數(shù)據(jù)的水深值波動相對于海底地形的變化較小，相比于傳統(tǒng)趨勢面濾波法和CUBE算法，保留了較多的正常水深點(diǎn)。造成上述差異的主要原因在于傳統(tǒng)趨勢面濾波法在復(fù)雜地區(qū)局部范圍內(nèi)所構(gòu)造的趨勢面不能合理反映真實(shí)海底地形且未考慮粗差數(shù)據(jù)對趨勢面構(gòu)造的影響，導(dǎo)致大量正常水深點(diǎn)誤判為粗差點(diǎn)而被濾除；而CUBE算法是在設(shè)定的格網(wǎng)節(jié)點(diǎn)下選擇水深，并對測區(qū)格網(wǎng)節(jié)點(diǎn)處進(jìn)行水深及其相關(guān)誤差估計，最終擬合的CUBE曲面與實(shí)際海底地形曲面存在一定的差異，將會導(dǎo)致一部分正常點(diǎn)被濾除。

(5) 在V區(qū)域中，由各算法濾波后的測深數(shù)據(jù)點(diǎn)數(shù)可以看出，改進(jìn)算法對突變地形的保留效果最好，CUBE算法次之，傳統(tǒng)趨勢面濾波法最差。對比圖12和圖14，CUBE算法和傳統(tǒng)趨勢面濾波法濾除的突變地形點(diǎn)主要為突變地形邊界點(diǎn)和突變地形邊緣以內(nèi)的部分局部連續(xù)點(diǎn)。造成這一現(xiàn)象的主要原因在于傳統(tǒng)趨勢面濾波法和CUBE算法是在假設(shè)海底地形連續(xù)的前提下進(jìn)行的，均無法有效識別突變地形點(diǎn)群。傳統(tǒng)趨勢面濾波法將突變地形點(diǎn)群視為海底連續(xù)地形的一部分，濾除掉與所構(gòu)造趨勢面之間偏差較大的水深點(diǎn)，導(dǎo)致出現(xiàn)突變地形邊界點(diǎn)和突變地形邊緣以內(nèi)的部分局部連續(xù)點(diǎn)被剔除的結(jié)果；CUBE算法運(yùn)用水深信息局部相關(guān)性原則為每一個水深節(jié)點(diǎn)確定水深估值，對于突變地形邊緣以內(nèi)的局部連續(xù)點(diǎn)，由于其局部連續(xù)性，得到最大限度的保留，但突變地形邊界點(diǎn)具有粗差單點(diǎn)突變的特性，由其鄰近水深節(jié)點(diǎn)構(gòu)成的CUBE曲面與突變地形邊界點(diǎn)仍存在較大差異，最終導(dǎo)致出現(xiàn)突變地形邊界點(diǎn)和少量的突變地形局部連續(xù)點(diǎn)被剔除的結(jié)果；而改進(jìn)算法采用自然鄰點(diǎn)影響域作為擬合范圍，由于其是海底地形表達(dá)的最小局部范圍，較好地保留了突變地形上的局部連續(xù)點(diǎn)，通過引入突變地形邊界點(diǎn)判斷準(zhǔn)則，有效識別了突變地形的邊界點(diǎn)，較好地刻畫了海底突變地形的實(shí)際形狀。

4 結(jié) 論

通過理論分析及試驗(yàn)比對得結(jié)論如下：

(1) 改進(jìn)算法和CUBE算法的粗差點(diǎn)識別能力優(yōu)于傳統(tǒng)趨勢面濾波法，特別是在海底地形復(fù)雜區(qū)域或粗差點(diǎn)成簇連續(xù)出現(xiàn)的區(qū)域，改進(jìn)算法可有效識別并剔除粗差點(diǎn)。

(2) 傳統(tǒng)趨勢面濾波法和CUBE算法在濾波過程中均有較多正常點(diǎn)被剔除掉，而改進(jìn)算法則對這些正常水深點(diǎn)進(jìn)行了保留，提高了海底地形表達(dá)的精度。

(3) 傳統(tǒng)趨勢面濾波法和CUBE算法由于在假設(shè)海底地形連續(xù)的前提下進(jìn)行的，均無法有效識別突變地形點(diǎn)群，導(dǎo)致突變地形部分局部連續(xù)點(diǎn)和全部邊界點(diǎn)被不合理剔除，而改進(jìn)算法采用自然鄰點(diǎn)影響域作為擬合范圍和引入突變地形邊界點(diǎn)判斷準(zhǔn)則，有效識別了突變地形點(diǎn)群，保證了海底突變地形的完整性。

本文對傳統(tǒng)趨勢面濾波法的改進(jìn)主要集中于海底實(shí)際趨勢面的構(gòu)建，而對粗差判斷準(zhǔn)則未進(jìn)行討論，仍然采用2倍(或3倍)中誤差判斷準(zhǔn)則作為改進(jìn)算法的粗差判斷準(zhǔn)則，這可能會導(dǎo)致部分正常水深點(diǎn)被不合理刪除和一些較小粗差被保留下來，下一步可針對粗差判斷準(zhǔn)則進(jìn)行研究及改進(jìn)。此外，改進(jìn)算法與CUBE算法在效率方面相差較大，下一步可在自然鄰點(diǎn)快速搜索、數(shù)據(jù)分塊處理和算法并行運(yùn)算等方面進(jìn)行研究。

[1] 劉雁春, 肖付民, 暴景陽, 等. 海道測量學(xué)概論[M]. 北京: 測繪出版社, 2006: 132-138. LIU Yanchun, XIAO Fumin, BAO Jingyang, et al. Introduction to Hydrogrophy[M]. Beijing: Surveying and Mapping Press, 2006: 132-138.

[2] 黃謨濤, 翟國君, 王瑞, 等. 海洋測量異常數(shù)據(jù)的檢測[J]. 測繪學(xué)報, 1999, 28(3): 269-276. DOI: 10.3321/j.issn:1001-1595.1999.03.015. HUANG Motao, ZHAI Guojun, WANG Rui, et al. The Detection of Abnormal Data in Marine Survey[J]. Acta Geodaetica et Cartographica Sinica, 1999, 28(3): 269-277. DOI: 10.3321/j.issn:1001-1595.1999.03.015.

[3] 陽凡林, 劉經(jīng)南, 趙建虎. 多波束測深數(shù)據(jù)的異常檢測和濾波[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2004, 29(1): 80-83. YANG Fanlin, LIU Jingnan, ZHAO Jianhu. Detecting Outliers and Filtering Noises in Multi-beam Data[J]. Geomatics and Information Science of Wuhan University, 2004, 29(1): 80-83.

[4] 趙建虎, 劉經(jīng)南. 多波束測深及圖像數(shù)據(jù)處理[M]. 武漢: 武漢大學(xué)出版社, 2008: 205-210. ZHAO Jianhu, LIU Jingnan. Multi-beam Sounding and Image Data Processing[M]. Wuhan: Wuhan University Press, 2008: 205-210.

[5] WARE C, KNIGHT W, WELLS D. Memory Intensive Statistical Algorithms for Multibeam Bathymetric Data[J]. Computers & Geosciences, 1991, 17(7): 985-993.

[6] 朱慶, 李德仁. 多波束測深數(shù)據(jù)的誤差分析與處理[J]. 武漢測繪科技大學(xué)學(xué)報, 1998, 23(1): 1-4, 46. ZHU Qing, LI Deren. Error Analysis and Processing of Multibeam Soundings[J]. Journal of Wuhan Technical University of Surveying and Mapping, 1998, 23(1): 1-4, 46.

[7] CALDER B R, MAYER L A. Automatic Processing of High-rate, High-density Multibeam Echosounder Data[J]. Geochemistry Geophysics Geosystems, 2003, 4(6): 1048.

[9] MALLACE D, ROBERTSON P. Alternative Use of CUBE: How to Fit a Square Peg in a Round Hole[C]∥Proceedings of US Hydrographic Conference. Norfolk, Virginia:[s.m.]，2007.

[10] LCDR Aluizio Maciel de Oliveira Junior, CDR Izabel King Jeck. Multibeam Processing for Nautical Charts (Using CUBE and “Surface Filter” to Enhance Multibeam Processing)[J]. International Hydrographic Review, 2009(2): 61-71.

[11] 王德剛, 葉燦銀. CUBE算法及其在多波束數(shù)據(jù)處理中的應(yīng)用[J]. 海洋學(xué)研究, 2008, 26(2): 82-87. WAND Degang, YE Canyin. The Theory of CUBE Algorithm and Its Application in the Processing of Multi-beam Data[J]. Journal of Marine Sciences, 2008, 26(2): 82-87.

[12] 黃謨濤, 翟國君, 柴洪洲, 等. 檢測多波束測深異常數(shù)據(jù)的CUBE算法模型解析[J]. 海洋測繪, 2011, 31(4): 1-4. HUANG Motao, ZHAI Guojun, CHAI Hongzhou, et al. Analysis on the Mathematical Models of CUBE Algorithm for the Detection of Abnormal Data in Multibeam Echosounding[J]. Hydrographic Surveying and Charting, 2011, 31(4): 1-4.

[13] 李新娜, 歸慶明, 許阿裴. 基于識別變量的粗差探測Bayes方法[J]. 測繪學(xué)報, 2008, 37(3): 355-360, 366. DOI: 10.3321/j.issn:1001-1595.2008.03.015. LI Xinna, GUI Qingming, XU Apei. Bayesian Method for Detection of Gross Errors Based on Classification Variables[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(3): 355-360, 366. DOI: 10.3321/j.issn:1001-1595.2008.03.015.

[14] 黃賢源, 隋立芬, 翟國君, 等. 利用Bayes估計進(jìn)行多波束測深異常數(shù)據(jù)探測[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2010, 35(2): 168-171. HUANG Xianyuan, SUI Lifen, ZHAI Guojun, et al. Outliers Detection of Multi-beam Data Based on Bayes Estimation[J]. Geomatics and Information Science of Wuhan University, 2010, 35(2): 168-171.

[15] 黃賢源, 翟國君, 黃謨濤, 等. 利用水深不確定度探測測深異常值的方法[J]. 測繪科學(xué)技術(shù)學(xué)報, 2011, 28(1): 70-74, 78. HUANG Xianyuan, ZHAI Guojun, HUANG Motao, et al. The Approach on Detecting Outliers of Multi-beam Data by Uncertainty[J]. Journal of Geomatics Science and Technology, 2011, 28(1): 70-74, 78.

[16] 董江, 任立生. 基于趨勢面的多波束測深數(shù)據(jù)濾波方法[J]. 海洋測繪, 2007, 27(6): 25-28. DONG Jiang, REN Lisheng. Filter of MBS Sounding Data Based on Trend Surface[J]. Hydrographic Surveying and Charting, 2007, 27(6): 25-28.

[17] 王海棟, 柴洪洲, 宋國大, 等. 多波束測深趨勢面系數(shù)的主成分估計[J]. 海洋測繪, 2009, 29(5): 5-7. WANG Haidong, CHAI Hongzhou, SONG Guoda, et al. Principal Components Estimation of Trend Surface Coefficients in Multibeam Bathymetry[J]. Hydrographic Surveying and Charting, 2009, 29(5): 5-7.

[18] 王海棟, 柴洪洲, 翟天增, 等. 多波束測深異常的兩種趨勢面檢測算法比較[J]. 海洋通報, 2010, 29(2): 182-186. WANG Haidong, CHAI Hongzhou, ZHAI Tianzeng, et al. Comparison of Two Trend Surface Detection Algorithms of Multibeam Bathymetry Outlier[J]. Marine Science Bulletin, 2010, 29(2): 182-186.

[19] 李志林, 朱慶. 數(shù)字高程模型[M]. 2版. 武漢: 武漢大學(xué)出版社, 2008: 17-19, 112-133. LI Zhilin, ZHU Qing. Digital Elevation Model[M]. 2nd ed. Wuhan: Wuhan University Press, 2008: 17-19, 112-133.

[20] YANG Xunnian, WANG Guozhao. Planar Point Set Fairing and Fitting by Arc Splines[J]. Computer-Aided Design, 2001, 33(1): 35-43.

[21] 陳曦. 反求工程中基于點(diǎn)云的特征挖掘技術(shù)研究[D]. 杭州: 浙江大學(xué), 2005: 38-39. CHEN Xi. Study on Feature Mining Technology Based on Point Cloud in Reverse Engineering[D]. Hangzhou: Zhejiang University, 2005: 38-39.

[22] 湯國安, 劉學(xué)軍, 閭國年. 數(shù)字高程模型及地學(xué)分析的原理與方法[M]. 北京: 科學(xué)出版社, 2005: 63-65. TANG Guoan, LIU Xuejun, Lü Guonian. Principle and Method of Digital Elevation Model and Geological Analysis[M]. Beijing: Science Press, 2005: 63-65.

[23] 王幼寧, 劉繼志. 微分幾何講義[M]. 北京: 北京師范大學(xué)出版社, 2011: 116-117. WANG Youning, LIU Jizhi. Differential Geometry Notes[M]. Beijing: Beijing Normal University Press, 2011: 116-117.

[24] 陳軍, 趙仁亮, 喬朝飛. 基于Voronoi圖的GIS空間分析研究[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2003, 28(S1): 32-37. CHEN Jun, ZHAO Renliang, QIAO Chaofei. Voronoi Diagram-based GIS Spatial Analysis[J]. Geomatics and Information Science of Wuhan University, 2003, 28(S1): 32-37.

[25] 董箭, 彭認(rèn)燦, 鄭義東, 等. 局部動態(tài)最優(yōu)Voronoi圖的NNI算法及其在格網(wǎng)數(shù)字水深模型中的應(yīng)用[J]. 測繪學(xué)報, 2013, 42(2): 284-289, 303. DONG Jian, PENG Rencan, ZHENG Yidong, et al. An algorithm of Natural Neighbor Interpolation Based on local Dynamic Optimal Voronoi Diagram and Its Application in Grid Digital Depth Model[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(2): 284-289, 303.

[26] 郭艷軍, 潘懋, 燕飛, 等. 自然鄰點(diǎn)插值方法在三維地質(zhì)建模中的應(yīng)用[J]. 解放軍理工大學(xué)學(xué)報(自然科學(xué)版), 2009, 10(6): 650-655. GUO Yanjun, PAN Mao, YAN Fei, et al. Application of Natural Neighbor Interpolation Method in Three-Dimensional Geological Modeling[J]. Journal of PLA University of Science and Technology (Natural Science Edition), 2009, 10(6): 650-655.

[27] 吳英姿. 多波束測深系統(tǒng)地形跟蹤與數(shù)據(jù)處理技術(shù)研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2001: 59-79. WU Yingzi. A Study on Multi-Beam Sounding System Seafloor Tracking & Data Processing Techniques[D]. Harbin: Harbin Engineering University, 2001: 59-79.

[28] 蔡永昌, 朱合華. 基于局部搜索算法的自然鄰接點(diǎn)方法[J]. 力學(xué)學(xué)報, 2004, 36(5): 623-628. CAI Yongchang, ZHU Hehua. Natural Neighbour Method Based on the Algorithm of Local Search[J]. Acta Mechanica Sinica, 2004, 36(5): 623-628.

[29] 董箭, 彭認(rèn)燦, 鄭義東. 利用局部動態(tài)最優(yōu)Delaunay三角網(wǎng)改進(jìn)逐點(diǎn)內(nèi)插算法[J]. 武漢大學(xué)學(xué)報(信息科學(xué)版), 2013, 38(5): 613-617. DONG Jian, PENG Rencan, ZHENG Yidong. An Improved Algorithm of Point-by-point Interpolation by Using Local Dynamic Optimal Delaunay Triangulation Network[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 613-617.

[30] 呂志平, 魏子卿, 李軍, 等. 我國CGCS2000高精度坐標(biāo)轉(zhuǎn)換格網(wǎng)模型的建立[J]. 測繪學(xué)報, 2013, 42(6): 791-797. Lü Zhiping, WEI Ziqing, LI Jun, et al. The Establishment of High Precise Coordinate Transformation Grid Model of CGCS2000[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(6): 791-797.