CN110632557A - 一种声发射源定位方法及系统 - Google Patents

Abstract

本发明公开了一种声发射源定位方法及系统，在监测系统中布置n+1个不共面的声发射传感器，其中n≥5；将其中一个声发射传感器作为参考传感器S₀，其坐标记为(x₀,y₀,z₀)；另外n个声发射传感器S_i的坐标记为(x_i,y_i,z_i)，i＝1,2,…,n；记录n+1个声发射传感器接收到声发射信号的时间；将非线性的声发射源距离时间方程化为线性方程，得到线性方程组，以各方程残差的权重平方和最小为目标，构建目标函数，并嵌入声发射参数之间的相关关系，得到声发射源坐标(x,y,z)。本发明易于实现，定位精度和计算效率高。

Description

一种声发射源定位方法及系统

技术领域

本发明涉一种声发射源定位方法及系统。

背景技术

声发射定位技术在岩石力学领域有着广泛的应用，精确确定声发射源位置对岩体损伤与裂纹扩展机理的研究以及岩爆发生位置的确定是至关重要的。为了提高声发射源定位精度，学者们提出了许多定位方法，但这些方法大多需要预先给定介质波速，大量实验和工程实践表明测量得到的介质波速往往与真实波速之间存在较大的偏差。比如在岩石力学测试过程中，随着加载过程中裂纹、空隙的不断扩展与贯通，介质平均波速也将不断地发生明显变化，特别是在测试后期，介质波速与预先测定的波速之间将存在巨大的误差。此外，在动态采矿开采过程中，由于巷道掘进、爆破振动等的影响，介质平均波速也将实时发生变化，并且预先爆破的地点不一定恰巧在真实岩爆发生的位置，因此预先测量的波速将存在较大的误差并最终影响定位结果的精度。

介于波速测量误差对定位精度的严重影响，近年来学者们提出了未知波速体系下声发射源定位方法，可以分为迭代和解析定位方法两种。由于无需预先测定波速，这些方法一定程度上提高了声发射源定位的精度。然而，他们仍然受到诸多因素的影响，其中迭代定位方法由于需要预先给定较为精确的迭代初值，从而大大限制了其应用，如果初值给定偏差较大，该方法很容易造成迭代不收敛或者收敛速度慢的现象。解析定位方法虽然可以避免上述问题，但仅适用于固定几何形状的声发射源定位问题，传感器需放置在特定位置，且定位精确并不高。

为解决上述问题，有必要提供一种能进一步提高定位精度的声发射源定位方法。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提出一种声发射源定位方法及系统，易于实现，定位精度和计算效率高。

本发明所提供的技术方案为：

一种声发射源定位方法，包括以下步骤：

在监测系统中布置n+1个不共面的声发射传感器(n+1个声发射传感器只需不全部共面，布置位置任意)，其中n≥5；

将其中一个声发射传感器作为参考传感器S₀，其坐标记为(x₀,y₀,z₀)；另外n个声发射传感器S_i的坐标记为(x_i,y_i,z_i)，i＝1,2,…,n；

记录n+1个声发射传感器接收到声发射信号的时间；将第i个声发射传感器S_i接收到声发射信号的时间与声发射信号的触发时间之差(声发射源到声发射传感器S_i的走时)记为t_i，i＝0,1,2,…,n；

构建以下线性方程组：

其中，L_i＝x_i ²+y_i ²+z_i ²-x₀ ²-y₀ ²-z₀ ²，(x,y,z)为声发射源的坐标，V和K为附加变量，V＝v²，K＝t₀V，v为介质平均波速；Δt_i＝t_i-t₀；

以线性方程组中各方程残差

的平方和最小为目标，构建目标函数；

求解目标函数，得到声发射源的坐标(x,y,z)，将其作为定位结果。

进一步地，所述目标函数包括第一目标函数，第一目标函数为：

其中，ε为误差向量，

待求解的变量

求解第一目标函数，得到关于θ的普通最小二乘解θ⁽⁰⁾：

θ⁽⁰⁾＝(A^TA)^-1A^TL

其中，矩阵右上角符号T和-1分别表示对矩阵转置以及对矩阵求逆；

