0 引　言

管道作为一种高效的输送手段，在船舶系统中占据至关重要的地位，船舶管道的正常运行对于船舶的航行安全和操作效率至关重要。随着管道的日益老化和潜在灾害风险的存在，加之管理不善，管道破损和泄漏事件频发，不仅造成资源浪费，还可能危及公共安全^[1,2]。因此，开展及时而有效的管道泄漏检测和定位对于确保船舶的安全运行至关重要。目前的研究主要集中在基于管道压力/流量数据的泄漏检测方法，以及利用光纤和声发射传感技术的方法^[3]。在船舶的实际应用中，基于声学技术的泄漏检测和定位方法的优势突出，因为船舶管道的工作环境恶劣且背景噪声较大，而基于声学技术的泄漏检测方法能够在这种环境下提供高敏感度和准确性的检测，因此该技术已成为管道泄漏定位的主流方法之一。

声学泄漏检测定位技术可以根据传感器的布设位置分为管内和管外检测定位技术。管内检测技术涉及将传感器布置在管道内部，直接接触管内介质，通常需要暂停管道服务以实施检测，且成本较高^[4]。管外检测技术可进一步分为接触式和非接触式方法，如将传感器安装在阀门或管道的暴露点，或者使用传声器阵列从远处捕捉泄漏声信号^[5]。接触式方法在有可接入点的情况下，可以提供较好的检测和定位效果，但受限于检测距离，并且船舶内部空间有限同时管道分布复杂，导致人工在船舶环境下难以直接接触到管道，这使得接触式方法在实际的船舶管道泄漏检测中存在较大挑战。而基于传声器阵列的非接触式管道泄漏定位技术，由于其高精度和非侵入性的特点，可以在不影响船舶正常运行的情况下有效检测管道泄漏，从而提高船舶运行的安全性和可靠性，在此应用领域展现出重要且前景广阔的应用潜力^[6]。

基于传声器阵列的泄漏定位技术主要分为3类方法：基于最大输出功率的波束形成技术，基于高分辨率谱估计的定位技术和基于时延估计的定位技术^[7]。后者因其较低的计算复杂性和高准确率而被广泛应用于定位任务。此技术基于声源信号到达不同传感器的时间差的估计，主要包括广义互相关（Generalized Cross Correlation，GCC）和自适应最小均方（Least Mean Square，LMS）时延算法^[8]。随后，结合传声器阵列的几何拓扑特征，设计相应的定位算法。郑晓亮等^[9]提出通过延迟求和法对附加传感器信号的时延进行估计，并通过构建权向量对线性阵列的结果进行加权，有效降低了互相关定位方法的平均误差。在此基础上，郑晓亮等^[10]考虑了基于线性阵列的两步定位方法，并通过使用六元阵列改进了定位的稳定性和抗干扰能力。张梅等^[11]通过运用四点传感器阵列结合小波去噪技术，在广义互相关的基础上计算时延，实现了误差小于1 m的定位精度。王强等^[12]通过使用十字形四元传感器阵列并考虑了阵列布局和泄漏位置对时延估计以及定位结果的影响，发现通过增加阵列的孔径和布设间距可以显著减小误差。由此可见，利用传感器阵列进行管道泄漏的定位具有明确的应用潜力。然而，该方法的效果主要受限于时延估计的精度，因此，当前研究的重点集中在提高时延估计精度上。

为了解决现存的管道泄漏定位难题，本研究提出了基于传感器阵列技术与变分模态分解（Variational Mode Decomposition，VMD）的管道泄漏定位技术。船舶的实际运行环境中存在复杂的噪声干扰，而这种定位技术能够显著提升管道泄漏检测的效率和准确性，具体来说，首先利用传感器阵列捕获管道泄漏产生的声信号，随后应用VMD技术对收集的信号进行多重分解，提取出关键的泄漏信号，同时筛除背景噪声等干扰，以增加时延估计的精度。同时，采用修正的权函数来处理信号之间的互功率谱密度函数，以进一步提升GCC的时延估计精度。最后，结合传感器阵列构型建立定位算法，以估计泄漏源的位置，从而实现有效且精确的管道泄漏定位。

1 VMD及广义互相关原理

图1为提出的管道泄漏定位算法的具体流程。

1.1 VMD

