输气管道泄漏的波达时差交叉定位方法*

王 强 郑晓亮 薛 生 袁宏永 付 明

(1 安徽理工大学电气与信息工程学院 淮南 232001)

(2 安徽理工大学能源与安全学院 淮南 232001)

(3 清华大学合肥公共安全研究院 合肥 230601)

0 引言

随着燃气需求量的逐年提高，由输气管道泄漏带来的潜在风险不可避免地随之增加。因此，对管道泄漏的安全监测已经成为公共安全领域的重要内容。目前常用的大范围监测方法[1-4]可实现管道泄漏位置的粗略定位，但还需要较大范围的搜索定位才能准确找到泄漏位置，整个过程工作量大、耗时长、不利于快速排查修复漏点。泄漏发生后，管内高压气体快速喷出并与泄漏孔摩擦振动产生泄漏声源，声源位置与泄漏位置高度重合，对该声源进行定位即可实现对管道泄漏的精确溯源。实际工况下必须对泄漏位置进行三维定位即同时获取其方位和深度信息，才具有应用价值。

基于传感器阵列的声源定位方法主要包括：(1)高分辨率谱估计法[5-6]，可实现对波达方向的超分辨率估计，计算量较大；(2)波束形成法[7-9]，理论上可对近场声源进行三维定位，但计算量大幅提高，丁浩等[10]提出一种可实现近场声源三维定位的波束形成方法，对2 m 以内的3000 Hz 以上中高频声源定位误差低于10%；(3)波达时差(Time difference of arrival,TDOA)法[11-15]，可实现对近场声源的三维定位且计算量小，杨祥清等[16]提出一种基于球形差值的随机梯度下降算法，实现了对三维空间点声源的精确定位。对管道泄漏声波信号特性的研究表明泄漏声波的能量集中在低频波段[17-18]，低频声波较小的传播衰减以及较强的障碍穿透能力有利于信号的检测与定位。根据近场声源判据r＜2d2/λ(r为声源与阵列距离，d为阵元间距，λ为波长)，由于低频声波波长λ较长，且实际定位中阵元间距d不能过大，因此较远距离的低频声源很难满足近场判据，其更趋近于远场声源。而平面阵列只能观测远场声源的空间方位信息，因此基于平面阵列的定位方法无法实现对远场声源的三维定位。综上所述，3 种方法中TDOA 法虽然能以较低的计算量进行近场声源三维定位，但不适用于对较远距离低频泄漏声源的三维定位。

在雷达定位中，常使用多站测向交叉[19-20]的方法，对不同位置基站所观测信源方向进行交叉以计算信源距离。为了将TDOA法应用到管道泄漏三维定位中，本文提出一种基于TDOA的交叉定位方法，将阵列设置在泄漏区域不同位置并获取两组空间方位信息，对其进行交叉求取空间交点从而完成定位。

1 管道泄漏定位原理

TDOA 法的基本原理是通过延时估计算法估计声波信号到达不同传感器阵元的时差，再结合阵列与声源的几何关系联立方程解出声源坐标。由泄漏源所产生的声波可分为近场球面波和远场平面波，理论上平面阵列可实现对近场声源的三维定位和远场声源的空间方位定位(相当于二维定位)。图1为平面阵列与两种声源的位置关系，当传感器阵列与声源足够近时，声波沿不同方向传播至4 个传感器，设声波到达传感器1～3相对于参考传感器的延时分别为τ1、τ2、τ3。结合延时τ1、τ2、τ3以及传感器与声源的位置关系，建立3 个不相关的方程即可解出近场声源坐标。

图1 平面阵列与声源位置关系示意图Fig.1 Relation of position between plane array and acoustic sources

