羅進華，楊修偉，朱凌霄，朱培民

（1.中海油田服務(wù)股份有限公司，天津300451；2.中國地質(zhì)大學(xué)（武漢）地球物理與空間信息學(xué)院，湖北 武漢430074）

側(cè)掃聲吶作為探測海底地貌反射特征的一種有效手段，具有分辨率高、覆蓋面積大、采集效率高、對海底地貌直觀成像等優(yōu)勢，其被越來越多的應(yīng)用于海洋測繪［1－2］、海洋工程［3－4］、海洋地質(zhì)勘探［5－6］、海洋生態(tài)與環(huán)境保護［7－8］和海底目標(biāo)探測與定位［9］等多個方面的研究。側(cè)掃聲吶原始圖像為時間序列記錄，且由于聲吶拖曳速度、姿態(tài)等多種因素影響，使得海底目標(biāo)的成像發(fā)生移位和變形，進而使相鄰條帶圖像不能完美鑲嵌拼接。為還原采樣點的真實平面位置，恢復(fù)海底地貌真實形態(tài)，需要對側(cè)掃聲吶圖像進行幾何校正［10－11］，以最大限度真實地對海底目標(biāo)成像。經(jīng)過幾何校正后的聲吶采樣數(shù)據(jù)，需要通過重采樣獲得規(guī)則的網(wǎng)格數(shù)據(jù)，以生成全覆蓋的海底地貌圖。由于側(cè)掃聲吶數(shù)據(jù)量通常較大，重采樣時選擇計算效率高的插值算法十分關(guān)鍵。

插值技術(shù)已普遍應(yīng)用于醫(yī)學(xué)、動畫、數(shù)字高程模型［12］、遙感［13］和氣象［14－15］等領(lǐng)域，但專門針對側(cè)掃聲吶數(shù)據(jù)的特點而展開的插值算法研究非常少。目前普遍采取直接填充數(shù)據(jù)空隙的方式［11］，但其存在較大弊端。由于經(jīng)過幾何校正后的采樣點存在單條聲吶掃描線（ping）呈直線分布、ping與ping之間旋轉(zhuǎn)角度不一且距離不等的特征，本文針對側(cè)掃聲吶數(shù)據(jù)的這些特征，采取了一種特殊的方法，即基于扭曲網(wǎng)格三角剖分有限元插值技術(shù)的重采樣方法，對側(cè)掃聲吶圖像進行重新構(gòu)建。在選擇合適重采樣間隔的基礎(chǔ)上，該方法比其他常用的插值方法計算效率高，而且能夠很好地保留原始數(shù)據(jù)的信息。

1 重采樣間隔的選擇

側(cè)掃聲吶數(shù)據(jù)重采樣的精度與原始數(shù)據(jù)的分辨率、數(shù)據(jù)分布密度、插值采樣間隔等有關(guān)［16－17］，因此需要對側(cè)掃聲吶數(shù)據(jù)進行深入分析，以確定合適的重采樣間隔。

1.1 側(cè)掃聲吶圖像分辨率

側(cè)掃聲吶數(shù)據(jù)分辨率分為橫向（垂直航跡）分辨率和縱向（沿航跡）分辨率。

1.1.1 橫向（垂直航跡）分辨率

側(cè)掃聲吶橫向（垂直航跡）分辨率（Δx）取決于脈寬的大?。▓D1）。理論上講，兩個目標(biāo)的最小距離應(yīng)至少等于脈寬的一半［18］，才能分辨出來。但由于波束是傾斜狀態(tài)，因此越靠近聲吶正下方，“足跡”越長，分辨率與傾斜角（β）的余弦有關(guān)（圖2），實際分辨率總小于理論值，且從中心向兩側(cè)逐漸增高。橫向分辨率Δx計算公式如下：

式中，c為聲速；τ為脈沖持續(xù)時間；h為聲吶距海底的高度；R為斜距（slant range）。

