Method for Bayesian inversion of sound source distance by using ice layer bending wave time-frequency spectrum

Abstract

The invention belongs to the field of polar acoustic detection, and provides a method for carrying out Bayesian inversion on a sound source distance by utilizing an ice layer bending wave time-frequency spectrum, which comprises the following steps of: acquiring a bending wave signal propagated in an ice layer by utilizing a sensor arranged in the ice layer; the collected bending wave signals are filtered to improve the signal-to-noise ratio, the filtered bending wave signals s (t) are subjected to time-frequency analysis to obtain the time-frequency spectrum of bending waves, and t represents time; calculating a phase velocity frequency dispersion curve c (f) based on the six acoustic parameters and a bending wave frequency dispersion formula, wherein f represents frequency; generating a pseudo transmission pulse w (t); generating a receiving signal; performing Bayesian inversion; determining a final output solution; the probability density of the parameters is obtained based on the iteration process, and the average value of the probability density is the expectation of the inversion parameters according to the traversal theorem. The method is low in cost; convergence is fast, global optimization can avoid falling into a local solution, and inversion parameters are few; and the requirement of high-cost ocean survey is reduced.

Description

Technical Field

[0001] The invention belongs to the field of polar acoustic detection and relates to a passive ranging method, in particular to a method for inverting the distance of a sound source by using the Bayesian inversion of the time spectrum of ice bending waves. Background Art

[0002] The Arctic region has undergone significant environmental changes in recent years. Multi-year ice is no longer common in the Arctic, instead Arctic ice melts in the summer and refreezes each winter. The changing Arctic environment has prompted new research into acoustic detection, identification, and tracking of anthropogenic sources. Understanding the propagation of sound in first-year ice is important, particularly in nearshore and shallow water areas where anthropogenic activity near shore is expected to increase, perhaps through Arctic shipping, natural resource exploration, and tourism.

[0003] Understanding how to detect, identify, and track anthropogenic sources in these environments, in the absence of large arrays, has a direct impact on military defense capabilities in the Arctic. Rapid advances in technology and methods now allow ice-surface seismometers to be deployed in polar conditions to continuously record seismic wavefields at a sampling rate of 500 Hz for more than 30 days, and the next generation of products is expected to be able to record for several months while continuously transmitting data via satellite communications. Having sensors that can detect and track data sources can improve people's perception awareness, strategic planning, and response capabilities. Sensor technologies that typically provide these capabilities consist of large sensor arrays and require intensive data processing. Reducing the number of sensors and miniaturizing packaging will increase versatility and rapid deployment capabilities in near-shore Arctic applications.

[0004] Research has found that there are bending waves with dispersion effects in the ice that can propagate over long distances, and the time-frequency spectrum of the bending waves is closely related to the propagation distance. Using the bending wave signals recorded by ice surface sensors to estimate the distance to the sound source has broad application prospects. Summary of the invention

[0005] In view of the above-mentioned existing research and technology, the technical problem to be solved by the present invention is to provide an ice acoustic detection method suitable for polar nearshore fixed ice areas, which utilizes the bending waves widely present in the ice layer.

[0006] In order to solve the above technical problems, the present invention provides a method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves, comprising the following steps:

[0007] Step 1: Use sensors installed in the ice layer to collect bending wave signals propagating in the ice layer; the sensors here include but are not limited to seismic detectors, accelerometers and other instruments that can collect ice surface particle vibrations.

[0008] Step 2: Filter the collected bending wave signal to improve the signal-to-noise ratio, and perform time-frequency analysis on the filtered bending wave signal s(t) to obtain the time-frequency spectrum of the bending wave, where t represents time. The filtering method here includes performing broadband filtering from a frequency domain perspective to filter out noise signals outside the bending wave bandwidth, and performing filtering from a time domain perspective to set the noise outside the bending wave duration to zero.

[0009] Step 3: Calculate the phase velocity dispersion curve c(f) based on the six acoustic parameters and the bending wave dispersion formula, where f represents the frequency; the six parameters include ice thickness d, density ρ, longitudinal wave velocity c s , shear wave velocity c t , and the speed of sound in water c L and density ρ L The parameters are based on historical data or field measurements with measuring instruments. The default value can be set to ρ = 910 kg / m 3 , c s =3800m / s,c t =1900m / s,c L =1450m / s,ρ L =1000kg / m 3 The dispersion formulas for flexural waves include the following two:

