0 引　言

光纤水听器阵列凭借其高灵敏度、抗干扰能力强、大规模组阵能力及高可靠性等，成为现代海洋监测和水下探测的重要技术手段^{[1 − 4]}。在近岸浅水区域，大规模的光纤水听器通常采用掩埋式布放，通过铠装光缆将感知的声学信号传输至岸基监测站，实现对水下目标的长期、隐蔽监测^[4]。

与传统的浮标式或拖曳式阵列相比，海底掩埋式部署进一步减少了外部环境对阵列的干扰和破坏，提高了系统的稳定性和使用寿命。然而，在实际埋设布放施工过程中，由于阵列尺度较大，受海底地形、沉积物特性、布放船机动性以及海流海浪变化等因素影响，往往导致实际阵元埋设位置与理论预期位置存在偏差^[5]。若不及时予以校准而采用预定阵形的波束形成算法，将导致波束分裂、测向不准确，甚至算法性能恶化失效^[6]。因此，有必要采取有效的措施对阵列的阵元位置进行精确校准。

尽管目前存在许多阵列位置校准方法，但能被切实用于工程实际的却很少见。基于传感器直接测量方法^[7]，可实现阵列位置的实时监测，但这种方法的阵位校准精度受传感器数量和自身精度的限制，增加了系统工程实施复杂度，降低了系统的可靠性。自校准方法不需要知道精确的信源方位，比有源校准方法对环境要求更低，但实际应用中由于阵列误差的复杂性，该方法难以保证参数估计的唯一性，且算法对初始值和门限的选取较为敏感，容易陷入局部最优解^[8]。有源校准方法是利用方位已知的辅助声源对阵列误差进行校准^{[5 − 7,9 − 12]}。该方法可以通过恰当选取符合外场试验要求的校准源来提高阵列校准精度。因此，有源校准方法具有较好的环境适应性和较高的校准精度。目前，经海试验证较为有效的有源校准方法主要分为两类：一类是基于时延估计的有源校准方法^[7,10]，通过布设2个或多个成一定张角的校准声源，利用它们估计各阵元信号的时延关系，利用声源和接收阵的GPS位置确定各阵元的位置。该方法工程实施简单、计算代价小，但校准精度受时延估计算法的性能影响较大。一类是基于子空间分解的有源校准方法^[5,6,11]，通过构建目标函数并利用粒子群等优化算法求解最优阵列位置。这类方法对校准源的信号特性要求相对较高，如果校准源的方位信息不准确或信号质量不理想，可能导致校准精度下降。此外，上述两类方法的理论框架主要建立在认为辅助声源为窄带或远场平面波条件下，而对于工作于浅海环境低频掩埋式大尺度阵列而言，远场校准面临水声信道多途干扰和辅助声源信号远距离传播信噪比不足问题^[9]。尽管文献[5,12]提出了近场有源校准解决方案，但算法在真实海洋环境下有效性和工程适用性还有待进一步验证。

针对浅海大尺度海底掩埋式光纤水听器阵列位置校准的工程需求实际，本文提出一种融合时延估计和遗传优化算法的三辅助源阵列位置校准方法。仿真和海试结果表明，本文提出的方法具有较理想的校准精度和环境适应性。

1 阵列位置校准模型

如图1所示，假定有$ N $个阵元组成的均匀线列阵埋设于海床上，阵元间距为$ d $ m。不失一般性，这里仅考虑二维$ {X}{Y} $平面内的阵列位置校准问题，认为其阵列位置仅在水平面上有一个缓慢的随机扰动。

假设每个阵元的预期埋设位置为 $ {p}_{0n}= \left({x}_{0n}, {y}_{0n}\right), $ $ n=\mathrm{1,2},\dots ,N $，各阵元的位置误差为$ {\delta }_{n}=\left(\Delta {x}_{n},\Delta {y}_{n}\right) $，选取第一个阵元$ {p}_{1}=\left({x}_{1},{y}_{1}\right) $作为参考阵元，并假定位置已知，其所在位置设为坐标轴的原点，建立坐标系。则阵列结构可以表示为：