圖1 側(cè)掃聲吶水中傳播的脈沖Fig.1 Acoustic pulse of sidescan sonar in water

圖2 垂直航跡方向分辨率Fig.2 Across－track resolution

1.1.2 縱向（沿航跡）分辨率

側(cè)掃聲吶縱向（沿航跡）分辨率（Δy）取決于聲吶的水平波束寬度（圖3）。波束寬度隨與聲吶的距離增加而逐漸變寬，因此，縱向分辨率從中心向兩側(cè)逐漸降低?？v向分辨率Δy計算公式如下：

式中，R為斜距；θh為水平波束角。

圖3 沿航跡方向分辨率Fig.3 Along－track resolution

1.2 數(shù)據(jù)分布特征

幾何校正后的采樣點數(shù)據(jù)分布是不規(guī)則且疏密不均的，但數(shù)據(jù)分布又不是完全雜亂無章的，縱向和橫向分布特征均可根據(jù)聲吶高度等信息對變化規(guī)律進行估計。

1.2.1 橫向（垂直航跡）分布特征

聲吶數(shù)據(jù)橫向變化是最大斜距（SRmax）以及高度（h）的函數(shù)，采樣點間距（dx）滿足以下表達式：

式中，xi為第i個采樣點與采集中心的距離。

在測量過程中，每ping的最大斜距（SRmax），采樣點數(shù)（nx），聲吶高度（h）均為定值，可以求得橫向采樣間距從中心向兩側(cè)逐漸減小，即形成中間采樣點稀疏，向兩側(cè)逐漸密集的分布規(guī)律。

1.2.2 縱向（沿航跡）分布特征

縱向采樣間隔（dy）表示為聲吶行進速度（v）與采樣時間間隔（Δt）的乘積，如公式（6），因此最小速度和最大速度分別決定了縱向采樣間隔的最小值和最大值。

1.2.3 姿態(tài)對分布特征的影響

當(dāng)今水下姿態(tài)傳感器精度越來越高，特別是在深拖、水下機器人（Remote Operated Vehicle，ROV）和自主式水下潛器（Autonomous Underwater Vehicle，AUV）等水下載體搭載有高精度姿態(tài)儀測量的情況下，需要考慮聲吶姿態(tài)對采樣點位置的搖影響。聲吶儀縱搖（Pitch）變化會引起采樣點間距的變化，聲吶儀艏向（Heading）變化會造成單ping呈直線分布、ping與ping之間存在旋轉(zhuǎn)角度的現(xiàn)象。聲吶儀橫搖（Roll）會造成ping沿與聲吶儀垂直行進方向的偏移。

1.3 重采樣間隔的選擇原則

重采樣間隔的選擇對圖像質(zhì)量有直接影響，大的采樣間隔會使圖像信息丟失，分辨率下降，而過小的采樣間隔則會大大增加計算量，不僅對圖像質(zhì)量無大的影響，而且往往致使處理后的聲吶文件非常大，應(yīng)用受到限制。

當(dāng)分辨率小于采樣間隔時，兩個采樣點之間存在未測量區(qū)域，需要通過插值獲得屬性值。此時，選擇圖像的最高分辨率作為重采樣間隔，可以在保證圖像分辨率的基礎(chǔ)上，獲得未測區(qū)域的屬性值。

當(dāng)分辨率大于采樣間隔時，由于不同的散射值平均，產(chǎn)生平滑效應(yīng)而丟失了小尺寸細(xì)節(jié)信息，此時選擇最小采樣間隔作為重采樣間隔，可以最大限度的保持原始圖像的信息量，減少原始數(shù)據(jù)信息的損失。

由于艏向的變化會引起每ping數(shù)據(jù)以拖魚所在處為原點的整體旋轉(zhuǎn)，與重采樣坐標(biāo)軸形成一定的夾角。選擇重采樣間隔時，應(yīng)將垂直航跡和沿航跡方向的重采樣間隔投影到相應(yīng)坐標(biāo)軸上，獲得沿坐標(biāo)軸方向的重采樣間隔。

