CN103135094B - 基于bfgs拟牛顿法的信号源定位方法 - Google Patents

Abstract

本发明公开一种基于BFGS拟牛顿法的信号源定位方法，主要解决现有定位方法对定位引入中间误差、定位结果模糊、实用性不强、计算量大且定位结果不稳定的问题。本发明实现的具体步骤是：1、获得到达时间差的测量值与达到增益比的测量值；2、利用真实值等于测量值减去测量误差的关系，建立联合定位误差方程；3、利用BFGS拟牛顿法获得定位误差方程的解；4、将定位方程的解加上参考监测节点的位置坐标值，得到目标源的位置坐标值。本发明可实现稳定及高精度定位，低信噪比的环境下也适用，可实用性强。

Description

基于BFGS拟牛顿法的信号源定位方法

技术领域

本发明属于通信技术领域，更进一步涉及无线通信技术和信号处理技术领域中基于信号的到达时间差(Time Difference ofArrival，TDOA)与到达增益比的(Gain Ratios ofArrival，GROA)的BFGS(Broyden-Fletcher-Goldfarb-Shanno)拟牛顿法信号源定位方法。本发明可用于无线传感器网络或蜂窝网中，实现对目标信号源的定位。

背景技术

随着通信技术与信息技术的迅速发展，信号源定位技术作为信号处理领域的重要研究内容，已在电磁频谱监测、雷达、传感器网络、无线通信等领域得到广泛应用，研究高精度与实时性的定位方法具有重要的意义。基于到达时间差的定位方法以其定位精度高、算法计算复杂性低、易于实现等诸多的优势而受到越来越多的重视。

西安电子科技大学提出的专利申请“基于监测节点圆周分布的TDOA定位方法”(申请号：201010598701.1，公开号：102026370A)中公开了一种基于监测节点呈圆周分布的TDOA定位方法。该方法根据已知的到达时间差的测量值和监测节点位置坐标，利用Chan算法多次计算得到多个目标节点的坐标，将圆周半径的0.2倍作为门限值，对于门限内定位坐标进行统计平均，将统计平均的得到的坐标作为泰勒算法的初始值。该定位方法对于监测节点呈圆周分布且目标节点处于距离圆心0.2倍的半径范围内的情况下定位精度有所提高。但是，该方法仍然存在的不足是，在实际环境中的实用性不强，在实际环境中监测节点是随机布设的，大多数情况下很难满足监测节点按圆周分布且目标节点处于距离圆心0.2倍的半径范围内的要求。算法中所涉及的Chan定位算法是一种基于到达时间差的定位法，该方法由于具有较高的定位精度而得到了广泛的应用。Chan定位算法采用了两步加权最小二乘算法，但为了使定位方程组线性化，在第一步加权最小二乘中引入了辅助变量，并假定辅助变量与信号源的位置是无关的(实际上两者之间是有关的)，这一假设使得求解过程引入了中间误差；在第二步利用加权最小二乘法时，所得的定位结果具有模糊性。

朱亚坤、冯立杰发表的“基于牛顿迭代搜索法的多节点协同振源定位研究”(传感器技术学报，2009年9月，第24卷，第9期)文章公开一种定位方法。该方法首先用最小二乘法的定位结果作为基本牛顿法的初始值，然后通过基本牛顿法来求解非线性定位方程组获得目标源的位置，定位精度有所提高。但是，该方法仍然存在的不足是，该方法对于可能存在海塞矩阵奇异的情况而不能确定后继点。当海塞矩阵非奇异时，也未必能保证海塞矩阵是正定的，由此导致算法失效。当海塞矩阵满足正定时，虽然算法有效，但是，需要计算海塞矩阵的逆使的计算量较大，尤其是当矩阵的阶数较高时，计算量的增加更为明显。