$ \mathit{p}=\left[\begin{matrix} p_{1}\\ {p}_{2}\\ \vdots\\ {p}_{N}\end{matrix}\right]=\left[\begin{matrix}\begin{matrix}{x}_{1}& {y}_{1}\\ {x}_{2}& {y}_{2}\end{matrix}\\ \begin{matrix}\vdots&\vdots\\ {x}_{N}& {y}_{N}\end{matrix}\end{matrix}\right]=\left[\begin{matrix}\begin{matrix}{x}_{01}+\Delta {x}_{1}& {y}_{01}+\Delta {y}_{1}\\ {x}_{02}+\Delta {x}_{2}& {y}_{02}+\Delta {y}_{2}\end{matrix}\\ \begin{matrix}\vdots&\vdots\\ {x}_{0N}+\Delta {x}_{N}& {y}_{0N}+\Delta {y}_{N}\end{matrix}\end{matrix}\right]。$

(1)

众所周知采用远场辅助源校准方法要求声源与阵列的距离$ R $满足^[11]：$ R\geqslant 2{L}^{2}/\lambda $，$ L $为阵列长度，$ \lambda $为辅助声源信号的波长。而对于工作于浅海环境的低频掩埋式大尺度阵列而言，远场校准距离$ R $达到公里级，导致存在水声信道多途干扰强、辅助声源发射声源级不足等现实问题。因此，本文提出首先采用近场双辅助源校准方法估计阵列的初始位置。选取2个（或多个）在空间和时间上均不相干的宽频带辅助声源$ {s}_{m}\left(t\right), m \geqslant 2 $，辅助声源的位置位于阵列的菲涅耳区域，声源以球面波形式向外辐射声波，则阵列接收到的信号具有距离相关的传播特性，在时刻$ t $对应第$ m $个声源，阵元$ n $的接收信号可以表示为：

$ {y}_{n}\left(t\right)={\alpha }_{mn}{s}_{m}\left(t-{\tau }_{mn}\right)+{w}_{n}\left(t\right)。$

(2)

式中：$ {s}_{m}\left(t\right) $为第$ m $个辅助声源的发射信号；$ {\alpha }_{mn}={r}_{m1}/ {r}_{mn} $为信号幅度，随传播距离衰减；$ {r}_{nm}= ||{p}_{n}-{s}_{m}||$为声源$ m $到阵元$ n $的距离；$ {\tau }_{mn}={r}_{nm}/c $声源$ m $到阵元$ n $的传播时间，$ c $为声波海水中的传播速度；$ {w}_{n}\left(t\right) $为均值为0，方差为$ {\sigma }^{2} $的高斯白噪声。

这里以双辅助源校准方案为例，假定声源$ {s}_{1} $和$ {s}_{2} $分别位于$ \left({x}_{s1},{y}_{s1}\right) $和$ \left({x}_{s2},{y}_{s2}\right) $，则在近场条件下，声源$ m $到阵元$ n $的信号传播时间满足：

$ {\tau }_{mn} = \frac{{r}_{mn}}{c} = \frac{||{p}_{n}-{s}_{m}||}{c} = \frac{\sqrt{{\left({x}_{n}-{x}_{sm}\right)}^{2}+{\left({y}_{n}-{y}_{sm}\right)}^{2}}}{c}。$

(3)

对于任意阵元$ n(n\geqslant 2) $，根据式(3)可建立如下2个球面方程：

$ {\left({x}_{n}-{x}_{s1}\right)}^{2}+{\left({y}_{n}-{y}_{s1}\right)}^{2}={\left(c\cdot {\tau }_{1n}\right)}^{2}，$

(4)

$ {\left({x}_{n}-{x}_{s2}\right)}^{2}+{\left({y}_{n}-{y}_{s2}\right)}^{2}={\left(c\cdot {\tau }_{2n}\right)}^{2}。$

(5)

求解这2个二元二次方程组，不难得到阵元$ n $的真实坐标值^[13]，即