为提高定位范围，无法始终保持传感器阵列与泄漏位置足够近。一方面，由远场判据r＞2d2/λ可知，将泄漏声源视为远场声源更有利于缩小阵列孔径；另一方面，管道泄漏所产生低频声波的波长更长，这进一步加剧了泄漏声源趋向远场声源的趋势。因此在实际工况下，将泄漏源作为远场声源更具有可操作性。当泄漏源被视为远场声源时，决定延时值的变量变为声源相对于阵列的空间方位而不再是三维坐标。设传感器1～3 与参考传感器的距离分别为r1、r2、r3，两传感器连线与x轴夹角为θ1、θ2、θ3。以传感器1为例，图2为远场平面声波到达传感器1和参考传感器的波程差示意图，θ为声源方位角，φ为仰角。

图2 波程差示意图Fig.2 Illustration of wave path difference

由图2 中几何关系可知，远场声源到达传感器1～3与参考传感器的波程差可表示为

延时τ1、τ2、τ3由算法估计得到，结合声速c可得3个包含变量θ、φ的独立方程：

求解式(2)即可得出空间方位解。为进一步实现对远场泄漏声源的三维定位，本文提出一种基于TDOA法的交叉定位方法，其步骤为：(1)将传感器阵列先后设置在泄漏区域的两个不同位置；(2)通过TDOA 法进行泄漏声源空间定向；(3)结合两组空间方位，使用交叉定位法获取泄漏声源的三维坐标，完成泄漏定位。

2 延时估计

空间定向步骤的关键是延时估计。广义互相关法从信号的相关分析基础上发展而来，首先计算两个传感器所接收相关信号的互功率谱，对互功率谱进行加权处理后经傅里叶反变换即可得到两个信号的互相关函数。设不同位置的两个传感器采集到的相关信号分别为x1(t)和x2(t)，求两个信号的互功率谱得

其中，W(ω)表示加权函数，取1 时不进行任何加权，此时R12(τ)表示基本互相关法的傅里叶变换形式。目前常用的加权函数有PHAT、Roth、SCOT、Echart、ML 等[21]，其中Echart 加权函数需要噪声的先验知识，因此本文只对基本互相关法和PHAT、Roth、SCOT、ML 加权广义互相关法进行时延估计效果对比。

在实验室环境下(无强噪声干扰)模拟管道泄漏，所使用管道为DN50 钢管，内压0.6 MPa，泄漏孔形状为圆形，孔径1.5 mm。分别在距离泄漏孔0.4 m 和2 m 处设置声波传感器，图3为时域连续信号的波形和频谱图。相较于距离泄漏孔0.4 m 处信号，2 m 处时域信号衰减幅度达到一个数量级，但500 Hz以下频段衰减较慢。考虑到所使用声波传感器频响范围下限为8 Hz，而高频声波衰减较快且穿透障碍物的能力较差，不利于提高传感器阵列有效定位范围，因此本文仅截取信号的10～500 Hz分量进行分析和定位。

使用基本互相关法和多种加权广义互相关法对两个传感器信号的10～500 Hz 分量进行延时估计，保持室内温度为23°C，管道内压和泄漏孔规格不变，改变信号到达两传感器的距离差，每组距离进行5 次实验。由表1 估计结果可得，PHAT、Roth、SCOT 以及ML加权广义互相关法的估计结果基本保持为0，基本互相关法的估计结果则随着距离的增大而增大。使用15 组数据的基本互相关分析结果计算声速，得到声速估计均值为349.5 m/s，与一个标准大气压和15°C 条件下的340 m/s 空气介质声速典型值较为接近。分析本文条件下加权广义互相关法估计效果不佳的原因为：本文采用低通滤波截取信号10～500 Hz 分量，导致信号互功率谱G12(ω)高频部分能量大幅降低。傅里叶反变换实际是对频域信号再次进行傅里叶变换，根据式(4)对G12(ω)进行傅里叶反变换后，其高频部分(可视为幅度接近0 的时域直流信号)本不应该被反映在反变换后的互相关系数上(可视为G12(ω)的频谱)。而多种加权函数[21]的分母均包含信号的互功率谱或自功率谱，对G12(ω)的高频部分起到了放大作用，则反变换后G12(ω)的频谱出现了直流分量，即对应互相关系数在0 时刻出现了峰值。为保证延时估计精度，采样率必须足够高，同时为避免高频信号干扰，有必要对信号进行低通滤波。因此加权广义互相关法不适合本文高采样率低频信号的延时估计，故采用基本互相关法进行延时估计并完成定位。