根据θ⁽⁰⁾中的前三个元素确定声发射源的坐标(x,y,z)，将其作为初步定位结果。

进一步地，所述目标函数包括第二目标函数，第二目标函数为：

其中，W为权重矩阵，W≈(V²BQB)^-1，矩阵Q＝diag(0.5)+0.5，矩阵

diag(·)函数表示生成一个以括号内元素为主对角线元素，不在主对角线上元素全为0的n阶方阵；{·}^o表示变量{·}的不含有噪音干扰的值，并且

表示θ⁽⁰⁾的第j个元素，j＝1,2,3,4；

求解第二目标函数，得到关于θ的无约束加权最小二乘解，记为θ⁽¹⁾：

θ⁽¹⁾＝(A^TWA)^-1A^TWL；

根据θ⁽¹⁾的前三个元素确定声发射源的坐标(x,y,z)，将其作为定位结果。

进一步地，所述目标函数包括第三目标函数，第三目标函数为：

其中，ε′是一个表示θ⁽¹⁾不准确度的变量，

θ′为与声发射参数相关的一个待求参数，

表示θ⁽¹⁾的第j个元素，j＝1,2,3,4,5；W′为权重矩阵，W′＝(4B′C_ΔθB′)^-1，并且

其中x^o、y^o、z^o、K^o和V^o分别用

和

近似；C_Δθ为声发射参数的协方差矩阵，C_Δθ＝(A^TWA)^-1；

求解第三目标函数，得到关于θ′的加权最小二乘解：

θ′＝(A′^TW′A′)^-1A′^TW′L′

根据θ′求得修正后的声发射参数θ⁽²⁾：

其中，表示θ⁽²⁾的第j个元素，θ_j′表示θ′的第j个元素，j＝1,2,3,4,5；

根据θ⁽²⁾的前三个元素确定声发射源的坐标(x,y,z)的最优解，将其作为最终的定位结果。

本发明还公开了一种声发射源定位系统，包括数据处理模块；数据处理模块采用上述的定位方法，基于各个声发射传感器的坐标和它们接收到声发射信号的时间计算声发射源坐标(x,y,z)，实现声发射源的定位。

进一步地，还包括在监测系统中布置的n+1个不共面的声发射传感器，其中n≥5。

---------------------------------------------------------------------

以下对上述方法原理进行进一步说明。

1.将非线性方程化为线性

设声发射源位置为(x,y,z)，该声发射源被n+1个传感器所包围，他们的坐标为S_i(x_i,y_i,z_i),i＝0,1,…,n。声发射源的控制方程为：

公式(1)为非线性方程组，其中v为介质平均波速，t₀为声发射源到参考传感器S₀的走时；Δt_i为传感器S_i(i＝1,2,…,n)与参考传感器S₀之间接收到声发射信号的时间差，即到时差，Δt_i＝t_i-t₀，t_i为声发射源到传感器S_i的走时，i＝1,2,…,n，当i＝0时，Δt_i＝0。

方程两边同时乘以v，然后平方，得到：

(x_i-x)²+(y_i-y)²+(z_i-z)²＝v²(Δt_i+t₀)²,i＝0,1,…,n (2)

将i>1的方程减去i＝1的方程，可得：

其中，L_i＝x_i ²+y_i ²+z_i ²-x₀ ²-y₀ ²-z₀ ²，V＝v²，K＝t₀V，n≥6。附加变量V，K使得(3)式成为一个线性方程组。但是x,y,z,V,K五个声发射参数仍然受到公式(2)中i＝0的方程约束，即：

接下来将对声发射源坐标进行求解，由于TDOA定位测量值误差的影响(假设误差服从高斯分布)，公式(3)左右两侧不可能完全相等，它们之间的差值用ε来表示，并写作矩阵形式：

ε＝Aθ-L (5)

其中,

虽然可以通过上式直接求得关于θ的普通最小二乘解：

θ⁽⁰⁾＝arg minε^Tε＝(A^TA)^-1A^TL (6)

但是方程组中声发射参数(x,y,z,V,K)中五个未知参数仍然受到公式(4)的约束。

