CN112285650B - 异常tdoa存在下未知波速声发射源定位方法、系统及存储介质 - Google Patents

Abstract

本发明公开了一种异常TDOA存在下未知波速声发射源定位方法、系统及存储介质，其中方法包括：利用TDOA测量值和各传感器坐标构建非线性双曲线方程并化为线性方程组；求解得到普通最小二乘解；通过迭代权重估计获取方程残差；对方程残差取绝对值得到绝对残差并排序；计算各绝对残差的间距、全局距离和局部距离来获取全局距离和局部距离判据；确定临界残差值，大于临界残差值即为异常残差，异常残差对应的TDOA为异常TDOA，将其排除；利用剩余正常的TDOA测量值和相应传感器坐标重新构造TDOA线性方程组，并获得声发射源坐标。该方案具有更加鲁棒的定位精度，体现在可排除异常TDOA，降低其对定位精度的影响；无需预先确定波速，消除了波速误差对定位精度的影响。

Description

异常TDOA存在下未知波速声发射源定位方法、系统及存储 介质

技术领域

本发明涉及结构健康监测领域，尤其涉及一种异常TDOA存在下的未知波速声发射源定位方法、系统及存储介质。

背景技术

声发射源定位技术作为一种重要的无损监测技术被广泛应用于矿山岩爆灾害预警、结构健康监测、材料损伤和破坏机理等关键课题的研究。一个高精度的声发射源定位方法是至关重要的，然而定位精度易受异常TDOA和波速误差的影响。现有大多数定位方法在不存在TDOA(到达时间差)异常值的情况下尚可取得较为理想的声发射源定位结果。但现实工程环境往往是非常复杂的，声发射传感器易被干扰信号等错误触发，进而产生异常信号。而现有自动到时拾取方法虽然可以减少了人工拾取到时的主观性，但仍无法彻底消除异常TDOA的产生，特别是对于高噪音信号。另外，通道串扰、声发射信号微弱、信号初至不明显、传感器位置错误以及传感器与采集设备故障等均会导致异常TDOA的产生，进而造成定位结果的严重误差。此外，波速也是影响声发射源定位精度的一个重要因素，现有源定位方法大都需要预先准确测定的波速。而在真实工程环境下，波速受外部条件的影响很大，预测波速与真实波速之间往往存在较大的偏差并将引起严重的定位误差。因此，一种抗异常TDOA干扰且不受波速测量误差影响的声发射源定位方法值得进一步研究。

发明内容

本发明提供了一种异常TDOA存在下未知波速声发射源定位方法、系统及存储介质，以解决现有的声发射源定位方法在存在异常TDOA时难以得到精确定位结果的问题。

第一方面，提供了一种异常TDOA存在下未知波速声发射源定位方法，包括：

从多个传感器中选取一个传感器作为参考传感器，利用TDOA测量值和各传感器坐标构建非线性双曲线方程，并通过引入两个中间变量V、K，将非线性方程组化为线性方程组；其中，V＝v²,K＝vr₀，v表示声发射信号的传播速度，r₀表示声发射源到参考传感器的距离；

设置线性方程组中各方程权重为1，对线性方程组进行求解得到普通最小二乘解；根据普通最小二乘解计算方程残差；

根据该方程残差估计方程权重，然后根据方程权重对线性方程组中各方程进行加权并计算权重最小二乘解；根据新的定位结果重新计算方程残差；如此反复直到方程残差的变化小于预设精度时停止迭代；最后一次迭代得到的方程残差即为最终稳定的方程残差；

对最终稳定的方程残差取绝对值得到绝对残差，并由小到大进行排序；计算各绝对残差的间距、全局距离和局部距离来获取全局距离和局部距离判据；当某一绝对残差的间距同时满足全局距离和局部距离判据时，该间距即为分割正常残差与异常残差的间距；该间距对应的绝对残差为临界残差值，大于该绝对残差的值即为异常残差，异常残差对应的TDOA为异常TDOA，并将其排除；

利用剩余正常的TDOA测量值和相应的传感器坐标重新构造TDOA线性方程组，并获得不受异常TDOA影响的声发射源坐标。