中国科学院计算机技术研究所提出的专利申请“一种无线传感器网络的节点定位方法”(申请号：200810103124.7，公开号：101251592A)中公开了一种无线传感器网络的节点定位方法。该定位方法首先为无线传感器网络各局部构建局部相对坐标，然后将各个局部相对坐标进行融合，得到所有节点的全局相对坐标，并使用位置信息已知的信标节点，把全局坐标转换为全局绝对坐标，获得全局绝对初始坐标后，进行节点定位迭代求精而获得定位目标的位置，定位精度有所提高。该方法仍然存在的不足是，该定位方法中，获得DFP拟牛顿法初始坐标值的过程比较复杂，并且DFP拟牛顿法并不是当前变尺度拟牛顿法中数值稳定性最好的一个方法。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于BFGS拟牛顿法的信号源定位方法，减小了定位的计算量，提高了定位精度。

实现本发明的技术思路是，针对已有的定位方法中存在对监测节点分布的限制、引入中间误差、模糊解及计算量较大的问题，对空间监测节点的布设不设任何限制，建立定位误差方程组时，直接利用了辅助变量与目标源位置的关系，在获得定位误差方程的解时，采用了BFGS拟牛顿法，用最小二乘法得到的估计值作为BFGS拟牛顿的迭代初始值，将定位误差方程的解与参考监测节点位置坐标值之和作为未知目标源的位置。

本发明的具体实现步骤包括如下：

(1)获得测量值：

1a)在空间分布的无线传感器节点中任选一个作为参考监测节点，其余的无线传感器节点作为辅助监测节点；

1b)将参考监测节点与辅助监测节点接收信号的互相关函数最大值所对应的相关运算滞后时间作为到达时间差的测量值；

1c)按照下式计算到达增益比的测量值；

g＝R(ξ)/(V(ξ)-U(ξ))

其中，g表示到达增益比的测量值，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，V(ξ)表示参考监测节点接收信号的自相关函数，U(ξ)表示参考监测节点处噪声的自相关函数，ξ表示两个接收信号之间相关运算滞后时间。

(2)根据真实值等于测量值减去测量误差的关系，得到定位误差方程：

$\mathrm{\δ}=b-G\left[\begin{array}{c}v\\ \sqrt{v^{T}v}\end{array}\right]$

其中，δ表示定位测量误差向量，b表示定位误差方程的常数向量；G表示定位误差方程系数矩阵，v表示目标源的位置坐标值与参考监测节点位置坐标值的差。

(3)获得定位误差方程的解：

3a)将定位误差方程最小二乘估计结果的前三项作为BFGS拟牛顿法定位的迭代初始值，设置一个大于零的迭代终止允许误差值，最大迭代次数为1000次；

3b)判断目标代价函数一阶导在迭代初始值处的值是否小于迭代终止允许误差值，若小于，则将迭代初始值作为迭代终止估计值，执行步骤3g)，否则，执行步骤3c)；

3c)将海塞矩阵逆的近似矩阵初始化为单位矩阵，迭代次数初始值设定为0；

3d)利用BFGS拟牛顿迭代公式，计算目标源位置坐标值与参考监测节点位置坐标值之差的迭代估计值；

3e)判断目标代价函数在迭代估计值处的值是否小于迭代终止允许误差值，若小于，则将此次迭代估计值作为迭代终止估计值，执行步骤3g)，否则，执行步骤3f)；

3f)判断迭代次数是否等于最大迭代次数，若等于，则将迭代初始值作为迭代终止估计值，执行步骤3g)，否则，将迭代次数加1，执行骤3d)；

3g)终止迭代。

(4)确定目标源位置：

将迭代终止估计值与参考监测节点位置坐标值相加，其和作为未知目标源的位置。

本发明与现有技术相比，具有以下优点：

第一、本发明由于在建立定位误差方程组时，直接利用了辅助变量与目标源位置的关系，克服了现有技术中Chan定位技术存在的中间误差及模糊解缺点，使得本发明具有定位精度高的优点。

第二、本发明由于在获得定位误差方程的解时，采用了BFGS拟牛顿法，克服了现有技术中基本牛顿迭代定位技术存在海塞矩阵可能是奇异的或非正定的情况而导致算法失效的缺点与计算海塞矩阵逆所需计算量大的缺点，使得本发明具有计算量小、实时性强的优点。

第三、本发明由于用最小二乘法得到的估计值作为BFGS拟牛顿的迭代初始值，克服了现有技术中无线传感器网络定位技术中获得变尺度拟牛顿法迭代初始值的过程较复杂的缺点，并使用目前变尺度拟牛顿法中数值稳定性最好的BFGS拟牛顿法取代DFP拟牛顿法，使得本发明具有实现简单、定位结果更稳定的优点。

第四、本发明由于对空间监测节点的布设不设任何限制，克服了现有技术中基于监测节点圆周分布的到达时间差定位技术对监测节点按圆周分布且目标节点处于距离圆心0.2倍的半径范围内限制的缺点，使得本发明具有使用范围更广的优点。

附图说明

图1为本发明的流程图；

图2为本发明与基本牛顿法、Chan算法定位精度的仿真对比图；

图3为本发明与基本牛顿法、Chan算法定位偏差的仿真对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参考图1，本发明实现的具体步骤如下：

步骤1，获得测量值。

在空间分布的M个无线传感器节点中任选一个作为参考监测节点，其余的无线传感器节点作为辅助监测节点，其中M是大于3的节点个数，M取值的多少是根据用户对定位精度的要求、设备的代价等综合考虑而确定的。

将参考监测节点的接收信号与辅助监测节点的接收信号做互相关运算，得到参考监测节点与辅助监测节点接收信号的互相关函数如下：

R(ξ)＝αβP(ξ-D)+Q(ξ)

其中，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，ξ表示两个接收信号之间相关运算滞后时间，α和β分别表示未知目标信号到达参考监测节点和辅助监测节点i后的衰减系数；P(ξ-D)表示未知目标信号的自相关函数，D表示未知目标信号到达参考监测节点与辅助监测节点之间的到达时间差，“-”为延迟符号；Q(ξ)表示参考监测节点与辅助监测节点噪声的互相关函数。

考虑噪声之间是否相关，分以下两种情况分别对参考监测节点与辅助监测节点接收信号互相关函数做简化或修正。

当噪声之间互不相关时，则有Q(ξ)＝0，故参考监测节点的接收信号与辅助监测节点的接收信号的互相关函数可简化为：

R(ξ)＝αβP(ξ-D)

其中，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，ξ表示两个接收信号之间相关运算滞后时间，α和β分别表示未知目标信号到达参考监测节点和辅助监测节点后的衰减系数；P(ξ-D)表示未知目标信号的自相关函数，D表示两个无线传感器节点之间的到达时间差，“-”为延迟符号。

由自相关函数的性质：

|P(ξ-D)|≤P(0)

则，互相关函数R(ξ)在ξ＝D处取得最大值。互相关函数取得最大值时的相关滞后时间ξ便是到达时间差D的测量值。

当噪声之间相关时，Q(ξ)≠0，R(ξ)在ξ＝D处不一定取得峰值，为了得到的测量值更准确，对互相关函数进行加权处理修正，使得R(ξ)在ξ＝D处能够取得峰值，同时根据维纳辛钦定理，信号的互相关函数等于互谱函数的傅里叶变换，则按下式计算参考监测节点与辅助监测节点接收信号的互相关函数：

$R\left(\mathrm{\ξ}\right)=\underset{-\∞}{\overset{\∞}{\∫}}S\left(f\right)\mathrm{\Φ}\left(f\right)e^{j2\mathrm{\πf\ξ}}\mathrm{df}$

其中，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，ξ表示两个接收信号之间相关运算滞后时间，S(f)表示参考监测节点与辅助监测节点接收信号的互谱函数，Φ(f)表示加权函数，f表示互谱函数的频率。

根据参考监测节点与辅助监测节点接收信号的互相关函数最大值所对应的相关运算滞后时间作为到达时间差的测量值，则按照下式计算到一个达时间差的测量值：

d＝argmaxR(ξ)

其中d表示到达时间差的测量值，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，ξ表示两个接收信号之间相关运算滞后时间，argmax(·)符号表示取函数(·)最大值所对应的自变量；

将参考监测节点接收信号做自相关运算，得到参考监测节点接收信号的自相关函数如下：

V(ξ)＝α²X(ξ)+U(ξ)

其中，V(ξ)表示参考监测节点接收信号的自相关函数，ξ表示两个接收信号之间相关运算滞后时间，α表示未知目标信号到达参考监测节点后的衰减系数；X(ξ)表示未知目标信号的自相关函数；U(ξ)表示参考监测节点处噪声的自相关函数。

按照下式计算到达增益比的测量值：

g＝R(ξ)/(V(ξ)-U(ξ))

其中，g表示到达增益比的测量值，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，V(ξ)表示参考监测节点接收信号的自相关函数，U(ξ)表示参考监测节点处噪声的自相关函数，ξ表示两个接收信号之间相关运算滞后时间。

步骤2，得到定位误差方程。

根据距离等于速度乘以时间，则到达距离差的测量值可按照下式计算：

e_i＝cd_i，i＝2，3，...，M

其中，e_i表示目标源到达辅助监测节点i与参考节点的到达距离差的测量值，c表示信号在空间中的传播速度，d_i表示目标源到达辅助监测节点i与参考监测节点的到达时间差的测量值，M表示监测节点的数目。

按照下式计算到达距离差的真实值：

$e_i^o=\sqrt{{(u-s_i)}^T(u-s_i)}-\sqrt{v^{T}v},i=\mathrm{2,3},...,M$

其中表示目标源到达辅助监测节点i与参考监测节点的到达距离差的真实值，u表示目标源的位置坐标值，s_i表示辅助监测节点i的三维位置坐标值，v表示目标源的三维位置坐标值与参考监测节点的三维坐标值之差，M表示监测节点的数目。

根据到达时间差的真实值等于到达时间差的测量值减去到达时间差测量误差的关系，得到到达时间差定位误差方程如下：

${\mathrm{\ξ}}_d(i-1)={e_i}^2-{(s_i-s)}^T(s_i-s)+2{(s_i-s)}^{T}v+{2e}_i\sqrt{v^{T}v},i=\mathrm{2,3},...,M$

其中ξ_d(i-1)表示第i-1个方程的到达时间差随机误差，e_i表示目标源到达辅助监测节点i与参考监测节点的到达距离差的测量值，s_i表示辅助监测节点i的三维位置坐标值，s表示参考监测节点的三维位置坐标值，v表示目标源的三维位置坐标值与参考监测节点的三维坐标值之差，M表示监测节点的数目。

根据信号的衰减与距离成反比，故目标源到达辅助监测节点i相对于到达参考监测节点的到达增益比的真实值可按照下式计算：

$y_i=\frac{\sqrt{{(u-s_i)}^T(u-s_i)}}{\sqrt{v^{T}v}},i=\mathrm{2,3},...,M$

其中y_i表示目标源到达辅助监测节点i相对于到达参考监测节点的到达增益比的真实值，u表示目标源的三维位置坐标值，s_i表示辅助监测节点i的三维位置坐标值，v表示目标源的三维位置坐标值与参考监测节点的三维坐标值之差，M表示监测节点的数日。

根据到达增益比的真实值等于到达增益比的测量值减去到达增益比测量误差的关系，得到的到达增益比的定位误差方程如下：

${\mathrm{\ξ}}_g(i-1)=e_i-(g_i-1)\sqrt{v^{T}v},i=\mathrm{2,3},...,M$

其中ξ_g(i-1)表示第i-1个方程的到达增益比随机误差，e_i表示目标源到达辅助监测节点i与参考监测节点的到达距离差的测量值，g_i表示目标源到达辅助监测节点i相对于到达参考监测节点的到达增益比的测量值，v表示目标源的三维位置坐标值与参考监测节点的三维坐标值之差，M表示监测节点的数目。

将到达时间差定位误差方程与到达增益比定位误差方程联列在一起，得到定位误差方程，整理成矩阵形式为：

$\mathrm{\δ}=b-G\left[\begin{array}{c}v\\ \sqrt{v^{T}v}\end{array}\right]$

其中，δ表示定位测量误差向量，b表示定位误差方程的常数向量；G表示定位误差方程系数矩阵，v表示目标源的位置坐标值与参考监测节点位置坐标值的差。

步骤3，获得定位误差方程的解。

3a)对定位误差方程按照下式计算得到最小二乘估计值：

$\widehat{\mathrm{\θ}}={\left(G^{T}G\right)}^{-1}{G}^{T}b$

其中向量的前三个元素表示目标源的位置坐标与参考监测节点位置坐标的差值的估计值，向量的第四个元素表示目标源到参考监测节点距离的估计值，G表示定位误差方程系数矩阵，b表示定位误差方程的常数向量。

3b)将向量的前三个元素作为BFGS拟牛顿法定位的迭代初始值m₀，设置一个大于零的迭代终止允许误差值，最大迭代次数为1000次。

3c)按照下式计算目标代价函数一阶导在迭代初始值m₀处的值，并判断该值是否小于迭代终止允许误差值，若小于，则将迭代初始值作为迭代终止估计值，执行步骤3h)，否则，执行步骤3d)：

$J\left(m_0\right)=(A^T+\sqrt{n_0^T{m}_0}\·m_0{B}^T)(W^{-1}+W^{-T})({\mathrm{Am}}_0+B\left(\sqrt{m_0^T{m}_0}\right)-b)$

其中，J(m₀)表示目标代价函数一阶导在迭代初始值处的值，m₀表示迭代初始值，A表示定位误差方程系数矩阵前三列元素构成的矩阵，B表示定位误差方程系数矩阵第四列元素构成的向量，W表示定位误差方程的加权系数矩阵，“T”表示共轭转置，“-1”表示矩阵的逆，b表示定位误差方程的常数向量。

3d)将海塞矩阵逆的近似矩阵初始化为单位矩阵，迭代次数k的初始值设定为0。

3e)利用如下BFGS拟牛顿迭代公式，计算目标源位置坐标值与参考监测节点位置坐标值之差的迭代估计值：

m_k+1＝m_k+λd

其中，m_k+1表示沿搜索方向第k+1次的迭代估计值，m_k表示沿搜索方向第k次的迭代估计值，λ表示搜索步长，d表示搜索方向。

3f)按照下式计算目标代价函数在迭代估计值处的值，并判断该值是否小于迭代终止允许误差值，若小于，则将此次迭代估计值作为迭代终止估计值，执行步骤3h)，否则，执行步骤3e)：

$J\left(m_{k+1}\right)=(A^T+\sqrt{m_{k+1}^T{m}_{k+1}}\·m_{k+1}{B}^T)(W^{-1}+W^{-T})({\mathrm{Am}}_{k+1}+B\left(\sqrt{m_{k+1}^T{m}_{k+1}}\right)-b)$

其中，J(m_k+1)表示目标代价函数一阶导在迭代估计值处的值，m_k+1表示搜索方向第k+1次的迭代估计值，A表示定位误差方程系数矩阵前三列元素构成的矩阵，B表示定位误差方程系数矩阵第四列元素构成的向量，W表示定位误差方程的加权系数矩阵，“T”表示共轭转置，“-1”表示矩阵的逆，b表示定位误差方程的常数向量。

3g)判断迭代次数k是否等于最大迭代次数1000，若等于，则将迭代初始值作为迭代终止估计值，执行步骤3h)，否则，将迭代次数加1，执行步骤3e)。

3h)终止迭代。

步骤4，确定目标源位置。

将迭代终止估计值与参考监测节点位置坐标值相加，其和作为未知目标源的位置。

下面结合附图2、附图3对本发明的效果做进一步的描述：

一、仿真条件

本发明的仿真均是在以下条件下进行的：到达时间差的测量值与到达增益比的测量值均是由真实值加上零均值的高斯随机误差生成，参与定位的监测节点数为M＝6，其中参考监测节点的位置坐标值为(0，100，100)，其他辅助监测节点的位置坐标值分别为(-5，500，5)、(-100，-100，-100)、(350，200，100)、(125，140，145)与(400，0，0)，目标源的位置坐标值为(500，650，550)，单位均为m，空间传播系数等于1，仿真信噪比的变化范围为-20dB到15dB，置定位循环次数为10000次。

二、仿真内容与结果

仿真1：用本发明与基本牛顿法、Chan算法在仿真条件下，分别在不同信噪比下，对目标节点进行10000次定位，计算定位结果与真实位置之间的误差值。分别在不同信噪比下，计算理论最小误差值。在图2中，以三角标示的曲线表示本发明的定位误差性能曲线，以星号标示的曲线表示基本牛顿发的定位误差性能曲线，以方块标示的曲线表示Chan算法的定位误差性能曲线，以圆圈标示的曲线表示克拉美罗下界曲线。

由图2可见，在不同的信噪比下，三种方法的定位误差性能曲线与克拉美罗界相比，应用本发明的定位误差性能曲线更接近克拉美罗下界曲线，说明本发明的定位误差最小，由此对目标节点的定位精度就更高。将本发明与Chan算法的定位误差性能曲线与克拉美罗下界曲线相比，在信噪比较高时，两种定位方法均可达到克拉美罗下界，即能达到理论最优定位，但在信噪比较低时，本发明的定位误差值小于Chan算法的定位误差值，本发明的优势就突显出来，可见本发明抗噪声性能较好，不仅可适用于高信噪比的环境，而且适合低信噪比的环境，在真实环境下的实用性强。比较可见，应用本发明的定位方法定位的精度更好，可实用性强，定位性能有进一步的提高。

仿真2：用本发明与基本牛顿法、Chan算法在仿真条件下，分别在不同信噪比下，对目标节点进行10000次定位，计算定位结果与真实位置之间的平均定位偏差值。在图3中，以三角标示的曲线表示本发明的定位偏差性能，以星号标示的曲线表示基本牛顿发的定位偏差性能曲线，以方块标示的曲线表示Chan算法的定位偏差性能曲线。

由图3可见，在不同的信噪比下，三种方法的定位偏差性能曲线相比，在信噪比较高时，本发明与Chan算法的定位偏差性能曲线几乎重合且接近于零，但在信噪比较低时，本发明的平均定位偏差值小于Chan算法的平均定位偏差值，本发明的定位稳定性更好。而基本牛顿法则无论在信噪比较低还是信噪比较高，平均定位偏差值都比较大，这是因为基本牛顿法因不能保证海塞矩阵是正定的，而导致收敛的结果可能偏离了真实的目标位置，定位误差较大，偏差也较大。比较可见，本发明抗噪声性能较好，对于低信噪比的环境也适用，本发明的平均定位偏差值最小，对目标节点的定位稳定性就更好，定位性能有进一步的提高。

综合上述仿真结果，应用本发明的定位算法，能达到较好的定位精度，且定位的稳定性较好，是一种有效的定位方法。

Claims (1)

1.一种基于BFGS拟牛顿法的信号源定位方法，其具体实现步骤如下：

(1)获得测量值：

1a)在空间分布的无线传感器节点中任选一个作为参考监测节点，其余的无线传感器节点作为辅助监测节点；

1b)将参考监测节点与辅助监测节点接收信号的互相关函数最大值所对应的相关运算滞后时间作为到达时间差的测量值；

1c)按照下式计算到达增益比的测量值；

g＝R(ξ)/(V(ξ)-U(ξ))

其中，g表示到达增益比的测量值，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，V(ξ)表示参考监测节点接收信号的自相关函数，U(ξ)表示参考监测节点处噪声的自相关函数，ξ表示两个接收信号之间相关运算滞后时间；

所述参考监测节点与辅助监测节点接收信号的互相关函数等于如下：

R(ξ)＝αβP(ξ-D)+Q(ξ)

其中，R(ξ)表示参考监测节点与辅助监测节点接收信号的互相关函数，ξ表示两个接收信号之间相关运算滞后时间，α和β分别表示未知目标信号到达参考监测节点和辅助监测节点后的衰减系数；P(ξ-D)表示未知目标信号的自相关函数，D表示未知目标信号到达参考监测节点与辅助监测节点之间的到达时间差，“-”为延迟符号；Q(ξ)表示参考监测节点与辅助监测节点噪声的互相关函数；

所述参考监测节点接收信号的自相关函数如下：

V(ξ)＝α²X(ξ)+U(ξ)

