三角剖分有限元法在声呐图像重构中的应用

罗进华，杨修伟，朱凌霄，朱培民

（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）取决于脉宽的大小（图1）。理论上讲，两个目标的最小距离应至少等于脉宽的一半［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）。波束宽度随与声呐的距离增加而逐渐变宽，因此，纵向分辨率从中心向两侧逐渐降低。纵向分辨率Δ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.