CN102829925A - 平面八元阵列电站锅炉承压管泄漏定位方法 - Google Patents

Abstract

本发明公开了属于火力发电电站锅炉泄漏监测技术领域的一种平面八元阵列电站锅炉承压管泄漏定位方法。依托现有的声波测温系统中开发出基于平面八元阵列的8个测点水冷壁泄漏定位系统。通过最小均方算法的自适应滤波器法得到时间迟延估计，并采用基于蜂窝网络的到达时间差定位TDOA定位系统中的改进成三维空间算法的CHAN算法实现锅炉承压管的泄漏的精确位置定位。本次阵列布置方式是参考了现有的声波测温技术的布置方式，即能够同时进行测温和泄漏信号的定位，本发明提供一种有效的、精确的锅炉承压管泄漏定位方法。

Description

平面八元阵列电站锅炉承压管泄漏定位方法

技术领域

本发明属于火力发电电站锅炉泄漏监测技术领域，特别涉及一种平面八元阵列电站锅炉承压管泄漏定位方法。具体说是涉及火力发电机组的锅炉内水冷壁、过热器、再热器和省煤器受热面管道泄漏的精确定位技术。

背景技术

炉膛管道泄漏声学监测法的基本原理为：在锅炉合适位置（例如观火孔）布置传声器，以接收水冷壁泄漏信号。当水冷壁发生泄漏时，管内的高温高压蒸汽从裂缝或破口喷射出来，产生宽频噪声，并沿空气和金属管道传播。然而，锅炉正常运行时会产生较强的背景噪声。传声器采集现场所有声音数据，并通过前置放大器转变为电信号，电信号远距离传输到控制室主机。监测系统采集炉膛声音信号，并将数据经计算机软件根据数学模型进行计算分析，监测背景噪声和异常泄漏信号，判断是否存在泄漏，发出报警，并经定位算法，进行精确定位

目前，国内外炉管泄漏检测装置采用大量布置测点覆盖锅炉受热面的方式，在滤掉锅炉背景噪声后检测泄漏声的声压级大小和泄漏声的频谱特征从而诊断泄漏是否发生。若某测点发生报警，则确定泄漏源位于该测点为圆心，10米为半径的半球空间内，所以装置的主要作用还是判断泄漏并确定泄漏的受热面，并不能定位到具体的管排上。对于泄漏孔径1～4mm的情况，检修任务往往需要花费大量的人力和时间，其面临的重大技术难题是解决泄漏源的精确定位问题。

发明内容

本发明的目的是提供一种平面八元阵列电站锅炉承压管泄漏定位方法，其特征在于，将八个传声器布置在电站锅炉炉膛内一个横截水平平面内组成八元阵列，接受泄漏信号，经过对声信号的互相关处理得出声音到达不同位置传声器的时间差，再根据目标与传声器位置之间的几何关系确定出其泄漏位置；具体定位步骤如下：

步骤一：定义传声器i为第i个传声器，i为1-8的正整数；S为泄漏点，各传声器坐标分别为第1传声器 (a,0,0)，第2传声器 (a1,0,0)，第3传声器 (a2,b,0)，第4传声器 (a2,b1,0)，第5传声器 (a1,b2,0)，第6传声器 (a,b2,0)，第7传声器 (0,b1,0)，第8传声器 (0,b,0)，若泄漏点S坐标为(x, y, z)，泄漏点S到传声器i与传声器j的时间迟延用τ_ij表示。

得到用于确定泄漏声源的双曲面定位方程组为：

${\mathrm{c\τ}}_{\mathrm{ij}}=\left|\right|{\mathrm{\ΔS}}_i|-|{\mathrm{\ΔS}}_j\left|\right|$ ；

式中：c为炉膛声波的传播速度；ΔS_i为泄漏声源到i传声器的距离；ΔS_j为泄漏声源到j传声器的距离；τ_ij为TDOA（到达时间差定位），用最小均方（LMS）自适应滤波器法求出；

步步骤二：通过最小均方（LMS）算法的自适应滤波器法计算各路信号之间的时间延迟效应；