[0010] The first dispersion formula satisfies:

[0011]

[0012] The symbols are defined as angular frequency ω = 2πf, horizontal wave number k = ω / c, longitudinal wave number k of ice s =ω / c s , the shear wave number k of ice t =ω / c t , the wave number of water k L =ω / c L , the vertical wave number of the ice longitudinal wave Ice shear wave vertical wave number By solving the root of the equation, we can get the phase velocity c at different frequencies f, that is, c(f). L is the water density, ρ is the ice density, and d is the ice thickness. r1 represents the shear-bending coupling term, and r2 is the interface response term. It usually appears as a bending wave between 1 Hz and 100 Hz, and as a Scholte interface wave above 1000 Hz. This method is used to describe high-frequency bending wave signals.

[0013] The second dispersion formula satisfies:

[0014]

[0015] The first kind of consistency is defined by Poisson's ratio v = [0.5(c s / c t ) 2 -1] / [(c s / c t ) 2 -1] and Young's modulus The bending stiffness of ice is calculated as D = Ed 3 / 12(1-v 2 ), substitute it into the dispersion formula and solve for the phase velocity c at different frequencies f, i.e. c(f). Usually very low frequency long period waves below 1 Hz appear as bending gravity waves, and low frequencies of 1 Hz-100 Hz are classic bending wave modes. This method is used to describe low-frequency bending wave signals. Because the signal frequency bands of bending waves are different, the method for calculating the dispersion curve should be selected according to the actual frequency band. A simple principle is that when the frequency of the bending wave is greater than 100 Hz, the first method should be selected, and when the frequency of the bending wave is less than 10 Hz, the second method should be selected; both methods can be used for other frequency bands.

[0016] Step 4: Generate a pseudo-transmitted pulse w(t); It should be noted that because the excitation mechanism of the sound source is complex and diverse, it affects the polarization, directionality and signal amplitude of the sound wave propagation, which will cause the received signal to be completely different from the explosion sound. If the elastic wave encounters an ice crack or ice ridge during propagation, the amplitude information will also change. These complex factors are usually not fully understood. Therefore, in order to reduce the impact of these factors, the simulation generated signal needs to be improved as follows: Generate a pulse signal h(t) with the same time length as s(t), and its spectrum is expressed as Using the spectrum amplitude of the bending wave signal s(t) to replace the spectrum amplitude of the simulated signal will produce a pseudo-transmitted pulse w(t) that better describes the real signal. If the spectrum of the received signal s(t) is expressed as Then the spectrum of the pseudo-transmitted pulse is The time domain signal w(t) is obtained through inverse Fourier transform.

[0017] Step 5: Generate the received signal; the acoustic propagation distance R and a time shift parameter T generated due to the unknown excitation time are two parameters of Bayesian inversion, and these two parameters are closely related to the received signal. Based on the phase velocity dispersion curve c(f) obtained in step 3, the pseudo-transmitted pulse spectrum W(f) obtained in step 4, and the propagation distance R and time shift parameter T set in step 5, the received signal r is generated. 1 (t) is expressed as:

[0018]

[0019] Where real means taking the real part of the data.1 (t) is cyclically shifted by T to obtain the final received signal r(t). The cyclic shift is expressed as:

[0020] r((t+T)%T max )=r 1 (t)

[0021] Tmax is the duration of the bending wave signal s(t), and % represents the remainder operation.

[0022] Step 6: Perform Bayesian inversion; based on the time-frequency analysis of the bending wave, define the cost function between the generated received signal r(t) and the true signal s(t) as the mean square value of the time-frequency spectrum error:

[0023]

[0024] Where A represents the time-frequency spectrum of the real bending wave signal in step 2, B(X) represents the time-frequency spectrum of the generated signal, and X represents the parameters to be inverted: the sound propagation distance R and a time shift parameter T generated due to the unknown excitation time. When the sound source distance and time shift correspond to the true values, the cost function has a minimum value. According to this description, there is a one-to-one correspondence between the model parameters and the global minimum of the cost function.

