基于Lebenberg_Marquardt的声发射源定位算法

Sound Source Localization Algorithm Based on Levenberg_Marquardt

摘要: 声发射是一种寻找材料破裂点或地震发声源位置非常有效的方法，在本次实验中介绍了声发射源定位中时差定位法的基本原理，将定位问题转化为多元非线性方程组求解问题。根据问题特点，将传统的Geiger定位算法与Lebenberg_Marquardt算法相结合，通过改进Geiger算法的迭代方向与迭代步，极大地减少了算法的时间复杂度，提高了算法的精确度。数值仿真算例表明：该方法能有效地解决了迭代时雅可比矩阵与其转置矩阵接近奇异时d_k过大的问题，保证了算法的收敛性并且提高了迭代算法的收敛速度。

Abstract: Acoustic emission is an extremely effective method for locating material fracture points or seismic sound sources. In this experiment, the basic principle of time difference localization in acoustic emission source localization is introduced, which transforms the localization problem into a multi-dimensional nonlinear equation solving problem. Based on the characteristics of the problem, the traditional Geiger localization algorithm is combined with the Levenberg_Marquardt algorithm by improving the iteration direction and step of the Geiger algorithm, greatly reducing the time complexity of the algorithm and improving its accuracy. Numerical simulation examples demonstrate that this method effectively solves the problem of excessive d_k when the Jacobian matrix and its transpose matrix are close to singular during iteration, ensuring the convergence of the algorithm and improving the convergence speed of the iterative algorithm.

1. 引言

声发射(acoustic emission, AE) [1] [2] [3] [4] 作为一种研究材料破裂重要手段，在岩石力学应用领域应用十分广泛。如何精确地找到材料破裂点或地震发声源的位置，是声发射中的一个重要问题。在传统的声发射问题中，比较常用的算法有：最小二乘法、Bayesian定位方法、相对定位法、Geiger定位方法和单纯形定位方法等。

目前，国内外学者针对声发射定位方法开展了大量的研究工作，比较常用的有最小二乘法 [5] ，Geiger法 [6] ，单纯形法 [7] 以及牛顿法 [8] 等。其中，Geiger法是一个经典的定位技术，其中使用到地震到时(即地震波到达不同地震站的时间)进行震源位置的估计。林峰 [9] 等研究人员提出了一种联合定位方法，将线性定位方法和Geiger定位方法相结合，有效解决了线性定位精度不高的问题，并对常规Geiger定位方法的初值选择进行了优化。同时，Li Nan [10] 等人提出了一种基于单纯形法的非线性微震震源定位方法，显著提高了震源定位的精度。此外，李会义 [11] 等人利用牛顿线性迭代法成功求解定位方程组，实现了微震定位的可靠监测。传统的Geiger方法可能会使用迭代重新定权最小二乘方法进行求解。该过程中，可能会对方程或模型的某些部分进行简化或假设以便计算。在这种迭代类中常会遇到迭代步中雅可比矩阵与其转置矩阵的乘积接近于奇异的情况，在这种情况下会直接导致迭代步过大，从而使误差过大。

针对上述问题，本文引进了Levenberg_Marquardt算法 [12] 。作为一种结合了梯度下降和高斯–牛顿方法的优化技术，可以用来提高Geiger定位方法的效度。利用LM算法，可以更准确地进行震源位置的求解，因为它在数据不是完全线性的情况下仍然有效，并且提供了一种智能化的方式来调整步长以平衡收敛速度和稳定性。不仅提高了算法本身的鲁棒性，而且很大程度上降低了算法的时间复杂度。

2. 声发射源三维时差定位方法

2.1. Geiger算法

Geiger算法是声发射源定位中常用的一种定位算法。该算法基于声发射信号在不同传感器之间的传播时间差(time difference of arrival，TDOA)来确定声源的位置。Geiger算法的基本思想是通过将声源位置的坐标表示为未知变量，并将TDOA表示为已知常数，将定位问题转化为一个非线性方程组求解问题。该方程组的解即为声源的位置。

算法的核心是迭代过程，其中使用的最小二乘法用于拟合实测和计算得到的TDOA值，以更新声源位置。在每次迭代中，根据当前的声源位置估计值和TDOA信息，计算残差，并构建目标函数。然后，使用数值优化算法(如Levenberg_Marquardt算法)求解目标函数，并得到更新的声源位置估计值。重复这一过程，直到迭代收敛或达到预定的停止条件。Geiger算法在声发射源定位中具有较高的定位精度和稳定性，并且在实践中得到广泛应用。通过结合其他改进算法，如Levenberg_Marquardt算法，可以进一步提高算法的性能和收敛速度。

Geiger定位算法是从一个给定的初始点 $\theta \left(x,y,z,t,v\right)$ 通过迭代而接近最终结果。每一次迭代都基于最小二乘法计算一个修正向量 $\Delta \theta \left(\Delta x,\Delta y,\Delta z,\Delta t,\Delta v\right)$ ，把 $\Delta \theta$ 向量加到上次迭代的结果(试验点)上，得到一个新的试验点，然后判断这个新试验点是否满足要求；如果满足要求，此点坐标即所求震源位置，如果不满足，则继续迭代。