綜上所述，重采樣間隔應(yīng)在綜合考慮圖像分辨率及原始采樣點分布特征的基礎(chǔ)上，選擇最高分辨率和最小采樣間隔中的較小值，并投影到相應(yīng)的重采樣坐標(biāo)軸上，獲取重采樣間隔。

2 插值方法的選擇

2.1 常規(guī)插值算法的適用性

常用的空間插值算法包括反距離加權(quán)法、徑向基函數(shù)法、克里格法、三角剖分線性插值法等［19］。根據(jù)插值過程中確定對插值點有影響的點的方式，可將以上方法大致分為2種。一種是通過搜索半徑確定對插值點有影響的采樣點，計算每個點的權(quán)重，求得插值點屬性。由于側(cè)掃聲吶圖像采樣點分布不規(guī)則，固定的搜索半徑?jīng)]有考慮采樣點分布密度變化的影響，小搜索半徑使采樣點稀疏的區(qū)域出現(xiàn)空白，大搜索半徑常常使采樣點密集的區(qū)域產(chǎn)生過度的平滑效應(yīng)。另一種是通過將離散的節(jié)點剖分為互不重疊單元，在每個單元內(nèi)建立插值函數(shù)，進行插值，如普通三角剖分線性插值法，該方法插值精度高，能很好地適應(yīng)不同采樣點疏密的變化。但由于側(cè)掃聲吶圖像數(shù)據(jù)量大，使用普通的三角剖分法計算效率較低，難以滿足實際生產(chǎn)的需求。

2.2 扭曲網(wǎng)格三角剖分有限元插值方法

普通的三角剖分針對任意散亂點，通過計算形成滿足一定條件的三角網(wǎng)，計算量大。側(cè)掃聲吶圖像幾何校正后采樣點分布雖然不規(guī)則，但并非完全雜亂無章，校正后的數(shù)據(jù)網(wǎng)格可以認(rèn)為是原數(shù)據(jù)規(guī)則網(wǎng)格的扭曲變形（圖4）。本文提出的基于扭曲網(wǎng)格的三角剖分法，利用了網(wǎng)格扭曲變形后的拓?fù)潢P(guān)系，直接將四邊形分成三角形，計算量大大減小。使用該方法對幾何校正后的側(cè)掃聲吶采樣數(shù)據(jù)進行三角剖分，可以快速生成三角網(wǎng)，節(jié)省了普通的離散點三角剖分所需的時間。

圖4 幾何校正前后采樣點分布變化示意圖Fig.4 Distribution of sampling points before and after post geometric correction

2.3 重疊區(qū)域的處理方法

側(cè)掃聲吶數(shù)據(jù)采集過程中由于姿態(tài)變化的影響，常常造成采樣點區(qū)域重疊的情況。對于這種情況一般有兩種處理方式：一種把重疊區(qū)域當(dāng)成同一區(qū)域的多次測量，因此對位于多個單元內(nèi)的插值點通過平均的方式進行求取；另一種是采取覆蓋的方式。本文采取了前一種方法。

3 扭曲網(wǎng)格三角剖分有限元插值的基本原理

3.1 扭曲網(wǎng)格三角剖分

三角剖分具有最大化最小角特性，是最接近于規(guī)則化的三角網(wǎng)［20］。對扭曲的四邊形網(wǎng)格單元進行三角剖分，同樣可以遵循該原則。任意凸四邊形單元都有2種剖分為三角形單元的方式（圖5a和圖5b），剖分時選擇使凸四邊形剖分后的6個內(nèi)角中最小角較大的剖分方式（圖5b）；凹四邊形則對凹角頂點與相對頂點連接，剖分為2個三角形（圖5c）。

圖5 四邊形三角剖分Fig.5 Triangulation of quadrilateral

3.2 面積坐標(biāo)

三角形中任意一點P與其3個角點相連形成3個子三角形，面積分別為Ai，Aj，Ak，大三角形面積為A（圖6）。P點的位置可以由3個比值來確定，即P（Li，Lj，Lk）。其中，Li＝Ai／A，Lj＝Aj／A，Lk＝Ak／A。

圖6 三角形單元的面積坐標(biāo)Fig.6 Area coordinates of triangular element

3.3 插值函數(shù)

三角形有限元插值函數(shù)的基本形式為

式中，d1，d2，d3為3個節(jié)點屬性值；d為插值點的屬性；Ni為形函數(shù)，且滿足。形函數(shù)Ni可以用三角形的面積坐標(biāo)表示［21］：Ni＝Li（i＝1，2，3）。

4 實驗設(shè)計

本文以美國Triton公司公布給Triton Isis用戶的示例側(cè)掃聲吶數(shù)據(jù)（圖7）為實驗對象，對扭曲網(wǎng)格三角剖分有限元插值算法在側(cè)掃聲吶圖像幾何校正重采樣中的應(yīng)用效果進行測試。該聲吶信號頻率100kHz，脈沖寬度0.1ms，脈沖間隔0.12s，水中聲速使用1 500m／s，水平波束角為1°，最大斜距為90m，橫向（垂直航跡）采樣點為4 096個，縱向（沿航跡）采樣點為2 009個。

圖7 原始聲吶數(shù)據(jù)Fig.7 Raw sonar data

4.1 側(cè)掃聲吶數(shù)據(jù)插值前處理

上述側(cè)掃聲吶數(shù)據(jù)測量過程中拖體的航跡和三維姿態(tài)均存在不同程度的跳躍，幾何校正前需要對數(shù)據(jù)進行預(yù)處理［22－23］。側(cè)掃聲吶文件記錄的導(dǎo)航數(shù)據(jù)由于更新率比聲吶儀的采樣速率低，導(dǎo)致多ping共用同一坐標(biāo)，如圖8中200ping數(shù)據(jù)只有18組有效坐標(biāo)（圖9）。本文使用多項式插值得到不重疊的各ping坐標(biāo)。對于跳躍的姿態(tài)數(shù)據(jù)，本文采取了中值濾波和均值平滑相結(jié)合的處理方式。預(yù)處理前后拖魚的坐標(biāo)和艏向分別見圖9和圖10所示。預(yù)處理后，拖體速度和姿態(tài)最大值、最小值統(tǒng)計見表1。根據(jù)拖魚高度對聲吶數(shù)據(jù)進行幾何校正，幾何校正后的聲吶數(shù)據(jù)如圖8所示，該圖坐標(biāo)已經(jīng)從時間域轉(zhuǎn)換為空間域，消除了各種變形和位移因素的影響。

圖8 幾何校正后聲吶數(shù)據(jù)Fig.8 Sonar data after geometric correction

圖9 預(yù)處理前后的坐標(biāo)Fig.9 Navigation before and after pre－processing

圖10 預(yù)處理前后的拖魚艏向Fig.10 Heading before and after pre－processing

表1 預(yù)處理后姿態(tài)數(shù)據(jù)統(tǒng)計結(jié)果Table 1 Statistics of attitude after pre－processing

4.2 插值間隔的選取

根據(jù)聲吶水平波束角、斜距、脈沖寬度、水中聲速和聲吶距海底的高度等信息及公式（1）～（3），可得縱向采樣間距為0.182～1.571m，橫向采樣最小間距為0.075 5m，向中心逐漸增大。

根據(jù)聲吶速度、脈沖間隔和采樣點個數(shù)等信息及式（4）～（6），可以獲得縱向采樣點分布間隔為0.141 4～0.190 8m，橫向采樣間隔為0.044 2～0.806 0m。

為了不損失原數(shù)據(jù)的信息，用最小的網(wǎng)格作為重采樣間隔，橫向重采樣間隔選擇0.044 2m，縱向重采樣間隔選擇0.141 4m，分別選擇東西向和南北向作為x和y軸進行采樣，通過投影可得，x方向和y方向重采樣間隔分別選擇0.033和0.105m。

5 結(jié)果分析

使用0.033和0.105m作為重采樣間隔，對扭曲網(wǎng)格三角剖分有限元插值算法進行重采樣實驗，并將所得結(jié)果與反距離加權(quán)法、克里格法和三角剖分線性插值法所得結(jié)果從定性及定量角度進行比較與評價。

5.1 定性評價

通過對搜索橢圓x方向半軸長為0.10m，y方向半軸長為0.15m時反距離加權(quán)法和克里格法重采樣圖像分析的結(jié)果表明，當(dāng)選擇較小的搜索半徑時，會造成很多區(qū)域空白（圖11），不能獲得完整的區(qū)域數(shù)據(jù)。如果加大搜索半徑（搜索橢圓x方向半軸為1.80m，y方向半軸為0.20m），則又因密集區(qū)參與插值的點太多，而造成平滑效應(yīng)，降低了圖像分辨率（圖12a和圖12b）。

普通三角剖分線性插值法（圖12c）與扭曲網(wǎng)格三角剖分有限元法（圖12d）重采樣均可以很好地保留原始圖像的信息，具有較高的分辨率。

圖11 重采樣結(jié)果（搜索橢圓半徑：橫向0.10m，縱向0.15m）Fig.11 Resampling results（setting the radius across－track and radius along－track of the searching ellipse to 0.10mand 0.15m）

圖12 重采樣結(jié)果（搜索橢圓半徑：橫向1.80m，縱向0.20m）Fig.12 Resampling results（setting the radius across－track and radius along－track of the searching ellipse to 1.80mand 0.20m）

5.2 定量分析

5.2.1 重采樣精度評價

本文選擇標(biāo)準(zhǔn)差（σ）、偏度（S）、峰度（K）、信息熵（En）為指標(biāo)，對不同插值方法的重采樣結(jié)果進行對比分析與評價。設(shè)N為采樣點總數(shù)，xi為第i個采樣點的屬性值，ˉx為采樣點屬性平均值。

標(biāo)準(zhǔn)差反映圖像采樣點相對于平均值的離散情況，在一定程度上反映圖像信息量的大小，標(biāo)準(zhǔn)差越大，反映圖像信息量越豐富。標(biāo)準(zhǔn)差計算公式為

偏度描述了信息分布偏離對稱分布的程度，偏度越大，數(shù)據(jù)的不對稱性越強；峰度則描述了信息分布的陡峭程度，峰度越大，數(shù)據(jù)密度函數(shù)曲線越陡。因此，偏度和峰度值越大，說明圖像越偏離正態(tài)分布，所含信息量越大。偏度和峰度計算公式分別為

信息熵值是圖像所具有信息量的度量，是測量數(shù)據(jù)屬性分布隨機性的特征參數(shù)，表征了圖像中紋理的復(fù)雜程度。數(shù)據(jù)的屬性越均勻，熵值越小；反之，數(shù)據(jù)的屬性越復(fù)雜，則熵值越大。對數(shù)據(jù)屬性為［0，M］范圍的聲吶數(shù)據(jù)進行統(tǒng)計，nk為第k級屬性的數(shù)據(jù)個數(shù)，n為數(shù)據(jù)總數(shù)。信息熵計算公式為

式中，Pk＝nk／n（k＝0，…，M）。

根據(jù)上述評價指標(biāo)及其計算方法，得到各指標(biāo)的評價結(jié)果（表2）。扭曲網(wǎng)格三角剖分有限元插值法重構(gòu)后的標(biāo)準(zhǔn)差、偏度、峰度和信息熵均高于其他3種方法，表明該方法信息量豐富，重采樣圖像分辨率高。