首先将时间延迟效应进行自适应滤波器的参数转化；即将信号之间的延时效应看做信号通过一个传递函数为H(f)=e^-j2tD的滤波器的效应，源信号s(n)到s(n-D)的延迟等效为H(f)的输出时刻源信号的表达式（n表示时间变量，n-D也为时间变量；s(n)表示n时刻的源信号表达式， s(n-D)表示n-D时刻的源信号表达式）：

$s(n-D)=s\left(n\right)*h\left(n\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h\left(m\right)s(n-m)$

式中，h(m)为H(f)的单位脉冲响应函数。

其次，FIR滤波器参数的自适应实现，频域等价形式：

$H\left(f\right)=\frac{G_{x1x2}\left(f\right)}{G_{x1x1}\left(f\right)}$

变换到时域，有：

$h\left(m\right)=F^{(-1)}\left[\frac{G_{x1x2}\left(f\right)}{G_{x1x1}\left(f\right)}\right]$

式中，G_x1x2(f)和G_x1x1(f)分别为两信号的自功率谱和互功率谱。

利用自适应滤波器，得到时延估计值。当自适应过程收敛时，h(m)取最大值时的m值，即为时延估计值。

步骤三：泄漏信号S的初始位置为(x, y, z)，已知第i个传声器的位置为(x_i,y_i,z_i)，和泄漏点S之间的距离为R_t，那么：

${R_i}^2=K_i-2X_{i}x-2Y_{i}y-2Y_{i}z+x^2+y^2+z^2---\left(2\right)$

其中，K_i=X_i ²+Y_i ²+Z_i ²,i=1,2,3...，令R_i,1表示泄漏点S与传声器i和泄漏点S与传声器1的距离差，则：

$R_{i,1}={\mathrm{c\τ}}_{i,1}=R_i-R_1---\left(3\right)$

其中，c为当地声速，τ_i,1为TDOA测量值。将方程（2）先进行线性化处理，则

${R_i}^2={(R_{i,1}+R_1)}^2={R_{i,1}}^2+2R_{i,1}{R}_1+{R_1}^2---\left(4\right)$

取i=1，则式（3）为：

${R_1}^2=K_1-2X_{1}x-2Y_{1}y-2Z_{1}z+x^2+y^2+z^2---\left(5\right)$

令式(3)与式(4)相减得：

${R_{i,1}}^2+2R_{i,1}{R}_1=K_i-2X_{i,1}x-2Y_{i,1}y-2Z_{i,1}z-K_1---\left(6\right)$

式中，X_i,1=X_i-X₁，Y_i,1=Y_i-Y₁，Z_i,1=Z_i-Z₁。若视x，y，z，R₁为未知数，则式（6）可以视为线性方程组，求解该方程组可得泄漏信号S的坐标位置；

若要进行三维空间定位，则最少需要4只传声器，当有且只有4个传声器时，可测得3个TDOA值。在此假定R₁为已知，则S位置(x，y，z)由式错误!未找到引用源。展开成下式：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=-{\left[\begin{array}{ccc}{X}_{\mathrm{2,1}}& Y_{\mathrm{2,1}}& Z_{\mathrm{2,1}}\\ X_{\mathrm{3,1}}& Y_{\mathrm{3,1}}& Z_{\mathrm{3,1}}\\ X_{4,,1}& Y_{4,1}& Z_{\mathrm{4,1}}\end{array}\right]}^{-1}\{\left[\begin{array}{c}{R}_{\mathrm{2,1}}\\ R_{\mathrm{3,1}}\\ R_{\mathrm{4,1}}\end{array}\right]R_1+\frac{1}{2}\left[\begin{array}{c}{R_{\mathrm{2,1}}}^2-K_2+K_1\\ {R_{\mathrm{3,1}}}^2-K_3+K_1\\ {R_{\mathrm{4,1}}}^2-K_4+K_1\end{array}\right]\}---\left(7\right)$

式中，K_i=X_i ²+Y_i ²+Z_i ²,i=1,2,3...

将式（5-6）代人式（5-1），取i=1，可得有关R₁的二次方程，将所得正根代回式（5-6），即得泄漏点S的估计位置；根据炉膛内运行环境等先验信息来消除式（5-6）中的模糊性。

所述炉膛截面处测点分布超过5个时，CHAN算法可利用冗余TDOA值，因此该方法能尽量减小某一传声器的损坏造成的测量误差；对于非线性TDOA方程组第一步应将其转换为线性方程组，即可得到一个初始定位值，第二步采用加权最小二乘(WLS)算法，将第一步所得初始定位值附加变量等已知约束条件进行第二次WLS估计，从而改进了定位结果；

假设z=[z_p ^T,R₁]^T为未知量，其中Z_p=[x,y,z]^T，从式（5）得出以Z为变量的线性方程组：h=GZ。则泄漏点S对应的误差矢量为：

式中，

$h=\left[\begin{array}{c}{R_{\mathrm{2,1}}}^2-K_2+K_1\\ {R_{\mathrm{3,1}}}^2-K3+K1\\ .\\ .\\ .\\ {R_{M,1}}^2-K_M+K_1\end{array}\right]$ $G=-\left[\begin{array}{cccc}{X}_{\mathrm{2,1}}& Y_{\mathrm{2,1}}& Z_{\mathrm{2,1}}& R_{\mathrm{2,1}}\\ X_{\mathrm{3,1}}& Y_{\mathrm{3,1}}& Z_{\mathrm{3,1}}& R_{\mathrm{3,1}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ X_{M,1}& Y_{M,1}& Z_{M,1}& R_{M,1}\end{array}\right]$ ，

若令无噪声时{Δ}的表达形式为{Δ}⁰，那么τ_i,j=τ_i,j ⁰+n_i,j,R_i,1=R_i,1 ⁰+cn_i,1；同时有R_i ⁰=R_i,1 ⁰+R₁ ⁰，由此可得噪声的误差矢量为：

$\mathrm{\ψ}=\mathrm{cBm}+0.5c^{2}n\⊗n\≈\mathrm{cBn}---\left(9\right)$

式中，B=diag{R₂ ⁰,R₃ ⁰,...R_M ⁰}，

代表schur乘积。

误差矢量的协方差矩阵可表示为：

式中，Q为TDOA协方差矩阵；

先令Z间元素相互独立，然后运用加权最小二乘(WLS)算法对式（7）进行处理，由此得Z的估计：

$\begin{array}{c}\mathrm{Za}=\arg \min \left\{{(h-\mathrm{GZ})}^T{\mathrm{\ψ}}^{-1}(h-\mathrm{GZ})\right\}\\ ={\left(G^T{\mathrm{\ψ}}^{-1}G\right)}^{-1}\left(G^T{\mathrm{\ψ}}^{-1}h\right)\end{array}---\left(11\right)$

由于B=diag{R₂ ⁰,R₃ ⁰,...R_M ⁰}中包含未知泄漏点S的位置信息，所以ψ是未知量，上式无法给出独立解，应将结果进一步近似。

式(10)中的协方差矩阵ψ用协方差矩阵Q近似代替得：

$\mathrm{Za}={\left(G^T{Q}^{-1}G\right)}^{-1}\left(G^T{Q}^{-1}h\right)---\left(12\right)$

B矩阵可由式（12）中的初始解计算得出，经过计算式（10）得出第一次WLS的估计值。进行重复迭代，可改进泄漏点估计位置。假定z=z⁰+ΔZ，ΔZ及z协方差矩阵为：

$\mathrm{\ΔZ}=c{\left(G^{0T}{\mathrm{\ψ}}^{-1}{G}^0\right)}^{-1}{G}^{0T}{\mathrm{\ψ}}^{-1}\mathrm{Bn}---\left(13\right)$

$\mathrm{cov}\left(Z\right)=E\left[{\mathrm{\ΔZ\ΔZ}}^T\right]={\left(G^{0T}{\mathrm{\ψ}}^{-1}{G}^0\right)}^{-1}---\left(14\right)$

以上计算过程中前提是假设元素x，y，z和R相互独立，但是实际测量中这几个元素之间存在一定关系，应用这一关系将能得出更准确的泄漏信号估计值。随机矢量Z的均值即为为实际值，其元素可表示为：

$Z_1=x^0+e_1,Z_2=y^0+e_2,Z_3=z^0+e_3,Z_4={R_1}^0+e_4---\left(15\right)$

其中e₁，e₂，e₃，e₄表示Z的估计误差。新的误差矢量ψ＇为：

${\mathrm{\ψ}}^{\′}=h^{\′}-G^{\′}{Z}^{\′}---\left(16\right)$

其中， $h^{\′}=\left[\begin{array}{c}{(Z_1-X_1)}^2\\ {(Z_2-Y_1)}^2\\ {(Z_3-Z_1)}^2\\ {\left({Z_1}^2\right)}^2\end{array}\right]$ ， $G^{\′}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 1& 1& 1\end{array}\right]$ ， $Z^{\′}=\left[\begin{array}{c}{(x-X_1)}^2\\ {(y-Y_1)}^2\\ {(z-Z_1)}^2\end{array}\right]$

将式（13）代人式（14）中得出ψ，当误差e_i，较小时ψ的协方差矩阵为：

$\mathrm{\ψ}=E\left[{\mathrm{\ψ\ψ}}^T\right]=4B^{\′}\mathrm{cov}\left(Z\right)B^{\′}---\left(17\right)$

$B^{\′}=\mathrm{diag}\{x^0-X_1,y^0-Y_1,z^0-Z_1,{R_1}^0\}---\left(18\right)$

其中ψ为高斯分布，Z＇的ML估计为：

$Z^{\′}={\left(G{\mathrm{\ψ}}^{-1}G\right)}^{-1}\left(G{\mathrm{\ψ}}^{-1}h\right)---\left(19\right)$

泄漏点的最终定位计算结果为：

${Z_p}^{\′}=\sqrt{Z^{\′}}+\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]---\left(20\right)$

