应用于声纳的高性能扇形插值方法研究

席伟光，童大为，费玉杰，徐天柏

（海鹰企业集团有限责任公司，江苏无锡，214115）

0 引言

高频多波束图像声纳因其卓越的图像质量，在海洋工程、环境监测等民用领域以及反水雷、反潜等军用领域都得到了广泛的应用[1]。

高频多波束图像声纳的探测页面是一个以声基阵[2]开角为圆心角的扇形区域，只有采用扇形显示方式，才能最大限度地还原水下环境的真实面貌，因此需要将信号处理输出的二维矩阵图像转换成扇形图像。声纳作为一种专用设备，集成度较高，配备的软硬件资源有限，很难适配通用扇形插值软件[3]，且大部分通用软件的扇形插值功能非专为声纳设计，在图像质量和实时性[4]方面无法满足使用要求。

RTheta数字扫描变换算法[5]的插值效果逼真，能满足声纳图像质量要求，但其复杂的计算过程无法满足声纳的实时性要求，故优化算法是本文研究的重要内容。

本文以高频多波束声纳图像为研究对象，以R-Theta数字扫描变换算法为基础建立扇形插值的数学模型；使用“查表法”[6]“图形对称法”“通用化设计”优化数学模型；使用面向对象的程序设计方法完成了软件设计；最后通过声纳设备的试验结果，说明方法的有效性。

1 扇形插值原理

本文讨论的扇形插值原理基于R-Theta数字扫描变换算法，它是一种二维线性插值算法[7]，该算法最先应用于B超仪，其原理是利用极坐标与直角坐标的转换关系，将B超探头扫描获得的波束数据[8]映射到矩形的显示区域中。

高频多波束图像声纳与B超仪的工作原理类似。本文以某型前视声纳为例进行说明，该声纳的探测开角为100度，探测扇面内均匀分布了200个波束，波束宽度[9]为0.5°，每个波束上等间隔分布600个样本点，假设扇形的右侧半径为角度参考线，即θ=0°，由此可求得任意波束上的任意样本点的径向距离R和波束偏移角（波束与参考线的夹角）θ。

算法的目的是把如图1所示的扇形波束图映射到存贮器矩阵中。为了说明上的方便，将图1中的局部存贮空间进行放大，局部放大图见图2，每个存贮单位（+号所在位置）中应填入合适的像素值。

图1 扇形波束图

图2 扇形图像的局部放大图

对于处在相邻两个波束之间的所有像素值，使用下列3个方程来插值。

(x,y)是扇形显示区域内的任意一点，s(*,*)是(x,y)相邻的四个样本点的值。而z(i)和z(i+1)是计算f(x,y)的过渡值，Rerr是径向距离误差，θerr是方位角相对误差，由下式确定：

其中R'为点(x,y)的径向距离，θ'为(x,y)与角度参考线的夹角，ΔR为同一波束上相邻采样点(远端)的径向距离差，Δθ为波束宽度。

综上，将(x,y)相邻两个波束上的4个样本点s(i,j)、s(i,j+1)、s(i+1,j)、s(i+1,j+1)代入公式⑴⑵⑶，可得：

由公式⑹可求得显示区域扇形范围内任意一点的像素值。

2 方法设计

2.1 建立数学模型

结合某型前视声纳，以扇形的圆心为原点建立直角坐标系，如图3所示。

图3 建立直角坐标系的扇形波束图

在直角坐标系中，矩形显示区中的任意一点都可以用(x,y)表示。根据扇形插值原理：

直角坐标系内的任意一点(x,y)可以转化为极坐标形式(R',θ')。

θ'的取值范围为根据图2设s(i+1,j+1)对应点的极坐标为(R,θ)，则s(i,j+1)对应的极坐标为(R,θ-Δθ)，(i,j)对应的极坐标为(R-ΔR,θ-Δθ)，(i+1,j)对应的极坐标为(R-ΔR,θ)，根据图2中(x,y)和s(i+1,j+1)的位置关系易推出：

将⑺⑻⑼⑽代入⑾⑿可知s(i+1,j+1)的极坐标可由(x,y)确定，同理f(x,y)周围的其余三个样本点的极坐标均可由(x,y)确定。

