CN117331024A - 基于全集到时差测量和残差估计的微震/声发射源定位方法和定位系统 - Google Patents

Abstract

本发明公开了一种基于全集到时差测量和残差估计的微震/声发射源定位方法。首先，根据N(N‑1)个全集到时差测量和N个传感器坐标构建未知波速的控制方程，并引入两个中间变量将其线性化；其次，根据对到时差测量误差的计算，得到线性到时差方程的残差，并估计权值；最后，通过引入正交投影算子分离中间变量，求解加权线性方程，得到加权线性最小二乘解，实现微震/声发射源定位。本发明的目的在于消减现有定位方法中预先波速测定和到时差的测量误差对定位精度的影响，该方法易于实现、定位精度高。

Description

基于全集到时差测量和残差估计的微震/声发射源定位方法 和定位系统

技术领域

本发明涉及无损检测技术领域，尤其涉及一种基于全集到时差测量和残差估计的微震/声发射源定位方法。

背景技术

微震/声发射源定位技术作为一种常用的无损监测方法，目前已广泛应用于航空、隧道、采矿、结构工程等领域。微震/声发射定位技术不仅可用于危险源预测以及风险评估，还可用于结构损伤机理的研究。为了提高微震/声发射源定位精度，学者们提出了许多定位方法，现有的微震/声发射源定位方法大多依赖于预先测定的波速。然而，由于工程环境和地质结构的复杂性，几乎不可能准确地测定波速。具体而言，随着工程任务的推进，工程环境的动态变化会引起应力和岩体结构的变化，从而导致区域平均波速实时变化。此外，由于地质结构的不均匀性和各向异性，不同方向的波速存在差异，在工程中几乎不可能测量出不同方向上所有路径的波速。若波速测量的偏差较大，将导致微震/声发射源定位误差增大。申请号为201910787751.5的发明专利提出了一种未知波速体系下声发射源线性定位方法和系统，该方法在一个三维监测系统中放置n个传感器，记录传感器坐标以及到时数据；通过声发射源坐标公式确定声发射源位置，并得到最优附加变量值，从而计算出介质波速和触发时刻。然而，由于现场环境噪声的影响，现有这种方法中到时差的测量误差几乎是不可避免的，通过控制方程的累计传递作用和平方放大效应，也会导致较大的定位误差。因此，一种消减波速和测量误差影响的新型微震/声发射源定位方法，仍需进一步的研究。

发明内容

针对现有技术的不足，本发明提供了一种基于全集到时差测量和残差估计的微震/声发射源定位方法，以消减现有定位方法中预先波速测定和到时差的测量误差对定位精度的影响，该方法易于实现、定位精度高。

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

一种基于全集到时差测量和残差估计的微震/声发射源定位方法，包括以下步骤：

步骤1、根据N(N-1)个全集到时差测量和N(N≥6)个传感器坐标构建未知波速的控制方程，并引入两个中间变量将其线性化；

构建未知波速的控制方程如下：

D_j-D_i＝vt_i,j,i≠j,i,j＝1,2,…,N

其中，N表示传感器的个数；D_i,D_j分别表示微震/声发射源与声发射传感器S_i,S_j之间的距离， (x,y,z)为微震/声发射源的坐标，(x_i,y_i,z_i)为传感器S_i的坐标，(x_j,y_j,z_j)为传感器S_j的坐标；t_i,j表示声发射传感器S_i接收到声发射信号的时间t_i与声发射传感器S_j接收到声发射信号的时间t_j之差，t_i,j＝t_j-t_i；v为未知的介质波速。

引入中间变量线性化后的控制方程如下：

其中，x_i,j＝x_j-x_i,y_i,j＝y_j-y_i，z_i,j＝z_j-z_i； V和R_i为中间变量，V＝v²，/>

步骤2、根据到时差测量误差计算线性到时差方程的残差，并据此估计方程权值；

到时差方程的残差表示如下：

其中，其中{·}^o表示变量{·}不含噪音干扰的值，n_i,j表示到时差测量中的噪音项，假设其服从高斯随机分布。

忽略平方项，将ε_k表示为矩阵形式为：

ε＝2VPn

其中，diag(·)表示一个以括号内元素为主对角线元素，不在主对角线上的元素全为0的N阶方阵；n为由n_i,j组成的列向量；/>其中，

表示θ⁽¹⁾中的第m个元素，m＝1,2,3。

θ⁽¹⁾为上述线性方程组关于x,y,z的最小二乘解：

θ⁽¹⁾＝(A^TA)^-1A^T(L-B(B^TP^⊥B)^-1B^TP^⊥L)

其中，上标-1表示矩阵的逆；上标T表示矩阵或向量的转置；待求解变量 P^⊥＝I-A(A^TA)^-1A^T，其中I为主对角线元素均为1的单位矩阵。

步骤3、通过引入正交投影算子分离中间变量，求解加权线性方程，得到微震/声发射源的加权线性最小二乘解。

用于分离中间变量的正交投影算子为：

P^⊥′＝I-B(B^TWB)^-1B^TW

其中，W为权重矩阵；W＝ψ^-1；ψ表示方程残差ε的协方差矩阵，具体表达式为E为数学期望求解符号，/>是关于n_i,j的协方差矩阵，为简单起见，假设到达时间的源信号是高斯随机过程，并且所有n_i,j的信噪比都是相同的，则

进一步地，所述线性方程的加权最小二乘解为：

θ⁽²⁾＝(A^TP^⊥′WP^⊥′A)^-1A^TP^⊥′WP^⊥′L

其中，待求解变量即为微震/声发射源定位结果。

另一方面，提供一种基于全集到时差测量和残差估计的微震/声发射源定位系统，包括数据处理模块；数据处理模块采用上述的微震/声发射源定位方法，基于N个传感器坐标和它们接收到声发射信号的时间数据计算微震/声发射源坐标，实现其定位。

进一步地，所述的一种基于全集到时差测量和残差估计的微震/声发射源定位系统，还包括在监测系统内随机布置的至少6个不全共面的传感器。

为了便于进一步理解本发明的技术方案，下面对上述定位方法原理进行进一步说明：

用(x_i,y_i,z_i)(i＝1,2,…,N)表示声发射传感器S_i的坐标，(x_j,y_j,z_j)(j＝1,2,…,N)表示声发射传感器S_j的坐标，(x,y,z)表示待求的微震/声发射源坐标,N表示传感器个数。用t_i表示声发射传感器S_i接收到声发射信号的时间，t_j表示声发射传感器S_j接收到声发射信号的时间。

首先，根据传感器S_i和S_j的坐标和其接收到声发射信号的时间建立其与声发射源之间的控制方程：

D_j-D_i＝vt_i,j,i≠j,i,j＝1,2,…,N (1)

其中，D_i,D_j分别表示微震/声发射源与声发射传感器S_i,S_j之间的距离，t_i,j表示声发射传感器S_i接收到声发射信号的时间t_i与声发射传感器S_j接收到声发射信号的时间t_j之差，t_i,j＝t_j-t_i；v为未知的介质波速。

为了线性化上述双曲方程，将式(1)平方并简化为：

其中，x_i,j＝x_j-x_i,y_i,j＝y_j-y_i，z_i,j＝z_j-z_i； V和R_i为中间变量，V＝v²，/>

将上述线性到时差方程(2)重新整理，可表示为Aθ+Bδ＝L的形式。其中，

然而，由于到时差测量误差的影响，方程两边是不相等的。方程残差的矩阵形式可表示为:

ε＝L-Aθ-Bδ (3)

将式(3)乘以正交算子P^⊥，去除有关震源的相关参数，用中间变量δ表示的新方程为：

ε′＝P^⊥(L-Bδ) (4)

其中，P^⊥＝I-A(A^TA)^-1A^T。

通过最小化线性方程(4)，得到中间变量δ的最小二乘解为：

其中，{v}^-1表示矩阵{·}的逆；上标T表示矩阵或向量的转置。

将式(5)代入式(3)，得到微震/声发射源坐标θ的最小二乘解为：

θ⁽¹⁾＝(A^TA)^-1A^T(L-B(B^TP^⊥B)^-1B^TP^⊥L) (6)接下来，考虑到实际工程场景中环境噪音的影响，到时差测量值t_i,j可表示为：

其中，{·}^o表示变量{·}不含噪音干扰的值，n_i,j为到时差测量中的噪音项，假设其服从高斯随机分布。

使用式(7)将式(2)中的t_i,j表示为得到线性方程中包含噪音项的残差为：

其中，

忽略式(8)中的平方项，可将ε_k表示为矩阵形式：

ε＝2VPn (9)

其中，diag(·)表示一个以括号内元素为主对角线元素，不在主对角线上的元素全为0的N阶方阵；n为由n_i,j组成的列向量；在实际工程中，参数V和P中的/>为未知量，上述线性方程的最小二乘解可用于近似x^o,y^o,z^o和/>则/> 其中，

表示θ⁽¹⁾中的第m个元素，m＝1,2,3。

在估计方程残差后，得到线性方程组的权重矩阵为：

W＝ψ^-1 (10)

其中，ψ为方程残差ε的协方差矩阵，具体表达式为E为数学期望求解符号，/>是关于n_i,j的协方差矩阵，可以用功率谱来表示，为简单起见，假设到达时间的源信号是高斯随机过程，并且所有n_i,j的信噪比都是相同的，则/>

然后，通过引入正交投影算子P^⊥′分离中间变量，求得线性方程的加权最小二乘解为：

θ⁽²⁾＝(A^TP^⊥′WP^⊥′A)^-1A^TP^⊥′WP^⊥′L (11)

其中，P^⊥′＝I-B(B^TWB)^-1B^TW，待求解变量即为最终微震/声发射源定位结果。

有益效果

本发明提出了一种基于全集到时差测量和残差估计的微震/声发射源定位方法，该方法易于实现、定位精度高。具体体现在：

(1)该方法使用冗余到时差测量作为新的约束，进一步提高了定位精度；

(2)该方法对到时差方程中的噪音项进行了估计，提高了定位方法对环境噪音所造成的随机误差的容忍性；

(3)该方法得到了微震/声发射源坐标唯一解析解，避免了计算复杂、不收敛的问题；

(4)该方法无需选择参考传感器，避免了到时误差对定位结果的影响；

(5)该方法可以动态反演平均波速，避免了波速测量误差对定位精度的影响。

附图说明

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

图2为本发明实施例所述方法与其它定位方法定位结果对比图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚，下面结合附图和具体实施例对本发明进行详细描述。

如图1所示，本实施例提供一种基于全集到时差测量和残差估计的微震/声发射源定位方法，包括以下步骤：

步骤1、根据N(N-1)个全集到时差测量和N(N≥6)个传感器坐标构建未知波速的控制方程，并引入两个中间变量将其线性化；

构建未知波速的控制方程如下：

D_j-D_i＝vt_i,j,i≠j,i,j＝1,2,…,N

其中，N表示传感器的个数；D_i,D_j分别表示微震/声发射源与声发射传感器S_i,S_j之间的距离， (x,y,z)为微震/声发射源的坐标，(x_i,y_i,z_i)为传感器S_i的坐标，(x_j,y_j,z_j)为传感器S_j的坐标；t_i,j表示声发射传感器S_i接收到声发射信号的时间t_i与声发射传感器S_j接收到声发射信号的时间t_j之差，t_i,j＝t_j-t_i；v为未知的介质波速。

引入中间变量线性化后的控制方程如下：

其中，x_i,j＝x_j-x_i,y_i,j＝y_j-y_i，z_i,j＝z_j-z_i； V和R_i为中间变量，V＝v²，/>

步骤2、根据到时差测量误差计算线性到时差方程的残差，并据此估计方程权值；

到时差方程的残差表示如下：

其中，其中{·}^o表示变量{·}不含噪音干扰的值，n_i,j表示到时差测量中的噪音项，假设其服从高斯随机分布。

忽略平方项，将ε_k表示为矩阵形式为：

ε＝2VPn

其中，diag(·)表示一个以括号内元素为主对角线元素，不在主对角线上的元素全为0的N阶方阵；n为由n_i,j组成的列向量；/>其中，

表示θ⁽¹⁾中的第m个元素，m＝1,2,3。

θ⁽¹⁾为上述线性方程组关于x,y,z的最小二乘解：

θ⁽¹⁾＝(A^TA)^-1A^T(L-B(B^TP^⊥B)^-1B^TP^⊥L)

其中，上标-1表示矩阵的逆；上标T表示矩阵或向量的转置；待求解变量

P^⊥＝I-A(A^TA)^-1A^T，其中I为主对角线元素均为1的单位矩阵。

步骤3、通过引入正交投影算子分离中间变量，求解加权线性方程，得到微震/声发射源的加权线性最小二乘解。

用于分离中间变量的正交投影算子为：

P^⊥′＝I-B(B^TWB)^-1B^TW

其中，W为权重矩阵；W＝ψ^-1；ψ表示方程残差ε的协方差矩阵，具体表达式为E为数学期望求解符号，/>是关于n_i,j的协方差矩阵，为简单起见，假设到达时间的源信号是高斯随机过程，并且所有n_i,j的信噪比都是相同的，则

具体的，所述线性方程的加权最小二乘解为：

θ⁽²⁾＝(A^TP^⊥′WP^⊥′A)^-1A^TP^⊥′WP^⊥′L

其中，待求解变量即为微震/声发射源定位结果。

本实施例提供一种基于全集到时差测量和残差估计的微震/声发射源定位系统，包括数据处理模块；数据处理模块采用上述的微震/声发射源定位方法，基于N个传感器坐标和它们接收到声发射信号的时间数据计算微震/声发射源坐标，实现其定位。

具体的，所述的一种基于全集到时差测量和残差估计的微震/声发射源定位系统，还包括在监测系统内随机布置的至少6个不全共面的传感器。

试验例：