表2 圖12中4幅圖的統(tǒng)計分析結(jié)果Table 2 Statistics of the four sidescan sonar images shown in Fig.12

5.2.2 重采樣計算效率評價

通過對4種方法重采樣計算時間（表3）進行比較可知，基于扭曲網(wǎng)格三角剖分有限元插值技術(shù)的重采樣計算速度高于普通三角剖分線性插值法、反距離加權(quán)法插值法、克里格法插值法的重采樣速度，結(jié)果表明其適合實際應(yīng)用中大數(shù)據(jù)量的處理。

表3 側(cè)掃聲吶重采樣計算時間Table 3 Computing time of sidescan sonar image resampling

6 結(jié)語

基于側(cè)掃聲吶原始數(shù)據(jù)采樣間隔特點的分析，確定了橫向和縱向重采樣間隔，充分保留了原始數(shù)據(jù)的信息。

根據(jù)側(cè)掃聲吶數(shù)據(jù)幾何校正后的空間分布特征，采用了基于扭曲網(wǎng)格的三角剖分有限元插值的重采樣方法對其進行了重構(gòu)。該方法與常用的一些插值算法的對比表明，本文采用的方法能很好地適應(yīng)側(cè)掃聲吶數(shù)據(jù)排列的特點，重構(gòu)后不僅保留了數(shù)據(jù)的細(xì)節(jié)信息，而且節(jié)省了三角剖分的時間，具有較高的計算效率，適合實際應(yīng)用中大數(shù)據(jù)量的處理。

（References）：

［1］YANG F L，LIU J N，ZHAO J H，et al.Seabed texture classification using BP neutral network based on GA［J］.Science of Surveying and Mapping，2006，31（2）：111－114.陽凡林，劉經(jīng)南，趙建虎，等.基于遺傳算法的 BP網(wǎng)絡(luò)實現(xiàn)海底底質(zhì)分類［J］.測繪科學(xué)，2006，31（2）：111－114.

［2］CARMICHAEL D R，LINNETT L M，CLARKE S J，et al.Seabed classification through multifractal analysis of sidescan sonar imagery［J］.IEE Proceedings－Radar，Sonar and Navigation，1996，143（3）：140－148.

［3］GAUER R C，MCFADZEAN A，REID C.An automated sidescan sonar pipeline inspection system［C］∥Oceans＇99MTS／IEEE.Riding the Crest into the 21st Century.Houston：IEEE，1999，2：811－816.

［4］LAI X H，PAN G F，GOU Z K，et al.Study on application of sidescan sonar in submarine pipeline inspection［J］.The Ocean Engineering，2011，41（3）：117－121.來向華，潘國富，茍諍慷，等.側(cè)掃聲吶系統(tǒng)在海底管道檢測中應(yīng)用研究［J］.海洋工程，2011，41（3）：117－121.

［5］FUSI N，KENYON N H.Distribution of mud diapirism and other geological structures from long－range sidescan sonar（GLORIA）data，in the Eastern Mediterranean Sea［J］.Marine Geology，1996，132（1）：21－38.

［6］ZITTER T A C，HUGUEN C，WOODSIDE J M.Geology of mud volcanoes in the eastern Mediterranean from combined sidescan sonar and submersible surveys［J］.Deep Sea Research Part I：Oceanographic Research Papers，2005，52（3）：457－475.

［7］BROWN C J，COOPER K M，MEADOWS W J，et al.Small－scale mapping of sea－bed assemblages in the eastern English Channel using sidescan sonar and remote sampling techniques［J］.Estuarine，Coastal and Shelf Science，2002，54（2）：263－278.

［8］COCHRANE G R，LAFFERTY K D.Use of acoustic classification of sidescan sonar data for mapping benthic habitat in the Northern Channel Islands，California［J］.Continental Shelf Research，2002，22（5）：683－690.

