考慮自然鄰點影響域的多波束測深數(shù)據(jù)趨勢面濾波改進算法

張志衡，彭認燦，黃文騫,3，董 箭，劉國輝

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

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

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

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

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

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

(1)

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

在假設局部海底地形復雜度適中且變化不大的前提下，局部曲面擬合范圍Q可以待檢測點q(xq,yq,zq)為中心且包含20～30個水深點的區(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)的平面坐標(xi,yi)和水深值zi擬合出趨勢面后，趨勢面濾波的異常值判斷準則為[4]

(4)

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

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

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

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

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

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

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

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

2 趨勢面濾波改進算法

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

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

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

(5)

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

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

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

(6)

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

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

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

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

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

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

(7)

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

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

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

(8)

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

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

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

B=TA+Δ

(9)

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

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

(10)

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

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

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

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

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

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

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

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

2.3 突變地形邊界點判斷準則

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

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

針對這一問題，本文在上述趨勢面濾波改進算法中加入突變地形邊界點判斷準則，避免了突變地形邊界點被不合理剔除。如圖5所示，在上述迭代趨勢面濾波方法的濾波過程中，假設N0點集的自然鄰點影響域在正常的海底地形上，N1點集的自然鄰點影響域為正常地形到突變地形的過渡區(qū)域，P1點為突變地形邊界點，N2點集的自然鄰點影響域在突變海底地形上，P2點為突變地形上的局部連續(xù)點，所提突變地形邊界點的判斷準則如下：

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

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

圖5 突變地形邊界點判斷準則示意圖Fig.5 The boundary point of mutation terrain detection

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

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

3.1 實測算例

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

表1 實測算例試驗結(jié)果分析

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

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

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

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

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

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

3.2 模擬算例

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

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

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

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

表2 模擬算例實驗結(jié)果分析

表3 I～IV區(qū)域粗差檢測正確率

(3) 對于復雜海底處出現(xiàn)的單個和多個成簇出現(xiàn)的粗差點(III和IV區(qū)域)，由于傳統(tǒng)趨勢面濾波法在復雜地區(qū)局部范圍內(nèi)所構(gòu)造的趨勢面不能合理反映真實海底地形，導致多波束測深數(shù)據(jù)中的粗差點剔除不徹底，尤其在復雜海底處多個成簇出現(xià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 改進算法濾波后的模擬數(shù)據(jù)海底地形表面及等深線圖Fig.13 The submarine surface and contour map from the simulative data after the improved filtering

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

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

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

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

4 結(jié) 論

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

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

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

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

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

[1] 劉雁春, 肖付民, 暴景陽, 等. 海道測量學概論[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]. 測繪學報, 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]. 武漢大學學報(信息科學版), 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]. 武漢: 武漢大學出版社, 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]. 武漢測繪科技大學學報, 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ù)處理中的應用[J]. 海洋學研究, 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]. 測繪學報, 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估計進行多波束測深異常數(shù)據(jù)探測[J]. 武漢大學學報(信息科學版), 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]. 測繪科學技術學報, 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版. 武漢: 武漢大學出版社, 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] 陳曦. 反求工程中基于點云的特征挖掘技術研究[D]. 杭州: 浙江大學, 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] 湯國安, 劉學軍, 閭國年. 數(shù)字高程模型及地學分析的原理與方法[M]. 北京: 科學出版社, 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]. 北京: 北京師范大學出版社, 2011: 116-117. WANG Youning, LIU Jizhi. Differential Geometry Notes[M]. Beijing: Beijing Normal University Press, 2011: 116-117.

[24] 陳軍, 趙仁亮, 喬朝飛. 基于Voronoi圖的GIS空間分析研究[J]. 武漢大學學報(信息科學版), 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] 董箭, 彭認燦, 鄭義東, 等. 局部動態(tài)最優(yōu)Voronoi圖的NNI算法及其在格網(wǎng)數(shù)字水深模型中的應用[J]. 測繪學報, 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] 郭艷軍, 潘懋, 燕飛, 等. 自然鄰點插值方法在三維地質(zhì)建模中的應用[J]. 解放軍理工大學學報(自然科學版), 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ù)處理技術研究[D]. 哈爾濱: 哈爾濱工程大學, 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] 蔡永昌, 朱合華. 基于局部搜索算法的自然鄰接點方法[J]. 力學學報, 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] 董箭, 彭認燦, 鄭義東. 利用局部動態(tài)最優(yōu)Delaunay三角網(wǎng)改進逐點內(nèi)插算法[J]. 武漢大學學報(信息科學版), 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高精度坐標轉(zhuǎn)換格網(wǎng)模型的建立[J]. 測繪學報, 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.