$ {{x}_{n}=\frac{{{(x}_{s2}^{2}-{x}_{s1}^{2})+{(y}_{s2}^{2}-{y}_{s1}^{2})+\left(c\cdot {\tau }_{1n}\right)}^{2}-{\left(c\cdot {\tau }_{2n}\right)}^{2}-2{y}_{n}({y}_{s2}-{y}_{s1})}{2\left({x}_{s2}-{x}_{s1}\right)}。} $

(6)

记$ \gamma ={{(x}_{s2}^{2}-{x}_{s1}^{2})+{(y}_{s2}^{2}-{y}_{s1}^{2})+\left(c\cdot {\tau }_{1n}\right)}^{2}-{\left(c\cdot {\tau }_{2n}\right)}^{2} $，将式(6)代入式(4)，经化简后得：

$ {y}_{n}=\frac{-b\pm \sqrt{{b}^{2}-4ac}}{2a}，$

(7)

其中，

$ {\left\{\begin{aligned}&a = 4{\left({y}_{s2}-{y}_{s1}\right)}^{2}+4{({x}_{s2}-{x}_{s1})}^{2}，\\ &b = 8\left( {x}_{s2} - {x}_{s1}\right) \left( {y}_{s2} - {y}_{s1}\right) {x}_{s1} - 8{({x}_{s2} - {x}_{s1})}^{2}{y}_{s1} - 4\gamma ({y}_{s2} - {y}_{s1})，\\ & c = 4{\left( {x}_{s2} - {x}_{s1}\right) }^{2}\left( {x}_{s1}^{2} + {y}_{s1}^{2} - {\left( c\cdot {\tau }_{1n}\right) }^{2}\right) - 4{x}_{s1}\left( {x}_{s2} - {x}_{s1}\right) \gamma + {\gamma }^{2}。\end{aligned}\right. }$

(8)

式(8)的求解计算非常简单，这里$ {y}_{n} $存在双解问题，根据工程实际中辅助声源的布设位置进行取舍。由式(6)和式(7)可看出，在阵列位置校准中，需要通过估计多路时延$ {\tau }_{nm} $来反推出阵元相对位置。考虑到本文研究中的辅助声源和掩埋式光纤水听器阵列之间没有建立严格的时间同步机制，难以采用传统的同步钟测距方法直接求取时延值，这里选用间接时延估计方法予以解决。

如前所述，选取阵元1作为参考阵元，则辅助源到阵元1时延表示为：

$ {\tau }_{m1} = \frac{{r}_{m1}}{c} = \frac{||{p}_{1}-{s}_{m}||}{c}=\frac{\sqrt{{\left({x}_{1}-{x}_{sm}\right)}^{2}+{\left({y}_{1}-{y}_{sm}\right)}^{2}}}{c}。$

(9)

对应声源$ m $发声工况，令每个阵元$ n(n\geqslant 2) $相对于参考阵元1的相对时延为$ {\tau }_{n1}^{m} $，则有：

$ {{\tau }_{mn}=\tau }_{n1}^{m}+{\tau }_{m1}。$

(10)

由此可见，只要估计出每个阵元$ n $与参考阵元1的相对时延，即可解算出阵列的几何位置。

求解上述时延估计问题的一种经典方法是采用广义互相关法^[14]，基本思想是利用两接收信号的广义相关函数来估计时间延迟，并根据不同最优准则，对输入信号进行加权处理，从而提高时延估计性能。

根据式(11)，广义相关函数$ {\widehat{\psi }}_{{x}_{1}{x}_{2}}\left(\tau \right) $在$ \tau $取值范围内取最大值的位置对应于最优时延估计值，即：

$ {\widehat{\tau }}_{GCC}=\underset{\tau }{\mathrm{argmax}}{\widehat{\psi }}_{{x}_{1}{x}_{2}}\left(\tau \right)。$

(11)

将式(11)所获得的时延估计结果分别代入式(6)和式(7)，可以得到阵列的初始位置，并作为后续寻优的种群初始化位置向量。

2 基于遗传算法和MUSIC的阵列位置寻优方法 2.1 遗传算法的基本原理

遗传算法（Genetic Algorithm, GA）是一种基于自然选择和遗传机制的全局优化方法。遗传算法具有自适应性强、全局搜索能力强、可处理非线性和复杂优化问题等特点，其核心步骤包括^[15]：初始化种群、适应度函数、选择操作、交叉操作、变异操作及迭代与终止条件。通过重复执行上述操作，直到达到最大迭代次数或适应度函数收敛终止条件，寻优结束后输出当前最优解。标准遗传算法采用固定的参数设置，这极大地限制和束缚了算法的性能。为搜索到最优解，这里引入自适应概念实现动态自适应进化过程，基本思想是^[16]：首先计算适应度值，采用轮盘赌法进行选择操作，并保存个体的适应度值，按照适应度分配交叉概率和变异概率，进行交叉变异操作。

2.2 基于MUSIC算法的适应度函数构造

遗传算法中的适应度函数是判别阵列位置布局优化好坏的关键标准。在优化阵元位置的过程中，适应度函数决定了进化的方向，为了优化得到全局最优的阵列位置，这里通过在阵列的远场引入第3个位置精确已知的辅助声源发射单频窄带信号，并选用工程应用较为广泛的MUSIC方位估计算法构造适应度函数^[17]。MUSIC算法对单频窄带信号具有良好的估计性能，尤其适用于窄带远场信号的方位估计。

对于一个远场方位估计问题，阵列流形向量中第$ n $个元素可以表示为：

$ {a\left( {{\theta _d},{\delta _n}} \right) = {\text{exp}}\left\{ {{{j}}\dfrac{{2{\text{π}} }}{\lambda }\left[ {\left( x_n^{{\text{est}}} + \Delta {x_n}){\text{sin}}{\theta _d} + (y_n^{{\text{est}}} + \Delta {y_n}){\text{cos}}{\theta _d} \right) } \right]} \right\}。}$