其中，V(ξ)表示参考监测节点接收信号的自相关函数，ξ表示两个接收信号之间相关运算滞后时间，α表示未知目标信号到达参考监测节点后的衰减系数；X(ξ)表示未知目标信号的自相关函数；U(ξ)表示参考监测节点处噪声的自相关函数；

(2)根据真实值等于测量值减去测量误差的关系，得到定位误差方程：

$\mathrm{\δ}=b-G\begin{array}{}\left[\begin{array}{c}V\\ \sqrt{v^{T}v}\end{array}\right]\end{array}$

其中，δ表示定位测量误差向量，b表示定位误差方程的常数向量；G表示定位误差方程系数矩阵，v表示目标源的位置坐标值与参考监测节点位置坐标值的差；

(3)获得定位误差方程的解：

3a)对定位误差方程按照下式计算得到最小二乘估计值：

$\widehat{\mathrm{\θ}}={\left(G^{T}G\right)}^{-1}{G}^{T}b$

其中，向量的前三个元素表示目标源的位置坐标与参考监测节点位置坐标的差值的估计值，向量的第四个元素表示目标源到参考监测节点距离的估计值，G表示定位误差方程系数矩阵，b表示定位误差方程的常数向量；

3b)将向量的前三个元素作为BFGS拟牛顿法定位的迭代初始值m₀，设置一个大于零的迭代终止允许误差值，最大迭代次数为1000次；

3c)按照下式，计算目标代价函数一阶导在迭代初始值m₀处的值，并判断该值是否小于迭代终止允许误差值，若小于，则将迭代初始值作为迭代终止估计值，执行步骤3h)，否则，执行步骤3d)：

$J\left(m_0\right)=(A^T+\sqrt{m_0^T{m}_0}\·m_0{B}^T)(W^{-1}+W^{-T})({\mathrm{Am}}_0+B\left(\sqrt{m_0^T{m}_0}\right)-b)$

其中，J(m₀)表示目标代价函数一阶导在迭代初始值处的值，m₀表示迭代初始值，A表示定位误差方程系数矩阵前三列元素构成的矩阵，B表示定位误差方程系数矩阵第四列元素构成的向量，W表示定位误差方程的加权系数矩阵，“T”表示共轭转置，“-1”表示矩阵的逆，b表示定位误差方程的常数向量；

3d)将海塞矩阵逆的近似矩阵初始化为单位矩阵，迭代次数初始值设定为0；

3e)利用如下BFGS拟牛顿迭代公式，计算目标源位置坐标值与参考监测节点位置坐标值之差的迭代估计值：

m_k+1＝m_k+λd

其中，m_k+1表示沿搜索方向第k+1次的迭代估计值，m_k表示沿搜索方向第k次的迭代估计值，λ表示搜索步长，d表示搜索方向；

3f)按照下式计算目标代价函数在迭代估计值处的值，并判断该值是否小于迭代终止允许误差值，若小于，则将此次迭代估计值作为迭代终止估计值，执行步骤3h)，否则，执行步骤3e)：

$J\left(m_{k+1}\right)=(A^T+\sqrt{m_{K+1}^T{m}_{k+1}}\·m_{k+1}{B}^T)(W^{-1}+W^{-T})({\mathrm{Am}}_{k+1}+B\left(\sqrt{m_{k+1}^T{m}_{k+1}}\right)-b)$

其中，J(m_k+1)表示目标代价函数一阶导在迭代估计值处的值，m_k+1表示搜索方向第k+1次的迭代估计值，A表示定位误差方程系数矩阵前三列元素构成的矩阵，B表示定位误差方程系数矩阵第四列元素构成的向量，W表示定位误差方程的加权系数矩阵，“T”表示共轭转置，“-1”表示矩阵的逆，b表示定位误差方程的常数向量；

3g)判断迭代次数是否等于最大迭代次数，若等于，则将迭代初始值作为迭代终止估计值，执行步骤3h)，否则，将迭代次数加1，执行骤3e)；

3h)终止迭代；

(4)确定目标源位置：

将迭代终止估计值与参考监测节点位置坐标值相加，其和作为未知目标源的位置。