[0025] In order to solve the problem of inverting the sound propagation distance and time shift parameters, the present invention adopts Bayesian reasoning, and the basic principle is expressed as follows:

[0026]

[0027] Where P(X|A) is the posterior probability, which represents the probability of estimating the parameter given the true time-frequency spectrum. P(A|X) is the likelihood function, assuming that the measurement error is uncorrelated and random, so it is modeled with a normal distribution with zero mean of the cost function value:

[0028]

[0029] P(X) is the prior probability. It is assumed that the model parameters R and T are uniformly distributed in a limited range. The propagation distance R is generally set within 1 km. The time shift parameter T ranges from -Tmax to Tmax, where Tmax is the duration of the bending wave signal s(t). The initial parameters are set to R = 100m and T = 0. P(A) is the marginal likelihood function, which is essentially a normalization factor that is the same for all probabilities. Simulated annealing Bayesian global optimization is used, and the σ in the cost function uses the simulated annealing temperature value. The initial temperature is set to H 0 =0.4, the final temperature is H N =0.0002, H k =H 0 e -αk, k = 1, 2, ... N, Compared with the logarithmic and inverse proportional temperature drop curves, the negative exponential type ensures that the simulated annealing algorithm has a faster convergence speed. The number of iterations is set to 20,000.

[0030] Step 7: Determine the final output solution. The probability density of the parameters is obtained based on the iterative process. According to the ergodic theorem, the average value of the probability density is the expectation of the inversion parameters.

[0031] Beneficial effects of the present invention: The present invention proposes a method for inverting the sound source distance using the Bayesian time-frequency spectrum of ice bending waves. The bending wave signals widely present in the ice layer can be collected using a single ice surface sensor, and the miniaturized equipment is easy to deploy and has low cost. The Bayesian inversion of the sound source distance has the advantages of fast convergence, global optimization to avoid falling into a local solution, and fewer inversion parameters. No detailed prior knowledge is required, reducing the need for high-cost ocean surveys. Finally, the effectiveness of the method was verified through field tests. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 is the inversion flow chart of the present invention;

[0033] Figure 2 The method is to collect bending wave signals in the present invention;

[0034] Figure 3 is the filtered noise-reduced bending wave signal of the present invention;

[0035] Figure 4 is the calculation result of the bending wave dispersion curve in the present invention;

[0036] Figure 5 The present invention generates and transmits a pseudo pulse;

[0037] Figure 6 is the generation of a received signal in the present invention;

[0038] Figure 7 is the Bayesian inversion process in the present invention;

[0039] Figure 8 is the probability density distribution in the present invention;

[0040] Fig. 9 This is the comparison and verification between the real signal and the generated signal in the present invention. DETAILED DESCRIPTION

[0041] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0042] like Figure 1 As shown, the present invention provides a method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves, comprising the following steps:

[0043] like Figure 2 As shown, in step 1, the bending wave signal with a propagation distance of 100m is collected in the field ice test. The bending wave signal propagating in the ice layer is collected by using a sensor arranged in the ice layer; the sensor here includes but is not limited to various instruments capable of collecting ice surface particle vibrations, such as a seismic detector and an accelerometer.

[0044] like Figure 3 As shown, step 2 filters the collected bending wave signal to improve the signal-to-noise ratio, performs time-frequency analysis on the filtered bending wave signal s(t), and obtains the time-frequency spectrum of the bending wave, where t represents time; the filtering method here includes performing broadband filtering from the frequency domain perspective to filter out the noise signal outside the bending wave bandwidth, and filtering from the time domain perspective to set the noise outside the bending wave duration to zero. Here, the maximum filtering frequency is set to 400 Hz and the minimum frequency is set to 50 Hz; and the noise signal after 0.5 s is set to zero to improve the signal-to-noise ratio.