［9］REED S，PETILLOT Y，BELL J.Automated approach to classification of mine－like objects in sidescan sonar using highlight and shadow information［J］.IEE Proceedings：Radar，Sonar and Navigation，IET，2004，151（1）：48－56.

［10］CHAVEZ J P S，ISBRECHT J A，GALANIS P，et al.Processing，mosaicking and management of the Monterey Bay digital sidescansonar images［J］.Marine Geology，2002，181（1）：305－315.

［11］ZHANG J B，PAN G F，DING W F.Research on Side－scan Sonar Images＇Correction［J］.Technical Acoustics，2009，28（6）：44－47.張濟博，潘國富，丁維鳳.側(cè)掃聲吶圖像改正研究［J］.聲學(xué)技術(shù)，2009，28（6）：44－47.

［12］LU H X，LIU X J，WANG Y J，et al.Noise error analysis of slope algorithms based on grid DEM derived from interpolation［J］.Acta Geodaetica et Catographica Sinica，2012，11（6）：926－932.盧華興，劉學(xué)軍，王永君，等.插值條件下格網(wǎng)DEM 坡度計算模型的噪聲誤差分析［J］.測繪學(xué)報，2012，11（6）：926－932.

［13］LOU X L，HUANG W G，ZHOU C B，et al.A method for fast resampling of remote sensing imagery［J］.Journal of Remote Sensing，2002，6（2）：96－101.樓琇林，黃韋艮，周長寶，等.遙感圖像數(shù)據(jù)重采樣的一種快速算法［J］.遙感學(xué)報，2002，6（2）：96－101.

［14］LIN Z H，MO X G，LI H X，et al.Comparison of three spatial interpolation methods for climate variable in China［J］.Acta Geographica Sinica，2002，57（1）：47－56.林忠輝，莫興國，李宏軒，等.中國陸地區(qū)域氣象要素的空間插值［J］.地理學(xué)報，2002，57（1）：47－56.

［15］CHANG W Y，DAI X G，CHEN H W.A case study of geoestatistical interpolation to meteorological fields［J］.Chinese Journal of Geophysics，2004，47（6）：982－990.常文淵，戴新剛，陳洪武.地質(zhì)統(tǒng)計學(xué)在氣象要素場插值的實例研究［J］.地球物理學(xué)報，2004，47（6）：982－990.

［16］HERITAGE G L，MILAN D J，ANDREW R G，et al.Influence of survey strategy and interpolation model on DEM quality［J］.Geo－Morphology，2009，（112）：334－334.

［17］LIN Z L，ZHU Q.Digital elevation model［M］.2nd ed.Wuhan：Wuhan University Press，2003.李志林，朱慶.數(shù)字高程模型［M］.第二版.武漢：武漢大學(xué)出版社，2003.

［18］MAZEL C.Side scan sonar record interpretation（Training Manual）［M］.Salem，New Hampshire：Klein Associates，INC，1985：14－16.

［19］LI X，CHENG G D，LU L.Advance in Earth Science，2000，15（3）：260－265.李新，程國棟，盧玲.空間內(nèi)插方法比較［J］.地球科學(xué)進展，2000，15（3）：260－265.

［20］WU X B，WANG S X，XIAO C S.A new study of Delaunay triangulation creation［J］.Acta Geodaetica et Catographica Sinica，1999，28（1）：30－37.武曉波，王世新，肖春生.Delaunay三角網(wǎng)的生成算法研究［J］.測繪學(xué)報，1999，28（1）：30－37.

［21］COOK R D，MALKUS D S，PLESHA M E，et al.Concepts and applications of finite element analysis［M］.4th ed.USA：John Wiley＆Sons Ine.，2001.

［22］PHILIPPE B.The handbook of Sidescan Sonar［M］.USA：Springer，2009.

［23］LE BAS T P，MASON D C，MILLARD N C.TOBI image processing－the state of the art［J］.IEEE Journal of Oceanic Engineering，1995，20（1）：85－93.