或者

${Z_p}^{\′}=-\sqrt{Z^{\′}}+\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]---\left(21\right)$

定位估计的模糊性可通过系统中的先验信息进行排除。

本发明的有益效果是突破了当前锅炉四管泄漏检测中单点各自探测的思路，考虑各个单点传声器之间的相干信息，将传声器组成了八元阵列接收泄漏信号，经过对声信号的互相关处理得出声音到达不同位置传声器的时间差，再根据目标与传声器位置之间的几何关系确定其泄漏位置。同时，本次阵列布置方式是参考了现有的声波测温技术的布置方式，即能够同时进行测温和泄漏信号的定位，本发明提供一种有效的、精确的锅炉承压管泄漏定位方法。

附图说明

图1 为传声器八元阵列结构及其在炉膛的分布示意图。

图2 为锅炉承压管泄漏定位系统图。

图3 为现场试验吹灰器的泄漏信号波形图，a关闭吹灰器的状态和b开启吹灰器的状态。

具体实施方式

本发明采用平面八元阵列对电站锅炉承压管泄漏进行精确定位，其传声器的阵列结构及其在炉膛的分布如图1所示，测点布置可根据现场情况适当改动。图2为锅炉承压管泄漏精确定位系统图。传声器和前置放大器组合一体（cs_1……cs_8.），采用1/2英寸预极化驻极体测量传声器MP201（灵敏度50Mv/Pa）和ICP前置放大器MA201。前置放大器输出与信号调理器、DAQ板卡、PXI系统和主机串联；定位系统采用18V直流电源供电；前置放大器用于传声器的ICCP供电，并可对信号进行滤波与放大。调节档分三个：1倍、10倍、100倍，通过BNC接头与信号调理器连接，调理后的信号通过68针同轴接头电缆与数据采集DAQ板卡连接。采用labview软件和NI PXI-6133采集卡（将DAQ板卡采集的数据进行转换），每通道采样速率最高可达3MS/s。设置采样频率为