根据多个传感器获得的到达时间所建立的式(1)所给出的固定方程的到声发射源位置坐标。

$\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+{\left(z-z_i\right)}^2}=v\left(t-t_i\right)$ (1)

式中： $x_i,y_i,z_i$ 代表第i个传感器的坐标位置； $x,y,z$ 代表声发射源坐标；v代表声波波速； $t_i$ 代表P波到达第i个传感器的时间，t代表发声源发声的时间。

在(1)式中未知变量共有5个，理论上我们只要建立5个方程就可以求出声发射源位置坐标。或通过建立5个以上的方程得到超定方程组对声发射源进行定位。现在假设有n (n > 5)个传感器接收到有效信号： $\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+{\left(z-z_i\right)}^2}=v\left(t-t_i\right)$ ( $i=1,2,3,\cdots ,n$ )。即可建立需求解的非线性方程组。

我们令： $f_i=\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+{\left(z-z_i\right)}^2}-v\left(t-t_i\right)$ 。

则我们需求解的非线性方程组为：

$\{\begin{array}{c}{f}_1=0\\ f_2=0\\ ⋮\\ f_n=0\end{array}$

其中 $f_i\left(1,2,\cdots ,n\right)$ 为给定n维欧式空间 $R^n$ 中的实值函数。

引进向量记号，令

$F\left(l\right)=\left[\begin{array}{c}{f}_1\left(l\right)\\ f_2\left(l\right)\\ ⋮\\ f_n\left(l\right)\end{array}\right]$ $l=\left[\begin{array}{l}x\\ y\\ z\\ t\\ v\end{array}\right]$

则上式可写成： $F\left(l\right)=\text{0}$ 。

若存在 $l^{*}\in R$ ，使得 $F\left(l^{\text{*}}\right)=0$ ，则称 $l^{*}$ 是方程的解。

考虑非线性方程组：

$F\left(l\right)=\text{0}$

其中 $F\left(l\right):R^n\to R^m$ 是连续可微的及 $F{\left(l\right)}^{\prime }$ 连续可微。

设 $J\left(l\right)$ 是 $F{\left(l\right)}^{\prime }$ 的Jacobian矩阵

$J\left(l\right)=\left[\begin{array}{ccccc}\frac{\partial f_1\left(l\right)}{x}& \frac{\partial f_1\left(l\right)}{y}& \frac{\partial f_1\left(l\right)}{z}& \frac{\partial f_1\left(l\right)}{t}& \frac{\partial f_1\left(l\right)}{v}\\ \frac{\partial f_2\left(l\right)}{x}& \frac{\partial f_2\left(l\right)}{y}& \frac{\partial f_2\left(l\right)}{z}& \frac{\partial f_2\left(l\right)}{t}& \frac{\partial f_2\left(l\right)}{v}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{\partial f_n\left(l\right)}{x}& \frac{\partial f_n\left(l\right)}{y}& \frac{\partial f_n\left(l\right)}{z}& \frac{\partial f_n\left(l\right)}{t}& \frac{\partial f_n\left(l\right)}{v}\end{array}\right]=\left[\begin{array}{ccccc}\frac{x-x_1}{R_1}& \frac{y-y_1}{R_1}& \frac{z-z_1}{R_1}& -v& t_1-t\\ \frac{x-x_2}{R_2}& \frac{y-y_2}{R_2}& \frac{z-z_2}{R_2}& -v& t_2-t\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ \frac{x-x_n}{R_n}& \frac{y-y_n}{R_n}& \frac{z-z_n}{R_n}& -v& t_n-t\end{array}\right]$

利用泰勒公式，我们有：

$F\left(l\right)=F\left(l_0\right)+J\left(l_0\right)\left(l-l_0\right)+o\left({\left(l-l_0\right)}^2\right)$ (2)

一般 $l_0$ 充分接近 $l^{*}$ ，略去上式的高阶小量，因而可用线性方程组

$F\left(l\right)=F\left(l_0\right)+J\left(l_0\right)\left(l-l_0\right)=\text{0}$ (3)

代替(1-1)，设(1-3)的解为 $l_1$ ，则

$l_1=l_0-{\left[J\left(l_{\text{0}}\right)\right]}^{-1}F\left(l_{\text{0}}\right)$ (4)

$l_{k+1}=l_k-{\left[J\left(l_k\right)\right]}^{-1}F\left(l_k\right)$ (5)

计算迭代步为：

$d_k=-J_k^{-1}{F}_k$ (6)

其中 $F_k=F\left(l_k\right)$ ， $J_k=F^{\prime }\left(x_k\right)$ 。

我们将问题(1-1)的求解转化为求解

$\underset{x\in R^4}{\min }\frac{1}{2}{\Vert F\left(l\right)\Vert }^2$ (7)

可知Geiger算法是牛顿法在非线性最小二乘法问题上的应用，其搜索方向为：

$d_k=-{\left(J_k^{\text{T}}{J}_k\right)}^{-1}{J}_k^{\text{T}}{F}_k$ (8)

故问题(1-1)的求解转化为：

$l_{k+1}=l_k-{\left[J_k^{\text{T}}{J}_k\right]}^{-1}{J}_k^{\text{T}}F\left(l_k\right)$ (9)

2.2. Levenberg_Marquardt算法

LM算法是一种非线性优化方法，处于牛顿法和梯度下降法之间，对于过参数化问题不敏感，能有效处理冗余参数问题，从而降低代价函数陷入局部极小值的可能性。因此，LM算法在计算机视觉等领域有着广泛的应用。

在上述的最小二乘问题中，由于我们所省略的二阶项对Geiger算法在收敛存在上有一定的影响。为了得到更加有效的方法，在运用Geiger算法时加上搜索策略即：

$x_{k+1}=x_k+{\alpha }_k{d}_k$ (10)

其中 ${\alpha }_k$ 为搜索因子， ${\alpha }_k$ 利用Wolfe-Powell准则确定

$\{\begin{array}{l}{\Vert F\left(l_k+{\alpha }_k{d}_k\right)\Vert }^2\le {\Vert F\left(l_k\right)\Vert }^2+\sigma {\alpha }_k{F}^{\text{T}}\left(l\right)d_k\\ {\Vert F\left(l_k+{\alpha }_k{d}_k\right)\Vert }^2\ge {\Vert F\left(l_k\right)\Vert }^2+\left(1-\sigma \right){\alpha }_k{F}^{\text{T}}\left(l\right)d_k\end{array}$ (11)

其中 $\sigma \in \left(0,\frac{1}{2}\right)$ 。

Geiger算法对于绝大多数的三维声发射问题中的最小二乘是全局收敛的，但对于某些问题仍然可能收敛很慢。且进行迭代求解时当迭代步中 $J_k^{\text{T}}{J}_k$ 接近于奇异时 $d_k$ 过大，而且当 $J_k^{\text{T}}{J}_k$ 是奇异时，其迭代步是没有意义的。

因此使用LM算法对Geiger算法进行修正，对于当前迭代点LM算法搜索方向为:

$d_k\left({\mu }_k\right)=-{\left(J_k^{\text{T}}{J}_k+{\mu }_{k}I\right)}^{-1}{J}_k^{\text{T}}{F}_k$ (12)

其中 ${\mu }_k$ 为正参数，LM是对G-N迭代步的修正。 ${\mu }_k$ 是为了防止 $J_k^{\text{T}}{J}_k$ 接近奇异时 $d_k$ 过大。而且当 $J_k^{\text{T}}{J}_k$ 是奇异时。迭代步是没有定义的。LM克服了迭代对矩阵 $J\left(l\right)$ 必须满秩的要求。

在得到修正向量 $d_k$ 后，以 $x_{k+1}=x_k+{\alpha }_k{d}_k$ 为新的实验点继续迭代，直到满足误差要求。

3. 基于Levenberg_Marquardt的声发射源定位算法

3.1. 数值仿真

数据来源为康玉梅等所做的实验仿真。为了定位声发射源在试件中的位置，需要根据虚拟声发射源位置、传感器坐标和材料波速生成一组模拟到时数据，并通过处理这组数据实现声发射源的定位。为了考察定位结果的重复性，在声发射源位置生成了三组带有随机误差的到时数据。由于一定量的到时误差与传感器布局有关，因此会影响定位结果的准确性。

假设试件尺寸为100 mm × 100 mm × 200 mm，模拟试样的传感器布置和模拟声发射源的位置如图1所示。以点A为原点，依此可确定传感器位置坐标，详见表1。试验中的破裂点为点M (30, 40, 50)和点C (100, 100, 0)，破裂时间为100秒。

Figure 1. Schematic diagram of the location of sensors and acoustic emission sources

图1. 传感器和声发射源的位置示意图

Table 1. Coordinates of the locations of various sensors

表1. 各个传感器的位置坐标

表2为各个传感器的真实到时以及含随机误差的模拟到时数据。根据不同的时差定位方法，可以计算出不同的声发射源位置。

Table 2. Actual arrival times of each sensor and data on arrival times with three sets of random errors

表2. 各个传感器的真实到时和三组随机误差到时数据

3.2. 算法的实现

针对有8个传感器的实验环境中，分别使用最小二乘法、Geiger法以及修正算法进行计算实验。表3为计算出的各个定位点的位置、计算误差以及迭代次数。从表3中还可以看出，Levenberg_Marquardt算法对Geiger算法进行修正后明显减少了迭代步数，使收敛速度更快，这对于声发射实时监测的海量数据的处理是非常有意义的。

Table 3. Actual arrival times of each sensor and data on arrival times with three sets of random errors

表3. 各个传感器的真实到时和三组随机误差到时数据

4. 结语

本文主要介绍了三维声发射定位的原理，以及以Gergir算法为基础，使用Levenberg_Marquardt算法对其迭代步进行了修正，一方面不仅仅在精度上相较于普通的定位算法得到了提高，而且相较于传统最小二乘法的迭代速度也得到了提高，另一方面解决了迭代中容易出现的迭代步过大的问题。很好地降低了算法的时间复杂度，提高了算法本身的鲁棒性。

NOTES

^*通讯作者。