进一步地，所述线性方程组为：

Aθ＝b

其中，

x、y、z分别表示声发射源坐标的三个分量，x₀、y₀、z₀分别表示参考传感器坐标的三个分量，x_i、y_i、z_i分别表示第i个传感器坐标的三个分量，t_i,0表示第i个传感器与参考传感器之间的TDOA测量值，且i＝1,2,…,n，n为传感器中非参考传感器的个数。

进一步地，所述普通最小二乘解为：

θ＝(A^TA)^-1A^Tb。

进一步地，所述方程残差为：

ε＝b-Aθ

其中，

ε_i表示线性方程组中第i个方程对应的方程残差，i＝1,2,…,n。

进一步地，所述根据方程残差估计方程权重的具体公式为：

其中，w_i表示线性方程组中第i个方程对应的方程权重，λ_m为绝对残差的中位数，并设置λ_m的最低下限为0.05max(|ε_i|)。

进一步地，所述权重最小二乘解为：

θ＝(A^TWA)^-1A^TWb

其中，W＝diag(w₁,w₂,…,w_n)，diag(·)表示生成一个以括号中元素为主对角线的对角矩阵。

进一步地，所述各绝对残差的间距为：

d_j＝|ε_j|-|ε_j-1|,j＝2,3,…,n

全局距离为：

其中，

局部距离为：

其中，

j＝3,4,…,n

k₂＝2。

进一步地，所述全局距离和局部距离判据为：

其中，k₂＝2，k₁的数值是根据在正态分布的方程残差集合中发现异常残差值的概率为P₀＝0.15来确定的。

第二方面，提供了一种异常TDOA存在下未知波速声发射源定位系统，包括数据处理模块；数据处理模块采用如上所述的异常TDOA存在下未知波速声发射源定位方法进行声发射源定位。

第三方面，提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序适于被处理器加载并执行如上所述的异常TDOA存在下未知波速声发射源定位方法。

为了便于进一步理解本发明的技术方案，下面具体叙述上述计算公式的推导过程。

假设一个三维定位系统包含n+1个传感器s_i(x_i,y_i,z_i)(i＝0,1,…,n)，声发射源位于o(x,y,z)，系统中声发射信号的传播速度为v。假设第一个传感器s₀最先接收到信号，将其视为参考传感器。那么传感器s_i与参考传感器s₀之间的非线TDOA方程为

其中，t_i,0表示传感器s_i与s₀之间的TDOA测量值，i＝1,2,…,n。

对公式(1)进行移项、平方并化简，得到

其中，V＝v²，

b_i＝(x_i-x₀)²+(y_i-y₀)²+(z_i-z₀)²，并且i＝1,2,…,n。

将公式(2)用矩阵表示为

Aθ＝b (3)

其中，

并且

在真实背景噪音下，TDOA测量值不可避免地包含误差，因此公式(3)是不一致的，在考虑方程权重相等(均为1)的情况下，我们可以求得声发射源坐标的普通最小二乘解为

θ＝(A^TA)^-1A^Tb (4)

在得到声发射源坐标之后，我们可以得到方程的初步残差为：

ε＝b-Aθ (5)

其中，

根据该残差我们可以得到方程的权重估计。我们直接将残差平方的倒数作为权重，但为了避免残差趋于零时，相应方程将被赋予过高的权重出现过拟合现象。我们将获得的权重设置一个上限，即残差下限λ_m

其中，λ_m为残差取绝对值后的中位数。此时将有半数方程权重相等，而当多个较小残差存在时，依然会出现过拟合现象，故对λ_m设置一个最低下限为0.05max(|ε_i|)。

然后根据获得的权重对方程(3)进行加权，得到权重最小二乘解

θ＝(A^TWA)^-1A^TWb (7)

其中，

根据新的定位结果重新计算方程残差；如此反复直到方程残差的变化小于预设精度时停止迭代；最后一次迭代得到的残差即为最终稳定的方程残差。接下来，我们根据获得的残差，来识别和排除异常TDOA。

由于异常TDOA会引起对应大的(异常)方程残差，因此我们只需确定异常残差即可确定异常TDOA。为确定异常残差，找到区分正常残差与异常残差之间的临界残差λ_L是关键，大于λ_L的残差即可确定为异常残差

发明人发现，当异常值存在时，正常残差往往较为集中，而异常残差则会远离正常残差集中的区域。根据这个特定点我们可以识别出异常残差。注意在异常TDOA检测过程中，需将方程残差取绝对值，以减小残差间的平均距离，一定程度上可以消除对异常TDOA的误判，尤其是当正负残差均匀分布且残差个数n较小时。

接下来将详细介绍如何识别出异常残差，首先将方程残差取绝对值得到的绝对残差并按照升序编号，计算相邻两残差之间的间隔距离d_j：

d_j＝|ε_j|-|ε_j-1|,j＝2,3,…,n (9)

接下来的讨论均基于绝对残差而言，为方便描述简称绝对残差为残差。如果残差中存在异常值，那么其与正常残差之间会存在一较大的距离间隔。为定量描述这一距离间隔，我们引入全局距离

和局部距离

如果距离d_j同时大于

与

的某一倍数那么d_j即为分割正常残差和异常残差的间距，如公式(10)：

其中k₂＝2，k₁的数值是根据在正态分布的残差集合中发现异常残差的概率为P₀＝0.15来确定的。如当残差个数为n＝20时，可以计算出k₁＝9.75。

为d_j之前的距离的加权平均，由下式得到：

为各距离间隔的权重，该权重随着残差序号的增大而减小，并且与残差个数n相关，

由下式得到：

其中，

为权重

的和，目的是对权重进行归一化处理。

局部距离

的处理方式与

相似，只是权值衰减速度更快，以突出其局部性

式中，

如果距离间隔d_j满足局部距离与全局距离判据，如公式(10)。那么d_j就是正常残差与异常残差之间的距离间隔。d_j对应的残差|ε_j|即为临界残差值λ_L，所有比λ_L更大的残差都将被视为异常残差，异常残差对应的TDOA测量值即为异常TDOA。为保证识别出TDOA异常值的个数不超过正常值的个数，我们可以仅分析j>0.5n的间距d_j。此外，由于

的权值下降速度快，故

的大小仅由d_j前边少数几个距离间隔所决定。而当这些距离都非常小、趋于0时，会高估距离的局部变化，基于这种情况

应调整为：

至此，公式(10)中的所有残差都得以确定。通过公式(10)可以识别并排除异常TDOA。然后我们利用剩余正常的TDOA数据并根据公式(4)重新计算声发射源坐标的普通最小二乘解，该解不受异常TDOA的影响。

有益效果

本发明提出了一种异常TDOA存在下未知波速声发射源定位方法、系统及存储介质，具有更加鲁棒的定位精度。具体体现在：(1)该方案可以识别和排除异常TDOA，降低了其对定位精度的影响；(2)该方案不需要预先确定波速，消除了波速误差对定位精度的影响；(3)该方案无需选择迭代初值，避免了迭代不收敛或收敛速度慢的问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的异常TDOA存在下未知波速声发射源定位方法的流程图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。

实施例1

本实施例提供了一种异常TDOA存在下未知波速声发射源定位方法，包括：

从多个传感器中选取一个传感器作为参考传感器，利用TDOA测量值和各传感器坐标构建非线性双曲线方程，并通过引入两个中间变量V、K，将非线性方程组化为线性方程组；其中，V＝v²,K＝vr₀，v表示声发射信号的传播速度，r₀表示声发射源到参考传感器的距离；

设置线性方程组中各方程权重为1，对线性方程组进行求解得到普通最小二乘解；根据普通最小二乘解计算方程残差；

根据该方程残差估计方程权重，然后根据方程权重对线性方程组中各方程进行加权并计算权重最小二乘解；根据新的定位结果重新计算方程残差；如此反复直到方程残差的变化小于预设精度时停止迭代；最后一次迭代得到的方程残差即为最终稳定的方程残差；其中预设精度根据实际情况选择，如本实施例中预设精度取10^-4；