图3 信号波形和频谱Fig.3 Wave form and spectrum of signal

表1 多种互相关法延时估计结果Table 1 Results of delay estimation using multiple methods

3 交叉定位

交叉定位是本文所提出定位方法的关键。由于TDOA空间定向必然存在误差，空间中的两条直线几乎不可能相交。为解决不相交直线的交点求取问题，将两条空间直线投影到多个平面上并在平面上求取两条直线的交点，再将多个平面交点聚焦到空间中某一点上，从而完成伪交点的求取。图4为基于投影法的空间直线伪交点求取过程原理图。泄漏点g位于xOy面下方，坐标为(x0,y0,z0)。阵列p1、p2位于xOy平面，坐标分别为(x1,y1,0)、(x2,y2,0)。利用TDOA法可获得泄漏点g相对于阵列p1、p2的两组空间方位角(θ1,φ1)、(θ2,φ2)，并根据该方位角形成两条空间直线l1、l2。为方便计算，将φ1、φ2定义为直线l1、l2与z轴负方向的夹角，θ1、θ2仍为直线l1、l2在xOy面投影与x轴正方向夹角。由于误差的存在，图4 中直线l1、l2均不与泄漏点g相交。将直线l1、l2投影到xOz面得交点s1，投影到xOy面得交点s2，过交点s1做平行于xOy面的平面C，最后将s2投影到平面C即为所求伪交点s。

图4 投影法原理图Fig.4 Schematic diagram of projection method

由p1(x1,y1,0)、p2(x2,y2,0)以及(θ1,φ1)、(θ2,φ2)得直线l1、l2：

分别将直线l1、l2投影到xOz面和xOy面得直线：

由式(7)即可得经投影后再聚焦的伪交点坐标s(x,y,z)。

4 管道泄漏定位实验

4.1 实验平台

根据本文所提出的基于TDOA 法的管道泄漏交叉定位方法，如图5所示搭建定位实验平台。该实验平台包含以下部分：

(1)气体压缩装置：包含空气压缩机、储气罐、减压阀以及球阀若干，当储气罐充满并关闭空压机时，气体压缩装置可向模拟测试管道提供最高1 MPa的持续稳定压力输出；

(2)泄漏管道：包含一根长度为6 m 的DN50镀锌钢管，预留圆形泄漏孔直径为1.5 mm；

(3)定位阵列：包含4 个声波传感器和信号放大器、采集仪和配套分析软件以及传感器阵列支架。

图5 定位实验平台Fig.5 Experimental platform for determining leakage localization

使用上述实验平台进行定位实验，调节减压阀保持管道内压为0.6 MPa，关闭空压机以避免其产生干扰噪声；采集仪采样率为3 MHz，采样时间为10 s；4 个声波传感器按十字交叉排列，图6为传感器阵列的实物图和示意图。由图3所示信号时域波形可知泄漏信号为持续稳定信号，保持管内压力、管道规格和泄漏孔规格等条件不变，则两次定位的信源条件可视为不变。传感器阵列采用支架固定，实验室地面平坦且支架高度固定，改变位置后可保证阵列高度、阵元相对位置等条件基本不变。按图4所示以阵列位置连线为y轴建立空间直角坐标系，保持阵列平面与xOz面重合(即阵列平面保持竖直)。预先设定阵列两次布放位置的坐标，计算阵列与泄漏点的距离和相对角度，并按照此相对位置布放阵列。因此，两次采样的实验条件、阵列安装、阵列相对位置等参数得以保证。综上所述，实验步骤为：储气罐充满气后关闭空压机，调节减压阀为管道提供0.6 MPa 压力输出；使用传感器阵列先后在两个不同位置对泄漏信号进行采集；导出不同位置所采集信号，使用本文定位算法完成定位。

图6 传感器阵列Fig.6 Sensor array

4.2 声源模型对比

