三角剖分有限元法在聲吶圖像重構中的應用

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

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

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

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

1 重采樣間隔的選擇

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

1.1 側掃聲吶圖像分辨率

側掃聲吶數據分辨率分為橫向（垂直航跡）分辨率和縱向（沿航跡）分辨率。

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

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

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

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

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

1.1.2 縱向（沿航跡）分辨率

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

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

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

1.2 數據分布特征

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

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

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

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

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

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

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

1.2.3 姿態對分布特征的影響

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

1.3 重采樣間隔的選擇原則

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

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

當分辨率大于采樣間隔時，由于不同的散射值平均，產生平滑效應而丟失了小尺寸細節信息，此時選擇最小采樣間隔作為重采樣間隔，可以最大限度的保持原始圖像的信息量，減少原始數據信息的損失。

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

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

2 插值方法的選擇

2.1 常規插值算法的適用性

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

2.2 扭曲網格三角剖分有限元插值方法

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

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

2.3 重疊區域的處理方法

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

3 扭曲網格三角剖分有限元插值的基本原理

3.1 扭曲網格三角剖分

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

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

3.2 面積坐標

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

圖6 三角形單元的面積坐標Fig.6 Area coordinates of triangular element

3.3 插值函數

三角形有限元插值函數的基本形式為

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

4 實驗設計

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

圖7 原始聲吶數據Fig.7 Raw sonar data

4.1 側掃聲吶數據插值前處理

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

圖8 幾何校正后聲吶數據Fig.8 Sonar data after geometric correction

圖9 預處理前后的坐標Fig.9 Navigation before and after pre－processing

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

表1 預處理后姿態數據統計結果Table 1 Statistics of attitude after pre－processing

4.2 插值間隔的選取

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

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

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

5 結果分析

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

5.1 定性評價

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

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

圖11 重采樣結果（搜索橢圓半徑：橫向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 重采樣結果（搜索橢圓半徑：橫向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 重采樣精度評價

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

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

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

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

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

根據上述評價指標及其計算方法，得到各指標的評價結果（表2）。扭曲網格三角剖分有限元插值法重構后的標準差、偏度、峰度和信息熵均高于其他3種方法，表明該方法信息量豐富，重采樣圖像分辨率高。

表2 圖12中4幅圖的統計分析結果Table 2 Statistics of the four sidescan sonar images shown in Fig.12

5.2.2 重采樣計算效率評價

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

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

6 結語

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

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

（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.陽凡林，劉經南，趙建虎，等.基于遺傳算法的 BP網絡實現海底底質分類［J］.測繪科學，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.來向華，潘國富，茍諍慷，等.側掃聲吶系統在海底管道檢測中應用研究［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.張濟博，潘國富，丁維鳳.側掃聲吶圖像改正研究［J］.聲學技術，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.盧華興，劉學軍，王永君，等.插值條件下格網DEM 坡度計算模型的噪聲誤差分析［J］.測繪學報，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.樓琇林，黃韋艮，周長寶，等.遙感圖像數據重采樣的一種快速算法［J］.遙感學報，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.林忠輝，莫興國，李宏軒，等.中國陸地區域氣象要素的空間插值［J］.地理學報，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.常文淵，戴新剛，陳洪武.地質統計學在氣象要素場插值的實例研究［J］.地球物理學報，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.李志林，朱慶.數字高程模型［M］.第二版.武漢：武漢大學出版社，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.李新，程國棟，盧玲.空間內插方法比較［J］.地球科學進展，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三角網的生成算法研究［J］.測繪學報，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.