。（分子中的S，表示采样点sample；分母中的s，为时间秒。总体意义为每秒多少个采样点）

我们在国内安装了声学测温系统的300MW电站锅炉机组上进行了冷态、热态定位实验，利用声学测温系统的空气吹扫装置，压缩空气来自电厂仪用空气，压力0.6MPa，在吹扫管出口安装喷嘴，测温层平面上3m处的观火孔上，将喷嘴深入炉膛0.3m，模仿泄漏噪声，进行冷态定位研究。测量得到的泄漏点位置为（1,0.3,3），（14.4,0.3,3），单轴方向上定位结果的最大平均误差为0.301m，符合工程应用要求。

热态实验中，我们利用短时吹灰器吹灰时发出的很大的喷流噪声，泄漏点为前墙左侧传声器上方2.7m处的短处，其坐标位置为（3.5,0.2,2.7），坐标为测量值。热态定位误差较大，单轴的平均绝对误差达到半米，但是对于炉膛大空间来说，定位到1m以内的区域泄漏，仍然具有应用价值。

Claims (2)

1.一种平面八元阵列电站锅炉承压管泄漏定位方法，其特征在于，将八个传声器布置在电站锅炉炉膛内一个横截水平平面内组成八元阵列，使用传声器收集锅炉内泄漏信号，通过对八个传声器收集的声信号进行处理，得出声音到达不同位置传声器的时间差，再根据目标与传声器位置之间的几何关系确定出其泄漏位置；具体定位步骤如下：