假设一个尺寸为300×300×300mm的立方体定位系统，16个传感器的坐标分别为(0,0,0),(300,0,0),(0,300,0),(0,0,300),(300,0,300),(300,300,300),(0,300,300),(150,0,150),(300,150,150),(150,300,150),(0,150,150),(150,150,0),(150,150,300),(0,0,150),(300,300,150)(单位：mm)。同时，在传感器阵列内部随机生成100个声发射源，用于验证定位精度。在进行声发射源定位过程中波速为未知量。本次试验通过模拟的方法产生一组到时数据，在到时数据中添加4μs的高斯随机误差，模拟环境噪音对定位精度的影响。

利用上述实施例中所述步骤和公式进行计算，将该定位方法与LLSUV、USBM、NIUV及CAS四种定位方法进行比较，根据图2可以看出通过上述实施例所述方法进行定位的100个声发射源定位结果的平均绝对距离误差最小，说明本发明提供的技术方案可行，且定位精度更高。

以上所述的仅是本发明的实施例，方案中公知的具体结构及特性等常识在此未作过多描述。应当指出，对于本领域的技术人员来说，在不脱离本发明结构的前提下，还可以作出若干变形和改进，这些也应该视为本发明的保护范围，这些都不会影响本发明实施的效果和本发明的实用性。

Claims (4)

1.一种基于全集到时差测量和残差估计的微震/声发射源定位方法，其特征在于，包括以下步骤：

步骤1、根据N(N-1)个全集到时差测量和N(N≥6)个传感器坐标构建未知波速的控制方程，并引入两个中间变量将其线性化；

构建未知波速的控制方程如下：

D_j-D_i＝vt_i,j,i≠j,i,j＝1,2,…,N

其中，N表示传感器的个数；D_i,D_j分别表示微震/声发射源与声发射传感器S_i,S_j之间的距离， (x,y,z)为微震/声发射源的坐标，(x_i,y_i,z_i)为传感器S_i的坐标，(x_j,y_j,z_j)为传感器S_j的坐标；t_i,j表示声发射传感器S_i接收到声发射信号的时间t_i与声发射传感器S_j接收到声发射信号的时间t_j之差，t_i,j＝t_j-t_i；v为未知的介质波速。

引入中间变量线性化后的控制方程如下：

其中，x_i,j＝x_j-x_i,y_i,j＝y_j-y_i，z_i,j＝z_j-z_i； V和R_i为中间变量，V＝v²，/>

步骤2、根据到时差测量误差计算线性到时差方程的残差，并据此估计方程权值；

到时差方程的残差表示如下：

其中，其中{·}^o表示变量{·}不含噪音干扰的值，n_i,j表示到时差测量中的噪音项，假设其服从高斯随机分布。

忽略平方项，将ε_k表示为矩阵形式为：

ε＝2VPn

其中，diag(·)表示一个以括号内元素为主对角线元素，不在主对角线上的元素全为0的N阶方阵；n为由n_i,j组成的列向量；/>其中，

表示θ⁽¹⁾中的第m个元素，m＝1,2,3。

θ⁽¹⁾为上述线性方程组关于x,y,z的最小二乘解：

θ⁽¹⁾＝(A^TA)^-1A^T(L-B(B^TP^⊥B)^-1B^TP^⊥L)

其中，上标-1表示矩阵的逆；上标T表示矩阵或向量的转置；待求解变量

P^⊥＝I-A(A^TA)^-1A^T，其中I为主对角线元素均为1的单位矩阵。

步骤3、通过引入正交投影算子分离中间变量，求解加权线性方程，得到微震/声发射源的加权线性最小二乘解。

用于分离中间变量的正交投影算子为：

P^⊥′＝I-B(B^TWB)^-1B^TW

其中，W为权重矩阵；W＝ψ^-1；ψ表示方程残差ε的协方差矩阵，具体表达式为E为数学期望求解符号，/>是关于n_i,j的协方差矩阵，为简单起见，假设到达时间的源信号是高斯随机过程，并且所有n_i,j的信噪比都是相同的，则

2.根据权利要求1所述的一种基于全集到时差测量和残差估计的微震/声发射源定位方法，其特征在于，所述线性方程的加权最小二乘解为：

θ⁽²⁾＝(A^TP^⊥′WP^⊥′A)^-1A^TP^⊥′WP^⊥′L

其中，待求解变量即为微震/声发射源定位结果。

3.一种基于全集到时差测量和残差估计的微震/声发射源定位系统，其特征在于：包括数据处理模块；数据处理模块采用上述的微震/声发射源定位方法，基于N个传感器坐标和它们接收到声发射信号的时间数据计算微震/声发射源坐标。

4.根据权利要求8所述的一种基于全集到时差测量和残差估计的微震/声发射源定位系统，其特征在于，还包括在监测系统内随机布置的至少6个不全共面的传感器。