Sound velocity profile inversion method based on inverted multi-beam echo sounder

Abstract

A sound velocity profile inversion method based on an inverted multi-beam echometer. Said method comprises the following steps: mounting, in an inverted manner, the multi-beam echometer on an underwater submerged buoy or fixing same to a water bottom, transmitting a beam to a water surface by means of a transmitting transducer array, and receiving an echo signal by means of a receiving transducer array; the multi-beam echometer obtaining an angle of arrival and arrival time of an echo according to the received echo signal; solving an EOF according to sound velocity profile prior information, and obtaining a dimension reduction primary function description method of the sound velocity profile; in combination with the EOF, a ray tracing algorithm, a surface sound velocity and multi-beam data, establishing an optimization model; according to the established optimization model, the arrival time and the angle of arrival of the received echo and the surface sound velocity, using an optimization algorithm to obtain an estimation result of the sound velocity profile of a measurement area; and further, calculating a water temperature profile of the measurement area by using an estimated value of the sound velocity profile. Said method can rapidly and accurately track fluctuations of the sound velocity profile and the temperature profile.

Claims

1. A sound velocity profile inversion method based on an inverted multi-beam echo sounder, comprising following steps: (1) mounting the multi-beam echo sounder on an underwater submerged buoy or fixing the multi-beam echo sounder to a water bottom in an inverted manner; transmitting a beam to a water surface by a transmitting transducer array; and receiving an echo signal by a receiving transducer array; (2) obtaining, by the multi-beam echo sounder, an angle of arrival of an echo and a time of arrival of the echo according to the received echo signal; (3) obtaining a description method of a dimension reduction primary function of the sound velocity profile, by solving an Empirical Orthogonal Function (EOF) according to sound velocity profile prior information, wherein a sound velocity profile at a time t is expressed as: $\begin{array}{cc}\begin{array}{c}c\left(z,t\right)=c_0\left(z\right)+\Delta c\left(z,t\right)\\ =c_0\left(z\right)+\underset{i=1}{\overset{I}{.Math.}}{\alpha }_i\left(t\right){\psi }_i\left(z\right)\end{array}& \left(1\right)\end{array}$ where in the formula (1), z denotes a depth, c.sub.0(z) denotes an average sound velocity at the depth z, and Δc(z t) denotes a disturbance quantity; i denotes an order of the EOF, 1≤i≤I, I denotes the total number of orders of the EOF, α.sub.i(t) denotes a coefficient of the i.sup.th-order EOF, and ψ.sub.i(z) denotes the i.sup.th-order EOF; (4) establishing an optimization model in combination with the EOF, a sound ray tracing algorithm, a surface sound velocity and multi-beam data, $\begin{array}{cc}\begin{array}{c}\hat{\alpha }\left(t\right)=\arg \underset{\alpha }{\min }\left({.Math.\mathrm{To}{a}_{meas}\left(t\right)-To{a}_{est}\left(t\right).Math.}_2+K{.Math.c_{0meas}\left(t\right)-c_{0est}\left(t\right).Math.}_2\right)\\ \mathrm{To}{a}_{est}\left(t\right)=\mathscr{H}\left[\alpha ,\mathrm{Do}{a}_{meas}\left(t\right)\right]\end{array}& \left(2\right)\end{array}$ where in the formula (2), K denotes a weight coefficient, Toa.sub.meas(t) denotes actual measurement data of a time of arrival of each beam at a time t, Toa.sub.est(t) denotes an estimated value of the time of arrival of each beam at the time t, c.sub.0meas(t) denotes measurement data of the sound velocity at a position of the receiving transducer array of the multi-beam echo sounder at the time t, and c.sub.0est(t) denotes an estimated value of the sound velocity at the position of the receiving transducer array of the multi-beam echo sounder at the time t; H denotes a measurement operator, which performs a calculation according to the sound ray tracing algorithm, and Doa.sub.meas(t) denotes a measurement value of a direction of arrival of each beam at the time t; and (5) obtaining an inversion value of an EOF coefficient in an optimization algorithm, according to the optimization model established in step (4), the time of arrival of the echo and the angle of arrival of the echo received in step (2), and a measurement value of the sound velocity at the position of the receiving transducer array; and further obtaining an estimation result of the sound velocity profile by using the formula (1).