步骤一：定义传声器i为第i个传声器，i为1-8的正整数；S为泄漏点，各传声器坐标分别为第1传声器 (a,0,0)，第2传声器 (a1,0,0)，第3传声器 (a2,b,0)，第4传声器 (a2,b1,0)，第5传声器 (a1,b2,0)，第6传声器 (a,b2,0)，第7传声器 (0,b1,0)，第8传声器 (0,b,0)，若泄漏点S坐标为(x, y, z)，泄漏点S到传声器i与传声器j的时间迟延用τ_ij表示；其中，a,a1,a2,b,b1,b2分别为所述平面声阵列的长度和宽度；

得到用于确定泄漏声源的双曲面定位方程组为：

${\mathrm{c\τ}}_{\mathrm{ij}}=\left|\right|{\mathrm{\ΔS}}_i|-|{\mathrm{\ΔS}}_j\left|\right|;$

式中：c为炉膛声波的传播速度；ΔS_i为泄漏声源到i传声器的距离；ΔS_j为泄漏声源到j传声器的距离；τ_ij为到达时间差定位TDOA，用最小均方LMS自适应滤波器法求出；

步骤二：通过LMS算法的自适应滤波器法计算各路信号之间的时间延迟效应；

首先将时间延迟效应进行自适应滤波器的参数转化；即将信号之间的延时效应看做信号通过一个传递函数为H(f)=e^-j2tD的滤波器的效应，源信号s(n)到s(n-D)的延迟等效为H(f)的输出：其中n表示时间变量，n-D也为时间变量。  s(n)表示n时刻源信号的表达式    s(n-D)表示n-D时刻的源信号表达式

$s(n-D)=s\left(n\right)*h\left(n\right)=\underset{m=-\∞}{\overset{\∞}{\mathrm{\Σ}}}h\left(m\right)s(n-m)$

式中，h(m)为H(f)的单位脉冲响应函数。

其次，FIR滤波器参数的自适应实现，频域等价形式：

$H\left(f\right)=\frac{G_{x1x2}\left(f\right)}{G_{x1x1}\left(f\right)}$

变换到时域，有：

$h\left(m\right)=F^{(-1)}\left[\frac{G_{x1x2}\left(f\right)}{G_{x1x1}\left(f\right)}\right]$

式中，G_x1x2(f)和G_x1x1(f)分别为两信号的自功率谱和互功率谱。

利用自适应滤波器，得到时延估计值。当自适应过程收敛时，h(m)取最大值时的m值，即为时延估计值。

步骤三：泄漏信号S的初始位置为(x, y, z)，已知第i个传声器的位置为(x_i,y_i,z_i)，和泄漏点S之间的距离为R_i，那么：

${R_i}^2=K_i-2X_{i}x-2Y_{i}y-2Y_{i}z+x^2+y^2+z^2---\left(2\right)$

其中，K_i=X_i ²+Y_i ²+Z_i ²,i=1,2,3...，令R_i,1表示泄漏点S与传声器i和泄漏点S与传声器1的距离差，则：

$R_{i,1}={\mathrm{c\τ}}_{i,1}=R_i-R_1---\left(3\right)$

其中，c为当地声速，τ_i,1为TDOA测量值。将方程（2）先进行线性化处理，则

${R_i}^2={(R_{i,1}+R_1)}^2={R_{i,1}}^2+2R_{i,1}{R}_1+{R_1}^2---\left(4\right)$