[0045] like Figure 4 As shown in Figure 1, the phase velocity dispersion curve c(f) is calculated based on six acoustic parameters and the bending wave dispersion formula, where f represents the frequency. The six parameters include ice thickness d, density ρ, longitudinal wave velocity c s , shear wave velocity c t , and the speed of sound in water c L and density ρ L The five parameters are based on historical data or field measurements with measuring instruments. The default value can be set to ρ = 910 kg / m 3 , c s =3800m / s,c t =1900m / s,c L =1450m / s,ρ L =1000kg / m 3 The dispersion formulas for flexural waves include the following two:

[0046] The first dispersion formula satisfies:

[0047]

[0048] The symbols are defined as angular frequency ω = 2πf, horizontal wave number k = ω / c, longitudinal wave number k of ice s =ω / c s , the shear wave number k of ice t =ω / c t , the wave number of water k L =ω / c L , the vertical wave number of the ice longitudinal wave Ice shear wave vertical wave number By solving the root of the equation, we can get the phase velocity c at different frequencies f, that is, c(f). L is the water density, ρ is the ice density, and d is the ice thickness. r1 represents the shear-bending coupling term, and r2 is the interface response term. Figure 4 As shown by the solid line, it usually appears as a bending wave between 1 Hz and 100 Hz, and as a Scholte interface wave above 1000 Hz. This method is used to describe high-frequency bending wave signals.

[0049] The second dispersion formula satisfies:

[0050]

[0051] The first kind of consistency is defined by Poisson's ratio v = [0.5(c s / c t ) 2 -1] / [(c s / c t ) 2 -1] and Young's modulus The bending stiffness of ice is calculated as D = Ed 3 / 12(1-v 2 ), substitute it into the dispersion formula, and solve to get the phase velocity c at different frequencies f, i.e. c(f). Figure 4 As shown by the dotted line, usually very low frequency long period waves below 1 Hz appear as flexural gravity waves, and low frequencies of 1 Hz-100 Hz are classic flexural wave modes. This method is used to describe flexural wave signals in the low frequency band.

[0052] Because the signal frequency bands of bending waves are different, the method for calculating the dispersion curve should be selected according to the actual frequency band. A simple principle is that when the frequency of the bending wave is greater than 100Hz, the first method should be selected, and when the frequency of the bending wave is less than 10Hz, the second method should be selected; both methods can be used for other frequency bands. Since the maximum frequency of the bending wave collected in the experiment is 400Hz, the first method is selected to calculate the dispersion curve of the ice bending wave.

[0053] like Figure 5 As shown, step 4 generates a pseudo-transmitted pulse w(t). It should be noted that because the excitation mechanism of the sound source is complex and diverse, it affects the polarization, directionality and signal amplitude of the sound wave propagation, which will cause the received signal to be completely different from the explosion sound. If the elastic wave encounters an ice crack or ice ridge during propagation, the amplitude information will also change. These complex factors are usually not fully grasped. Therefore, in order to reduce the impact of these factors, the simulation generated signal needs to be improved as follows: generate a pulse signal h(t) with the same time length as s(t), and its spectrum is expressed as Using the spectrum amplitude of the bending wave signal s(t) to replace the spectrum amplitude of the simulated signal will produce a pseudo-transmitted pulse w(t) that better describes the real signal. If the spectrum of the received signal s(t) is expressed as Then the spectrum of the pseudo-transmitted pulse is The time domain signal w(t) is obtained through inverse Fourier transform.

[0054] like Figure 6 As shown, step 5 generates the received signal. The acoustic propagation distance R and a time shift parameter T generated due to the unknown excitation time are two parameters of Bayesian inversion, and these two parameters are closely related to the received signal. According to the phase velocity dispersion curve c(f) obtained in step 3 and the pseudo-transmitted pulse spectrum W(f)=F[w(t)] obtained in step 4, F represents the Fourier transform operation, the time domain signal w(t) is transformed into the spectrum W(f), and the propagation distance R and time shift parameter T set in step 5, the received signal r is generated. 1 (t) is expressed as:

[0055]

[0056] Where real means taking the real part of the data. 1 (t) is cyclically shifted by T to obtain the final received signal r(t). The cyclic shift is expressed as:

[0057] r((t+T)%T max )=r 1 (t)

[0058] Tmax is the duration of the bending wave signal s(t), and % represents the remainder operation.

[0059] like Figure 7 As shown, step 6 performs Bayesian inversion; based on the time-frequency analysis of the bending wave, the cost function between the generated received signal r(t) and the true signal s(t) is defined as the mean square value of the time-frequency spectrum error:

[0060]

[0061] Where A represents the time-frequency spectrum of the real bending wave signal in step 2, B(X) represents the time-frequency spectrum of the generated signal, and X represents the parameters to be inverted: the sound propagation distance R and a time shift parameter T generated due to the unknown excitation time. When the sound source distance and time shift correspond to the true values, the cost function has a minimum value. According to this description, there is a one-to-one correspondence between the model parameters and the global minimum of the cost function.

[0062] In order to solve the problem of inverting the sound propagation distance and time shift parameters, the present invention adopts Bayesian reasoning, and the basic principle is expressed as follows:

[0063]

[0064] Where P(X|A) is the posterior probability, which represents the probability of estimating the parameter given the true time-frequency spectrum. P(A|X) is the likelihood function, assuming that the measurement error is uncorrelated and random, so it is modeled with a normal distribution with zero mean of the cost function value:

[0065]

[0066] P(X) is the prior probability. It is assumed that the model parameters R and T are uniformly distributed in a limited range. The propagation distance R is generally set within 1 km. The time shift parameter T ranges from -Tmax to Tmax, where Tmax is the duration of the bending wave signal s(t). The initial parameters are set to R = 100m and T = 0. P(A) is the marginal likelihood function, which is essentially a normalization factor that is the same for all probabilities. Simulated annealing (SA) Bayesian global optimization is used, and the σ in the cost function uses the temperature value of simulated annealing. The initial temperature is set to H 0 =0.4, the final temperature is H N =0.0002, H k =H 0 e -αk , k = 1, 2, ..., N, Compared with the logarithmic and inverse proportional temperature drop curves, the negative exponential type ensures that the simulated annealing algorithm has a faster convergence speed. The number of iterations is set to 20,000.

[0067] like Figure 8 As shown, step 7 obtains the probability density of the parameters based on the Bayesian iterative process. According to the ergodic theorem, the average value of the probability density is the expectation of the inversion parameter, and the final output solution is determined to obtain the final estimated distance of about 107m, which is consistent with the actual propagation distance. Fig. 9 As shown in the figure, a waveform comparison between the generated signal calculated based on the inversion parameters and the real bending wave signal is given, and it is found that the waveform inversion results are relatively consistent.

[0068] The measured data results show that the method of the present invention can effectively locate the distance of the sound source through the accelerometer on the ice. Although the accuracy and robustness of the present invention have been demonstrated by the ice test data, suppressing the mismatch between the generated data and the real data is crucial to reducing the uncertainty of the results. Therefore, future work should focus on two main directions. The first is to improve the signal-to-noise ratio, such as obtaining the best results through automatic selection and denoising, such as machine learning-based methods. The second is to use an effective forward model that can explain the local changes in the properties of ice while keeping the computational cost at an acceptable level. Another interesting point is that lower-frequency elastic wave propagation can propagate to a farther range, such as sea ice deformation and rupture has passive advantages and does not require active sources. These low-frequency sound sources are widespread. This will open the way for tomography including sea ice thickness and elastic parameters. In addition, the horizontal component may provide ideas for the inversion of ice thickness and elastic properties by including waveforms of the other two basic modes (transverse waves and longitudinal plate waves) in the cost function.

[0069] It should be understood that the application of the present invention is not limited to the above examples. For ordinary technicians in this field, improvements or changes can be made based on the above description. All these improvements and changes should fall within the scope of protection of the claims attached to the present invention.

Claims