接下来，我们将讨论如何求解线性方程组(5)的加权最小二乘解，以及如何嵌入公式(4)的约束关系，从而得到更加精确的未知波速下声发射源坐标。

a)不考虑声发射参数相关约束的加权最小二乘解

首先，我们不考虑公式(4)中表述的约束关系。此时，使得线性方程组(5)中各方程残差加权平方和最小的无约束解为：

θ⁽¹⁾＝arg minε^TWε＝(A^TWA)^-1A^TWL (7)

其中，W为权重矩阵，本发明将线性方程组中各方程残差的协方差矩阵作为权重矩阵。接下来将说明如何对权重进行估计。

将变量{·}的不含有噪音干扰的值定义为{·}^o，TDOA测量值可以记为：

其中，Δe_i为服从高斯分布的TDOA定位误差，其矢量表达式为Δe。将该式代入公式(5)，并利用关系Δt_i＝t_i-t₀求得残差ε的值为：

ε＝VBΔe+0.5VΔe·Δe (9)

其中，

diag(·)函数表示生成一个以括号内元素为主对角线元素的方阵。因此，权重矩阵可以写成W＝[E(εε^T)]^-1，E(·)为期望，上标-1表示矩阵的逆。事实上Δe_i总是远小于

因此可以忽略(9)式等号右边的平方项，ε将变成高斯随机矢量，此时权重矩阵可近似表示为：

W≈(V²BQB)^-1 (10)

其中矩阵Q是TDOA定位测量值的协方差矩阵，其对角线元素为1，其他元素为0.5，如公式(11)所示。然而，矩阵B是不能直接获得的，需要在初始定位结果的基础上进行估计。初始定位的执行，一般假设各公式权重相等并根据式(6)计算声发射参数的普通最小二乘解θ⁽⁰⁾。

得到初步定位结果θ⁽⁰⁾后，到达时间

可由下式计算得到，进而得到矩阵B的估计值。

其中，

表示θ⁽⁰⁾中第j个元素，j＝1,2,3,4。

b)利用声发射参数的相关约束进行修正

接下来，通过考虑声发射参数x,y,z,V,K在公式(4)之间的约束关系，进一步推导一个更加精确的修正解。由于测量误差的存在，式(7)的估计值将会偏离真实值，因此可以表示为：

θ⁽¹⁾＝θ^o+Δθ (13)

声发射参数θ⁽¹⁾的偏差Δθ和协方差矩阵C_Δθ可以用扰动法近似求解：

C_Δθ＝E{ΔθΔθ^T}＝(A^TWA)^-1 (15)

其中，

由于矢量ε是服从零均值高斯分布的，因此Δθ也服从均值为零的高斯分布，矢量θ⁽¹⁾则是一个均值以真值为中心的无偏估计。

用θ⁽¹⁾的前三项分别减去x₀,y₀和z₀，然后求平方，θ⁽¹⁾的第四项保持不变，第五项平方后除以V^o，可以得到一个新的关于声发射参数的方程组：

ε′＝L′-A′θ′ (17)

其中，

表示θ⁽¹⁾的第j个元素，j＝1,2,3,4,5。ε′是一个表示θ⁽¹⁾不准确度的变量。将公式(13)代入公式(17)，省略二次项及分母扰动项得到：

将公式(18)写作矩阵形式：

ε′≈2B′Δθ (19)

其中，

由于偏差Δθ很小，该近似是有效的。由于Δθ服从均值为零的高斯分布，因此ε′也近似服从零均值高斯分布，相应的权重矩阵W′可以写成：

W′＝Ε[ε′ε′^T]^-1＝(4B′C_ΔθB′)^-1 (20)

在计算过程中x^o、y^o、z^o、K^o和V^o可用公式(7)的计算结果，即

和

近似代替。

对新的线性方程组(17)最小化权重残差平方和，可得到θ′的加权最小二乘估计为：

θ′＝arg minε′^TW′ε′＝(A′^TW′A′)^-1A′^TW′L′ (21)

根据θ′便可求得修正后的声发射参数θ⁽²⁾：