取i=1，则式（3）为：

${R_1}^2=K_1-2X_{1}x-2Y_{1}y-2Z_{1}z+x^2+y^2+z^2---\left(5\right)$

令式(3)与式(4)相减得：

${R_{i,1}}^2+2R_{i,1}{R}_1=K_i-2X_{i,1}x-2Y_{i,1}y-2Z_{i,1}z-K_1---\left(6\right)$

式中，X_i,1=X_i-X₁，Y_i,1=Y_i-Y₁，Z_i,1=Z_i-Z₁。若视x，y，z，R₁为未知数，则式（6）可以视为线性方程组，求解该方程组可得泄漏信号S的坐标位置；

若要进行三维空间定位，则最少需要4只传声器，当有且只有4个传声器时，可测得3个TDOA值。在此假定R₁为已知，则S位置(x，y，z)由式错误!未找到引用源。展开成下式：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=-{\left[\begin{array}{ccc}{X}_{\mathrm{2,1}}& Y_{\mathrm{2,1}}& Z_{\mathrm{2,1}}\\ X_{\mathrm{3,1}}& Y_{\mathrm{3,1}}& Z_{\mathrm{3,1}}\\ X_{4,,1}& Y_{4,1}& Z_{\mathrm{4,1}}\end{array}\right]}^{-1}\{\left[\begin{array}{c}{R}_{\mathrm{2,1}}\\ R_{\mathrm{3,1}}\\ R_{\mathrm{4,1}}\end{array}\right]R_1+\frac{1}{2}\left[\begin{array}{c}{R_{\mathrm{2,1}}}^2-K_2+K_1\\ {R_{\mathrm{3,1}}}^2-K_3+K_1\\ {R_{\mathrm{4,1}}}^2-K_4+K_1\end{array}\right]\}---\left(7\right)$

式中，K_i=X_i ²+Y_i ²+Z_i ²,i=1,2,3...

将式（5-6）代人式（5-1），取i=1，可得有关R₁的二次方程，将所得正根代回式（5-6），即得泄漏点S的估计位置；根据炉膛内运行环境等先验信息来消除式（5-6）中的模糊性。

2.根据权利要求1所述平面八元阵列电站锅炉承压管泄漏定位方法，其特征在于，所述炉膛截面处测点分布超过5个时，CHAN算法可利用冗余TDOA值，因此该方法能尽量减小某一传声器的损坏造成的测量误差；对于非线性TDOA方程组，第一步应将其转换为线性方程组，即得到一个初始定位值，第二步采用加权最小二乘WLS算法，将第一步所得初始定位值附加变量等已知约束条件进行第二次WLS估计，从而改进了定位结果；

假设z=[z_p ^T,R₁]^T为未知量，其中Z_p=[x,y,z]^T，从式（5）得出以Z为变量的线性方程组：h=GZ。则泄漏点S对应的误差矢量为：

式中，

$h=\left[\begin{array}{c}{R_{\mathrm{2,1}}}^2-K_2+K_1\\ {R_{\mathrm{3,1}}}^2-K3+K1\\ .\\ .\\ .\\ {R_{M,1}}^2-K_M+K_1\end{array}\right]$ $G=-\left[\begin{array}{cccc}{X}_{\mathrm{2,1}}& Y_{\mathrm{2,1}}& Z_{\mathrm{2,1}}& R_{\mathrm{2,1}}\\ X_{\mathrm{3,1}}& Y_{\mathrm{3,1}}& Z_{\mathrm{3,1}}& R_{\mathrm{3,1}}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ X_{M,1}& Y_{M,1}& Z_{M,1}& R_{M,1}\end{array}\right]$ ，

若令无噪声时{Δ}的表达形式为{Δ}⁰，那么τ_i,j=τ_i,j ⁰+n_i,j,R_i,1=R_i,1 ⁰+cn_i,1；同时有R_i ⁰=R_i,1 ⁰+R₁ ⁰，由此可得噪声的误差矢量为：

$\mathrm{\ψ}=\mathrm{cBm}+0.5c^{2}n\⊗n\≈\mathrm{cBn}---\left(9\right)$

式中，B=diag{R₂ ⁰,R₃ ⁰,...R_M ⁰}，