阵列第一次布放的位置p1坐标为(1.90 m,0,0)，第二次布放位置p2坐标为(-1.90 m,0,0)，泄漏点g的坐标为(0,3.10 m,-0.84 m)，则声源与两阵列的距离r均为3.73 m。根据本文第2 节多次延时估计结果得到的声速平均值349.5 m/s，以最高频率500 Hz 计算所截取信号最小波长λ为787.8 mm。如图6(b)阵列示意图所示可得传感器1和传感器4 的间距(最大阵元间距)d为54.4 cm，由判据r＞2d2/λ计算得到远场声源与阵列之间距离的临界值为0.75 m，因此本文实验条件下泄漏声源满足远场声源判据。根据本文第3 节对空间方位角(θ,φ)的定义，计算得泄漏点g相对于布放位置p1、p2的空间方位角坐标为(121.5°,77.07°)和(58.5°,77.07°)。将传感器1 设为参考传感器，传感器2、3、4 相对于参考传感器延时τ2-1、τ3-1、τ4-1的估计结果如表2所示。根据阵列中各传感器和泄漏点的位置，计算出远近场声源的理论延时值并在表2 中一同给出。对比可知，实验条件下的泄漏声源更趋近于远场声源。如图2所示以阵列平面为xOy面建立空间直角坐标系，根据图6(b)所示传感器位置关系，将延时估计结果带入式(2)得到两组空间方位坐标。再将空间方位计算结果转换至图4所示坐标系，得到转换后的坐标为(122.0°,78.1°)和(58.1°,77.0°)，与理论值较为接近。将阵列位置坐标以及方位角带入式(7)得伪交点s的坐标为(-0.004 m,3.047 m,-0.791 m)，相较于与泄漏点g(0,3.10 m,-0.84 m)的误差为0.07 m。将延时估计结果带入近场声源模型进行空间定位，此时无法得出有效解，进一步说明实验条件下泄漏源为远场声源。

表2 延时估计值与理论值对比Table 2 Comparison between estimated and theoretical values of time delay

4.3 定位实验

上述定位实验将传感器阵列对称放置在泄漏点两侧，但实际泄漏点位置未知，无法保证传感器对称放置。为模拟实际场景，仍以阵列两次布放的位置连线作为y轴、连线中点作为原点建立空间直角坐标系，但泄漏点分布放置。本文仅考虑阵列前向泄漏点的定位，因此实际泄漏点的分布范围应为xOz面前半空间。考虑到管道多以直线铺设，泄漏点高度变化有限，并且xOz面前半空间关于yOz面对称、阵列位置关于原点对称，讨论yOz面右半空间的定位情况即可推广至左半空间。因此，定位实验划定泄漏点区域为x、y坐标取[0,10 m]范围、z坐标取-0.84 m。图7为阵列位置与泄漏位置的俯视图，阵列两次布放的位置分别为(M1,N1)、(M2,N2)和(M3,N3)，A1～C3为9个泄漏位置(每次实验仅一个位置发生泄漏)，各点x和y坐标如图7所示，阵列位置z坐标为0，泄漏位置z坐标为-0.84 m。为了探究阵列孔径对定位精度的影响，还需要改变阵列模型孔径参数进行实验。记图6(b)所示阵列模型孔径为孔径1，进一步缩小孔径使传感器1、3 与十字交点的距离为35 cm，传感器2、4与交点距离为15 cm，记为孔径2；传感器1、3 与交点的距离为20 cm，传感器2、4 与交点距离为15 cm，记为孔径3。

图7 阵列与泄漏位置俯视图Fig.7 Position of arrays and leakage points

使用3 种孔径的阵列，在3 组阵列布放位置下分别对9 个泄漏位置进行定位实验，每个条件下进行5 次实验并计算平均定位误差。表3为定位结果，对于最大的孔径1，阵列布放位置为M1、N1(即阵列间距2 m)时，A 组泄漏点定位精度相对较高，B、C两组明显下降；阵列布放位置为M2、N2(即阵列间距6 m)时，整体定位误差有所降低；阵列布放位置为M3、N3(即阵列间距10 m)时，A、B 两组泄漏点定位误差明显降低，C 组有所下降但仍较高。同样条件下，A 组定位精度高于B、C 两组，A 组内距离原点最远的A1点定位误差最大。随着阵元间距的减小，定位误差呈增大趋势，且孔径3 下误差明显增大。分析原因为