2. The sound velocity profile inversion method based on the inverted multi-beam echo sounder according to claim 1, wherein in step (5), the optimization algorithm comprises a simulated annealing algorithm, a hill climbing algorithm, a particle swarm optimization algorithm, a genetic algorithm, and an ant colony algorithm.

3. A measurement system of an inverted multi-beam echo sounder, wherein the sound velocity profile inversion method according to claim 1 is applied to the measurement system.

4. The measurement system of an inverted multi-beam echo sounder according to claim 3, wherein in step (5), the optimization algorithm comprises a simulated annealing algorithm, a hill climbing algorithm, a particle swarm optimization algorithm, a genetic algorithm, and an ant colony algorithm.

5. A water temperature profile inversion method based on an inverted multi-beam echo sounder, comprising following steps: (1) mounting the multi-beam echo sounder on an underwater submerged buoy or fixing the multi-beam echo sounder to a water bottom in an inverted manner, transmitting a beam to a water surface by a transmitting transducer array, and receiving an echo signal by a receiving transducer array; (2) obtaining, by the multi-beam echo sounder, an angle of arrival of an echo and a time of arrival of the echo according to the received echo signal; (3) obtaining a description method of a dimension reduction primary function of the sound velocity profile, by solving an Empirical Orthogonal Function (EOF) according to sound velocity profile prior information, wherein a sound velocity profile at a time t is expressed as: $\begin{array}{cc}\begin{array}{c}c\left(z,t\right)=c_0\left(z\right)+\Delta c\left(z,t\right)\\ =c_0\left(z\right)+\underset{i=1}{\overset{I}{.Math.}}{\alpha }_i\left(t\right){\psi }_i\left(z\right)\end{array}& \left(3\right)\end{array}$ where in the formula (3), z denotes a depth; c.sub.0(z) denotes an average sound velocity at the depth z, and Δc(z,t) denotes a disturbance quantity; i denotes an order of the EOF, 1≤i≤I, I denotes the total number of orders of the EOF, α.sub.i(t) denotes a coefficient of the i.sup.th-order EOF, and ψ.sub.i(z) denotes the i.sup.th-order EOF; (4) establishing an optimization model in combination with the EOF, a sound ray tracing algorithm, a surface sound velocity and multi-beam data: $\begin{array}{cc}\begin{array}{c}\hat{\alpha }\left(t\right)=\arg \underset{\alpha }{\min }\left({.Math.\mathrm{To}{a}_{meas}\left(t\right)-To{a}_{est}\left(t\right).Math.}_2+K{.Math.c_{0meas}\left(t\right)-c_{0est}\left(t\right).Math.}_2\right)\\ \mathrm{To}{a}_{est}\left(t\right)=\mathscr{H}\left[\alpha ,\mathrm{Do}{a}_{meas}\left(t\right)\right]\end{array}& \left(4\right)\end{array}$ where in the formula (4), K denotes a weight coefficient, taking a value in [0,1]; Toa.sub.meas(t) denotes actual measurement data of a time of arrival of each beam at a time t, Toa.sub.est(t), denotes an estimated value of a time of arrival of each beam at the time 1, c.sub.0meas(t) denotes measurement data of the sound velocity at a position of the receiving transducer array of the multi-beam echo sounder at the time t, and c.sub.0est(t) denotes an estimated value of the sound velocity at the position of the receiving transducer array of the multi-beam echo sounder at the time t; H denotes a measurement operator, which performs a calculation according to the sound ray tracing algorithm, and Doa.sub.meas(t) denotes a measurement value of a direction of arrival of each beam at the time t; (5) obtaining an inversion value of an EOF coefficient in an optimization algorithm, according to the optimization model established in step (4), the time of arrival of the echo and the angle of arrival of the echo received in step (2), and a measurement value of the sound velocity at the position of the receiving transducer array; and obtaining an estimation result of the sound velocity profile by using the formula (3); and (6) calculating a water temperature profile of a measurement area according to an estimation value of the sound velocity profile in step (5).

6. The water temperature profile inversion method based on the inverted multi-beam echo sounder according to claim 5, wherein in step (5), the optimization algorithm comprises a simulated annealing algorithm, a hill climbing algorithm, a particle swarm optimization algorithm, a genetic algorithm, and an ant colony algorithm.

7. The water temperature profile inversion method based on the inverted multi-beam echo sounder according to claim 5, wherein in step (6), the water temperature profile of the measurement area is calculated by using a Mackenzie formula for a sound velocity, a Chen-Millero-Li formula for the sound velocity, a Dell Grosso formula for the sound velocity, a Leroy formula for the sound velocity, a Wilson accurate formula for the sound velocity, a simplified Wilson formula for the sound velocity, or a simplified EM Hierarchical formula for the sound velocity.

8. A measurement system of an inverted multi-beam echo sounder, wherein the water temperature profile inversion method according to claim 5 is applied to the measurement system.

9. The measurement system of an inverted multi-beam echo sounder according to claim 8, wherein in step (5), the optimization algorithm comprises a simulated annealing algorithm, a hill climbing algorithm, a particle swarm optimization algorithm, a genetic algorithm, and an ant colony algorithm.

10. The measurement system of an inverted multi-beam echo sounder according to claim 8, wherein in step (6), the water temperature profile of the measurement area is calculated by using a Mackenzie formula for a sound velocity, a Chen-Millero-Li formula for the sound velocity, a Dell Grosso formula for the sound velocity, a Leroy formula for the sound velocity, a Wilson accurate formula for the sound velocity, a simplified Wilson formula for the sound velocity, or a simplified EM Hierarchical formula for the sound velocity.

Description

BRIEF DESCRIPTION OF DRAWINGS

(1) FIG. 1 is a structural block diagram of a multi-beam echo sounder;

(2) FIG. 2 is a schematic flowchart of an inversion method according to the present disclosure;

(3) FIG. 3 is a schematic diagram of a sound ray tracing algorithm for constant sound velocity;

(4) FIG. 4 is a comparison diagram of estimation results and true values of a sound velocity profile obtained with an inversion method according to the present disclosure; and

(5) FIG. 5 is a comparison diagram of estimation results and true values of a water temperature profile obtained with an inversion method according to the present disclosure.

DESCRIPTION OF EMBODIMENTS

(6) The present disclosure is specifically described in further detail as follows with reference to the accompanying drawings and embodiments.

(7) As shown in FIG. 1, a multi-beam echo sounder includes the followings.

(8) A transmitting circuit is configured to generate a transmission electrical signal to be transmitted.

(9) A transmitting transducer array is configured to convert the generated transmission electrical signal into an acoustic signal for transmission.

(10) A receiving transducer array is configured to receive an echo reflected and scattered by water and a water surface, and convert an echo acoustic signal into an electrical signal.

(11) A receiving circuit is configured to pre-process such as amplify and filter the electrical signal obtained by a conversion of the receiving transducer array.

(12) An attitude sensor is configured to determine an attitude of the multi-beam echo sounder.

(13) A pressure sensor is configured to measure a pressure on the multi-beam echo sounder.

(14) A surface sound velocity probe is configured to determine a sound velocity at a mounting position of the multi-beam echo sounder.

(15) A system control and signal processing system is configured to implement high-speed collection and real-time processing of signals, and transmit collected data and processing results to a display control and post-processing platform via a network, wherein the collected data includes attitude sensor data, pressure sensor data and surface sound velocity data; and the processing results include an angle of arrival of the echo and a time of arrival of the echo.

(16) A display control and post-processing platform is configured to issue an instruction to the system control and signal processing system, and receive the processing results and collected data from the system control and signal processing system; invert a sound velocity profile, by using an optimization model established on the processing results and collected data, and further calculate a water temperature profile of a measurement area by using an estimated value of the sound velocity profile.

(17) A flow of the sound velocity profile inversion method based on an inverted multi-beam echo sounder according to the present disclosure is as shown in FIG. 2. The sound velocity profile inversion method includes following steps.

(18) (I) The multi-beam echo sounder is mounted on an underwater submerged buoy or fixed to a water bottom in an inverted manner; and is configured to transmit a beam to a water surface by a transmitting transducer array, and receive an echo signal by a receiving transducer array.