将⑺⑻⑼⑽⑾⑿代入⑷⑸可知Rerr和θerr由(x,y)确定。由⑹式可知，f(x,y)的值可由周围四点s(*,*)的值、Rerr和θerr确定，因此只需求得s(*,*)的值，即可得到f(x,y)的值。下文介绍了通过s(*,*)的极坐标获取s(*,*)的值的方法。

如图4所示，将待进行扇形插值的200×600（W×H）个样本点的声纳波束数据按顺序分布到内存中。

设图2中点s(i+1,j+1)，s(i,j+1)，s(i,j)，s(i+1,j)对应图4中点A、B、C、D；波束样本点j的取值范围为（0,599），波束号i的取值范围为（0,199）。

设内存中波束数据的起始地址为pBmembuf[0],则图4中点A、B、C、D对应的内存值如下：

图4 声纳波束数据内存分布图

A：pBmembuf[200(j + 1)+ 200 - i - 1];

B：pBmembuf[200(j + 1)+ 200 - i ];

C：pBmembuf[200j + 200 - i];

D：pBmembuf[200j + 200 - i - 1];

上文设s(i+1,j+1)对应的极坐标为(R,θ)，且由⑾⑿可知R是ΔR的整数倍，θ是Δθ的整数倍,结合⑺⑻⑼⑽可进一步推出：

由⒀⒁可知i和j可通过(x,y)计算获得，则点A、B、C、D的内存序号也可通过(x,y)计算获得。

由此可实现将直角坐标系的扇形波束图中的任意一点(x,y)周围的四点映射为波束域数据中对应的内存值，再根据公式⑹计算得到f(x,y)的值。

将图3扇形显示区内的像素点按此方法逐一赋值，即可完成扇形插值。

2.2 优化数学模型

在用软件实现R-Theta算法中，需要多次调用三角函数、浮点运算等大开销算子，为了节省运算资源，提高图像刷新率和系统实时性，需优化数学模型。由⒀⒁可知，计算任意一个f(x,y)的值，至少需要进行一次和一次计算，以图3显示区尺寸为800×650（W×H）为例，需分别进行52万次计算，采用查表法和图形对称法可大幅降低计算量。通过上述优化，在奔4处理器上测试扇形插值算法，150米量程的单图刷新时间从140ms降为80ms，性能明显提高。

本例列举的数学模型是根据某型前视声纳的专用需求设计的，不同声纳对扇形显示的需求不同，具体体现在扇形圆心角、扇形半径、显示区尺寸、显示量程等参数有差异。为了使数学模型更具通用性，应进行通用化设计。

2.2.1 查表法

查表法的原理是通过预处理[10]的方式提前计算复杂算法在一定范围的运算结果，并按规律储存到指定内存中，当程序需要调用算法时，查询并访问内存中的运算结果即可，因计算机访问内存地址的速度远高于执行算法本身[11]，故可达到降低计算量的目的。使用查表法的约束条件是：算法输入值的范围有限且可确定。

图3显示区尺寸为800×650（W×H），则x的取值范围为[-400,400]，y的取值范围为[0,650]，可通过以下方法将取值范围内的进行预处理。

初始化预处理结果存储数组：float m_sqrtTable [401][651]; float m_atanTable [401][651];

预处理伪代码如下：

2.2.2 图形对称法

图形对称法的原理是利用显示区左右两侧以y轴对称的特性，首先使用数学模型计算得到显示区一侧像素点相邻四点的内存序号以及Rerr和θerr的值，再利用轴对称性质，推算出另一侧对称点相邻四点的内存序号及Rerr和θerr的值，该方法可将复杂算法和的运算次数降低50%，大幅提升运算效率。

举例说明：图4中点A、B、C、D为(x,y)的四个相邻点，由轴对称性质可知，一定存在对称点A'、B'、C'、D'、(x',y')，已知点A对应的内存序号为200(j+1)+200-i-1，则很容易推出点A'的内存序号为200(j+1)+i，同理可推出点B'、C'、D'的内存序号，且R'err=Rerr；θ'err=θerr。然后根据公式⑹计算得到，f(x',y')的值。

