Next Article in Journal
Improving Soil Water Content and Surface Flux Estimation Based on Data Assimilation Technique
Next Article in Special Issue
Decadal Continuous Meteor-Radar Estimation of the Mesopause Gravity Wave Momentum Fluxes over Mohe: Capability Evaluation and Interannual Variation
Previous Article in Journal
CloudSatNet-1: FPGA-Based Hardware-Accelerated Quantized CNN for Satellite On-Board Cloud Coverage Classification
Previous Article in Special Issue
Features of Winter Stratosphere Small-Scale Disturbance during Sudden Stratospheric Warmings
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Infrasound Source Localization of Distributed Stations Using Sparse Bayesian Learning and Bayesian Information Fusion

1
College of Logistics Engineering, Shanghai Maritime University, Shanghai 201306, China
2
Institute of Vibration, Shock and Noise, State Key Laboratory of Mechanical System and Vibration, Shanghai Jiao Tong University, Shanghai 200240, China
3
Northwest Institute of Nuclear Technology, Xi’an 710024, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(13), 3181; https://doi.org/10.3390/rs14133181
Submission received: 20 April 2022 / Revised: 24 June 2022 / Accepted: 29 June 2022 / Published: 2 July 2022
(This article belongs to the Special Issue Infrasound, Acoustic-Gravity Waves, and Atmospheric Dynamics)

Abstract

:
The precise localization of the infrasound source is important for infrasound event monitoring. The localization of infrasound sources is influenced by the atmospheric propagation environment and infrasound measurement equipment in the large-scale global distribution of infrasound arrays. A distributed infrasound source localization method based on sparse Bayesian learning (SBL) and Bayesian information fusion is proposed to reduce the localization error. First, the arrival azimuth of the infrasound source is obtained based on the SBL algorithm. Then, the infrasound source localization result is obtained by the Bayesian information fusion algorithm. The localization error of the infrasound source can be reduced by this infrasound source method, which incorporates the uncertainty of the infrasound propagation environment and infrasound measurement equipment into the infrasound source localization results. The effectiveness of the proposed algorithm was validated using rocket motor explosion data from the Utah Test and Training Range (UTTR). The experimental results show that the arrival azimuth estimation error can be within 2° and the localization distance error is 3.5 km.

Graphical Abstract

1. Introduction

The infrasound is a sound wave characterized by low frequencies (below 20 Hz) and long wavelengths. The infrasound waves are mainly divided into two categories according to their sources: natural infrasound sources and artificial infrasound sources [1]. The natural infrasound sources are caused primarily by tsunamis, earthquakes, storms, etc. The artificial infrasound sources are mainly caused by nuclear explosions, other large explosions, and rocket launches. The infrasound is used by the comprehensive nuclear-test-ban treaty organization (CTBTO) to detect infrasound events in the atmospheric environment, such as explosions and nuclear tests. The global infrasound monitoring network IMS (infrasound monitoring system) was established by the CTBTO to monitor nuclear explosions, earthquakes, and atmospheric climate [1]. The infrasound source localization technology is mainly used to monitor volcanic eruptions [2,3,4,5], earthquakes [6,7], thunderstorms [1] and other geological disasters [8].
The localization of the infrasound source is obtained by fusing the results of the arrival azimuth estimation from multiple infrasound stations in the context of the large-scale global distribution [6]. The final infrasound source localization results are influenced by the accuracy of the arrival azimuth estimation of individual infrasound stations. A conventional estimation algorithm for infrasound arrival azimuth is the frequency-wave number analysis (FK analysis) [9]. The power spectral density distribution of the infrasound signal between different slownesses and azimuths is calculated, and the azimuth corresponding to the value with the highest power is the result of the arrival azimuth estimation. The FK analysis can only handle a short infrasound signal time window (a few seconds). This method deals with large time windows because the window may contain multiple interference phases on different slow vectors, making it difficult to accurately identify the arrival azimuth information of the infrasound source [9].
The sparse Bayesian learning (SBL) algorithm was proposed in a probability framework for solving regression and classification problems to obtain sparse solutions [10,11,12] and was then applied to the field of acoustic array signal processing for direction-of-arrival (DOA) estimation [13,14,15,16]. Many advantages are captured by the SBL algorithm. The covariances of the weights are estimated directly by the SBL algorithm, thus greatly reducing the number of parameters to be estimated. The sparsity is automatically estimated by the SBL algorithm without any user input [17]. Better performance can be shown by the SBL algorithm with fewer snapshots, compared to the multiple signal classification (MUSIC) algorithm [13]. It is also less sensitive to array geometry, which usually requires a uniform linear array (ULA). The incoming wave azimuth of infrasound is sparse for the angularly discretized infrasound source plane, so the SBL algorithm is suitable for estimating the incoming wave azimuth of infrasound sources with sparsity and thus obtaining sparse solutions. The potential to improve the accuracy of the infrasound source arrival azimuth estimation is possessed by the SBL algorithm based on the above advantages.
The infrasound stations are usually distributed and arranged in outdoor environments due to the strict requirements of the infrasound stations for outdoor environmental conditions. The arrival azimuth, rather than the exact location of the infrasound source, can be obtained by individual infrasound stations. The estimation error of the arrival azimuth of a single infrasound station can be reduced by the arrival azimuth estimation algorithm for infrasound. The next consideration is how to fuse the arrival azimuth estimations from multiple infrasound stations, the uncertainty of the propagation environment, and the uncertainty of the measurement equipment in order to obtain more accurate results of infrasound source localization. Accurate infrasound source localization results are difficult to obtain because of the following uncertainties [6]: (1) the uncertainty of the infrasound propagation environment, such as temperature, wind speed, and atmospheric pressure [18]; and (2) uncertainties in the infrasound measurement equipment, such as the size of the array aperture, the sampling frequency, and the wavelength of the infrasound station [19]. Multiple types of information can be fused by the Bayesian information fusion algorithm [19] based on the Bayesian statistical model. The posterior probability density function and Bayesian credibility contours are obtained by the algorithm through Bayesian inference. Therefore, the SBL algorithm and Bayesian information fusion algorithm are combined to improve the accuracy of the infrasound localization.
The purpose of this paper is to improve the accuracy of infrasound source localization in the presence of the atmospheric uncertainty of infrasound propagation and the uncertainty of infrasound measurement equipment. An infrasound source localization method based on the SBL algorithm and Bayesian information fusion is proposed in this paper. First, the infrasound source arrival azimuth estimation for the individual infrasound station is acquired by the SBL algorithm. Second, the arrival azimuth estimation results of multiple stations, the uncertainty of the propagation environment, and the uncertainty of the measurement equipment are fused by the Bayesian information fusion algorithm into the infrasound source localization results. The rocket motor explosion data obtained in the Utah Test and Training Range (UTTR) are used to verify the feasibility of the algorithm. The main contributions of this paper are as follows:
  • An infrasound source localization method is developed based on the SBL algorithm under the plane wave assumption of infrasound propagation.
  • The arrival azimuth estimation results, the uncertainty of the infrasound propagation environment, and the uncertainty of the measurement equipment are fused by the Bayesian information fusion algorithm to obtain the localization, making the infrasound source localization results robust to the actual atmosphere and measurement environment.
  • The validity and feasibility of the proposed method are verified using a rocket motor explosion at the Utah Test and Training Range (UTTR).
The rest of this paper is organized as follows. In Section 2.1, The tau-p model for the infrasound propagation and problem statement is given. In Section 2.2, the method for estimating the arrival azimuth of infrasound sources by one arrayed station based on sources by one arrayed station based on the SBL algorithm is introduced. In Section 2.3, the Bayesian information fusion for infrasound source localization by using multiple stations is introduced. In Section 2.4, the technical flowchart for the infrasound source localization algorithm is introduced. In Section 3, the validity and feasibility of the algorithm are verified using measured infrasound data. The SBL algorithm is compared with the conventional infrasound arrival azimuth estimation algorithm. In Section 4, the availability of the proposed infrasound source localization algorithm is discussed.

2. Materials and Methods

2.1. The Tau-p Model for the Infrasound Propagation and Problem Statement

2.1.1. The Tau-p Model for the Infrasound Propagation