(19) (II) The multi-beam echo sounder obtains an angle of arrival of an echo and a time of arrival of the echo according to the received echo signal.

(20) (III) A description method of a dimension reduction primary function of the sound velocity profile is obtained, by solving an EOF according to sound velocity profile prior information.

(21) The EOF may be calculated according to a sample set of known sound velocity profile in a measured sea area. It is assumed that there is sound velocity profile data actually measured M times, for example:

(22) $\begin{array}{cc}C=\left[\begin{array}{cccc}{c}_1\left(z_1\right)& c_2\left(z_1\right)& .Math.& c_M\left(z_1\right)\\ c_1\left(z_2\right)& c_2\left(z_2\right)& .Math.& c_M\left(z_2\right)\\ .Math.& .Math.& .Math.& .Math.\\ c_1\left(z_{N_z}\right)& c_2\left(z_{N_z}\right)& .Math.& c_M\left(z_{N_z}\right)\end{array}\right].& \left(5\right)\end{array}$

(23) In the formula (5), z.sub.i is a depth value. Then, a determined value of the sound velocity is:

(24) $\begin{array}{cc}{c}_0\left(z\right)=\frac{1}{M}\underset{m=1}{\overset{M}{.Math.}}{c}_m\left(z\right),& \left(6\right)\end{array}$

(25) A covariance matrix is

(26) $\begin{array}{cc}R=\frac{1}{M-1}\underset{m=1}{\overset{M}{.Math.}}{\left[c_m-c_0\right]\left[c_m-c_0\right]}^T.& \left(7\right)\end{array}$

(27) Finally, eigenvalue decomposition is performed on the covariance matrix to obtain an expansion coefficient {λ.sub.i}.sub.i=1.sup.I and a corresponding eigenvector {u.sub.i}.sub.i=1.sup.I. This set of orthogonal vectors corresponds to the EOF, namely {ψ.sub.i}.sub.i=1.sup.I.

(28) (IV) An optimization model is established in combination with the EOF, a sound ray tracing algorithm, a surface sound velocity and multi-beam data.

(29) The measurement operator H in the optimization model may be calculated according to the sound ray tracing algorithm.

(30) Because of the heterogeneity and anisotropy of seawater, sound waves may be refracted when they travel through different sound velocity layers, causing rays to bend, as shown in FIG. 3. In the case of inhomogeneous media, it is assumed that the sound velocity changes only in a vertical direction and is uniform in a horizontal direction. The sound ray propagation follows Snell's theorem:

(31) $\begin{array}{cc}\frac{\sin \theta }{c_0}=\frac{\sin \theta \left(z\right)}{c\left(z\right)}\equiv a.& \left(8\right)\end{array}$

(32) In the formula (8), c.sub.0 denotes a sound velocity at a position of a transducer array, referred to as a surface sound velocity; θ denotes an angle between a beam echo and the vertical direction when the beam echo reaches the receiving array; c(z) denotes a sound velocity at the depth z, and θ(z) denotes an angle between a tangent to the sound ray and the vertical direction at the depth z. It is assumed that a horizontal distance dr, a vertical distance dz and a path length ds of the sound ray propagation within a short period of time dt satisfy:

(33) $$\begin{gathered} \begin{array}{cc}dr=\tan \theta \left(z\right)\mathrm{dz}, \\ \mathrm{ds}=\frac{\mathrm{dz}}{\cos \theta \left(z\right)}, \\ \mathrm{dt}=\frac{\mathrm{ds}}{c\left(z\right)}=\frac{\mathrm{dz}}{c\left(z\right)\cos \theta \left(z\right)},& \left(9\right)\end{array} \end{gathered}$$

(34) The propagation time of the sound ray emitted by the transducer array and reaching at the depth z is:

(35) 0$\begin{array}{cc}t={\int }_{z_0}^z\mathrm{dt}={\int }_{z_0}^z\frac{1}{c\left(z\right)\sqrt{1-a^2{c}^2\left(z\right)}}\mathrm{dz},& \left(10\right)\end{array}$

(36) and the horizontal distance of the sound ray propagation is:

(37) $\begin{array}{cc}r={\int }_{z_0}^z\mathrm{dr}={\int }_{z_0}^z\frac{\mathrm{ac}\left(z\right)}{\sqrt{1-a^2{c}^2\left(z\right)}}\mathrm{dz}.& \left(11\right)\end{array}$

(38) In the formula (11), z.sub.0 is a depth of the transducer array.

(39) Because of the complexity and variability of ocean environment, it is difficult to represent the sound velocity profile c(z) as a simple continuous function. The whole water body is divided into several layers according to actually measured depths of measurement points of the sound velocity profile, as shown in FIG. 3. A sound ray tracing calculation is performed with an assumption of in-layer equal gradient.

(40) Firstly, the water body is divided into N layers in the depth direction. It is assumed that a depth and a sound velocity at a hierarchical position are respectively z=[z.sub.0, z.sub.1, . . . , z.sub.N].sup.T and c=[c.sub.0, c.sub.1, . . . , c.sub.N].sup.T, and it is assumed that the in-layer sound velocity has constant gradient, namely
c(z)=c.sub.i+g.sub.i(z−z.sub.i).  (12)

(41) In the formula (12), i is the number of layer (i=0, 2, . . . , N−1), c.sub.i is a sound velocity value at a boundary of the i.sup.th layer, z.sub.i is a depth of the i.sup.th layer, g.sub.i is a gradient of the sound velocity in the i.sup.th layer, and

(42) $g_i=\frac{c_{i+1}-c_i}{z_{i+1}-z_i}.$
It can be known according to the formula (10) that on the assumption of the in-layer constant gradient, the round trip propagation time of the sound ray in the i.sup.th layer may be expressed as:

(43) $\begin{array}{cc}{t}_i=\frac{2}{g_i}\ln \left[\frac{c_{i+1}\left(1+\sqrt{1-a^2{c}_i^2}\right)}{c_i\left(1+\sqrt{1-a^2{c}_{i+1}^2}\right)}\right].& \left(13\right)\end{array}$

(44) It can be known according to the formula (11) that a horizontal distance of sound ray propagation in the i.sup.th layer may be expressed as:

(45) $\begin{array}{cc}{r}_i=-\frac{1}{{\mathrm{ag}}_i}\left(\sqrt{1-a^2{c}_{i+1}^2}-\sqrt{1-a^2{c}_i^2}\right),& \left(14\right)\end{array}$

(46) In the whole depth range, relations between the total propagation time of sound rays and the depth and between the horizontal distance of sound rays and the depth may be expressed as:

(47) $\begin{array}{cc}t\left(z,\theta \right)=\{\begin{array}{c}\mathrm{for}{z}_0\le z\le z_1:\\ \frac{2}{g_0}\left(\ln \left[\frac{c\left(z\right)}{1+\sqrt{1-a^2{c\left(z\right)}^2}}\right]-\ln \left[\frac{c_0}{1+\sqrt{1-{\left(\sin \theta \right)}^2}}\right]\right)\\ \mathrm{for}{z}_i\le z<z_{i+1},i=1,2,.Math.,N-1:\\ \underset{n=0}{\overset{i-1}{.Math.}}{t}_n+\frac{2}{g_i}\left(\ln \left[\frac{c\left(z\right)}{1+\sqrt{1-a^2{c\left(z\right)}^2}}\right]-\ln \left[\frac{c_i}{1+\sqrt{1-(a^2{c}_i^2}}\right]\right)\end{array}& \left(15\right)\end{array}$$\begin{array}{cc}r\left(z,\theta \right)=\{\begin{array}{c}\mathrm{for}{z}_0\le z\le z_1:\\ \frac{-1}{a{g}_0}\sqrt{1-a^2{c\left(z\right)}^2}-\sqrt{1-{\left(\sin \theta \right)}^2)},\\ \mathrm{for}{z}_i\le z<z_{i+1},i=1,2,.Math.,N-1:\\ \underset{n=0}{\overset{i-1}{.Math.}}{r}_n+\frac{-1}{{\mathrm{ag}}_i}\left(\sqrt{1-a^2{c\left(z\right)}^2}-\sqrt{1-(a^2{c}_i^2}\right).\end{array}& \left(16\right)\end{array}$

(48) The measurement operator H in the formula (2) for the optimization model is obtained on the basis of the formula (15).

(49) (V) According to the established optimization model, data of an angle of arrival of the echo and a time of arrival of the echo as well as data of collected attitude, pressure and surface sound velocity are transmitted to the display control and post-processing platform through the network, and inversion of coefficients of the sound velocity profile primary function EOF is obtained in an optimization algorithm.

(50) In general, the optimization algorithm such as a simulated annealing algorithm, a hill climbing algorithm, a particle swarm optimization algorithm, a genetic algorithm or an ant colony algorithm can be used to estimate to-be-estimated parameters. FIG. 4 shows a simulation example of genetic algorithm in Matlab software. The number of orders of EOFs is 3. In the genetic algorithm, a lower limit of parameter estimation is [−100 −100 −100], an upper limit of parameter estimation is [100 100 100] the population number is 200, and the number of iterations is 100. As can be seen from FIG. 4, the sound velocity profile estimated by the method according to the present disclosure is very close to the real profile, and a mean square error of the whole profile is 0.9726 m/s, indicating that the method according to the present disclosure has a good estimation effect.

(51) (VI) A water temperature profile of a measurement area is calculated according to an estimated value of the sound velocity profile in step (5).

(52) The water temperature profile of the measurement area can be further calculated, according to the estimation result of the sound velocity profile in step (5) by using an empirical formula for sound velocity applicable to the measurement area. At present, generally accepted empirical formulas for sound velocity mainly includes a Mackenzie formula for the sound velocity, a Chen-Millero-Li formula for the sound velocity, a Dell Grosso formula for the sound velocity, a Leroy formula for the sound velocity, a Wilson accurate formula for the sound velocity, a simplified Wilson formula for the sound velocity, and an simplified EM Hierarchical formula for the sound velocity. The simplified EM Hierarchical formula for the sound velocity has a simple structure and high calculation accuracy, and can be applied to measurement areas with a wide range of seawater salinity, depth or temperature change. Thus, the simplified EM Hierarchical formula for the sound velocity is a preferred implementation method of the present disclosure.

(53) In the following, the water temperature profile is calculated by taking the simplified EM Hierarchical formula for the sound velocity as an example.

(54) The simplified EM Hierarchical formula for the sound velocity may be expressed as the following formulas.

(55) A calculation formula of a surface sound velocity is
c(T,0,S)=1449.05+T[4.57−T(0.0521−0.00023T)]+[1.33−T(0.0126−0.00009T)](S−35).  (17)

(56) When the depth of fresh water reaches 200 m and the depth of seawater reaches 1000 m, the sound velocity calculation formula is
c(T,z,S)=c(T,0,S)+16.5z.  (18)

(57) When the depth of fresh water reaches 2000 m and the depth of seawater reaches 11000 m, the sound velocity calculation formula is
c(T,z,S)=c(T,0,S)+z[16.3+z(0.22−0.003z√{square root over (T+2)})].  (19)

(58) When the depth is greater than 5000 m, a latitude correction should be taken into account in the calculation formula of the sound velocity:
c(T,z,S)=c(T,0,S)+z′[16.3+z′(0.22−0.003z′√{square root over (T+2)})].  (20)
z′=z(1−0.0026 cos(2φ)).  (21)

(59) In the formulas (17)-(21), T denotes a seawater temperature; z denotes a seawater depth; S denotes seawater salinity; p denotes a latitude of the calculated position; c(T,z,S) denotes a sound velocity when the seawater temperature is T, the seawater depth is z, and the seawater salinity is S; and the depth z is in units of km.

(60) A depth of a placement position of the measuring system can be determined according to data of the pressure sensor. According to the applicable range of the seawater depth of formulas (17)-(20), an appropriate formula can be selected to estimate the water temperature profile. In the simulation example, the system is placed at a depth of about 200 m, and thus the formula (17) is selected to estimate the water temperature profile of the measurement area, as shown in FIG. 5. As can be seen from FIG. 4, the water temperature profile estimated by the method according to the present disclosure is very close to the real profile, which indicates that the estimation effect of the method according to the present disclosure is good.

(61) The above embodiments are used to explain the present disclosure and but not intended to limit the present disclosure. Within the protection scope of the spirit and claims of the present disclosure, any modifications and changes made to the present disclosure all fall into the protection scope of the present disclosure.