对最终稳定的方程残差取绝对值得到绝对残差，并由小到大进行排序；计算各绝对残差的间距、全局距离和局部距离来获取全局距离和局部距离判据；当某一绝对残差的间距同时满足全局距离和局部距离判据时，该间距即为分割正常残差与异常残差的间距；该间距对应的绝对残差为临界残差值，大于该绝对残差的值即为异常残差，异常残差对应的TDOA为异常TDOA，并将其排除；

利用剩余正常的TDOA测量值和相应的传感器坐标重新构造TDOA线性方程组，并获得不受异常TDOA影响的声发射源坐标。

具体的，其中所述线性方程组为：

Aθ＝b

其中，

x、y、z分别表示声发射源坐标的三个分量，x₀、y₀、z₀分别表示参考传感器坐标的三个分量，x_i、y_i、z_i分别表示第i个传感器坐标的三个分量，t_i,0表示第i个传感器与参考传感器之间的TDOA测量值，且i＝1,2,…,n，n为传感器中非参考传感器的个数。

所述普通最小二乘解为：

θ＝(A^TA)^-1A^Tb。

所述方程残差为：

ε＝b-Aθ

其中，

ε_i表示线性方程组中第i个方程对应的方程残差，i＝1,2,…,n。

所述根据方程残差估计方程权重的具体公式为：

其中，w_i表示线性方程组中第i个方程对应的方程权重，λ_m为绝对残差的中位数，并设置λ_m的最低下限为0.05max(|ε_i|)。

所述权重最小二乘解为：

θ＝(A^TWA)^-1A^TWb

其中，W＝diag(w₁,w₂,…,w_n)，diag(·)表示生成一个以括号中元素为主对角线的对角矩阵。

所述各绝对残差的间距为：

d_j＝|ε_j|-|ε_j-1|,j＝2,3,…,n

全局距离为：

其中，

局部距离为：

其中，

j＝3,4,…,n

k₂＝2。

所述全局距离和局部距离判据为：

本实施例中，k₂＝2，k₁的数值是根据在正态分布的方程残差集合中发现异常残差值的概率为P₀＝0.15来确定的。如当残差个数为n＝20时，可以计算出k₁＝9.75。

本实施例提供的异常TDOA存在下未知波速声发射源定位方法的完整流程图如图1所示，I表示迭代次数。

实施例2

本实施例提供了一种异常TDOA存在下未知波速声发射源定位系统，包括数据处理模块；数据处理模块采用如实施例1所述的异常TDOA存在下未知波速声发射源定位方法进行声发射源定位。

实施例3

本实施例提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序适于被处理器加载并执行如实施例1所述的异常TDOA存在下未知波速声发射源定位方法。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

可以理解的是，上述各实施例中相同或相似部分可以相互参考，在一些实施例中未详细说明的内容可以参见其他实施例中相同或相似的内容。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。

实验验证：

假设一个边长为300mm的立方体定位系统，21个传感器分散且均匀地包围着该定位系统，他们的坐标分别为(150,150,300),(0,0,300),(0,150,300),(150,0,150),(0,0,150),(0,150,150),(300,0,300),(300,150,300),(300,150,150),(0,300,300),(150,300,150),(150,0,0),(0,300,150),(150,150,0),(0,0,0),(300,300,300),(300,300,150),(300,0,0),(150,300,0),(0,300,0),以及(300,300,0)，单位mm。传感器的具体编号按照接收到信号的先后顺序依次排序。真实的声发射源位于定位系统内，坐标为(88,266,211)(单位:mm)。假设该监测系统中声发射信号的真实传播速度为5000m/s。通过模拟的方式产生一组TDOA到时数据，并添加标准差为0.5×10^-6s的随机误差来模拟小的随机噪音误差。并向最大的两个TDOA数据中添加误差为±5×10^-6s的异常值来模拟异常TDOA数据。通过上述处理我们得到一组TDOA数据为：3.48,5.97,11.94,12.25,14.89,25.9,25.19,27.2,27.5,30.91,31.15,32.57,36.72,37.52,39.52,40.89,50.11,51.23,55.75,65.98,单位10^-6s。