VMD是一种用于多分量信号分解的算法，由Dragomiretskiy和Zosso提出。它旨在将一个复杂信号分解成一系列有限带宽，具有稀疏特性的固有模态函数（Intrinsic Mode Functions，IMFs），每个模态对应于信号中的一个基本频率成分^[13]。VMD具有较为完备的数学理论基础，本质上是一个自适应最优Weiner滤波器组，由于其优良的特性，已在较多领域取得了广泛的应用。

VMD的目标函数的构造分为3个步骤:

步骤1　对每个IMF进行Hilbert变换得到单边谱;

步骤2　利用指数函数调制每个IMF的中心频带转移到对应的基带；

步骤3　利用高斯平滑度和梯度平方准则，估算每个IMF的信号带宽。

基于此，可以构造相应的约束变分模型

$\left\{ \begin{array}{*{20}{l}}\min_{\left\{u_k\right\},\left\{\omega_k\right\}}\left\{ \sum\limits_k^{ } \left\| {\partial_t} \left[ \left( \delta(t) + \frac{\mathrm{i}}{\text{π}t} \right) * u_k(t)\right]\exp\left( -\text{i}\omega_kt\right) \right\| _2^2 \right\}。\\ \text{ s}\text{.t }\sum\limits_k^{ }u_k(t)=f(t)\end{array}\right.$

(1)

式中：${u_k}\left( t \right)$为第$k$固有模态函数；${\omega _k}$为该模态的中心频率；${\partial_t}$为函数对时间的偏导运算；$*$为卷积运算；i为虚数单位；$\delta \left( t \right)$为单位脉冲函数；$f\left( t \right)$为原信号。

引入惩罚因子$\alpha$和Lagrange乘子$\lambda$，将VMD的约束变分问题转化为非约束变分问题，可得：

$\begin{split} & L\left( {\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\},\lambda } \right)= \\& \alpha \sum\limits_k {\left\| {{\partial _{\text{t}}}\left[ {\left( {\delta (t) + \frac{{\text{i}}}{{\pi t}}} \right)*{u_k}(t)} \right]\exp \left( { - {\text{i}}{\omega _k}t} \right)} \right\|_2^2} + \\&\left\| {f(t) - \sum\limits_k {{u_k}} (t)} \right\|_2^2 + \left\langle {\lambda (t),f(t) - \sum\limits_k {{u_k}} (t)} \right\rangle 。\end{split}$

(2)

采用交替方向乘子法（Alternate Direction Method of Multiplers，ADMM）来解决上述非约束变分问题^[14]，通过不断迭更新${u_k}\left( t \right)$${\omega _k}$和$\lambda$。

因此，VMD将信号$f\left( t \right)$分解为$K$个IMF分量的步骤如下：

步骤1　初始化$\left\{ {\hat u_k^1} \right\}$，$\left\{ {\hat \omega _k^1} \right\}$和${\hat \lambda ^1}$，$n = 0$；

步骤2　迭代更新${u_k}$和${\omega _k}$

$\hat{u}_k^{n+1}(\omega)=\frac{\hat{f}(\omega)-\sum\limits_{i=1}^{k-1}\hat{u}_i^{n+1}(\omega)-\sum\limits_{i=k+1}^K\hat{u}_i^n(\omega)+\hat{\lambda}(\omega)/2}{1+2\alpha\left(\omega-\omega_k^n\right)^2}，$

(3)

$\omega_k^{n+1}=\frac{\int_0^{\infty}\omega\left|\hat{u}_k^{n+1}(\omega)\right|^2\text{d}\omega}{\int_0^{\infty}\left|\hat{u}_k^{n+1}(\omega)\right|^2\text{d}\omega}。$

(4)

式中：$\hat f(\omega )$、$\hat \lambda (\omega )$分别为$f\left( t \right)$和$\lambda$的傅里叶变换。

步骤3更新迭代$\lambda$

$\hat{\lambda}^{n+1}(\omega)=\hat{\lambda}^n(\omega)+\gamma\left(\hat{f}(\omega)-\sum\limits_{k=1}^K\hat{u}_k^{n+1}(\omega)\right)。$

(5)

式中：$\gamma$为对噪声的容许参数；

步骤4　重复步骤2和步骤3，直到满足终止条件

$\sum\limits_{k=1}^K\left( \left\| \hat{u}_k^{n+1}(\omega)-\hat{u}_k^n(\omega) \right\| _2^2/ \left\| \hat{u}_k^n(\omega) \right\| _2^2\right) < \varepsilon。$

(6)

式中：$\varepsilon$为预设精度。