(1)A 组泄漏点位于y轴上而B、C 两组偏离y轴，则A 组各点相对于两个阵列位置的方位角θ分别属于(0,90°)和(90°,180°)范围且关于θ=90°对称，B、C 两组各点两个方位角θ均小于90°且距离原点越远两方位角大小越接近，由式(7)中x、y坐标表达式可知计算A 组各点坐标时分母tanθ1-tanθ2的绝对值更大，此时分子值的上下浮动所引起的误差更小，反之计算B、C两组各点坐标时分母绝对值较小，则分子波动对误差影响增大；

(2)对于A 组各点，距离原点越远或阵列间距越小则方位角越接近90°，由正切函数变化率可知此时方位角θ的误差对定位误差的影响更大，因此表3 中同一条件下A1点的定位误差相对于A2、A3更大，阵列布放位置为M1、N1(阵列间距最小)时A组定位误差同样更大；

(3)对于B、C两组，阵列布放位置为M1、N1(阵列间距最小)时，各点两个方位角θ的差值较小，由原因(1)可知计算坐标时分子波动对误差影响较大，扩大阵列间距则误差得以降低；

(4)阵列孔径决定空间方位计算的分别率，孔径较小时信源方向的改变所引起的延时变化较小，不利于互相关法对延时进行准确估计，从而导致定位误差增大。

对比本文方法与现有基于TDOA 的声源空间定位方法。杨祥清等[16]使用SI-LMS算法在办公室环境进行语音信号定位，定位结果表明：SI-LMS 算法对距离阵列1 m 以内的近场声源定位误差小于0.05 m，对2～3 m 的中场声源定位误差小于0.1 m，对4 m 以上的远场声源定位误差则显著提高，最大误差可达到0.57 m。本文方法定位结果显示：阵列布放于M3、N3时，距离阵列或原点较近的泄漏点(A3、B3点)定位误差与SI-LMS 法基本一致，但距离超过4 m 且靠近y轴的A1、A2、B1和B2点定位误差仍保持较低，而偏离y轴的C1、C2和C3点定位误差则明显增加。相较于现有适用于近场声源三维定位的TDOA 法，本文方法对y轴附近较远距离的泄漏点定位精度有较明显提升。

表3 定位结果Table 3 Results of localization

综上所述，使用本文方法进行泄漏定位时，在有效信号检测范围内，应增大阵列布放间距；在不影响远场声源判据成立的前提下，应使用较大孔径阵列；当定位结果偏离y轴时应适当调整阵列位置，以保证泄漏点位于y轴附近，提高定位精度。

5 结论

为实现对管道泄漏远场声源的三维定位，本文使用交叉定位法对现有TDOA 定位方法进行了改进。将泄漏声源视为远场声源并使用TDOA 法进行不同位置的两次定向，提出空间不相交直线的伪交点求取方法，对两组空间方位进行交叉求取交点，从而完成定位。针对本文信号采样率高、频率低的特点，选取基本互相关法进行延时估计。建立管道泄漏定位实验平台，对阵列孔径、布放间距以及泄漏位置等因素对定位精度的影响进行了实验分析。定位结果表明：适当提高阵列间距、增大阵列孔径以及调整阵列指向，能够提高定位精度；与现有TDOA法相比，基于交叉定位改进的TDOA法对距离原点4 m 以上远场泄漏声源的定位精度有明显提升。

由于实验条件所限，本文采用阵列先后两次布放的方式进行定位，且阵元数和阵列孔径无法进一步扩大。实际定位中，应采用双阵列同时布放、采集的方式提高定阵列参数的一致性和可控性，同时应进一步探索更多阵元数、更复杂阵列模型以及更大阵列孔径对定位精度的有益影响。