根据上述实施案例中所述的步骤与方程式排除的异常值恰好是最大的两个TDOA测量值数据，这与模拟设置的异常TDOA完全吻合。进一步计算得到的声发射源坐标为(87.27,262.88,209.57)(单位:mm)，波速为4993.5m/s，这两个估计值与真实声源坐标(88,266,211)(单位:mm)和预设波速5000m/s吻合较好，因此本发明所提供的技术方案可以有效地识别并排除异常TDOA，具有较高的声发射源坐标与波速估计精度。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。

Claims (3)

1.一种异常TDOA存在下未知波速声发射源定位方法，其特征在于，包括：

从多个传感器中选取一个传感器作为参考传感器，利用TDOA测量值和各传感器坐标构建非线性双曲线方程，并通过引入两个中间变量V、K，将非线性方程组化为线性方程组；其中，V＝v²,K＝vr₀，v表示声发射信号的传播速度，r₀表示声发射源到参考传感器的距离；

设置线性方程组中各方程权重为1，对线性方程组进行求解得到普通最小二乘解；根据普通最小二乘解计算方程残差；

根据该方程残差估计方程权重，然后根据方程权重对线性方程组中各方程进行加权并计算权重最小二乘解；根据新的定位结果重新计算方程残差；如此反复直到方程残差的变化小于预设精度时停止迭代；最后一次迭代得到的方程残差即为最终稳定的方程残差；

对最终稳定的方程残差取绝对值得到绝对残差，并由小到大进行排序；计算各绝对残差的间距、全局距离和局部距离来获取全局距离和局部距离判据；当某一绝对残差的间距同时满足全局距离和局部距离判据时，该间距即为分割正常残差与异常残差的间距；该间距对应的绝对残差为临界残差值，大于该绝对残差的值即为异常残差，异常残差对应的TDOA为异常TDOA，并将其排除；

利用剩余正常的TDOA测量值和相应的传感器坐标重新构造TDOA线性方程组，并获得不受异常TDOA影响的声发射源坐标；

所述线性方程组为：

Aθ＝b

其中，

x、y、z分别表示声发射源坐标的三个分量，x₀、y₀、z₀分别表示参考传感器坐标的三个分量，x_i、y_i、z_i分别表示第i个传感器坐标的三个分量，t_i,0表示第i个传感器与参考传感器之间的TDOA测量值，且i＝1,2,…,n，n为传感器中非参考传感器的个数；

所述普通最小二乘解为：

θ＝(A^TA)^-1A^Tb；

所述方程残差为：

ε＝b-Aθ

其中，

ε_i表示线性方程组中第i个方程对应的方程残差，i＝1,2,…,n；

所述根据方程残差估计方程权重的具体公式为：

其中，w_i表示线性方程组中第i个方程对应的方程权重，λ_m为绝对残差的中位数，并设置λ_m的最低下限为0.05max(|ε_i|)；

所述权重最小二乘解为：

θ＝(A^TWA)^-1A^TWb

其中，W＝diag(w₁,w₂,…,w_n)，diag(·)表示生成一个以括号中元素为主对角线的对角矩阵；

所述各绝对残差的间距为：

d_j＝|ε_j|-|ε_j-1|,j＝2,3,…,n

全局距离为：

其中，

局部距离为：

其中，

j＝3,4,…,n

k₂＝2；

所述全局距离和局部距离判据为：

其中，k₂＝2，k₁的数值是根据在正态分布的方程残差集合中发现异常残差值的概率为P₀＝0.15来确定的。

2.一种异常TDOA存在下未知波速声发射源定位系统，其特征在于，包括数据处理模块；数据处理模块采用如权利要求1所述的异常TDOA存在下未知波速声发射源定位方法进行声发射源定位。

3.一种计算机可读存储介质，其存储有计算机程序，其特征在于，所述计算机程序适于被处理器加载并执行如权利要求1所述的异常TDOA存在下未知波速声发射源定位方法。

Images