代表schur乘积，

误差矢量的协方差矩阵可表示为：

式中，Q为TDOA协方差矩阵；

先令Z间元素相互独立，然后运用加权最小二乘(WLS)算法对式（7）进行处理，由此得Z的估计：

$\begin{array}{c}\mathrm{Za}=\arg \min \left\{{(h-\mathrm{GZ})}^T{\mathrm{\ψ}}^{-1}(h-\mathrm{GZ})\right\}\\ ={\left(G^T{\mathrm{\ψ}}^{-1}G\right)}^{-1}\left(G^T{\mathrm{\ψ}}^{-1}h\right)\end{array}---\left(11\right)$

由于B=diag{R₂ ⁰,R₃ ⁰,...R_M ⁰}中包含未知泄漏点S的位置信息，所以ψ是未知量，上式无法给出独立解，应将结果进一步近似；

式(10)中的协方差矩阵ψ用协方差矩阵Q近似代替得：

$\mathrm{Za}={\left(G^T{Q}^{-1}G\right)}^{-1}\left(G^T{Q}^{-1}h\right)---\left(12\right)$

B矩阵可由式（12）中的初始解计算得出，经过计算式（10）得出第一次WLS的估计值；进行重复迭代，可改进泄漏点估计位置；假定z=z⁰+ΔZ,，ΔZ及z协方差矩阵为：

$\mathrm{\ΔZ}=c{\left(G^{0T}{\mathrm{\ψ}}^{-1}{G}^0\right)}^{-1}{G}^{0T}{\mathrm{\ψ}}^{-1}\mathrm{Bn}---\left(13\right)$

$\mathrm{cov}\left(Z\right)=E\left[{\mathrm{\ΔZ\ΔZ}}^T\right]={\left(G^{0T}{\mathrm{\ψ}}^{-1}{G}^0\right)}^{-1}---\left(14\right)$

以上计算过程中前提是假设元素x，y，z和R相互独立，但是实际测量中这几个元素之间存在一定关系，应用这一关系将能得出更准确的泄漏信号估计值；随机矢量Z的均值即为为实际值，其元素可表示为：

$Z_1=x^0+e_1,Z_2=y^0+e_2,Z_3=z^0+e_3,Z_4={R_1}^0+e_4---\left(15\right)$

其中e₁，e₂，e₃，e₄表示Z的估计误差。新的误差矢量ψ＇为：

${\mathrm{\ψ}}^{\′}=h^{\′}-G^{\′}{Z}^{\′}---\left(16\right)$

其中， $h^{\′}=\left[\begin{array}{c}{(Z_1-X_1)}^2\\ {(Z_2-Y_1)}^2\\ {(Z_3-Z_1)}^2\\ {\left({Z_1}^2\right)}^2\end{array}\right]$ ， $G^{\′}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\\ 1& 1& 1\end{array}\right]$ ， $Z^{\′}=\left[\begin{array}{c}{(x-X_1)}^2\\ {(y-Y_1)}^2\\ {(z-Z_1)}^2\end{array}\right]$

将式（13）代人式（14）中得出ψ，当误差e_i，较小时ψ的协方差矩阵为：

$\mathrm{\ψ}=E\left[{\mathrm{\ψ\ψ}}^T\right]=4B^{\′}\mathrm{cov}\left(Z\right)B^{\′}---\left(17\right)$

$B^{\′}=\mathrm{diag}\{x^0-X_1,y^0-Y_1,z^0-Z_1,{R_1}^0\}---\left(18\right)$

其中ψ为高斯分布，Z＇的ML估计为：

$Z^{\′}={\left(G{\mathrm{\ψ}}^{-1}G\right)}^{-1}\left(G{\mathrm{\ψ}}^{-1}h\right)---\left(19\right)$

泄漏点的最终定位计算结果为：

${Z_p}^{\′}=\sqrt{Z^{\′}}+\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]---\left(20\right)$

或者

${Z_p}^{\′}=-\sqrt{Z^{\′}}+\left[\begin{array}{c}{X}_1\\ Y_1\\ Z_1\end{array}\right]---\left(21\right)$

定位估计的模糊性可通过系统中的先验信息进行排除。

Images