2.2.3 通用化设计

高频多波束图像声纳的探测开角、波束数、采样点数均为声纳固有参数，不同型号的声纳设备或有差异；而显示量程为声纳工作参数，可在声纳使用过程中根据按需设置。以上参数均为扇形插值的输入值，参与了扇形插值的计算过程，只有将这些参数对用户开放，才能使扇形插值算法更具通用性。

下文介绍了声纳固有参数通用化和显示量程参数通用化的详细设计方法。

2.2.3.1 声纳固有参数通用化

本节只讨论显示起始量程为0米的情况，非0情况在“显示量程参数通用化”中介绍。

探测开角决定了扇形圆心角，扇形圆心角结合波束数决定了Δθ，采样点数决定了扇形半径和ΔR，声纳图像显示区尺寸可通过以上参数估算获得。

在程序中，声纳固有参数的用户接口，可通过宏定义设定：

2.2.3.2 显示量程参数通用化

图像声纳的距离分辨率[12]较高，可达厘米级，假设某声纳的距离分辨率为2.5cm，则在200m量程下每个波束包含的采样点数为8000，远高于普通显示器的分辨率，为了实现满量程下的图像显示，一般会对样本点进行抽点处理，以输出符合图像显示区尺寸的图像数据，即牺牲了距离分辨率来满足全量程显示的需求；但也可通过缩小显示量程范围的方法，提升声纳图像的距离分辨率。

如图5中(1)所示，声纳全量程为0～200m，当目标处于50～100m范围内时，为了提高目标的距离分辨率，可将50～100m量程区间（阴影部分）的图像进行放大显示，50m量程定义为显示起始量程，100m量程定义为显示终止量程。

图5 显示量程放大图

一般有两种方法实现显示量程范围内的图像放大：

方法1：如图5中(2)所示，首先将显示量程内的样本点重新抽点处理，使之与扇形半径的像素点一致，然后将显示量程内的数据按照扇形的半径、圆心角进行插值，从图上可以发现，显示起始量程50m处的圆弧已缩小成了一个点，丢失了大量有效数据，越接近显示起始量程，丢失的有效数据元素越多，同时显示终止量程100m处的圆弧被拉伸了一倍，降低了方位分辨率。可见本方法与声纳实际应用场景有较大差异，容易造成目标变形。

方法2：如图5中(3)所示，同样，首先将显示量程内的样本点重新抽点处理，使之与原扇形半径的像素点一致，然后保持显示起始量程50m处的弧长和显示终止量程100m处的弧长不变，并将两圆弧间的环宽拉伸至与原扇形半径一致，如阴影部分所示，阴影部分是新生成扇形中的一部分，根据上文数学模型，只需知道扇形的圆心角、半径、波束数即可完成扇形插值。新扇形的波束数不变，圆心角β和半径R可通过以下方法计算获得：

r为原扇形半径，endrange为显示终止量程，startrange为显示起始量程，α为原扇形圆心角。由公式（15）可知endrange与startrange越接近新扇形的半径R越大，应用中可确定endrange与startrange间的最小差值，即声纳在最大距离分辨率时，显示区能够显示的最大量程范围。

方法2最大限度还原了声纳探测扇面，声纳目标变形小。

2.3 软件设计

软件设计基于C++编程语言，建立了扇形图像插值类：PolarImage。共由三个文件组成：PolarImage.h、PolarImage.cpp、UsrIf.h，PolarImage.h声明了PolarImage类，PolarImage.cpp实现了类成员函数，UsrIf.h为用户提供了声纳固有参数和显示尺寸的配置接口。

如图6所示为PolarImage类的结构组成。

从图6可知transform()函数是PolarImage类提供的唯一用户接口，其完整声明应为：static bool transform (unsigned char *InBuf,unsigned char *OutBuf, float startrange, float endrange)；InBuf为待插值的波束数据所在的内存地址，OutBuf为插值结果存放地址，startrange为显示起始量程，endrange为显示终止量程。

图6 PolarImage类结构图