The tau-p model [18,20] is commonly used to simulate the trajectory of infrasound propagation in the atmosphere. It is a forward model of infasound propagation, which is suitable for visualizing the process of infrasound propagation. The principle of tau-p model is as follows.
The ray parameter p in tau-p model is
p = sin ( ϕ ) c 0 1 + sin ( ϕ ) u 0 c 0 1 ,
where (1) c 0 is the static sound velocity at the receiving point, (2) u 0 is the horizontal wind speed along the infrasound propagation direction at the receiving point, and (3) ϕ is the elevation angle of the emission.
The range along the ray direction in a phase loop is R. The characteristic function is ψ , which is formulated as
R ( z , p ) = 2 z 0 z ( p ) ψ ( z , p ) p ( 1 u ( z ) p ) + u ( z ) c ( z ) 2 d z ,
ψ ( z , p ) = 1 c ( z ) 2 p 2 ( 1 u ( z ) p ) 2 1 / 2 ,
where (1) z 0 is lower limit of the integration (usually zero or the surface height of the infrasound source); (2) z ( p ) is the upper limit of the integration, i.e., the maximum height of the infrasound propagation; (3) ψ ( z , p ) is the characteristic function for which the maximum height z m a x of the infrasound trajectory can be obtained by rooting the function; (4) c ( z ) = γ k T m is the adiabatic sound speed; and (5) u ( z ) is the horizontal wind speed along the propagation direction.
The transverse offset Q of the infrasound propagation model is
Q ( z , p ) = z o z ( p ) 1 c ( z ) 2 ψ ( z , p ) v ( z ) d z ,
where v ( z ) represents the horizontal wind speed in the vertical propagation direction, and the root of ψ ( z , p ) represents the turning point of the propagating sound line.
The diagram of the infrasound propagation model is shown in Figure 1. The infrasound propagation model is derived from the classical Wentzel–Kramers–Brillouin ( WKB ) sound line theory, which suggests that the sound line will turn when the horizontal phase velocity V θ matches the background effective sound velocity c ( z ) + u ( z ) . The propagation trajectory of the tau-p model in a phase loop is shown in Figure 1. The tau-p is a forward model, which is not applicable to the infrasound source localization problem. For this reason, infrasound is assumed to be a plane wave in the following. The infrasound source localization problem is only concerned with the source of infrasound, not with the propagation trajectory of infrasound.

2.1.2. Problem Statement for the Infrasound Source Localization

The arrival azimuth information of infrasound can be obtained from a single infrasound station, while the location information of the infrasound source is not available. From the above tau-p model, it can be seen that the accuracy of the arrival azimuth estimation of a single station is difficult to be guaranteed due to the complexity of infrasound propagation. The arrival azimuth estimation results of individual infrasound stations greatly influence the final infrasound localization results. Accurate arrival azimuth estimations are the basis for the subsequent fusion of information from multiple infrasound stations to obtain accurate infrasound source localization results.
In addition, the infrasound stations are usually distributed in the field. The localization error of the infrasound source is also affected by the uncertainty of the infrasound propagation environment and the uncertainty of the infrasound measurement equipment.

2.2. Estimating the Arrival Azimuth of Infrasound Source Based on the SBL Algorithm

The infrasound forward propagation model, i.e., the tau-p model, is introduced in Section 2.1. The trajectory of infrasound propagation can be visualized by this model. However, this model is not applicable to the arrival azimuth estimation algorithm of infrasound sources. A plane wave assumption is adopted for infrasound waves due to the ultralong distance propagation of infrasound waves, and a signal model for infrasound propagation is constructed. The infrasound source localization algorithm based on the SBL algorithm is proposed on the basis of the constructed signal model of infrasound propagation.

2.2.1. Signal Model for Arrival Azimuth Estimation Based on Plane Wave Assumption

It is difficult to infer the infrasound source signal from the measured signal of the infrasound station and thus obtain the arrival azimuth of the infrasound source. The above problem can be solved by building the signal model in the Equation (5) to obtain the infrasound source arrival azimuth estimation results. In this model, infrasound is assumed to be a plane wave due to the far-field propagation characteristics of infrasound in the atmosphere. The schematic diagram of three-dimensional infrasound propagation under the assumption of plane waves is shown in Figure 2. P is the acoustic signal when the initial time is t and the propagation distance is r . The horizontal azimuth information α is contained in the unit normal to the wavefront u .
The infrasound source signal in Equation (5) is linked to the infrasound measurement signal by the signal model through a linear regression model. The measured infrasound signal is Fourier transformed to obtain a multi-snapshot of the measured data Y = y 1 , , y L C N × L (N is the number of infrasound sensor and L is the number of snapshots) since the infrasound signal monitored by the infrasound station is a non-smooth signal. Large time windows of infrasound signals can be processed by this infrasound source localization algorithm using a multi-snapshot approach.
Y = A X + N ,
where X = x 1 , , x L C M × L are the actual infrasound source amplitudes and x m l denotes the infrasound source amplitude at the m-th discretization angle under the plane wave assumption (e.g., θ m = 90 + m 1 M 180 ) and at the l-th snapshot; M is the number of discretizations in the plane of the infrasound incoming wave direction (the plane ranges from [−90°, 90°]) discretized by the angular resolution Δ θ ; and N = n 1 , , n L C N × L is additive noise, which is independent across the infrasound sensors and snapshots L. The additive noise N of each infrasound sensor is assumed to follow a zero-averaged circularly symmetric complex Gaussian. A = a 1 , , a M C N × M is the transfer matrix connecting the infrasound station and the infrasound source. The infrasound array steering vectors a m are contained in the transfer matrix for all hypothetical arrival azimuths of infrasound as columns. The form of a m is as follows:
a m = 1 , e j 2 π d λ sin θ m , , e j 2 π ( N 1 ) d λ sin θ m T for m = 1 , 2 , M ,
where d is the infrasound array aperture; λ is the infrasound wavelength; and θ m is the m-th discretized angle. The infrasound signal model is constructed under the following assumptions: (1) the infrasound source vector x l and the additive noise n l are independent for the l-th snapshot; and (2) the noise n l and the infrasound source vector x l are assumed to be Gaussian, independently and identically distributed across all snapshots for l = 1 , , L .
Solving Equation (5) is an underdetermined problem due to the relationship between M and N as M N . The infrasound source vector x l is K-sparse (K denotes the number of infrasound sources), usually K M . The l-th activity set is defined as
M l = m N x m l 0 .
and M l = M = m 1 , , m K is assumed to be constant throughout the snapshot l. A M C N × K is defined to contain only the K “active” columns of A , i.e., the arrival azimuth corresponding to the index of the K non-zero values is the arrival azimuth of the infrasound signal source.
The arrival azimuth estimation of the infrasound source X from the infrasound measurement signal Y is a very challenging problem for Equation (5). The existence, uniqueness, and stability of the solution of Equation (5) are difficult to guarantee since Equation (5) is underdetermined. The following part of sparse Bayesian learning is introduced in Section 2.2.2 to solve this underdetermined problem so that the arrival azimuth of the infrasound source is obtained.

2.2.2. The SBL Algorithm for Infrasound Source Arrival Azimuth Estimation