其中，

表示θ⁽²⁾的第j个元素，θ_j′表示θ′的第j个元素，j＝1,2,3,4,5。

θ⁽²⁾的前三个元素

即声发射源的坐标(x,y,z)的最优解，将其作为最终的定位结果。

该定位方法对求解公式进行了权重估计，并考虑了声发射参数之间的相关约束，进一步提高了声发射源定位精度。

上述方案针对未知波速体系下声发射源线性定位方法没有考虑声发射参数之间相关约束的问题，通过引入权重最小二乘原理进行解析求解，并嵌入声发射参数的相关约束进行修正，进一步提升了声发射源定位精度，具体操作可分为三个步骤：首先，将非线性的声发射源距离时间方程(非线性超定方程组)化为线性，然后引入权重最小二乘原理对线性方程组加权，并设置各方程初始权重相等，在不考虑声发射源参数相关约束的情况下计算声发射参数的普通最小二乘解θ⁽⁰⁾。其次，根据计算得到的声发射参数θ⁽⁰⁾进一步估计方程权重W，并得到无约束加权最小二成解θ⁽¹⁾。最后，利用声发射参数之间的相关约束构建新的线性方程组，并再次计算关于声发射源参数的加权最小二成解θ′进行修正，最终可以得到最优的声发射参数θ⁽²⁾。

有益效益：

本发明提供了一种声发射源定位方法及系统，其优势在于：

(1)该方法无需预先测定介质波速，可以实时对波速进行反演；由于可以实时反演波速，因此一定程度上避免了波速测量误差以及波速动态变化对定位结果的影响；

(2)该方法属于解析定位方法，无需选择迭代初值，避免了传统迭代方法计算不收敛以及运算效率低的问题，稳定性高，计算效率高；

(3)该方法对声发射传感器安装位置要求只需6个以上声发射传感器不同时共面，各声发射传感器无需安装在特定位置；

(4)该方法通过引入最小二乘原理可以综合利用所有声发射传感器的信息进行解析定位，提高了定位精度；

(5)该方法对线性方程组进行了权重估计，并在计算过程中考虑了声发射参数之间的相关约束，避免了现有解析定位方法没有考虑声发射参数之间的相关约束得到的定位坐标并不是最优的问题，进一步提高了声发射源定位的精确度。

附图说明

图1为本发明实施例定位方法流程图。

图2为本发明实施例声发射源三维定位示意图。

具体实施方法

以下结合附图和具体实施例对本发明进行进一步具体说明。

如图2所示为一个模拟的声发射定位系统，该系统由八个传感器S_i(x_i,y_i,z_i)(i＝0,1,…,7)组成，它们的坐标分别为(0,0,1),(100,0,0),(100,100,0),(0,100,0),(0,0,100),(100,0,100),(100,100,100),以及(0,100,100)(单位:mm)。预设一声发射源坐标位于坐标点O(40mm,50mm,60mm)，声发射监测系统中介质波速为v＝5000m/s，该参数在实际定位过程中被当作未知量，此处仅用于模拟产生到时数据，此外对每个到时数据添加标准差为0.1μs的误差扰动。通过模拟可以得到传感器S_i(i＝1,2,…,7)与参考传感器S₀之间的一组到时差数据，△t_i为2.3885,2.4431,0.4171,-2.0384,0.2654,0.3605以及-2.2143(单位:×10^-6s)。

以上述定位系统为例，对声发射源定位方法进行详细说明，具体步骤如下：1.首先根据传感器坐标S_i(x_i,y_i,z_i)以及到时差数据△t_i计算声发射源参数的普通最小二乘解θ⁽⁰⁾：

其中，

2.计算无约束加权最小二乘解θ⁽¹⁾：

其中，方程权重

Q＝diag(0.5)+0.5，

3.计算声发射源参数的最小二乘解θ′：

其中权重矩阵W′＝Ε[ε′ε′^T]^-1＝(4B′C_ΔθB′)^-1，并且

4.计算最优的未知波速体系下声发射源坐标θ⁽²⁾：

从上式可以得到，声发射源的定位结果(39.606mm,49.705mm,59.073mm)，与真实的声发射源坐标(40mm,50mm,60mm)非常接近，因此定位结果有效。

Claims (6)

1.一种声发射源定位方法，其特征在于，包括以下步骤：

在监测系统中布置n+1个不共面的声发射传感器，其中n≥5；

将其中一个声发射传感器作为参考传感器S₀，其坐标记为(x₀,y₀,z₀)；另外n个声发射传感器S_i的坐标记为(x_i,y_i,z_i)，i＝1,2,…,n；

记录n+1个声发射传感器接收到声发射信号的时间；将第i个声发射传感器S_i接收到声发射信号的时间与声发射信号的触发时间之差记为t_i，i＝0,1,2,…,n；

构建以下线性方程组：

其中，L_i＝x_i ²+y_i ²+z_i ²-x₀ ²-y₀ ²-z₀ ²，(x,y,z)为声发射源的坐标，V和K为附加变量，V＝v²，K＝t₀V，v为介质平均波速；Δt_i＝t_i-t₀；

以线性方程组中各方程残差的平方和最小为目标，构建目标函数；

求解目标函数，得到声发射源的坐标(x,y,z)，将其作为初步定位结果。

2.根据权利要求1所述的声发射源定位方法，其特征在于，所述目标函数包括第一目标函数，第一目标函数为：

其中，ε为误差向量，

待求解的变量

求解第一目标函数，得到关于θ的普通最小二乘解θ⁽⁰⁾：

θ⁽⁰⁾＝(A^TA)^-1A^TL

其中，

矩阵右上角符号T和-1分别表示对矩阵转置以及对矩阵求逆；

根据θ⁽⁰⁾中的前三个元素确定声发射源的坐标(x,y,z)，将其作为定位结果。

3.根据权利要求2所述的声发射源定位方法，其特征在于，所述目标函数包括第二目标函数，第二目标函数为：

其中，W为方程权重，W≈(V²BQB)^-1，矩阵Q＝diag(0.5)+0.5，矩阵

diag(·)函数表示生成一个以括号内元素为主对角线元素，不在主对角线上元素全为0的n阶方阵；{·}^o表示变量{·}的不含有噪音干扰的值，并且

表示θ⁽⁰⁾的第j个元素，j＝1,2,3,4；

求解第二目标函数，得到关于θ的无约束加权最小二乘解，记为θ⁽¹⁾：

θ⁽¹⁾＝(A^TWA)^-1A^TWL；

根据θ⁽¹⁾的前三个元素确定声发射源的坐标(x,y,z)，将其作为定位结果。

4.根据权利要求3所述的声发射源定位方法，其特征在于，所述目标函数包括第三目标函数，第三目标函数为：

其中，ε′是一个表示θ⁽¹⁾不准确度的变量，ε′＝A′θ′-L′，

θ′为与声发射参数相关的一个待求参数，

表示θ⁽¹⁾的第j个元素，j＝1,2,3,4,5；W′为权重矩阵，W′＝(4B′C_ΔθB′)^-1，并且

其中x^o、y^o、z^o、K^o和V^o分别用

和

近似；C_Δθ为协方差矩阵，C_Δθ＝(A^TWA)^-1；

求解第三目标函数，得到关于θ′的加权最小二乘解：

θ′＝(A′^TW′A′)^-1A′^TW′L′

根据θ′求得修正后的声发射参数θ⁽²⁾：

其中，

表示θ⁽²⁾的第j个元素，θ_j′表示θ′的第j个元素，j＝1,2,3,4,5；

由θ⁽²⁾的前三个元素确定声发射源的坐标(x,y,z)的最优解，将其作为最终的定位结果。

5.一种声发射源定位系统，其特征在于，包括数据处理模块；数据处理模块采用权利要求1～4中任一项所述的方法，基于各个声发射传感器的坐标和它们接收到声发射信号的时间计算声发射源坐标(x,y,z)，实现声发射源的定位。

6.根据权利要求5所述的基于随机传感器网络的声发射源定位系统，其特征在于，还包括在监测系统中布置的n+1个不共面的声发射传感器，其中n≥5。

Images