如图7所示为声纳应用程序及PolarImage类成员之间调用关系。声纳应用程序利用UsrIf.h完成对静态成员变量的初始化。

图7 函数调用关系图

声纳应用程序需要进行扇形插值时调用transform()函数，该函数先对输入参数进行有效性检查，通过后调用instance()函数获取PolarImage类对象（若对象不存在，则instance()函数创建对象再返回，创建对象时调用构造函数PolarImage()对m_sqrtTable和m_atanTable表空间进行初始化），再利用获取的PolarImage类对象调用_transform()函数，该函数完成了扇形插值数学模型的代码实现，是核心函数。

声纳应用程序删除PolarImage类对象时调用析构函数～PolarImage()，该函数释放了m_sqrtTable和m_atan Table表空间。

3 试验效果

本文研究的扇形插值方法已成功应用于多个型号的声纳设备。

3.1 某型前视声纳应用效果图

通过水池试验[13]验证了扇形显示效果，试验水池为矩形，长宽约为20×6m，在声基础阵正前方约15m处放置了一个直径约60cm的环型铁块，正前方1m处放置了一个直径约2cm的铁球，分别采用矩形显示和扇形显示两种方式进行了对比。

从图8矩形显示图可见，原本矩形的水池壁，被显示成塔形，且越接近基阵目标变形越明显，1m处的小球尺寸已严重失真，远大于实际尺寸，15m处的环型铁块失真不明显。说明在矩形显示模式下，目标越近失真越明显。

图8 前视声纳矩形显示图

从图9扇形显示图可见，矩形的水池壁无变形，1m处的小球呈现为1个亮点，15m处的环型铁块与矩形图差别不大。可见扇形模式下，近距离目标无失真。

图9 前视声纳扇形显示图

经测试，在奔四处理器上完成图像扇形插值时间约为80ms；在i7 2.4G处理器上完成图像扇形插值时间约为30ms。均能满足前视声纳的实时性要求。

3.2 某型分类识别声纳应用效果

该型声纳的固定参数如下：

探测开角：30°；

波束数：512；

波束上样本点数：1024；

扇形图像显示区大小：540×1038（W×H）；

矩形图像显示区大小：512×1024（W×H）；

通过湖试[14]验证了扇形显示效果，试验目标为一根φ0.066m×64m的阵缆[15]，将阵缆悬浮（利用浮球）于水中，两端用绳子固定，使其姿态保持直线，左端距离声基阵40余米，右端距离声基阵60余米。

图10为全量程[16]（0～350m）下的矩形显示图和扇形显示图，图中白圈标注部分为目标阵缆，可见在全量程模式下声纳图像显示的范围较大，但距离分辨率较低；矩形显示模式下的阵缆图像被拉伸发生了变形，扇形显示模式下的阵缆目标图像较小，真实还原了声基阵的水中姿态。

图10 识别声纳全量程显示图

图11 为显示量程 (0～80m)下的矩形显示图和扇形显示图。通过减小显示量程范围提升了目标的距离分辨率；矩形显示模式下的阵缆图像被拉伸发生了一定程度的变形，呈现为弧形；扇形显示模式下能够还原基阵的真实姿态。

图11 识别声纳0～80m显示量程显示图

图12 为显示量程(30～80m)下的矩形显示图和扇形显示图。显示量程范围进一步缩小，距离分辨率得到更大提高；矩形显示模式下的阵缆图像被拉伸发生了一定程度的变形，呈现为弧形，扇形显示模式下基本能够还原基阵的真实姿态，但较之图11的扇形图像有了些许弯曲。

图12 识别声纳30～80显示量程图

经测试在TMS320C6678处理器[17]上完成图像扇形插值的时间约为50ms。

4 结语

本文以高频多波束声纳图像为研究对象，采用R-Theta数字扫描变换算法建立了扇形插值的数学模型；根据装备实际使用要求，提出三种优化模型的方法：（1）查表法；（2）图形对称法；（3）通用化设计；并针对数学模型完成了软件设计。通过两型声纳设备的试验验证结果表明，该扇形插值方法平台适应性强，且图像质量和实时性均能满足声纳使用要求，可在行业内进一步推广应用。