(12)

式中：$ {\theta }_{d} $为远场辅助源的真实方位；$ {x}^{\rm{est}} $和$ {y}^{\rm{est}} $分别为阵元位置初始坐标值；$ \mathrm{\Delta }{x}_{i} $和$ \mathrm{\Delta }{y}_{i} $分别为阵列扰动误差，取决于遗传算法中变量的取值范围。

不难得到，N个阵元接收到的第k次快拍数据向量为：

$ \mathit{Y}\left(k\right)=\mathit{a}\left({\theta }_{d},\mathit{\delta }\right)\mathit{S}\left(k\right)+\mathit{N}\left(k\right)。$

(13)

式中：$ \mathit{a}\left({\theta }_{d},\mathit{\delta }\right) $为$ N\times D $维阵列流形向量；$ \mathit{S}\left(k\right) $为$ D\times 1 $维信号向量；$ \mathit{N}\left(k\right) $分别为$ N\times 1 $维噪声向量；这里$ D=1 $。一般假设信号和噪声不相关，噪声为零均值、方差为$ {\sigma }^{2} $的加性高斯白噪声。阵列协方差矩阵为：

$ {\boldsymbol{R}}_{y}=\mathit{E}\left[\left[\mathit{Y}\left(k\right){\mathit{Y}}^{H}\left(k\right)\right]\right]=\mathit{a}\left({\theta }_{d},\mathit{\delta }\right){\mathit{R}}_{\mathit{s}}{\mathit{a}}^{H}\left({\theta }_{d},\mathit{\delta }\right)+{\sigma }^{2}\mathit{I}。$

(14)

式中：$ {\boldsymbol{R}}_{\mathit{s}}=\mathit{E}\left\{\mathit{S}{\mathit{S}}^{H}\right\} $为信号协方差矩阵。在实际应用中，阵列协方差矩阵可以采用一段快拍样本估计，即样本协方差矩阵为：

$ {\widehat{\boldsymbol{R}}}_{y}=\frac{1}{K}\sum _{k=1}^{K}\left[\mathit{Y}\left(k\right){\mathit{Y}}^{H}\left(k\right)\right]。$

(15)

式中：$ K $为快拍数。

通过对协方差矩阵$ {\widehat{\boldsymbol{R}}}_{y} $进行特征分解，分别得到由较大的$ D $个特征值与较小的$ N-D $个特征值对应的特征向量组成的信号子空间$ {\mathit{U}}_{s} $与噪声子空间$ {\mathit{U}}_{n} $，则基于MUSIC方位谱构造目标函数为:

$ {P}_{\rm music}\left({\theta }_{d},\mathit{\delta }\right)=\underset{{\theta }_{d}}{\mathbf{max}}\frac{1}{\left|{\mathit{a}}^{H}\left({\theta }_{d},\mathit{\delta }\right){\mathit{U}}_{n}{\mathit{U}}_{n}^{H}\mathit{a}\left({\theta }_{d},\mathit{\delta }\right)\right|}。$

(16)

方位估计的目的是最大化MUSIC方位谱，这等同于找出目标函数分母的最小值，即

$ {E}_{\rm music}\left({\theta }_{d},\mathit{\delta }\right) = \underset{{\theta }_{d}}{\mathrm{min}}\left|{\mathit{a}}^{H}\left({\theta }_{d},\mathit{\delta }\right){\mathit{U}}_{n}{\mathit{U}}_{n}^{H}\mathit{a}\left({\theta }_{d},\mathit{\delta }\right)\right|。$

(17)

这里$ {E}_{\rm music}\left({\theta }_{d},\mathit{\delta }\right) $作为个体的适应度函数。

为防止阵元位置超出物理约束范围，进一步采用最小化阵元间距误差约束来限制适应度函数无效解的产生。根据优化的阵元位置坐标$ ({x}_{n}^{\rm opt},{y}_{n}^{\rm opt}) $估计相邻阵元间距离，表示为：

$ {d}_{ij}^{\rm opt}=\sqrt{{({x}_{i}^{\rm opt}-{x}_{j}^{\rm opt})}^{2}+{({y}_{i}^{\rm opt}-{y}_{j}^{\rm opt})}^{2}} 。$

(18)

定义阵元间距误差适应度函数$ {E}_{d} $为：

$ {E}_{d}=\sum _{i=1}^{N-1}({d}^{\rm opt}-d{)}^{2}。$

(19)

最终适应度函数可以构造为：

$ Fit\left(\mathit{\delta }\right)={K}_{1}{E}_{\rm music}+{K}_{2}{E}_{d}。$

(20)

式中：$ {K}_{1} $和$ {K}_{2} $是权重参数，根据优化目标进行适当调整。上述优化过程可以选用Matlab自带的函数求解^[15]，适应度值越小，个体越优。随着迭代的进行，种群中的个体逐渐向最优解靠近，适应度值也会逐渐趋于稳定。最后，利用平均误差和均方根误差来评价阵列位置校准方法的性能。

2.3 阵列位置寻优实施步骤

基于遗传算法和MUSIC的阵列位置寻优方法的具体实施步骤为：

1）利用位于阵列近场位置已知、分时工作的双辅助源发射宽频带信号，分别采集阵列输出信号并作时延估计分析，利用式(6)和式(7)获取阵列的初始几何位置；

2）利用位于阵列远场位置已知的单一辅助源发射单频窄带信号，采集阵列输出信号，利用式(17)～式(20)构造适应度函数；

3）设置遗传算法的种群规模、变量数量、变量上下界、交叉率、变异率及迭代次数等相关参数，并进行阵元位置迭代寻优；

4）如达到最大迭代次数或适应度收敛，则将寻优后结果代入式(12)和式(16)中，并对比分析不同阵列位置校准阶段的方位谱，验证校准方法是否有效。

3 仿真分析与海试验证 3.1 仿真实验与结果分析

考虑二维平面内有一条20元均匀线列阵，阵元间距为2 m，以阵元1为参考阵元。结合前期浅海光纤水听器阵列埋设施工经验，这里假设阵列位置在XY平面内的随机扰动服从高斯分布，且误差最大值不超过$ \lambda /10 $。设近场分时工作的双辅助声源位置分别为（5，20）和（20，25），周期性发射500～2000 Hz宽频带线性调频信号，信噪比为10 dB；远场分时工作的单一辅助声源位置为（2 000，2 000），周期性发射800 Hz单频窄带信号, 信噪比为5 dB；信号采集时长为1 s，采样频率为62500 Hz。图2分别给出了阵列的几何布局和近远场辅助源的时域波形及频谱图。

由图2(a)可知，理想阵元位置均匀排列在一条直线上，真实阵元位置由于存在埋设施工误差，与理想位置存在随机偏移。设置遗传优化算法中群体规模为150个个体，变量数量为40，变量的上下界为$ [-0.1,0.1] $，交叉概率为0.8，迭代次数为500，权重参数$ {K}_{1} $和$ {K}_{2} $分别取1和0.01，其余参数选用默认值。仿真结果如图3所示。

与真实阵列位置相比，采用近场双辅助声源校准方法估计阵列的初始位置，其平均误差为0.016 m，均方根误差为0.022 m；引入远场单一辅助声源，采用基于遗传算法和MUSIC的阵列位置寻优后，得到的阵列最终位置平均误差为0.006 m，均方根误差为0.008 m。由图3可知，校准后的阵元位置更接近真实阵元位置。

由图4可知，算法在进化到150代时就开始收敛于真值附近。

由图5可知，通过基于遗传算法和MUSIC的阵列位置校准后，得到的MUSIC方位谱的谱峰更逼近真实位置输出谱峰，表明文中提出的寻优方法可在相对较大的搜索空间内实现全局优化，得到更为准确的阵列位置。

3.2 海试验证与结果分析

为了验证阵列位置校准方法的有效性和工程适用性，以大连海域临岸海底埋设的一条20元、2 m间距的光纤水听器阵列为平台，开展了阵列位置校准试验。试验以船载模拟声源为辅助声源，根据光纤阵列的埋设位置坐标预先设置就位点，声源吊放到水下5 m，利用差分北斗/GPS定位系统获取辅助声源的绝对位置坐标，并设计了试验船机动工况来验证校准后阵列对移动目标的跟踪能力，其时域波形及频谱图如图6所示。

经试验后复盘，近场分时工作的2个辅助声源就位点坐标分别为（E121°40.3155', N38°50.119'）和(E121°40.359', N38°50.101')，发射500～2000Hz宽频带线性调频信号，远场分时工作的单一辅助声源就位点坐标为(E121°42.338', N38°49.558')，周期性发射600 Hz单频窄带信号，信号时长截取1 s，采样频率为62500Hz。

这里设置遗传优化算法中相关参数与3.1节仿真实验保持一致。校准结果如图7所示。从图7(a)可知，阵列埋设过程由于受海流海浪、敷设设备和工艺限制等因素影响，会导致阵列埋设位置与理想预期存在一定偏差。从图7(b)可知，通过引入阵元间距误差适应度函数约束，收敛后的真实阵元位置关系满足相邻阵元间距不大于2 m的工程实际情况。

由图8可知，试验船从距离光纤水听器阵列较远（3.5 km左右）处以4 kn航速由远及近机动，尽管阵列的接收信噪比较低，但通过对400 Hz附近的一根窄带线谱进行跟踪处理，仍可以获取较为准确的目标方位信息，且与试验船安装的差分北斗/GPS定位系统记录的位置轨迹趋势基本一致，说明本文提出的算法具有一定的环境适应能力，较好地满足了工程实际需求。

4 结　语

本文提出了一种融合时延估计和遗传优化算法的三辅助源阵列位置校准方法，分别采用位于阵列近场和远场位置已知、分时工作的辅助源发射宽频带信号和单频窄带信号，利用经典时延估计获取阵列位置初始值，利用基于遗传算法和MUSIC算法的阵列位置寻优方法获取阵列的优化位置。仿真和实验结果表明，本文提出的方法较好克服了浅海环境单一阵列远场校准面临的水声信道强多途干扰及远距离声传播衰减导致的接收信噪比不足等问题，具有良好的校准精度和一定的环境适应性。实际海上校准试验中，3个辅助源完全可以利用船拖曳单一宽频带模拟声源替代，通过机动到预先设定的校准位置、装订不同形式的宽带/窄带信号实现校准实验，在工程实施方面具有较高的便利性，展现出良好的应用前景。