The Bayesian formulation is used to solve linear problems (Equation (5)). The estimated parameters are considered variables of some prior distribution by the SBL algorithm. Then, this prior distribution is determined using the available knowledge. Finally, the posterior probabilities of the unknown parameters are inferred accordingly by the Bayesian rule.
(1)
Stochastic likelihood for the SBL algorithm
The additive noise in the Equation (5) is assumed to be a complex Gaussian data likelihood. The likelihood can then be written as
Y X ; σ 2 = exp 1 σ 2 Y A X F 2 π σ 2 N L .
(2)
Prior (prior distribution) on the infrasound sources
The complex infrasound source amplitude x m l is assumed to be independent in all snapshots and all discrete infrasound source arrival azimuths, and it follows a circularly symmetric complex Gaussian with zero-mean and dependent discrete infrasound source arrival azimuth variance γ m , m = 1 , , M ,
x m l ; γ m = δ x m l , for γ m = 0 1 π γ m e x m l 2 / γ m , for γ m > 0 ,
[ X ; γ ] = l = 1 L m = 1 M x m l ; γ m = l = 1 L N C ( x m l ; 0 , Γ ) ,
where N C is the normal distribution sign, and the unknown covariance matrix Γ for the infrasound source vector x l is assumed to be diagonal
Γ = diag ( γ ) , γ = γ 1 , , γ M .
The infrasound source power is denoted as a diagonal element of Γ , i.e., the hyperparameter γ 0 . The variance γ m 0 means that the arrival azimuth of this discrete infrasound source has an infrasound signal. When the variance γ m = 0 , then x m l = 0 with probability 1. Therefore, the sparsity of the model (Equation (5)) is controlled with the hyperparameters γ , as rank ( Γ ) = K M and the covariance Γ is estimated finally by the SBL algorithm.
(3)
Posterior (posterior distribution) on the infrasound sources
The posterior probability density function is solved by the Bayesian formula on the basis of the likelihood and prior
X Y ; γ , σ 2 Y X ; σ 2 [ X ; γ ] Y ; γ , σ 2 ,
where Y ; γ , σ 2 is the evidence term, which can be ignored because for a given γ , σ 2 , which is the normalization factor. Thus, Equation (12) can be simplified to
X Y ; γ , σ 2 Y X ; σ 2 [ X ; γ ]
e tr X μ X H Σ x 1 X μ X π N det Σ x L = N C μ X , Σ x .
where tr means trace of square matrix, and ( · ) H means Hermitian transpose of the matrix. Their product Equation (13) is a Gaussian distribution with mean μ X in Equation (15) and covariance Σ x in Equation (16) since both the likelihood Equation (8) and the prior in Equation (10) are Gaussian distributions.
μ X = E X Y ; γ , σ 2 = Γ A H Σ y 1 Y ,
Σ x = E x l μ x l x l μ x l H Y ; γ , σ 2 = Γ Γ A H Σ y 1 A Γ ,
where the infrasound sensor data covariance Σ y
Σ y = E y l y l H = σ 2 I N + A Γ A H ,
Σ y 1 = σ 2 I N σ 2 A 1 σ 2 A H A + Γ 1 1 A H σ 2 = σ 2 I N σ 2 A Σ x A H σ 2
where I N is the identity matrix of order N and E means expected value. The posterior mean X ^ MAP can be obtained using the maximum a posteriori (MAP) method when γ and σ 2 are known.
X ^ MAP = Γ A H Σ y 1 Y = μ X .
(4)
Evidence on the infrasound sources
The γ and σ 2 in the Equations (16)–(19) are considered as hyperparameters. They are estimated by maximizing the evidence. The evidence that is considered as a constant in Equation (12) is the product of the integration of the likelihood part in the Equation (8) and the prior part in the Equation (10) over the complex infrasound source amplitude X .
Y ; γ , σ 2 = l = 1 L N C y l ; 0 , Σ y l = R 2 M L Y X ; σ 2 [ X ; γ ] d X ,
The logarithm of the evidence is
log Y ; γ , σ 2 l = 1 L y l H Σ y l 1 y l + log det Σ y l ,
where the infrasound array data sample covariance matrix (SCM) is defined by the following equation
S y = Y Y H / L .
(5)
Infrasound source power estimation using the SBL algorithm
The above is the Bayesian statistical modeling phase, and the following is the solution phase. The infrasound source power (i.e., γ m ) was estimated by the SBL algorithm following the type-II maximum likelihood (maximizing the evidence) in Equation (23) [14]. The hyperparameters γ ^ and σ ^ 2 can be obtained by the type-II maximum likelihood (maximizing the evidence). That is, the following equation:
γ ^ , σ ^ 2 = arg m a x γ 0 , σ 2 > 0 log Y ; γ , σ 2 ,
where arg m a x is an abbreviation for “arguments of the maxima”. The objective function in Equation (23) is differentiated to obtain a local maximum. The following are derivative relations:
Σ y l 1 γ m = Σ y l 1 Σ y l γ m Σ y l 1 = Σ y l 1 a m a m H Σ y l 1 ,
log det Σ y l γ m = tr Σ y l 1 Σ y l γ m = a m H Σ y l 1 a m ,
the derivative of Equation (21) is
log p Y ; γ , σ 2 γ m = l = 1 L y l H Σ y l 1 a m a m H Σ y l 1 y l a m H Σ y l 1 a m = l = 1 L y l H Σ y l 1 a m 2 l = 1 L a m H Σ y l 1 a m .
For the solution to Equation (21), the necessary conditions are set, i.e., log   p Y ; γ , σ 2 γ m = 0 . To obtain an iterative equation in γ m , the first term in the Equation (26) is multiplied by the factor γ m old γ m new , whereby
γ m old γ m new l = 1 L y l H Σ y l 1 a m 2 l = 1 L a m H Σ y l 1 a m = 0 .
Assuming γ m old and Σ y l are known from previous iteration or initialization, the fixed point iteration for γ m is obtained by rewriting Equation (27) (with similarities to [10,11,13,21]:
γ m new = γ m old l = 1 L y l H Σ y l 1 a m 2 l = 1 L a m H Σ y l 1 a m .
Usually, the convergence of such fixed-point iterations is not guaranteed [11]. A larger number of iterations (about 1000) is required to guarantee the convergence of the iteration parameters. In this paper, the maximum number of iterations of the parameter γ is set to 1000.
(6)
Infrasound source noise variance estimation (hyperparameter σ 2 )
The noise is also a part of the model (Equation (5)) because, in the SBL algorithm, the sharpness of the local maxima in the γ spectrum is controlled by the noise variance, i.e., the higher the noise level, the wider the peaks in the γ spectrum. By maximizing the evidence-based approach, similar to the approach used to estimate the γ , one can estimate the infasound source noise variance σ [16,22]. The sparsity of the matrix X ^ MAP is controlled by Γ . Thus, the active set M is defined as
M = m N γ m > 0 .
The infrasound sensors data covariance matrix is
Σ y = Σ n l + A M Γ M A M H ,
where Σ n l = σ 2 I N ; A M is the active infrasound array steering vectors and it consists of K columns of A , indexed by M ; the set M denotes the position of γ non-zero entries with cardinality K, and it is estimated by selecting the strongest K peak from γ ^ new ; and Γ M = diag γ M new is the covariance matrix of the K active infrasound sources. The Jaffer’s necessary condition at the optimal solution Γ M , σ 2 must be satisfied, i.e., the following relationship is satisfied:
A M H S y Σ y A M = 0 ,
where S y is the infrasound array data sample covariance matrix. The necessary condition of Jaffer is satisfied since arbitrary correlations between infrasound source signals are allowed. Substituting Equation (30) into Equation (31) gives
A M H S y Σ n l A M = A M H A M Γ M A M H A M ,
multiplying (32) from right and left with the A M + = A ˜ M H ˙ A M 1 A M H and A M + H , respectively, and subtracting S y from both sides yields [22]. The σ 2 estimate of the infrasound source noise variance is
σ 2 new = 1 N K tr I N A M A M + S y .

2.2.3. The Pseudocode for the SBL-Based Estimation Algorithm for Infrasound Arrival Azimuth

The SBL-based infrasound source arrival azimuth estimation algorithm is summarized in Algorithm 1. The structure diagram of the Bayesian hierarchical model is shown in Figure 3. The μ X in Equation (15) and Σ y in Equation (17) are updated iteratively using the current γ for a given observed infrasound sensors data Y . Then γ m for m = 1 , , M is updated via Equation (28) and then Equation (33) is uesd to estimate the noise variance σ 2 . The infrasound source is assumed to be single in this paper, assuming that only one infrasound event occurs at a given moment, i.e., K = 1 , for estimating the noise variance in the Equation (33). The algorithm performance was not found to be sensitive to the number of infrasound sources K in the simulations. The algorithm is made more flexible by this assumption since there is no need to know the true number of infrasound sources. The relative improvement in the total infrasound source power is measured by the convergence rate ϵ in the Equation (34), where · 1 is the 1 norm. The algorithm stops when ϵ ϵ m i n . In this paper, the ϵ m i n is set to 10 4 . The j m a x is the maximum number of iterations. In this paper, the j m a x is set to 1000. The output is the active set M in the Equation (29). The estimated value of the arrival azimuth θ is the angle corresponding to the maximum gamma value γ m a x .
ϵ = γ m new γ m old 1 / γ m old 1 ,
Algorithm 1: The SBL-based infrasound source orientation algorithm
Input: A C N × M , Y C N × L , K = 1 , σ 0 2 = 0.1 , γ 0 = 1 , ϵ m i n = 10 4 , j m a x = 1000 ;
Output: γ new , σ 2 new , θ ;
  1: initialize: j = 0 , σ 2 = σ 0 2 , γ = γ 0 ;
  2: while ϵ > ϵ m i n and j < j m a x ;
  3:       Compute: Σ y = σ 2 I N + A Γ A H
  4:       Compute: μ m = γ m a m H Σ y 1 Y
  5:        γ m new update for m using Equation ( 28 )
  6:        σ 2 new = 1 N K tr I N A M A M + S y
  7:        ϵ = γ new γ old 1 / γ old 1
  8:        θ is the angle corresponding to γ m a x
  9: end

2.3. Bayesian Information Fusion for Infrasound Source Localization

The arrival azimuth estimation results of individual infrasound stations can be derived from the SBL algorithm, and the accuracy of infrasound source localization can be greatly improved when the fusion of information from infrasound stations is considered. Multiple information is fused by the Bayesian information fusion algorithm based on Bayesian formulation [23,24]. In this paper, the azimuth, measurement error, and model error of the infrasound station are combined by a Bayesian information fusion algorithm. The final fused multi-information infrasound source localization results are obtained.

2.3.1. Posteriori Distribution for Infrasound Source Localization

The arrival azimuths of the n infrasound stations estimated by the SBL algorithm are expressed as θ = θ 1 , , θ n . The location of the i-th infrasound station can be represented by x i , y i , and this is the Cartesian coordinate system, which can be transformed from the geodetic coordinate system.
The location parameter m = x 0 , y 0 to be estimated is the candidate infrasound source location. The posterior probability density function of the location parameter is
[ m θ ] = c [ θ ] [ m ] [ θ m ] ,
where [ θ m ] is the likelihood function; c [ θ ] ensures that the likelihood function [ θ m ] integrates to unity; and [ m ] is the prior for the infrasound source location parameters.

2.3.2. Likelihood Function

The likelihood function can be derived from the product of the estimated components of the arrival azimuths of all infrasound stations, under the assumption that the detections of each station are independent of each other.
[ θ m ] i = 1 n Θ i θ i m .
The consistency between the observed arrival azimuth, and the selected model parameters is measured by the Θ i , for a given i-th infrasound station. It is assumed that the errors between the observed arrival azimuths and the selected model parameters conform to a Gaussian distribution. The arrival azimuth likelihood function measured for the ith station is
Θ i θ i m 1 2 π σ θ 2 exp 1 2 γ i σ θ 2 .
Let ( x i , y i ) represent the candidate infrasound source position to the i-th infrasound station, then the arrival azimuth estimation residuals are
γ i θ i arctan y i y 0 x i x 0 .
The total variances in the arrival azimuth estimation are denoted σ θ 2 , which, accounts for both infrasound measurement equipment σ θ , meas 2 and infrasound propagation model σ θ , mod 2 error contributions
σ θ 2 = σ θ , meas 2 + σ θ , mod 2 ,
where σ θ 2 can be estimated using historical infrasound events. Usually, σ θ = 3 . 5 is suggested by an analysis of the western United States earthquakes [6]. The σ θ = 3 . 5 is mainly determined by the combination of discretization used in the slowness plane, the array aperture, and the signal-to-noise ratio (SNR). The variance on the infrasound measurement equipment error σ θ , meas 2 is a function of the distance between the sensors of the infrasound station, sampling frequency, and wavelength. The variance on the infrasound propagation model error σ θ , mod 2 is derived from the variation of wind speed, wind direction, and temperature.

2.3.3. The Framework for Infrasound Source Localization

An infrasound source localization algorithm for distributed stations based on the SBL algorithm and Bayesian information fusion is proposed in the context of the large-scale global distribution of infrasound stations. The accuracy of the infrasound source can be improved by the proposed algorithm while incorporating the uncertainty of the infrasound measurement equipment and the infrasound propagation uncertainty into the infrasound source localization results. First, the infrasound signals recorded at individual stations are processed based on the SBL algorithm to obtain the azimuth information of the corresponding stations. Second, infrasound propagation in the atmosphere is susceptible to wind speed, temperature, and atmospheric density, and the effect is modeled as the model variance σ θ , mod . Third, the measurement error generated during infrasound measurement is modeled as the measurement variance σ θ , meas . Finally, the Bayesian information fusion is used to fuse the localization results, model errors, and measurement errors from multiple stations to obtain the localization results of infrasound sources. The schematic diagram for Bayesian information fusion algorithm is shown in Figure 4.

2.4. Technical Flow Chart of the Proposed Infrasound Source Localization Method

The actual measured infrasound data are first used to obtain the arrival azimuth information of each station using a sparse Bayesian learning algorithm. Then, the arrival azimuth information obtained from each station is imported into the Bayesian information fusion algorithm to obtain the localization results of the infrasound source. The technical flow chart of the infrasound source localization method proposed in this paper is shown in Figure 5.

3. Results

3.1. Experimetntal Setup

The SBL algorithm was tested by UTTR’s rocket motor explosion data (explosive weight 39,000 lb (approximately 17,690 kg)) [25,26,27] to verify the feasibility of the SBL algorithm. For the UTTR event, the real source location and the location information of the four observatories are shown in the Table 1.
Firstly, the range of reconstruction frequency is selected according to the spectrum of infrasound. In this paper, the frequency range of 1–2 Hz is chosen for signal reconstruction, and then the reconstruction signal corresponding to each reconstruction frequency is solved at an interval of 0.1 Hz. The reconstructed signal is substituted into the SBL algorithm to solve the gamma value γ corresponding to the reconstruction frequency. The estimated value of the arrival azimuth θ is the angle corresponding to γ m a x .
The infrasound data recorded at the four stations BGU, BRP, HWU, and WMU were used as the experimental data, which were sampled at a frequency of 100 Hz. The infrasound signals with a passband of 1–5 Hz were selected using a Butterworth bandpass filter. The array structure of the four infrasound stations is shown in Figure 6.
The time-domain waveforms of the four stations are shown in Figure 7. The interception of the infrasound signals recorded at each station is selected as follows for each station. The BGU station intercepts 20 s signal as well as the sensor selections BGU1 and BGU3, and the array aperture is 135 m; the BRP station intercepts the 30 s signal as well as the sensor selections BRP1 and BRP3, and array aperture is 150 m; the HWU station intercepts the 100 s signal as well as the sensor selections HWU1 and HWU3, and the array aperture is 140 m; and the WMU station intercepts the 15 s signal as well as the sensor selections WMU2 and WMU3, and array aperture is 150 m. The intercepted infrasound is substituted into the SBL algorithm to obtain the azimuth information of the four stations. The infrasound signal and the SBL algorithm parameter setting are shown in Table 2.
The main parameters of the sparse Bayesian learning algorithm are set as follows: the number of snapshots is set to 39, the angular resolution Δ θ is 0.1°, the convergence error ϵ m i n is 0.0001, the maximum number of iterations j m a x is 1000, and the sampling frequency is 100 Hz. The process of transforming the infrasound signal from the time domain to the frequency domain is shown in the Figure 8. The infrasound frequency domain signal is substituted into the sparse Bayesian learning algorithm for the solution, and the results are shown in Figure 9.

3.2. Arrival Azimuth Estimation Results

The result of the arrival azimuth estimation of the SBL algorithm is the angle value corresponding to the maximum gamma value γ m a x . The results of the SBL algorithm arrival azimuth estimation θ for the four infrasound stations BGU, BRP, HWU, and WMU are 27.2°, 44.5°, −38.4°, and 27.9°, respectively, as shown in Figure 9. The arrival azimuth estimations with the SBL algorithm for the four stations were converted to azimuths in the geographic coordinate system by considering the effect of the array tilt angles of the infrasound stations, which were 4.71°, 8.15°, −33.67°, and 73.54° for the BGU, BRP, HWU, and WMU stations, respectively. The array tilt angle is defined as the acute angle between a linear array of two array sensors and the latitude line in the geographic coordinate system. The azimuths obtained by the four stations are shown in Figure 10.
It can be seen that the converted arrival azimuths of the SBL algorithm results for the four infrasound stations BGU, BRP, HWU, and WMU are 32.91°, 307.35°, 252.07°, and 314.35°, respectively. The actual azimuths between the infrasound source and the infrasound stations were compared to obtain the arrival azimuth estimation errors of the SBL algorithm, as shown in Table 3.

3.3. Comparison of the Error of Arrival Azimuth Estimation between SBL Algorithm, FK Analysis, and BF Algorithm

The SBL-based arrival azimuth estimation algorithm and the conventional infrasound arrival azimuth estimation algorithm (FK analysis) are applied for comparison in terms of the accuracy of arrival azimuth estimation. The complete slowness vector of infasound in the plane-wave hypothesis can be measured by the FK analysis. Additionally, the infrasound power distribution at different slownesses is calculated by FK analysis [28,29]. For the same infrasound events, the arrival azimuth estimation results of the FK analysis are shown in Figure 11. The azimuths obtained from the FK analysis of the four stations BGU, BRP, HWU, and WMU are 29.05°, 312.13°, 247.83°, and 323.53°, respectively.
The beamforming [30] algorithm is used to process the infrasound signal for the same infrasound event. The delay-and-sum technique is utilized by the algorithm to enhance the coherent signal and suppress the non-coherent background noise to obtain the arrival azimuth information. The results of beamforming are shown in Figure 12. From Figure 12, it can be seen that the results of the beamforming algorithm for the BGU, BRP, HWU, and WMU stations are 23.8°, 14.1°, −41.3°, and 26.9°, respectively. The results of the beamforming algorithms are converted into arrival azimuths by considering the effect of array inclination. The estimated azimuth of arrival results for BGU, BRP, HWU, and WMU stations are 28.51°, 310.74°, 254.97°, and 313.36°, and the orientation errors are 3.48°, 3.23°, 4.43°, and 1.30°, respectively.
The arrival azimuth estimation errors of the three arrival azimuth estimation algorithms are shown in Table 4.
According to Table 4, the arrival azimuth estimation errors of all four stations of the SBL algorithm are within 2°. The error is mainly due to the influence of low-frequency noise (such as thunderstorm noise, ground noise, and the noise of the measurement station itself) during the long-distance propagation of infrasound. The superiority of the SBL algorithm over the FK analysis in terms of arrival azimuth estimation accuracy is shown in Table 4. Compared with beamforming, the estimation error of the arrival azimuth of SBL is smaller and the side flaps are smaller, as can be seen from Table 4 and Figure 12. Therefore, the arrival azimuth estimation results of the SBL algorithm are used as prior knowledge for the Bayesian information fusion algorithm to perform fusion localization in the following.

3.4. Bayesian Information Fusion to Obtain Localization Results

The arrival azimuth estimation results of the infrasound sources are affected by the uncertainty of the infrasound propagation and the infrasound measurement equipment. The uncertainty of the arrival azimuth estimation results is modeled by the Bayesian information fusion algorithm as the variance in a probabilistic model. The localization result of the infrasound source is the maximum value of the posterior probability density function. The credibility contour curve represents the possible area of the infrasound source. The localization results are shown in Figure 13.
Figure 13. The localization result of the infrasound source by fusing the SBL arrival azimuth estimations. (The three blue contour curves are credibility contour curves with Bayesian credibility values of 0.75, 0.90, and 0.95, respectively; the red pentagram indicates the location of the infrasound source; the blue hexagram indicates the distributed array localization results; and the enlarged scale is 1:10).
Figure 13. The localization result of the infrasound source by fusing the SBL arrival azimuth estimations. (The three blue contour curves are credibility contour curves with Bayesian credibility values of 0.75, 0.90, and 0.95, respectively; the red pentagram indicates the location of the infrasound source; the blue hexagram indicates the distributed array localization results; and the enlarged scale is 1:10).
Remotesensing 14 03181 g013
In Figure 13, it can be seen that the infrasound source localization result and ground truth location are (41.10 N, 112.89 W) and (41.13 N, 112.89 W), respectively. The distance between the infrasound source and the localization result is 3.4855 km. Three of the curves shown in Figure 13 are credibility contour curves with Bayesian credibility values of 0.75, 0.90, and 0.95, respectively, which are the possible regions of the infrasound source. From the above localization results, it can be seen that accurate infrasound source localization results can be obtained by the proposed infrasound source algorithm for infrasound sources within 250 km distance from the infrasound station. The arrival azimuth estimation results from the FK analysis are used to obtain the infrasound source localization results using the Bayesian information fusion algorithm, which are shown in Figure 14.
Compared with Figure 13, the results of infrasound source localization are the same; both are 41.10 N, 112.90 W. From Table 4, it can be seen that the estimated errors of the arrival azimuths of the FK analysis are 2.02°, 4.62°, 1.18°, and 8.87° larger than those of the SBL for the four stations BGU, BRP, HWU, and WMU, respectively. The area of the credibility contour corresponding to the Bayesian credibility of 0.75 is 918.6532 km 2 and 923.7934 km 2 in Figure 13 and Figure 14, respectively. The area of the credibility contour corresponding to the Bayesian credibility of 0.75 obtained by the SBL is 5.1402 km 2 smaller than that of the FK analysis, which means that the corresponding mode of the posterior probability density function is more concentrated (i.e., the posterior probability density function is visually sharper at the peak). The above results indicate that the proposed method is more reliable than the combination method of FK analysis and Bayesian information fusion.

4. Discussion

The performance of the proposed infrasound source localization algorithm is discussed in this section. The accuracy of the infrasound source localization is related to the number of stations involved in localization. The results of Bayesian information fusion localization for different number of stations are shown in Figure 15. As shown in Figure 15, the localization errors of two stations, three stations, and four stations are 78.3306 km, 13.0329 km, and 3.4855 km, respectively. It can be seen that the error of infrasound source localization is further reduced with the increase in the number of infrasound stations. The same conclusion can be obtained by performing Bayesian information fusion algorithm on the results of the FK analysis and beamforming.
It can be seen from Figure 16a that the four stations NOQ, BGU, WMU, and BRP are basically arranged in a straight line. By contrast, in Figure 16b, the azimuth of the HWU station with respect to the infrasound source is different from the azimuths of the BGU, WMU, and BRP stations. Therefore, the spatial distribution of the four stations with respect to the infrasound source in Figure 16b is more dispersed compared with that in Figure 16a. The localization error of the algorithm is also related to the arrangement of infrasound stations. The NOQ station is added in order to discuss whether the localization error of the proposed infrasound source localization algorithm is related to the arrangement of infrasound stations. The position of this infrasound station is (40.6526 N, 112.1180 W), and the arrival azimuths obtained by the SBL algorithm is 301.10°. When the station arrangement is more dispersed relative to the infrasound source, as shown in Figure 16a, the localization result is (41.10 N, 112.90 W). It can be seen that the arrangement of the stations has little influence on the localization error of the infrasound source localization algorithm proposed in this paper. The same conclusion can be obtained by performing the Bayesian information fusion algorithm on the results of the FK analysis and beamforming. In Figure 16, the area of the credibility contours corresponding to the Bayesian credibility value of 0.75 is 387.6142 km 2 (Figure 16a) and 993.1705 km 2 (Figure 16b), respectively. Therefore, the station arrangement has a certain influence on the credibility contour areas of the infrasound source localization result. The estimation error of the arrival azimuth of the HWU station is 0.99° larger than that of the NOQ station, but the localization error is the same for both arrangements. The impact on the final localization error can be reduced by the Bayesian information fusion algorithm when the estimation error of the arrival azimuth of a single station is 0.99° larger.
Next, we discuss the influence of the distance between the station and the infrasound source on the error of the localization results. In Figure 17, the sum of the distances from BGU, NOQ, and WMU stations to the infrasound source; the BRP, NOQ, and WMU stations to the infrasound source; and the BRP, HWU, and WMU stations to the infrasound source are 295.1944 km, 569.0806 km, and 624.5797 km, respectively. The errors of the infrasound source localization for these three cases are 3.4855 km, 23.0137 km, and 45.6028 km, respectively. As can be seen from Figure 17, the localization error of the algorithm increases as the stations are arranged farther away relative to the infrasound source. Therefore, it is necessary to select the stations that are arranged closer to the infrasound source for analysis to obtain accurate infrasound source localization results. The same conclusion can be obtained by performing the Bayesian information fusion algorithm on the results of the FK analysis and beamforming.

5. Conclusions

The precise localization of infrasound sources is a problem of considerable interest. The infrasound source localization can be affected by the atmospheric environment and infrasound measurement equipments, resulting in poor accuracy of infrasound source localization. In this paper, the SBL algorithm is applied to infrasound source arrival azimuth estimation for the first time. The uncertainties in the atmospheric environment and measurement equipments can be modeled as the variance in the Bayesian information fusion algorithm. The maximum value of posterior probability is used as the infrasound source localization result. Infrasound source credibility contours with Bayesian credibility values of 0.75, 0.90 and 0.95, respectively, are obtainable by the Bayesian information fusion algorithm. Subsequently, rocket motor explosion infrasound data obtained in UTTR were used to verify feasibility of the algorithm. The infrasound signals recorded at the four infrasound stations are processed by the SBL algorithm and the Bayesian information fusion algorithm. The analysis of the measured data shows that the arrival azimuth estimation angle error of the SBL algorithm is within 2°. The arrival azimuth estimation accuracy is significantly improved, compared with the FK analysis and beamforming. Finally, the Bayesian information fusion algorithm is used to fuse the information of multiple stations to obtain the localization results of infrasound sources. The analysis of the measured data shows that for the infrasound source within 250 km from the infrasound station, the error of the positioned infrasound source distance is within 3.5 km.
In summary, (i) the SBL algorithm is applied to infrasound source localization in this paper, which significantly reduces the estimation error of the arrival azimuth; (ii) the Bayesian information fusion algorithm is used to incorporate the uncertainty brought by the infrasound propagation environment and the measurement equipment into the infrasound source localization, making the localization results more robust.

Author Contributions

Conceptualization, R.W. and L.Y.; methodology, R.W., L.Y., X.Y. and C.Z.; formal analysis, X.Y.; data curation, X.Y.; software, X.Y., C.Z.; writing—original draft, X.Y., R.W., L.Y., C.Z., T.W. and X.Z.; writing—review and editing, R.W., L.Y., C.Z., T.W. and X.Z.; supervision, R.W. and L.Y.; and funding acquisition, R.W. and L.Y. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the Natural Science Foundation of Shanghai (Grant No. 21ZR1434100, Founder: Liang Yu) and the National Natural Science Foundation of China (Grant No. 12074254, Founder: Liang Yu; Grant No. 51505277, Founder: Ran Wang). This work was also supported by the Marine Interdisciplinary Program of Shanghai Jiao Tong University (Project No. SL2021MS009, Founder: Liang Yu).

Data Availability Statement

Actual measured infrasound data are available from Los Alamos National Laboratory (https://www.lanl.gov/org/ddste/aldcels/earth-environmental-sciences/geophysics/software/seismoacoustics/inframonitor.php, accessed on 15 May 2021).

Acknowledgments

We sincerely thank Los Alamos National Laboratory for providing the infrasound data of the UTTR rocket motor explosion experiment [25], which occurred on 4 June 2007, and the geographic location of the infrasound event (41.13 N, 112.895 W). The authors would like to thank Xun Wang (Beihang University) for his support and full discussion. The authors sincerely appreciate the advice and assistance of the anonymous reviewers.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
pray parameter in tau-p model
R ( z , p ) range along the ray direction
ψ ( z , p ) the ray characteristic function
Q ( z , p )                               transverse offset
u ( z ) the horizontal wind speed along the propagation direction
z m a x the maximum height of the infrasound trajectory
ϕ the launch elevation angle at the start of the infrasound trajectory
z (axis)infrasound propagation height
α angle of wave number vector with x-axis
β angle of wave number vector with z-axis
u wave number vector
P (function)the acoustic signal when the initial time is t
r (vector)propagation position vector
tinitial time
X infrasound source amplitudes
Lsnapshot
Y infrasound signal of L snapshots observed by N infrasound sensors
N additive noise
Nthe number of infrasound sensor
N C the normal distribution sign
A infrasound array steering vectors
Kthe number of infrasound sources
Mnumber of infrasound planar angle discretizations
PDFprobability density function
σ 2 complex Gaussian with noise variance
Γ the diagonal matrix formed by the diagonal elements γ
Y = y 1 , , y L C N × L the measured infrasound signal
Σ x covariance of the posterior distribution on the infrasound source
Σ y the infrasound sensor data covariance
S y the infrasound array data sample covariance matrix (SCM)
· F Frobenius norm
trtrace of square matrix
( · ) H Hermitian transpose of the matrix
Eexpected value
a m the infrasound array steering vector
· 1 1 norm
[ Y X ] stochastic likelihood for the SBL algorithm
[ X Y ] posterior on the infrasound sources
[ X ] prior on the infrasound sources
[ Y ] evidence on the infrasound sources
[ m θ ] posterior pdf for the Bayesian information fusion
c [ θ ] enables the integration of [ θ m ] to be uniform
[ m ] the prior pdf for the Bayesian information fusion
[ θ m ] the likelihood function for the Bayesian information fusion
γ m the variance of azimuth on the random variable x m l at the m-th
snapshot
θ = θ 1 , , θ n arrival azimuth estimation
m = x 0 , y 0 the candidate source location
σ θ , meas 2 the variance of the arrival azimuth estimation from the infrasound
measurement equipment
σ θ , mod 2 the variance of the arrival azimuth estimation from the infrasound
propagation model

References

  1. Le Pichon, A.; Blanc, E.; Hauchecorne, A. Infrasound Monitoring for Atmospheric Studies; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010. [Google Scholar]
  2. Freret-Lorgeril, V.; Bonadonna, C.; Corradini, S.; Donnadieu, F.; Guerrieri, L.; Lacanna, G.; Marzano, F.S.; Mereu, L.; Merucci, L.; Ripepe, M.; et al. Examples of Multi-Sensor Determination of Eruptive Source Parameters of Explosive Events at Mount Etna. Remote Sens. 2021, 13, 2097. [Google Scholar] [CrossRef]
  3. Cigna, F.; Tapete, D.; Lu, Z. Remote Sensing of Volcanic Processes and Risk. Remote Sens. 2020, 12, 2567. [Google Scholar] [CrossRef]
  4. Batubara, M.; Yamamoto, M.-y. Infrasound Observations of Atmospheric Disturbances Due to a Sequence of Explosive Eruptions at Mt. Shinmoedake in Japan on March 2018. Remote Sens. 2020, 12, 728. [Google Scholar] [CrossRef] [Green Version]
  5. De Angelis, S.; Diaz-Moreno, A.; Zuccarello, L. Recent Developments and Applications of Acoustic Infrasound to Monitor Volcanic Emissions. Remote Sens. 2019, 11, 1302. [Google Scholar] [CrossRef] [Green Version]
  6. Mutschlecner, J.P.; Whitaker, R.W. Infrasound from earthquakes. J. Geophys. Res. Atmos. 2005, 110. [Google Scholar] [CrossRef] [Green Version]
  7. Garces, M.; Pichon, A.L. Infrasound from Earthquakes, Tsunamis and Volcanoes. In Encyclopedia of Complexity and Systems Science; Springer: New York, NY, USA, 2011; pp. 663–679. [Google Scholar] [CrossRef]
  8. Laiolo, M.; Ripepe, M.; Cigolini, C.; Coppola, D.; Della Schiava, M.; Genco, R.; Innocenti, L.; Lacanna, G.; Marchetti, E.; Massimetti, F.; et al. Space- and Ground-Based Geophysical Data Tracking of Magma Migration in Shallow Feeding System of Mount Etna Volcano. Remote Sens. 2019, 11, 1182. [Google Scholar] [CrossRef] [Green Version]
  9. Rost, S.; Thomas, C. Array Seismology: Methods and Applications. Rev. Geophys. 2002, 40, 2-1–2-27. [Google Scholar] [CrossRef] [Green Version]
  10. Tipping, M.E. Sparse Bayesian learning and the relevance vector machine. J. Mach. Learn. Res. 2001, 1, 211–244. [Google Scholar]
  11. Wipf, D.; Rao, B. Sparse Bayesian learning for basis selection. IEEE Trans. Signal Process. 2004, 52, 2153–2164. [Google Scholar] [CrossRef]
  12. Zhang, Z.; Rao, B.D. Sparse Signal Recovery With Temporally Correlated Source Vectors Using Sparse Bayesian Learning. IEEE J. Sel. Top. Signal Process. 2011, 5, 912–926. [Google Scholar] [CrossRef] [Green Version]
  13. Gerstoft, P.; Mecklenbräuker, C.F.; Xenaki, A.; Nannuru, S. Multisnapshot Sparse Bayesian Learning for DOA. IEEE Signal Process. Lett. 2016, 23, 1469–1473. [Google Scholar] [CrossRef] [Green Version]
  14. Gemba, K.L.; Nannuru, S.; Gerstoft, P. Robust Ocean Acoustic Localization With Sparse Bayesian Learning. IEEE J. Sel. Top. Signal Process. 2019, 13, 49–60. [Google Scholar] [CrossRef]
  15. Nannuru, S.; Gemba, K.L.; Gerstoft, P.; Hodgkiss, W.S.; Mecklenbräuker, C.F. Sparse Bayesian learning with multiple dictionaries. Signal Process. 2019, 159, 159–170. [Google Scholar] [CrossRef] [Green Version]
  16. Liu, Z.M.; Huang, Z.T.; Zhou, Y.Y. An Efficient Maximum Likelihood Method for Direction-of-Arrival Estimation via Sparse Bayesian Learning. IEEE Trans. Wirel. Commun. 2012, 11, 1–11. [Google Scholar] [CrossRef]
  17. Chen, S.S.; Donoho, D.L.; Saunders, M.A. Atomic Decomposition by Basis Pursuit. SIAM Rev. 2001, 43, 129–159. [Google Scholar] [CrossRef] [Green Version]
  18. Garcés, M.A.; Hansen, R.A.; Lindquist, K.G. Traveltimes for infrasonic waves propagating in a stratified atmosphere. Geophys. J. Int. 1998, 135, 255–263. [Google Scholar] [CrossRef] [Green Version]
  19. Modrak, R.T.; Arrowsmith, S.J.; Anderson, D.N. A Bayesian framework for infrasound location. Geophys. J. Int. 2010, 181, 399–405. [Google Scholar] [CrossRef] [Green Version]
  20. Drob, D.P.; Garcés, M.; Hedlin, M.; Brachet, N. The temporal morphology of infrasound propagation. Pure Appl. Geophys. 2010, 167, 437–453. [Google Scholar] [CrossRef]
  21. Wipf, D.P.; Rao, B.D. An Empirical Bayesian Strategy for Solving the Simultaneous Sparse Approximation Problem. IEEE Trans. Signal Process. 2007, 55, 3704–3716. [Google Scholar] [CrossRef]
  22. Bohme, J. Source-parameter estimation by approximate maximum likelihood and nonlinear regression. IEEE J. Ocean. Eng. 1985, 10, 206–212. [Google Scholar] [CrossRef]
  23. Ji, S.; Xue, Y.; Carin, L. Bayesian Compressive Sensing. IEEE Trans. Signal Process. 2008, 56, 2346–2356. [Google Scholar] [CrossRef]
  24. Blom, P.S.; Marcillo, O.; Arrowsmith, S.J. Improved Bayesian Infrasonic Source Localization for regional infrasound. Geophys. J. Int. 2015, 203, 1682–1693. [Google Scholar] [CrossRef] [Green Version]
  25. Stump, B.; Burlacu, R.; Hayward, C.; Pankow, K.; Nava, S.; Bonner, J.; Hock, S.; Whiteman, D.; Fisher, A.; Kim, T.S. Seismic and Infrasound Energy Generation and Propagation at Local and Regional Distances Phase 1-Divine Strake Experiment; Southern Methodist University: Dallas, TX, USA, 2008. [Google Scholar]
  26. Arrowsmith, S.J.; Whitaker, R.; Taylor, S.R.; Burlacu, R.; Stump, B.; Hedlin, M.; Randall, G.; Hayward, C.; ReVelle, D. Regional monitoring of infrasound events using multiple arrays: Application to Utah and Washington State. Geophys. J. Int. 2008, 175, 291–300. [Google Scholar] [CrossRef] [Green Version]
  27. Stump, B.W.; Zhou, R.M.; Kim, T.S.; Chen, Y.T.; Yang, Z.X.; Herrmann, R.B.; Burlacu, R.; Hayward, C.; Pankow, K. Shear Velocity Structure in NE China and Characterization of Infrasound Wave Propagation in the 1–210 km Range; Southern Methodist University: Dallas, TX, USA, 2008. [Google Scholar]
  28. Capon, J.; Bolt, B. Signal processing and frequency-wavenumber spectrum analysis for a large aperture seismic array. In Methods in Computational Physics; Elsevier: Amsterdam, The Netherlands, 1973; Volume 13, pp. 1–59. [Google Scholar]
  29. Aki, K.; Richards, P.G. Quantative seismology: Theory and methods. Earth Sci. Rev. 1981, 17, 296–297. [Google Scholar]
  30. Shani-Kadmiel, S.; Averbuch, G.; Smets, P.; Assink, J.; Evers, L. The 2010 Haiti earthquake revisited: An acoustic intensity map from remote atmospheric infrasound observations. Earth Planet. Sci. Lett. 2021, 560, 116795. [Google Scholar] [CrossRef]
Figure 1. Propagation trajectory of the tau-p model in a phase loop, the azimuth angle θ , the angle of elevation ϕ , the altitude of infrasound propagation z, the infrasound propagation distance R of a phase loop in Equation (2), the maximum altitude z m a x that a phase loop can propagate, and the transverse offset Q of infrasound propagation in Equation (4) (the red arrow and red dot indicate the direction of infrasound propagation at the starting point and the location of the maximum elevation of infrasound propagation, respectively).
Figure 1. Propagation trajectory of the tau-p model in a phase loop, the azimuth angle θ , the angle of elevation ϕ , the altitude of infrasound propagation z, the infrasound propagation distance R of a phase loop in Equation (2), the maximum altitude z m a x that a phase loop can propagate, and the transverse offset Q of infrasound propagation in Equation (4) (the red arrow and red dot indicate the direction of infrasound propagation at the starting point and the location of the maximum elevation of infrasound propagation, respectively).
Remotesensing 14 03181 g001
Figure 2. Schematic diagram of infrasound three-dimensional propagation under the assumption of plane waves. (The z-axis is the infrasound propagation height, which corresponds to the z-axis in Figure 1. u is the wave number vector; c is the sound velocity; r is the position vector; P is the acoustic signal when the initial time is t; and α , β are the angle values).
Figure 2. Schematic diagram of infrasound three-dimensional propagation under the assumption of plane waves. (The z-axis is the infrasound propagation height, which corresponds to the z-axis in Figure 1. u is the wave number vector; c is the sound velocity; r is the position vector; P is the acoustic signal when the initial time is t; and α , β are the angle values).
Remotesensing 14 03181 g002
Figure 3. The structure diagram of Bayesian hierarchical model (blue rectangules denote the known parameters; red circles denote the unknown parameters). The initialized infrasound source power parameters γ 0 are used to update the parameters γ in the Equation (28), and the initialized infrasound source noise variance parameters σ 0 are used to update the parameters σ in the Equation (33).
Figure 3. The structure diagram of Bayesian hierarchical model (blue rectangules denote the known parameters; red circles denote the unknown parameters). The initialized infrasound source power parameters γ 0 are used to update the parameters γ in the Equation (28), and the initialized infrasound source noise variance parameters σ 0 are used to update the parameters σ in the Equation (33).
Remotesensing 14 03181 g003
Figure 4. The schematic diagram for Bayesian information fusion algorithm. (a) SBL:SBL-based infrasound signal processing of stations to obtain the azimuth of stations (azimuth = θ 1 , θ 2 , , θ N ); (b) Bayesian information fusion:information fusion of station-acquired azimuths using Bayesian information fusion to obtain localization results (red pentagrams indicate localized infrasound source locations, yellow triangles indicate infrasound stations, and the three blue contour lines are credibility curves with different Bayesian credibility values).
Figure 4. The schematic diagram for Bayesian information fusion algorithm. (a) SBL:SBL-based infrasound signal processing of stations to obtain the azimuth of stations (azimuth = θ 1 , θ 2 , , θ N ); (b) Bayesian information fusion:information fusion of station-acquired azimuths using Bayesian information fusion to obtain localization results (red pentagrams indicate localized infrasound source locations, yellow triangles indicate infrasound stations, and the three blue contour lines are credibility curves with different Bayesian credibility values).
Remotesensing 14 03181 g004
Figure 5. The technical flowchart of the infrasound source localization algorithm proposed in this paper. The first step is to estimate the arrival azimuth by the SBL algorithm by updating the hyperparameter gamma using Equation (28). The estimated value of the arrival azimuth is the angle corresponding to the γ m a x . The second step is the Bayesian information fusion for infrasound source localization. The posterior PDF maximum value is the infrasound source location m = ( x , y ) .
Figure 5. The technical flowchart of the infrasound source localization algorithm proposed in this paper. The first step is to estimate the arrival azimuth by the SBL algorithm by updating the hyperparameter gamma using Equation (28). The estimated value of the arrival azimuth is the angle corresponding to the γ m a x . The second step is the Bayesian information fusion for infrasound source localization. The posterior PDF maximum value is the infrasound source location m = ( x , y ) .
Remotesensing 14 03181 g005
Figure 6. The array structures of BGU, BRP, HWU, and WMU are shown in (ad), respectively.
Figure 6. The array structures of BGU, BRP, HWU, and WMU are shown in (ad), respectively.
Remotesensing 14 03181 g006
Figure 7. Time domain signal after 1-5 Hz filtering for four stations: (a) the total signal duration of the BGU station is 20 s; (b) the total signal duration of the BRP station is 30 s; (c) the total signal duration of the HWU station is 100 s; and (d) the total signal duration of the WMU station is 15 s.
Figure 7. Time domain signal after 1-5 Hz filtering for four stations: (a) the total signal duration of the BGU station is 20 s; (b) the total signal duration of the BRP station is 30 s; (c) the total signal duration of the HWU station is 100 s; and (d) the total signal duration of the WMU station is 15 s.
Remotesensing 14 03181 g007
Figure 8. The process of converting infrasound signals from the time domain to the frequency domain using the snapshot method. The variable d m is the observed position of the m-th sensor, and the variable L is the total number of snapshots.
Figure 8. The process of converting infrasound signals from the time domain to the frequency domain using the snapshot method. The variable d m is the observed position of the m-th sensor, and the variable L is the total number of snapshots.
Remotesensing 14 03181 g008
Figure 9. The results of the SBL algorithm for four infrasound stations (the angular resolution Δ θ is 0.1°): (a) the SBL algorithm arrival azimuth estimation result for the BGU station is 27.2°; (b) the SBL algorithm arrival azimuth estimation result for the BRP station is 44.5°; (c) the SBL algorithm arrival azimuth estimation result for the HWU station is −38.4°; and (d) the SBL algorithm arrival azimuth estimation result for the WMU station is 27.9°.
Figure 9. The results of the SBL algorithm for four infrasound stations (the angular resolution Δ θ is 0.1°): (a) the SBL algorithm arrival azimuth estimation result for the BGU station is 27.2°; (b) the SBL algorithm arrival azimuth estimation result for the BRP station is 44.5°; (c) the SBL algorithm arrival azimuth estimation result for the HWU station is −38.4°; and (d) the SBL algorithm arrival azimuth estimation result for the WMU station is 27.9°.
Remotesensing 14 03181 g009
Figure 10. The azimuths obtained from the BGU station, BRP station, HWU station, and WMU station are (a) 32.91°, (b) 307.35°, (c) 252. 07°, and (d) 314.35°, respectively (the array tilt angle is defined as the acute angle between a linear array of two array sensors and the latitude line in the geographic coordinate system).
Figure 10. The azimuths obtained from the BGU station, BRP station, HWU station, and WMU station are (a) 32.91°, (b) 307.35°, (c) 252. 07°, and (d) 314.35°, respectively (the array tilt angle is defined as the acute angle between a linear array of two array sensors and the latitude line in the geographic coordinate system).
Remotesensing 14 03181 g010
Figure 11. The azimuths from the FK analysis of the BGU, BRP, HWU, and WMU stations are (a) 29.05°, (b) 312.13°, (c) 247.83°, and (d) 323.53°, respectively.
Figure 11. The azimuths from the FK analysis of the BGU, BRP, HWU, and WMU stations are (a) 29.05°, (b) 312.13°, (c) 247.83°, and (d) 323.53°, respectively.
Remotesensing 14 03181 g011
Figure 12. The results from the beamforming of the BGU, BRP, HWU, and WMU stations are (a) 23.8°, (b) 41.1°, (c) −41.3°, and (d) 26.9°, respectively.
Figure 12. The results from the beamforming of the BGU, BRP, HWU, and WMU stations are (a) 23.8°, (b) 41.1°, (c) −41.3°, and (d) 26.9°, respectively.
Remotesensing 14 03181 g012
Figure 14. The localization result of the infrasound source by fusing the FK arrival azimuth estimations. (The three blue contour curves are credibility contour curves with Bayesian credibility values of 0.75, 0.90, and 0.95, respectively; the red pentagram indicates the location of the infrasound source; and the blue hexagram indicates the distributed array localization results).
Figure 14. The localization result of the infrasound source by fusing the FK arrival azimuth estimations. (The three blue contour curves are credibility contour curves with Bayesian credibility values of 0.75, 0.90, and 0.95, respectively; the red pentagram indicates the location of the infrasound source; and the blue hexagram indicates the distributed array localization results).
Remotesensing 14 03181 g014
Figure 15. Localization errors for different numbers of stations (arrival azimuth estimation results from the SBL algorithm); (a) 78.3306 km for two stations (HWU and BRP); (b) 13.0329 km for three stations (BGU, HWU, and BRP); and (c) 3.4855 km for four stations (BGU, HWU, BRP, and WMU).
Figure 15. Localization errors for different numbers of stations (arrival azimuth estimation results from the SBL algorithm); (a) 78.3306 km for two stations (HWU and BRP); (b) 13.0329 km for three stations (BGU, HWU, and BRP); and (c) 3.4855 km for four stations (BGU, HWU, BRP, and WMU).
Remotesensing 14 03181 g015
Figure 16. The influence of the spatial distribution of infrasound stations on the localization results (arrival azimuth estimation results from the SBL algorithm); (a) BGU, NOQ, WMU, and BRP; (b) BGU, HWU, WMU, and BRP.
Figure 16. The influence of the spatial distribution of infrasound stations on the localization results (arrival azimuth estimation results from the SBL algorithm); (a) BGU, NOQ, WMU, and BRP; (b) BGU, HWU, WMU, and BRP.
Remotesensing 14 03181 g016
Figure 17. The influence of the distance between the station and the infrasound source on the error of localization result (arrival azimuth estimation results from the SBL algorithm); (a) the sum of the distances from the BGU, NOQ, and WMU stations to the infrasound source is 295.1944 km; (b) the sum of the distances from the BRP, NOQ, and WMU stations to the infrasound source is 569.0806 km; and (c) the sum of the distances from BRP, HWU, and WMU stations to the infrasound source is 624.5797 km.
Figure 17. The influence of the distance between the station and the infrasound source on the error of localization result (arrival azimuth estimation results from the SBL algorithm); (a) the sum of the distances from the BGU, NOQ, and WMU stations to the infrasound source is 295.1944 km; (b) the sum of the distances from the BRP, NOQ, and WMU stations to the infrasound source is 569.0806 km; and (c) the sum of the distances from BRP, HWU, and WMU stations to the infrasound source is 624.5797 km.
Remotesensing 14 03181 g017
Table 1. UTTR source location and observations (See also in Figure 13).
Table 1. UTTR source location and observations (See also in Figure 13).
StationLocation
Source location41.1310°N, 112.8950°W
BGU40.9204°N, 112.0309°W
BRP39.4727°N, 110.7409°W
HWU41.6071°N, 111.5652°W
WMU40.0795°N, 111.8310°W
Table 2. Infrasound signal and the SBL algorithm parameter setting.
Table 2. Infrasound signal and the SBL algorithm parameter setting.
Station NumberSnapshotsArray ApertureAngular ResolutionSensor
BGU39135 m0.1BGU1, BGU3
BRP39150 m0.1BRP1, BRP3
HWU39140 m0.1HWU1, HWU3
WMU39150 m0.1WMU2, WMU3
Table 3. Arrival azimuth estimation error of the SBL algorithm for four stations.
Table 3. Arrival azimuth estimation error of the SBL algorithm for four stations.
Station NumberArrival Azimuth Estimation Results (°)Actual Azimuths (°)Estimation Error (°)
BGU32.9131.990.92
BRP307.35307.510.16
HWU252.07250.541.53
WMU314.35314.660.31
Table 4. Comparison of the error of arrival azimuth estimation between SBL algorithm, FK analysis, and beamforming algorithm.
Table 4. Comparison of the error of arrival azimuth estimation between SBL algorithm, FK analysis, and beamforming algorithm.
Station NumberSBL Arrival Azimuth Estimation Error (°)FK Arrival Azimuth Estimation Error (°)Beamforming Arrival Azimuth Estimation Error (°)
BGU0.922.943.48
BRP0.164.783.23
HWU1.532.714.43
WMU0.319.181.30
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Wang, R.; Yi, X.; Yu, L.; Zhang, C.; Wang, T.; Zhang, X. Infrasound Source Localization of Distributed Stations Using Sparse Bayesian Learning and Bayesian Information Fusion. Remote Sens. 2022, 14, 3181. https://doi.org/10.3390/rs14133181

AMA Style

Wang R, Yi X, Yu L, Zhang C, Wang T, Zhang X. Infrasound Source Localization of Distributed Stations Using Sparse Bayesian Learning and Bayesian Information Fusion. Remote Sensing. 2022; 14(13):3181. https://doi.org/10.3390/rs14133181

Chicago/Turabian Style

Wang, Ran, Xiaoquan Yi, Liang Yu, Chenyu Zhang, Tongdong Wang, and Xiaopeng Zhang. 2022. "Infrasound Source Localization of Distributed Stations Using Sparse Bayesian Learning and Bayesian Information Fusion" Remote Sensing 14, no. 13: 3181. https://doi.org/10.3390/rs14133181

APA Style

Wang, R., Yi, X., Yu, L., Zhang, C., Wang, T., & Zhang, X. (2022). Infrasound Source Localization of Distributed Stations Using Sparse Bayesian Learning and Bayesian Information Fusion. Remote Sensing, 14(13), 3181. https://doi.org/10.3390/rs14133181

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop