Underwater target sensing method based on Monte Carlo

Abstract

The invention relates to the field of ocean exploration, and discloses an underwater target sensing method based on Monte Carlo, which comprises the steps of data extraction and uncertainty evaluation thereof, sound wave travel time calculation, underwater moving target position modeling, position inversion and uncertainty evaluation. As an underwater target has strong concealment, the position of the underwater target is difficult to detect through visual observation and electromagnetic waves. However, the interaction between the underwater moving target and the water body generates a sound wave signal, and the position of the underwater moving target can be accurately detected based on sonar or ocean bottom seismograph detection data by using the characteristic that sound waves are almost not attenuated in the water body. Through the propagation time of sound waves in a water body, the position of an underwater target can be modeled based on a nonlinear Monte Carlo inversion method and a Bayesian formula, the position of the underwater target is obtained, and meanwhile, the uncertainty of the position of the underwater target is evaluated through probability statistics. The method can play an important role in underwater vehicle positioning, marine organism migration research, ocean exploration and seabed navigation.

Description

Technical Field

[0001] The present invention relates to the fields of underwater vehicle positioning, marine life migration research, ocean exploration, seabed navigation and ocean investigation, and in particular to an underwater target perception method based on Monte Carlo. Background Art

[0002] With the deepening of marine research and development, there is a growing need to locate and sense underwater targets. However, due to the rapid attenuation of light and electromagnetic waves in water due to the seawater layer, it is difficult to sense the position of underwater targets (such as submersibles, underwater robots, marine life, and ocean detection instruments) through visual and electromagnetic wave means. Traditionally, underwater targets are sensed through acoustic signals, taking advantage of the weak attenuation of acoustic signals in water.

[0003] At present, the mainstream acoustic signal underwater target perception uses the least squares method, gradient descent, simulated annealing and genetic algorithm based on the principle of three-point positioning. As for linear inversion methods such as least squares and gradient descent, in order to obtain gradient information, it is usually necessary to ignore high-order terms and linearize the nonlinear equations to solve them. Only approximate results are obtained, and the results are strongly dependent on the initial position. For underwater targets without any prior information, the results are easy to fall into local extreme values, and the obtained solutions are not optimal. Simulated annealing, genetic algorithms, etc. can avoid falling into local extreme values ​​by randomly sampling the entire model space, but there are still some shortcomings, such as the sampling method seriously affects the efficiency and cannot provide uncertainty assessment of the results. In addition, when the acoustic signal generated by the underwater target is weak, its signal-to-noise ratio is poor, which will have a strong impact on the reliability of the results. Conventional gradient inversion and nonlinear inversion cannot give the uncertainty of the results. As for some current solutions, in addition to the shortcomings in inversion technology, most of the solutions also strongly rely on advanced equipment (such as ultra-short baselines, low-frequency hydrophones, etc.) and high-precision observation data. As a result, the scheme is too complicated and difficult to implement, and it is troublesome in actual operation because it requires additional equipment layout and field working conditions. In addition, it is difficult to promote because of the high price of the equipment. Summary of the invention

[0004] In view of the shortcomings of the above-mentioned conventional methods, the purpose of the present invention is to provide an underwater target perception method based on Monte Carlo. Compared with the existing technology, it uses publicly available data, does not rely on external equipment assistance, and the data processing will perform signal-to-noise ratio analysis and multiple data screening. The results are more stable, reliable, and accurate, and can also provide the uncertainty of the results.

[0005] The present invention adopts the following technical solution:

[0006] A Monte Carlo-based underwater target perception method comprises the following steps:

[0007] Step 1: Detecting the abnormal arrival time of the sound wave from the sonar or seafloor seismometer, i.e., the receiver recording data;

[0008] Step 2: Control the signal-to-noise ratio of the detected abnormal information waveform and retain the abnormal sound wave data with a good signal-to-noise ratio;

[0009] Step 3: Perform coordinate transformation and mesh the model;

[0010] Step 4: Calculate the acoustic wave propagation time from each grid node in the model to the receiver, and use the consistency of the time of generating abnormal signals to obtain the grid node with the minimum propagation time difference;

[0011] Step 5: Use spatial multi-scale iteration to further search for the best estimated position of the underwater target;

[0012] Step 6: Calculate the propagation time using the best estimated position of the underwater target, and remove the information with large propagation time differences to ensure the stability and accuracy of subsequent calculation results;

[0013] Step 7: Further invert the exact position of the underwater target through the Monte Carlo algorithm.

[0014] Step 1: Extracting the acoustic wave signal from the data recorded by the receiver to the time T i (φ, θ, d), where T i represents the abnormal acoustic wave time recorded by the i-th receiver, φ, θ and d represent the longitude, latitude and water depth of the receiver respectively; the extracted abnormal acoustic wave time is assumed to satisfy the Gaussian distribution, and the probability density function is obtained by the mean and standard deviation of all extracted times:

[0015]

[0016] In formula (1), μ and σ represent the mean and standard deviation respectively; T is the extracted acoustic anomaly time; ρ is the probability density function.

[0017] Step 2: The waveform data of a period of time before and after the detected abnormal sound wave signal is used as noise and effective signal respectively, and the signal-to-noise ratio is calculated; the information whose signal-to-noise ratio does not meet the requirements is removed from the data set. The signal-to-noise ratio is calculated as shown in formula (2).

[0018]

[0019] In formula (2), R sn represents the calculated signal-to-noise ratio, N s and N n Represents the number of samples of signal and noise respectively, s i and n i They represent the amplitude of the signal and noise at the i-th point respectively, and It represents the average value of signal and noise during this period of time.

[0020] Step 3: According to the spherical assumption, the longitude, latitude and water depth of the receiver are transformed into (x, y, z) in rectangular coordinates, where the directions of x, y and z in the rectangular coordinates point from the center of the sphere to the east, south and north respectively; the radius of the receiver in the spherical coordinates is r=Rd, where R=6371.0 km is the radius of the earth; then, in the rectangular coordinate system, the model space is gridded with a step size of (Δx, Δy, Δz), and the propagation time from all grids to the receiver is calculated as shown in formula (3):

[0021]

[0022] Where x n ,y n and z n Indicates the position of the nth grid in rectangular coordinates; x i ,y i and z i represents the position of the i-th receiver in rectangular coordinates; v w It represents the speed of sound waves in water. represents the propagation time from the nth grid to the i-th receiver.

[0023] Step 4: Assume that the time when the underwater target and the water body generate the sound wave signal is T0(x0, y0, z0), wherein x0, y0 and z0 represent the position in rectangular coordinates when the underwater target generates the sound wave signal; Calculate the estimated acoustic signal generation time of the nth grid pair for the i-th receiver, and calculate the minimum value of the difference in the estimated acoustic signal generation time between each grid pair of receivers to evaluate the approximate position of the underwater target, that is, Δt n The position of the grid at the minimum is as shown in formula (4):

[0024]

[0025] Where N is the number of receivers; Δt n The sum of the time differences between the estimated acoustic wave signals of all receivers in the nth grid. If the nth grid happens to be the location of the underwater target, then Δt n Theoretically it should be equal to 0s.

[0026] Step 5: Change the value of Δt in step 4 to n The smallest grid is further subdivided, and the calculation in step 4 is repeated with the subdivided grid. n The calculation stops when ∈ is the threshold value. At this time, the corresponding grid position is the estimated optimal underwater target position.

[0027] Step 6: Using the estimated acoustic signal generation time of each receiver at the grid position obtained in step 5, obtain the probability density function of the Gaussian distribution, and remove the data falling outside [μ-1.5σ, μ+1.5σ], where μ and σ represent the mean and standard deviation respectively, to ensure data reliability.

[0028] The step 7: using the Monte Carlo inversion algorithm to further locate the underwater target position according to the grid position obtained in step 5, and using the probability distribution of each parameter to evaluate its uncertainty.

[0029] Beneficial effects of the present invention:

[0030] 1) Using spatial multi-scale iterative grid search, the estimated position of underwater targets can be automatically refined until the error reaches the preset accuracy.

[0031] 2) The underwater target position obtained by spatial multi-scale iteration is used as the initial value of the Metropolis-Hastings algorithm, which can obtain accurate results with less calculation amount.

[0032] 3) Monte Carlo random sampling in the model space can effectively avoid falling into local extreme values, and the results obtained are more reliable.

[0033] 4) The posterior probability is calculated by the Bayesian formula, which can effectively evaluate the uncertainty of underwater targets in longitude, latitude and water depth. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a schematic diagram of the process of the present invention.

[0035] Figure 2 It is a schematic diagram of converting longitude, latitude and water depth into rectangular coordinates.

[0036] Figure 3 It is a schematic diagram of rectangular coordinate grid division, where the fine grid is a further refinement of the nth grid in the coarse grid.

[0037] Figure 4 The distribution of 500,000 sampling points in the seafloor seismometer relocation. The triangles represent the locations of active source shots.

[0038] Figure 5 is the histogram distribution of the sampling positions of the seafloor seismometers in the X, Y and Z directions.

[0039] Figure 6 The distribution of 500,000 sampling points in underwater target relocation. The triangles represent the locations of surface buoys and the circles represent the locations of seafloor seismometers.

[0040] Figure 7 is the sampling position histogram distribution of underwater targets in the X, Y and Z directions. DETAILED DESCRIPTION

[0041] In order to further understand the technical content, features and functions of the present invention, the following examples are listed and described in detail with reference to the accompanying drawings.

[0042] Example 1

[0043] like Figure 1 The Monte Carlo-based underwater target perception method shown in the figure provides a technical solution for accurately and reliably perceiving underwater targets. After the ship reaches the deployment location, the seafloor seismometer is released on the deck and freely falls to the seafloor. However, the actual position of the instrument at the seafloor is very different from the position obtained by the ship-borne positioning system. Based on the perception of the detection instrument placed on the seafloor, this case uses the seafloor seismometer as an example for numerical simulation, and locates its seafloor position through this method. The specific steps are as follows:

[0044] Step 1: Pick up the arrival time of the direct P wave from the active source shot recorded by the seafloor seismometer and calculate its probability density function. This case uses a total of 9 active source shot points to simulate data (all meet the signal-to-noise ratio requirements).

[0045] Step 2: According to the location of the active source and the placement point, the geographic coordinates are converted into rectangular coordinates for subsequent calculation of the direct P-wave propagation time. Figure 2 It is a schematic diagram of converting longitude, latitude and water depth into rectangular coordinates. Figure 3 is a schematic diagram of the rectangular coordinate grid division, and the location of the earthquake source in the rectangular coordinate system. Figure 4 As shown in the triangle in the figure, the focal depth is set to 0 km.

[0046] Step 3: Determine the model range X: 0-20km, Y: 0-20km and Z: 0-20km according to the distribution range of the earthquake source and the deployment point. At the same time, grid the model at an appropriate spacing and calculate the P-wave propagation time from each grid to the earthquake source according to the seawater sound wave speed of 1.5km / s.

[0047] Step 4: The preliminary positioning result is obtained by minimizing the difference between the actual earthquake source explosion measurement time and the calculation time. There is still a certain gap with the actual (10km, 10km, 10km) (it should be noted here that the water depth of most sea basins is about 3-4km, and the error in the water depth direction in positioning is relatively large. In order to prove the reliability and accuracy of this method, the water depth used in this case is 10km, which is much greater than the actual water depth).

[0048] Step 5: Further refine the initial positioning grid and iterate twice according to step 4 to obtain a relatively accurate position (9km, 8km, 9km). There is still a certain gap.

[0049] Step 6: According to the threshold range of [μ-1.5σ, μ+1.5σ], information with poor signal-to-noise ratio is eliminated. In this implementation case, the selected signals are all within the range, so no data is eliminated.

[0050] Step 7: Use the Metropolis Hastings algorithm to perform inversion, taking the (9km, 8km, 9km) obtained in step 5 as the initial sampling point, and perform 500,000 samplings in the model space (such as Figure 4 The final position of the seafloor seismometer is (10±0.2km, 10±0.2km, 10±0.4km), which is consistent with the real position (as shown in Figure 5 shown).

[0051] Example 2

[0052] In order to further illustrate the positioning process of underwater targets (such as underwater suspended instruments, whales or underwater vehicles) in the present invention, this case takes the suspended instrument as an example to perform numerical simulation. Combined with the surface sonar and seabed seismograph data, the data processing and positioning process are described in detail.

[0053] Step 1: Extract P-wave arrival information from sonar and seafloor seismometer data (all meeting the signal-to-noise ratio requirement) and calculate their probability density function. Nine sonars float on the water surface and four seafloor seismometers are located 15 km below the seafloor (the water depth is greater than the actual seafloor depth to demonstrate the accuracy of this method in the depth direction).

[0054] Step 2: Based on the sonar and seafloor seismic positions (such as Figure 6 As shown in the figure, where the triangle represents the sonar position and the circle represents the seafloor seismometer position), its geographic coordinates and water depth are converted into rectangular coordinates to facilitate the solution of the P-wave propagation time.

[0055] Step 3: Determine the model range X: 0-20km, Y: 0-20km and Z: 0-20km according to the location of the sonar and seafloor seismometer, then mesh the model and calculate the propagation time from each grid to the sonar and seafloor seismometer based on the actual propagation speed of sound waves in water, 1.5km / s.

[0056] Step 4: Based on the time when the underwater target generates a signal as a constant, calculate the estimated time of each sonar and seafloor seismometer, and calculate the corresponding difference. Rough positioning is performed based on the principle of minimizing the difference. The result is still far from (10km, 10km, 10km).

[0057] Step 5: The grid obtained in step 4 is further refined and a relatively more accurate position (9km, 9km, 9km) is obtained through 2 iterations. However, the result is still not accurate enough.

[0058] Step 6: Based on the signal generation time estimated in step 4, the signals that do not meet the requirements are eliminated according to the threshold range. In this experiment, Gaussian noise is added with a standard deviation of 0.1s.

[0059] Step 7: Use the Metropolis Hastings algorithm to invert the signals that meet the requirements. After sampling the model space 500,000 times, the final positioning result is (10.0±0.1km, 10±0.1km, 10±0.1km), which is exactly the same as the preset (10km, 10km, 10km). Compared with the above case, its standard deviation is smaller, mainly due to the constraints from the seafloor seismograph data.

[0060] The above-described embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the invention. It should be pointed out that for ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the attached claims.

Claims

1. A Monte Carlo-based underwater target perception method, characterized in that: The following steps are involved: Step 1: Detecting the abnormal arrival time of the sound wave from the sonar or seafloor seismometer, i.e., the receiver recording data; Step 2: Control the signal-to-noise ratio of the detected abnormal information waveform and retain the abnormal sound wave data with a good signal-to-noise ratio; Step 3: Perform coordinate transformation and mesh the model; Step 4: Calculate the acoustic wave propagation time from each grid node in the model to the receiver, and use the consistency of the time of generating abnormal signals to obtain the grid node with the minimum propagation time difference; Step 5: Use spatial multi-scale iteration to further search for the best estimated position of the underwater target; Step 6: Calculate the propagation time using the best estimated position of the underwater target, and remove the information with large propagation time differences to ensure the stability and accuracy of subsequent calculation results; Step 7: Further invert the exact position of the underwater target through the Monte Carlo algorithm.

2. The Monte Carlo-based underwater target perception method according to claim 1, characterized in that: Step 1: Extracting the acoustic wave signal from the data recorded by the receiver to the time T i (φ,θ,d), where T i represents the abnormal acoustic wave time recorded by the i-th receiver, φ, θ and d represent the longitude, latitude and water depth of the receiver respectively; the extracted abnormal acoustic wave time is assumed to satisfy the Gaussian distribution, and the probability density function is obtained by the mean and standard deviation of all extracted times: In formula (1), μ and σ represent the mean and standard deviation respectively; T is the extracted acoustic anomaly time; ρ is the probability density function.

3. The Monte Carlo-based underwater target perception method according to claim 2, characterized in that: Step 2: The waveform data of a period of time before and after the detected abnormal sound wave signal is used as noise and effective signal respectively, and the signal-to-noise ratio is calculated; the information whose signal-to-noise ratio does not meet the requirements is removed from the data set. The signal-to-noise ratio is calculated as shown in formula (2). In formula (2), R sn represents the calculated signal-to-noise ratio, N s and N n Represents the number of samples of signal and noise respectively, s i and n i They represent the amplitude of the signal and noise at the i-th point respectively, and It represents the average value of signal and noise during this period of time.

4. The Monte Carlo-based underwater target perception method according to claim 3, characterized in that: Step 3: According to the spherical assumption, the longitude, latitude and water depth of the receiver are transformed into rectangular coordinates (x, y, z), where the directions of x, y and z in the rectangular coordinates point from the center of the sphere to the east, south and north respectively; the radius of the receiver in the spherical coordinates is r = Rd, where R = 6371.0 km is the radius of the earth; then, in the rectangular coordinate system, the model space is gridded with a step size of (Δx, Δy, Δz), and the propagation time from all grids to the receiver is calculated as shown in formula (3): Where x n ,y n and z n Indicates the position of the nth grid in rectangular coordinates; x i ,y i and z i represents the position of the i-th receiver in rectangular coordinates; v w It represents the speed of sound waves in water. represents the propagation time from the nth grid to the i-th receiver.

5. The Monte Carlo-based underwater target perception method according to claim 4, characterized in that: Step 4: Assume that the time when the underwater target and the water body generate the sound wave signal is T0(x0, y0, z0), where x0, y0 and z0 represent the position of the underwater target in rectangular coordinates when the sound wave signal is generated; Calculate the estimated acoustic signal generation time of the nth grid pair for the i-th receiver, and calculate the minimum value of the difference in the estimated acoustic signal generation time between each grid pair of receivers to evaluate the approximate position of the underwater target, that is, Δt n The position of the grid at the minimum is as shown in formula (4): Where N is the number of receivers; Δt n The sum of the time differences between the estimated acoustic wave signals of all receivers in the nth grid. If the nth grid happens to be the location of the underwater target, then Δt n Theoretically it should be equal to 0s.

6. The Monte Carlo-based underwater target perception method according to claim 5, characterized in that: Step 5: Change the value of Δt in step 4 to n The smallest grid is further subdivided, and the calculation in step 4 is repeated with the subdivided grid. n The calculation stops when <∈, where ∈ is the threshold. At this time, the corresponding grid position is the estimated optimal underwater target position.

7. The Monte Carlo-based underwater target perception method according to claim 6, characterized in that: Step 6: Using the estimated acoustic signal generation time of each receiver at the grid position obtained in step 5, obtain the probability density function of the Gaussian distribution, and remove the data falling outside [μ-1.5σ,μ+1.5σ], where μ and σ represent the mean and standard deviation respectively, to ensure data reliability.

8. The Monte Carlo-based underwater target perception method according to claim 1, characterized in that: The step 7: using the Monte Carlo inversion algorithm to further locate the underwater target position according to the grid position obtained in step 5, and using the probability distribution of each parameter to evaluate its uncertainty.

Images