1. A method for inverting the distance of sound sources using the Bayesian inversion of the time spectrum of ice bending waves, characterized in that: The following steps are involved: Step 1: Using sensors placed on the ice layer to collect bending wave signals propagating in the ice layer; Step 2: filtering the collected bending wave signal to improve the signal-to-noise ratio, and performing time-frequency analysis on the filtered bending wave signal s(t) to obtain the time-frequency spectrum of the bending wave, where t represents time; Step 3: Calculate the phase velocity dispersion curve c(f) based on the six acoustic parameters and the bending wave dispersion formula, where f represents the frequency; the six parameters include ice thickness d, density ρ, longitudinal wave velocity c s , shear wave velocity c t , and the speed of sound in water c L and density ρ L ; The parameters are based on historical data or field measurements with measuring instruments; the default value is set to ρ = 910 kg / m 3 , c s =3800m / s,c t =1900m / s,c L =1450m / s,ρ L =1000kg / m 3 ; Step 4: Generate a pseudo transmit pulse w(t); Step 5: Generate a receiving signal; Step 6: Perform Bayesian inversion; Step 7: Determine the final output solution; The probability density of the parameters is obtained based on the iterative process. According to the ergodic theorem, the average value of the probability density is the expectation of the inversion parameters.

2. The method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves according to claim 1 is characterized in that: The filtering method in step 2 includes performing broadband filtering from a frequency domain perspective to filter out noise signals outside the bending wave bandwidth, and performing filtering from a time domain perspective to set the noise outside the bending wave duration to zero.

3. The method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves according to claim 1 is characterized in that: The dispersion formulas for bending waves in step 3 include the following two: The first dispersion formula satisfies: The symbols are defined as angular frequency ω = 2πf, horizontal wave number k = ω / c, longitudinal wave number k of ice s =ω / c s , the shear wave number k of ice t =ω / c t , the wave number of water k L =ω / c L , the vertical wave number of the ice longitudinal wave Ice shear wave vertical wave number By solving the root of the equation, we can obtain the phase velocity c at different frequencies f, i.e. c(f); ρ L is the water density, ρ is the ice density, and d is the ice thickness. r1 represents the shear-bending coupling term, and r2 is the interface response term; The second dispersion formula satisfies: The first kind of consistency is defined by Poisson's ratio v = [0.5(c s / c t ) 2 -1] / [(c s / c t ) 2 -1] and Young's modulus The bending stiffness of ice is calculated as D = Ed 3 / 12(1-v 2 ), and substitute it into the dispersion formula to obtain the phase velocity c at different frequencies f, i.e. c(f).

4. The method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves according to claim 1 is characterized in that: In step 4, the simulation generated signal is improved as follows: a pulse signal h(t) with the same time length as s(t) is generated, and its spectrum is expressed as The spectrum amplitude of the simulated signal is replaced by the spectrum amplitude of the bending wave signal s(t), which will produce a pseudo-transmitted pulse w(t) that better describes the real signal; if the spectrum of the received signal s(t) is expressed as Then the spectrum of the pseudo-transmitted pulse is W(f) = M2(f) The time domain signal w(t) is obtained through inverse Fourier transform.

5. The method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves according to claim 1 is characterized in that: The specific method of step 5 is as follows: the acoustic propagation distance R and a time shift parameter T generated due to the unknown excitation time are two parameters of Bayesian inversion. According to the phase velocity dispersion curve c(f) obtained in step 3, the pseudo-transmitted pulse spectrum W(f) obtained in step 4, and the propagation distance R and time shift parameter T set in step 5, the received signal r1(t) is generated and expressed as: Where real represents taking the real part of the data; r1(t) is cyclically shifted by T to obtain the final received signal r(t), and the cyclic shift is expressed as: r((t+T)%T max )=r1(t) Tmax is the duration of the bending wave signal s(t), and % represents the remainder operation.

6. The method for inverting the sound source distance using the Bayesian inversion of the time spectrum of ice bending waves according to claim 1 is characterized in that: The specific method of step 6 is: based on the time-frequency analysis of the bending wave, the cost function between the generated received signal r(t) and the true signal s(t) is defined as the mean square value of the time-frequency spectrum error: Where A represents the time-frequency spectrum of the real bending wave signal in step 2, B(X) represents the time-frequency spectrum of the generated signal, and X represents the parameters to be inverted: the sound propagation distance R and a time shift parameter T generated due to the unknown excitation time; when the sound source distance and time shift correspond to the true values, the cost function has a minimum value.

Images