在采用VMD对采集的信号$f\left( t \right)$进行分解后，需要考虑选择含有效信息的IMF分量以重构信号，并减少噪声对时延估计的影响。考虑到有效的IMF分量理应包含原始信号的有效信息，并且这些信息与干扰噪声间应为非相关关系，本研究采用Pearson相关系数作为评判标准来选择有效的IMF分量，其计算公式为^[15]：

$\rho_{u_if}=\frac{Cov\left(u_i,f\right)}{\sigma_{u_i}\sigma_f}。$

(7)

式中：$Cov\left( {{u_i},{u_j}} \right)$为IMF分量${u_i}$和信号$f\left( t \right)$的协方差；${\sigma _{{u_i}}}$和${\sigma _f}$分别为${u_i}$和信号$f\left( t \right)$的标准差。

当IMF分量包含原始信号的有效成分时，它应当与原信号具有较高的相似度，从而具有较高的Pearson相关系数值。相反，若IMF分量主要为干扰噪声，则其与原始信号的相关系数将较低，反映出较小的Pearson相关系数值。通过设定一个相关系数阈值$\theta$，当IMF分量与原信号的Pearson相关系数超过此阈值时，将该IMF分量保留；低于阈值的IMF分量则被剔除。遍历所有分量并筛选后，所保留的IMF分量用于重构出去噪后的信号，并用于后续的时延估计。

1.2 广义互相关

互相关是进行时延估计中应用最为广泛的算法之一，该方法主要基于信号间的互相关函数来估计时延。由于传感器阵列所采集的信号源自同一声源，因此不同传感器接收的信号间存在一定的相关性。假设2个传感器所采集到的信号可以表示为^{[16 - 17]}：

$x_1\left(t\right)=\alpha_1s\left(t-\tau_1\right)+n_1\left(t\right)，$

(8)

$x_2\left(t\right)=\alpha_2s\left(t-\tau_2\right)+n_2\left(t\right)。$

(9)

式中：$s\left( t \right)$为声源信号；${\alpha _i}$为信号由声源到传感器的衰减系数；${n_i}\left( t \right)$为环境噪声；${\tau _i}$为声波从声源到第$i$个传声器的时间。

信号${x_1}\left( t \right)$和${x_2}\left( t \right)$的互相关函数${R_{{x_1}{x_2}}}\left( \tau \right)$可表示为：

$R_{x_1x_2}\left(\tau\right)=E\left[x_1\left(t\right)x_2\left(t-\tau\right)\right]。$

(10)

式中：$\tau$为信号${x_1}\left( t \right)$和${x_2}\left( t \right)$间的时延。

将式(8)和式(9)代入式(10)，并考虑到信号$s\left( t \right)$和背景噪声${n_i}\left( t \right)$间，以及背景噪声间的互不相关性，上式可写为：

$R_{x_1x_2}\left(\tau\right)=\alpha_1\alpha_2R_{ss}\left(\tau-\tau_{12}\right)。$

(11)

式中：${\tau _{12}} = {\tau _1} - {\tau _2}$；${R_{ss}}\left( \tau \right)$为声源信号$s\left( t \right)$的自相关函数。

由相关函数的性质可知，当$\tau = {\tau _{12}}$时，${R_{{x_1}{x_2}}}\left( \tau \right)$取得最大值，此时${R_{{x_1}{x_2}}}\left( \tau \right)$所对应的时延$\tau$即是两传感器采集的信号间的时延。

根据维纳辛钦定理，互相关函数与互功率谱间存在如下关系^[18]：

${R_{{x_1}{x_2}}}\left( \tau \right) = \int_{ - \infty }^\infty {{G_{{x_1}{x_2}}}\left( \omega \right){e^{{\text{i}}\omega \tau }}{\text{d}}} \omega 。$

(12)

实际环境中的噪声会影响信号的峰值，从而降低时延估计的准确性。一般可以信号进行频域加权处理，在互功率谱中抑制噪声，以便锐化峰值。通过对加权后的互功率谱进行傅里叶反变换，变可得到一个锐化的广义互相关函数：

${R_{{x_1}{x_2}}}\left( \tau \right) = \int_{ - \infty }^\infty {{\psi _{12}}\left( \omega \right){G_{{x_1}{x_2}}}\left( \omega \right){e^{{\text{i}}\omega \tau }}{\text{d}}} \omega 。$

(13)

式中：${\psi _{12}}\left( \omega \right)$为加权函数。

图2为利用广义互相关估计时延的程序框图^[19]。根据使用的权函数${\psi _{12}}\left( \omega \right)$的不同，对互相关函数的锐化效果也有差异。表1总结了常用的加权函数极其表达式和各自的特性。

实际使用中，PHAT和ML权函数具有良好的效果，且由于PHAT权函数形式的简洁性，取得了较为广泛的应用。然而，PHAT权函数对噪声较为敏感，在信噪比（Signal-to-Noise Ratio，SNR）较高的环境下有较高的准确度。但是，随着SNR的降低，噪声在互功率谱中的比重增加，可能导致伪峰的出现而干扰时延估计结果，进而导致错误的定位结果。因此，考虑引入一个与SNR相关的调节参数，通过预估信号的SNR来调整该参数，动态地控制信号与噪声频谱在权重比例。此外，鉴于信号能量低时，权函数的分母可能趋近于0，导致误差，在权函数分母中加入一个相干因子以稳定计算。据此，修正的权函数形式为：

$\psi_{12}\left(\omega\right)=\frac{1}{\left|G_{x_1x_2}\left(\omega\right)\right|^{\beta}+\left|\gamma\left(\omega\right)\right|^2}。$

(14)

式中：$\beta$为与信号SNR的相关参数，当SNR≤10时，可取0.5；10<SNR≤25时，可取0.75，当25<SNR时，可取0.9。

2 泄漏声源定位方法

在获取传感器间信号的时延之后，下一步是对管道泄漏声源进行精确定位。一般来说，不同的传感器阵列拓扑结构会影响定位的性能。本文以五元十字形阵列作为例，阐述管道泄漏声源定位的位置估计算法，该阵列的具体布局如图3所示。

假定5个传感器的坐标分别为${M_1}\left( {0,0,0} \right)$、${M_2}( - d, 0,0 )$、${M_3}\left( {0,d,0} \right)$、${M_4}\left( {d,0,0} \right)$和${M_5}\left( {0, - d,0} \right)$，副传感器与位于坐标原点的主传感器的距离为$d$。声源的坐标为$S\left( {x,y,z} \right)$，其方位角和俯仰角分别为$\varphi$和$\theta$，声源距主传感器的距离为${r_1}$，距各副传感器的距离分别为${r_2}$、${r_3}$、${r_4}$和${r_5}$。声源信号到达主传感器的时间表示为$t$，而到达不同副传感器的时间与到达主传感器的延迟用${\tau _{ij}}$表示。基于几何关系，可以建立以下方程组：

$\left\{ \begin{gathered} {x^2} + {y^2} + {z^2} = r_1^2 ，\\ {\left( {x + d} \right)^2} + {y^2} + {z^2} = r_2^2，\\ {x^2} + {\left( {y - d} \right)^2} + {z^2} = r_3^2，\\ {\left( {x - d} \right)^2} + {y^2} + {z^2} = r_4^2，\\ {x^2} + {\left( {y + d} \right)^2} + {z^2} = r_5^2。\\ \end{gathered} \right.$

(15)

利用声波的传播速度$c$，将声源距离各个传感器的距离表示为时延，可得：

$\left\{ \begin{gathered} {r_1} = ct，\\ {r_2} = c\left( {t + {\tau _{12}}} \right) ，\\ {r_3} = c\left( {t + {\tau _{13}}} \right)，\\ {r_4} = c\left( {t + {\tau _{14}}} \right) ，\\ {r_5} = c\left( {t + {\tau _{15}}} \right) 。\\ \end{gathered} \right.$

(16)

将声源S的笛卡尔坐标转换为极坐标：

$\left\{ \begin{gathered} x = {r_1}\sin \theta \cos \varphi ，\\ y = {r_1}\sin \theta \sin \varphi ，\\ z = {r_1}\cos \theta 。\\ \end{gathered} \right.$

(17)

联立上式，可以求得声源距主传感器的距离${r_1}$、方位角$\varphi$和俯仰角$\theta$的解析式为：

$r_1=ct=c\frac{\tau_{13}^2+\tau_{15}^2-\tau_{12}^2-\tau_{14}^2}{2\left(\tau_{12}+\tau_{14}-\tau_{13}-\tau_{15}\right)}，$

(18)

$\varphi=\arctan\frac{\left(2t+\tau_{13}+\tau_{15}\right)\left(\tau_{15}-\tau_{13}\right)}{\left(2t+\tau_{12}+\tau_{14}\right)\left(\tau_{12}-\tau_{14}\right)}，$

(19)

$\theta=\arcsin\frac{x\cos\varphi+x\tan\varphi\sin\varphi}{ct}。$

(20)

式中，$x = \left( {r_2^2 - r_4^2} \right)/4d$。

据此，可以利用基于VMD和GCC的技术估计各传感器与主传感器之间的时延结果，来确定声源相对于主传声器的距离、方位角以及俯仰角，并依据式(17)将结果转换至笛卡尔坐标系，以获得声源的定位结果。

3 实验与结果分析

为了验证本研究提出算法的有效性及精度，本文设计了实验以模拟管道泄漏信号，并使用Matlab软件对采集的数据进行分析。实验使用直径为DN40的钢管，通过在管道上开孔的方式模拟泄漏点。管内介质为压缩空气，压缩空气由空压机提供，实验中管道的压力为0.4 MPa，实验现场照片如图4所示。声源定位实验采用了5个传声器组成的阵列，阵列中传感器的坐标分别为(0, 0, 0)、(–0.2 m, 0, 0)、(0, 0.2 m, 0)、(0.2 m, 0, 0)和(0, –0.2 m, 0)。每个通道的采样频率为25.6 kHz，采样时长为10 s。图5为一个传感器采集的管道泄漏声信号的波形。为了测试所提算法在不同SNR条件下的性能，在信号中添加高斯白噪声以改变SNR的大小，使其以5 dB为步长由−10 dB变化至30 dB，并与经典的时延估计算法的结果进行对比分析。

图6为在SNR=−5 dB的条件下，本文所提出的方法与PHAT权函数所得到的时延估计结果的对比，其中VMD的分解层数设为$K = 10$，相关系数阈值设定为0.4。值得注意的是，在较高的SNR下，几种方法均能获得明显的主峰，并能够准确地估计时延。然而，在较低的SNR条件下（见图6），干扰噪声会导致伪峰的出现，真实时延的主峰受到伪峰的显著干扰，这使得GCC-PHAT方法难以区分真实的主峰，从而无法准确估计时延，进一步导致错误的管道泄漏定位结果。相对地，在低SNR条件下，本文提出的VMD-GCC-β-PHAT方法能有效抑制干扰峰值的影响，结果中依然可以清晰地识别出真实时延所对应的主峰，从而获取准确的时延估计。此外，可以观察到本文提出的算法不仅能够锐化互相关函数，提取出真实的时延，并且与GCC-PHAT相比，其对应的互相关系数峰值更大，表明处理后的信号相关性更强，更能显著突出主峰，因此增加了时延估计的精度和鲁棒性。

利用本文提出的算法计算得到的副传感器与主传感器之间的时延结果，并结合第3节描述的定位算法，对管道泄漏声源进行空间定位。如表2所示，可以观察到当传感器阵列与泄漏声源位于不同位置时，本文提出的定位方法均能够获得较为准确的定位结果。其中，最大单轴定位误差为0.033 m，在可接受的误差范围之内。值得一提的是，前述时延估计结果已经证实，基于VMD-GCC-β-PHAT的方法能够获得更精确的时延估计，这也表明基于该时延估计得到的空间定位结果同样具有更高的准确性。因此，本研究提出的结合传声器阵列和VMD的管道泄漏声源定位方法，能够有效提升时延估计的精度，从而增强管道空间位置的估计精度。

4 结　语

为解决船舶管道泄漏定位的难题，本文提出一种结合GCC和VMD技术的改进时延估计算法，旨在实现更为精确的声源定位。该算法通过VMD对采集的信号进行多重分解，随后依据相关性系数筛选出主要的IMF，以此来提取信号的主要成分并去除噪声。通过对GCC-PHAT算法进行了改进以考虑SNR的影响，以提高时延估计的精度。研究以五元十字形传感器阵列为例，阐述了泄漏声源定位算法。实验结果验证了所提算法的鲁棒性和有效性，即便在低SNR条件下，也能实现准确的时延估计。此外，不同位置关系的传感器阵列与声源之间的定位结果进一步表明了本文所提出的管道泄漏声源定位方法的有效性。最后本研究所提出的改进时延估计算法，不仅提高了船舶管道系统的运行安全性，还为船舶维护提供了先进的技术支持，在船舶管道泄漏检测的实际应用中拥有巨大的发展前景。