Automatic Inversion Method and System for Equivalent Sound Velocity Based on Multi-Beam Overlap Region Constraint

Abstract

The present invention belongs to the technical field of marine surveying and mapping, and discloses an automatic inversion method and system for equivalent sound velocity based on multi-beam overlap area constraint. The method uses the surface sound velocity of multi-beam data to track the beam point positions of the multi-beam original data, and obtains the overlapping area of the multi-beam strip; grids the overlapping area of the multi-beam strip, and uses the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtains the optimal equivalent sound velocity of each survey line; uses the corresponding equivalent sound velocity obtained by inversion for each survey line to re-track the beam point positions, weakens the sound velocity deviation existing in the beam point position tracking using the surface sound velocity, and weakens the curved terrain between the strips in the overlapping area, and obtains a continuous seabed terrain for the splicing of the adjacent strip overlapping areas. The present invention can ensure the accuracy and reliability of seabed terrain measurement and improve the speed of the field work for seabed terrain measurement.

Description

Technical Field

[0001] The present invention belongs to the technical field of marine surveying and mapping, and particularly relates to an automatic inversion method and system for equivalent sound speed based on multi-beam overlap area constraint. Background Art

[0002] The development and exploration of the ocean are increasing day by day. However, underwater topographic survey is a fundamental task in marine surveying and mapping. By using sound waves and measuring information such as sound speed, propagation time, attitude, water level, etc., it can obtain information such as water depth and position, which can provide services for the development of various marine projects and also provide basic data for the application of nautical charts and marine geographic information. It is a crucial fundamental task for people to develop and explore the ocean. Compared with the single-beam operation mode, the multi-beam sounding technology has been upgraded from linear measurement to planar coverage. It can obtain a large number of sounding point data in the vertical direction of the track line and efficiently measure the size, shape and undulation changes of underwater targets. This technology has the advantages of high efficiency, high precision, high resolution and full coverage, and is widely used in fields such as seabed resource exploration, waterway mapping, island and reef surveying, port monitoring and underwater archaeology, becoming the core means to obtain seabed topographic and geomorphic information from shallow sea to deep sea.

[0003] The multi-beam sounding system measures the water depth by emitting sound waves and receiving the signals reflected from the seabed. However, the propagation speed of sound waves in water changes with the water depth, which is mainly due to the uneven distribution of temperature and salinity in the water body. The change of sound speed in seawater will affect the propagation speed of sound waves, and thus affect the accuracy of the sounding results. Therefore, accurate sound speed profile information is crucial for converting the propagation time of sound waves into water depth. In some cases, such as the environment where fresh water and seawater are mixed in a certain area, the change of sound speed may be very drastic, resulting in significant deviations in the results of traditional methods that use a small number of measured sound speed profiles to replace the entire area, that is, the measured water depth is underestimated or overestimated in the outer part of the covered area.

[0004] Currently, there are towing systems for obtaining sound speed profiles, such as mobile profilers or by measuring seawater conductivity, temperature, pressure, etc. to obtain the measurement data of sound speed profiles. However, due to the risk of entanglement or stranding of the towing instrument in each deployment, they have not been widely used. Instead, sound speed profilers and CTDs are used for measurement. To perform such measurements, the ship needs to remain stationary. Therefore, obtaining sound speed profile measurement data becomes a time-consuming process, with poor practicability and low efficiency in actual engineering applications. Summary of the Invention

[0005] To overcome the problems existing in the related technologies, the disclosed embodiments of the present invention provide an automatic inversion method and system for equivalent sound speed based on multi-beam overlap area constraint. The present invention solves the problems of topographic distortion caused by inaccurate sound speed measurement and low efficiency of multi-beam field measurement.

[0006] The technical solution is as follows: An automatic inversion method for equivalent sound velocity based on multi-beam overlap area constraint, comprising the following steps:

[0007] S1. Use the surface sound velocity of multi-beam data to track the beam point positions of the multi-beam original data, and obtain the overlapping area of the multi-beam strip;

[0008] S2. Grid the overlapping area of the multi-beam strip, and use the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtain the optimal equivalent sound velocity of each survey line;

[0009] S3. For each survey line, re-track the beam point positions using the corresponding equivalent sound velocity obtained by inversion, weaken the sound velocity deviation existing in the beam point position tracking using the surface sound velocity, and weaken the topographic distortion between the strips in the overlapping area, so as to obtain the continuous real seabed topography of the adjacent strip overlapping area, thereby getting rid of the dependence on the measured sound velocity information for beam point position tracking.

[0010] In step S1, obtaining the overlapping area of the multi-beam strip includes:

[0011] S101. Parse the obtained multi-beam original data, and use the surface sound velocity information measured by the surface sound velocity instrument in the data or assume a sound velocity profile to perform ray tracing on the data, and obtain the initial water depth and position of all measurement points;

[0012] S102. Search for the overlapping area between each survey line along the direction perpendicular to the multi-beam course.

[0013] In step S101, the ray tracing includes ray tracing based on constant sound velocity within a layer and ray tracing based on constant gradient within a layer.

[0014] Furthermore, the ray tracing based on constant sound velocity within a layer includes:

[0015] The sound velocity propagates at a constant speed within the layer, and the depths of the upper and lower interfaces of the layer are respectively and and , the layer thickness is , , are respectively the beam incident angle and the sound velocity of the layer, according to Snell's law, , then the horizontal displacement and the propagation time of the beam within the layer are respectively:

[0016] (1)

[0017] The horizontal propagation distance of the beam through the entire water column and the propagation time are as follows:

[0018] (2)

[0019] In the formula, is the number of water layers;

[0020] If the incident angle is 0, that is, in the case of single-beam sounding, then there is:

[0021] ;

[0022] .

[0023] Furthermore, the ray tracing based on a constant gradient within the layer includes:

[0024] The sound speed changes with a constant gradient within the layer , and the actual propagation trajectory of the beam within the layer is a continuous arc segment with a radius of curvature of , and the expression is:

[0025] (3)

[0026] (4)

[0027] In the formula, is the sound speed of the next water layer;

[0028] The horizontal displacement experienced by the beam in this layer and the arc length

[0029] (5)

[0030] In the formula, is the incident angle of the next water layer;

[0031] According to Snell's law: , then the horizontal displacement of the ray within the th layer and the time

[0032] (6)

[0033] In the formula, is the harmonic mean sound speed within the layer ;

[0034] (7)

[0035] If the incident angle is 0, the propagation trajectory of the sound wave in each layer is a straight line, vertically downward, and the propagation time in each layer is:

[0036] (8).

[0037] In step S2, the genetic algorithm is used to minimize the water depth deviation between overlapping grids, including:

[0038] Input the overlapping area data of the gridded multi-beam strip into the genetic algorithm and use the genetic algorithm for optimization; the genetic algorithm randomly generates a population within a certain range, and then performs crossover and mutation operations to calculate the fitness function. The fitness function is:

[0039] ;

[0040] In the formula, is the equivalent sound speed to be obtained, representing the sound speed gradient corresponding to stripes in the area; is the water depth of each grid node, is the inverse distance weighted average water depth of all measurement points in the grid, is the total number of grids in the selected area, is in the area the total number of measurement points on the

[0041] Judge whether the average water depth deviation of the overlapping area calculated by the current fitness function is less than 1% of the average water depth. If it reaches, exit the loop and return the optimal equivalent sound speed of each survey line. If not, perform iterative calculation again.

[0042] Furthermore, the specific method for the genetic algorithm to randomly generate a population within a certain range, and then perform crossover and mutation operations to calculate the fitness function is as follows:

[0043] S201, according to the range of the initial value , randomly generate the initial value with the data volume , and call the generated initial value the initial population;

[0044] S202, calculate the function value of each seed value of the generated initial population respectively, screen the initial population according to the size of each seed value, and screen out the seeds with small function values and keep them; then perform mutation and crossover on the seeds with deviations exceeding the threshold according to the mutation probability and crossover probability to form a new population;

[0045] S203, recalculate the respective Values, genetically retain excellent seed values, replace inferior individuals, and then generate new population individuals;

[0046] S204. Perform mutation and crossover operations on the newly generated population individuals to produce a new population after mutation and crossover, and then calculate the ;

[0047] S205. Repeat steps S203 and S204 until the optimal value is found or the maximum number of iterations is reached and then exit, and output the best equivalent sound speed value.

[0048] In step S2, invert to obtain the optimal equivalent sound speed for each survey line, including:

[0049] The process of equivalent sound speed migration is:

[0050] (9)

[0051] In the formula, The sound speed gradient is a constant, is the depth, is the sound speed, is the bottom sound speed, is the surface sound speed, is the bottom water depth value, is the incident surface sound speed, is a constant;

[0052] The expression for the propagation of sound waves in water that conforms to Snell's law is:

[0053] (10)

[0054] In the formula, 、 are the upper and lower beam incident angles of a single layer respectively, is a constant obtained based on Snell;

[0055] According to the principle of the equivalent sound speed profile method, the propagation time is:

[0056] (11)

[0057] Combining formula (10) and formula (11), derive the propagation time as:

[0058] (12)

[0059] Since the propagation time is accurately measured and is a known quantity, therefore, convert formula (12) to obtain as:

[0060] (13)

[0061] In the formula, are all one-way propagation times;

[0062] According to the theory of ray acoustics, the propagation trajectory of a sound ray in a constant gradient layer is a circular arc, and the radius of the circular arc is:

[0063] (14)

[0064] Through derivation, the displacement of the sound ray in the depth direction and the displacement in the horizontal direction are respectively:

[0065] (15).

[0066] Another object of the present invention is to provide an automatic inversion system for equivalent sound speed based on multi-beam overlap area constraint. This system implements the automatic inversion method for equivalent sound speed based on multi-beam overlap area constraint. This system includes:

[0067] An overlap area acquisition module for multi-beam strips, which is used to track the beam point positions of the multi-beam original data by using the surface sound speed of the multi-beam data, and obtain the overlap area of the multi-beam strips;

[0068] An equivalent sound speed inversion module, which is used to grid the overlap area of the multi-beam strips, and use the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtain the optimal equivalent sound speed for each survey line;

[0069] A re-beam point position tracking module, which is used to re-track the beam point positions of each survey line by using the corresponding equivalent sound speed obtained by inversion, weaken the sound speed deviation existing in the beam point position tracking by using the surface sound speed, and weaken the topographic distortion between the strips in the overlap area, so as to obtain the continuous real seabed topography of the adjacent strip overlap area, thereby getting rid of the dependence of the beam point position tracking on the measured sound speed information.

[0070] Furthermore, the automatic inversion system for equivalent sound speed based on multi-beam overlap area constraint is carried on a computer device. This computer device includes: at least one processor, a memory, and a computer program stored in the memory and executable on the at least one processor. When the processor executes the computer program, it realizes the functions in the above-mentioned automatic inversion system for equivalent sound speed based on multi-beam overlap area constraint.

[0071] Combining all the above technical solutions, the beneficial effects of the present invention are:

[0072] The present invention uses the surface sound velocity of multi-beam data to perform beam point tracking on the multi-beam raw data to obtain the overlapping area of ​​the strips; the multi-beam strips are gridded, and the water depth deviation between the overlapping grids is minimized by using a genetic algorithm, and the optimal equivalent sound velocity of each survey line is inverted; the beam point tracking is re-performed on each survey line using the corresponding equivalent sound velocity obtained by inversion, thereby reducing the sound velocity deviation existing in the beam point tracking using the surface sound velocity, so that the crying and laughing face-shaped curved terrain between the strips in the overlapping area is weakened, and a more continuous seabed terrain is obtained in the overlapping area of ​​adjacent strips, thereby reducing the dependence of the beam point tracking on the measured sound velocity profile.

[0073] In actual measurements of multi-beam systems, adjacent strips will have a certain area of ​​overlap. The present invention will make full use of the overlapping areas between multi-beam measurement strips to re-track the sound line, thereby recalculating the measured depth. The present invention inverts the equivalent sound speed in seawater by minimizing the difference in water depth in the overlapping areas of different tracks as a constraint, thereby achieving re-tracking of the sound line. The advantage of this method is that it can use the redundant information in the multi-beam data to achieve accurate multi-beam measurement point position reduction without direct sound speed measurement. Through this method, the accuracy and reliability of seabed topography measurement can be guaranteed, and the speed of seabed topography field work can be improved.

[0074] The present invention can provide accurate sound velocity information for seabed detection by using the equivalent sound velocity automatic inversion method constrained by the multi-beam overlapping area, improve the efficiency of field acquisition, and provide more accurate basic data for scientific research and marine engineering. The present invention innovatively uses the multi-beam overlapping area as a constraint and combines the equivalent sound velocity principle to achieve the reduction of the sound velocity error of multi-beam bathymetry, get rid of the time-consuming and labor-intensive process of measuring dense sound velocity profiles in large-scale sea areas with drastic sound velocity changes, and objectively restore the real seabed topography. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] The accompanying drawings herein are incorporated in and constitute a part of the specification, illustrate embodiments consistent with the present disclosure, and together with the description, serve to explain the principles of the present disclosure;

[0076] Figure 1 It is a flow chart of an equivalent sound velocity automatic inversion method based on multi-beam overlap area constraint provided by an embodiment of the present invention;

[0077] Figure 2 It is a principle diagram of an equivalent sound velocity automatic inversion method based on multi-beam overlapping area constraint provided by an embodiment of the present invention;

[0078] Figure 3 It is a schematic diagram of the gridding of the survey lines in the present invention;

[0079] Figure 4Schematic diagram of ray tracing based on constant sound speed within a layer in the present invention;

[0080] Figure 5 Schematic diagram of ray tracing based on constant gradient within a layer in the present invention;

[0081] Figure 6 Schematic diagram of the area difference of the sound speed profile in the present invention;

[0082] Figure 7 Schematic diagram of the equivalent sound speed profile in the present invention;

[0083] Figure 8 Schematic diagram of the automatic inversion system of equivalent sound speed based on the constraint of the multi-beam overlapping area provided by the embodiment of the present invention;

[0084] Figure 9 Experimental diagram of the accuracy of the inversion algorithm provided by the embodiment of the present invention;

[0085] Figure 10 Histogram of water depth deviation in the overlapping area provided by the embodiment of the present invention;

[0086] In the figure: 1. Overlapping area acquisition module of multi-beam strips; 2. Equivalent sound speed inversion module; 3. Point position tracking module of heavy beams. Detailed implementation manners

[0087] To make the above objects, features, and advantages of the present invention more obvious and understandable, the following will describe the detailed implementation manners of the present invention in conjunction with the accompanying drawings. Many specific details are set forth in the following description to fully understand the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without departing from the connotation of the present invention. Therefore, the present invention is not limited by the specific implementations disclosed below.

[0088] The innovation of the present invention lies in: the present invention utilizes the overlapping area between multi-beam sounding strips, takes minimizing the water depth difference in the overlapping area of adjacent strips as a constraint, and inversely calculates the equivalent sound speed through a genetic algorithm to achieve point position tracking of beam positions without an actually measured sound speed profile, and its accuracy can reach the same level as that of ray tracing using a sound speed profile, thereby getting rid of the need for sound speed acquisition in ocean terrain survey.

[0089] Embodiment 1. In the underwater measurement part of the multi-beam sounding system of the present invention, a multi-beam sonar transducer installed on a survey ship is adopted. As the detection carrier moves, three-dimensional data of the detection target and the seabed topography are obtained. The positioning and navigation part combines GNSS and attitude measurement devices (including a compass) to obtain the accurate position, attitude and heading information of the survey ship. The surface sound velocity data is measured by a surface sound velocity instrument and is used for beam pointing correction of the multi-beam transducer. The water body sound velocity profile information is inversely obtained by the present invention and no prior sound velocity information is required.

[0090] As Figure 1 shown, the automatic equivalent sound velocity inversion method based on multi-beam overlap area constraint provided by the embodiment of the present invention includes the following steps:

[0091] S1. Use the surface sound velocity of the multi-beam data to perform beam point position tracking on the multi-beam original data, and obtain the overlapping area of the multi-beam strip;

[0092] S2. Grid the overlapping area of the multi-beam strip, and use the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtain the optimal equivalent sound velocity of each survey line;

[0093] S3. Use the corresponding equivalent sound velocity inversely obtained for each survey line to re-perform beam point position tracking, weaken the sound velocity deviation existing in beam point position tracking using the surface sound velocity, and weaken the topographic distortion between strips in the overlapping area, so as to obtain the continuous real seabed topography of the adjacent strip overlapping area splicing, thereby getting rid of the dependence of beam point position tracking on the measured sound velocity information.

[0094] Exemplarily, in step S1, using the surface sound velocity of the multi-beam data to perform beam point position tracking on the multi-beam original data and obtain the overlapping area of the multi-beam strip includes:

[0095] S101. Parse the obtained multi-beam original data, use the surface sound velocity information measured by the surface sound velocity instrument in the data or assume a sound velocity profile to perform sound ray tracking on the data, and obtain the initial water depth and position of all measurement points;

[0096] When obtaining the multi-beam original data, it is necessary to ensure that there is a certain overlapping area between adjacent survey lines.

[0097] S102. Search for the overlapping area between each survey line along the direction perpendicular to the multi-beam heading.

[0098] Obtain the bounding box of the survey line, traverse all pairs of survey lines, and check one by one whether the coverage areas of the two survey lines overlap. Use a geometric algorithm to judge whether the coverage areas of the two survey lines overlap. If the coverage areas of the two survey lines overlap, calculate the boundary of the overlapping part and save it.

[0099] Exemplarily, in step S2, the meshing of the overlapping region of the multi-beam strips includes:

[0100] Due to the influence of factors such as ship speed, attitude, and heading, the overlapping regions of the two closest Ping data in adjacent survey lines cannot completely coincide. Therefore, in order to evaluate the water depth consistency of the overlapping part, the overlapping region is meshed, and the overlapping region is divided into grids: according to the defined grid size, the bounding box of the overlapping region is divided into multiple grid cells. The size of each grid cell is the pre-set grid resolution.

[0101] For example, if the bounding box of the overlapping region is (x_min, y_min, x_max, y_max) and the size of the grid cell is (dx, dy), then the number of rows and columns of the grid can be calculated by the following formula:

[0102] Number of rows = (y_max - y_min) / dy

[0103] Number of columns = (x_max - x_min) / dx

[0104] By dividing the bounding box of the overlapping region into several small grids, representing them with the row and column indices of the grid, and calculating the points in the grid; finally, the meshing result is saved, and indexing is performed through the rows and columns of the grid, avoiding the process of matching points with the same name.

[0105] This can avoid the complex process of matching points with the same name and can speed up the calculation. The grid resolution also needs to consider the size of the measurement point density, and the size of the grid can be set according to the actual water depth value to ensure that there are at least two measurement points from different survey lines within one grid;

[0106] Exemplarily, in step S2, using the genetic algorithm to minimize the water depth deviation between overlapping grids and inversely obtaining the optimal equivalent sound speed of each survey line includes:

[0107] The data of the overlapping region of the meshed multi-beam strips is input into the genetic algorithm, and the genetic algorithm is used for optimization. The genetic algorithm randomly generates a certain number of populations within the range, and then performs crossover and mutation operations, and then calculates the fitness function. The fitness function is:

[0108] ;

[0109] In the formula, is the equivalent sound speed to be obtained, representing the sound speed gradient corresponding to stripes in the region; is the water depth of each grid node, is the weighted average water depth of all measurement points in the grid by inverse distance, is the total number of grids in the selected area, in the area the total number of measurement points on the

[0110] Judge the Whether the calculated average water depth deviation of the overlapping area is less than 1% of the average water depth. If it reaches, exit the loop and return the optimal equivalent sound speed of each survey line. If it doesn't reach, perform iterative calculation again.

[0111] Exemplarily, the genetic algorithm randomly generates a certain number of populations within a range, then performs crossover and mutation operations, and then calculates the fitness function. The specific method is as follows:

[0112] S201, according to the range of the initial value , randomly generate the initial value of the data volume , and call the generated initial value the initial population;

[0113] S202, calculate the function values of each seed value of the generated initial population respectively, screen the initial population according to the size of each seed value, and screen out seeds with small function values, sort the function values of the seed values, and select the first 80% for retention; sort the function values of the seed values, and the last 20% are considered seeds with large deviations. Then, for the seeds with deviations exceeding the threshold, perform mutation and crossover according to the mutation probability and crossover probability to form a new population; select the parent generation: select individuals with higher fitness as the parent generation according to the fitness function. The crossover operation generates new offspring by exchanging part of the genes of the parent generation. The goal of crossover is to combine the excellent characteristics of the parent generation to generate more potential offspring. Mutation operation: perform mutation operation on the offspring individuals, randomly changing some genes of the individuals. The purpose of mutation is to increase the diversity of the population and avoid the algorithm falling into local optimum.

[0114] S203, recalculate the respective values for the new population, genetically retain the excellent seed values, replace the inferior individuals, and then generate new population individuals;

[0115] S204, perform mutation and crossover operations on the newly generated population individuals to generate a new population after mutation and crossover, and then calculate the of the new population;

[0116] S205, repeat steps S203 and S204 until the optimal value is found or the maximum number of iterations is reached and then exit, and output the best equivalent sound speed value.

[0117] Exemplarily, in step S3, assume that the true surface sound speed of the water body is , the surface sound velocity measured by the surface sound velocimeter is , when using the measured surface sound velocity to calculate the beam steering angle, according to the principle of beam steering, since is less than the true value, so when the time delay of each array element is not changed, the angle recorded by the beam will be smaller. When the attitude angle and installation deviation angle are ignored, the beam incident angle is , therefore, the water depth calculation formula at constant sound velocity is as follows:

[0118] ;

[0119] In the formula, represents the water depth value, represents the propagation vector corresponding to the beam. For any beam, when the propagation vector is constant and decreases, the water depth value will increase. In this case, the water depth value calculated by the edge beam will increase. Since the central beam is not affected by beam steering, the measured underwater topography shows a "crying face" shape. If is greater than the true value, the situation will be the opposite, and the topography will show a "smiling face" shape.

[0120] When there is an error in the sound velocity profile, such as when is greater than the actual sound velocity value, according to the following formula, will be smaller, and at this time, the propagation angle passed to the next layer will also decrease accordingly. The calculation formula is as follows:

[0121] ;

[0122] In the formula, is the sound velocity value of the next layer. Since decreases, according to the formula, the water depth value will increase accordingly. Since the central beam is almost vertically incident, the central beam has no refraction and is not affected by it, and the water depth error is relatively small. Therefore, the topography will show a "crying face" shape. If is less than the actual sound velocity value, the topography will show a "smiling face" shape.

[0123] If there is no surface sound velocity error and only a sound velocity profile error exists, the Snell constant corresponding to each beam is the correct value. However, for a certain water layer , if the sound velocity in the layer is greater than the true sound velocity, it can be known that the acoustic wave incident angle It will increase, resulting in a smaller water depth value obtained by subsequent sound ray tracking. Since the central beam is vertically incident and does not refract, it is not affected by this, and the water depth error is relatively small, which will cause the terrain to appear in a "smiling face" shape. If When the value is less than the actual sound speed value, it will cause the terrain to appear in a "crying face" shape.

[0124] Use the optimal equivalent sound speed of each sounding line obtained in Step 2 to re - perform sound speed tracking, and obtain the seabed terrain after the inverted sound speed. Thus, the automatic inversion of the equivalent sound speed is completed, getting rid of the dependence on the measured sound speed profile, and reducing the "crying and smiling face" terrain of the seabed terrain caused by inaccurate sound speed measurement.

[0125] Exemplarily, the principle of the automatic inversion method of equivalent sound speed based on the multi - beam overlap area constraint is as Figure 2 shown.

[0126] Exemplarily, the theoretical basis related to the present invention includes:

[0127] (1) The principle of multi - beam sounding ray tracking.

[0128] Ray tracking is a sound speed correction method that uses the sound speed profile to stack the positions of sound rays layer by layer, so as to calculate the coordinates of the underwater projection point of the sound ray (also known as the beam footprint) in the ship body coordinate system. Ray tracking usually divides the water layer between two adjacent sound speed sampling points among the N + 1 sampling points of the sound speed profile into one layer. Then, the entire water column through which the sound ray propagates can be regarded as composed of N water layers stacked. If the vertical displacement and horizontal displacement of the sound ray in each layer are obtained, the vertical displacement and horizontal displacement of the beam through the entire water column can be obtained by superposition. The change of sound speed within the layer is generally divided into two cases: when it is assumed that the sound speed within the layer is a constant value, the propagation trajectory of the sound ray is a straight line, and the calculation process of ray tracking is relatively simple, but the sound speed will change suddenly at the junction of adjacent layers; when it is assumed that the sound speed within the layer changes with a constant gradient, the propagation trajectory of the sound ray is an arc, which is more in line with the actual change of the sound ray underwater.

[0129] (1.1) Ray tracking based on constant sound speed within the layer.

[0130] As Figure 3 shown, the sound speed propagates at a constant speed within the th layer. The depths at the upper and lower interfaces of the layer are and respectively, the layer thickness is , , , are the beam incident angle and sound speed of the th layer respectively. According to Snell's law, , then the horizontal displacement of the beam within the layer and the propagation time They are respectively:

[0131] (1)

[0132] The horizontal propagation distance of the beam through the entire water column and the propagation time are:

[0133] (2)

[0134] In the formula, is the number of water layers;

[0135] If the incident angle is 0, that is, in the case of single-beam sounding, then there is:

[0136] ;

[0137] .

[0138] (1.2)Ray tracing based on a constant gradient within the layer.

[0139] As Figure 4 shown, assuming that the sound speed changes with a constant gradient within the layer and the actual propagation trajectory of the beam within the layer is a continuous arc segment with a radius of curvature of , and the expression is:

[0140] (3)

[0141] (4)

[0142] In the formula, is the sound speed of the next water layer;

[0143] The horizontal displacement and the arc length experienced by the beam in this layer are:

[0144] (5)

[0145] In the formula, is the incident angle of the next water layer;

[0146] According to Snell's law: , then the horizontal displacement of the sound ray in the th layer and the time are:

[0147] (6)

[0148] In the formula, is the harmonic mean sound speed within the layer ;

[0149] (7)

[0150] If the incident angle is 0, the propagation trajectory of the sound wave within each layer is a straight line, vertically downward, and the propagation time within each layer is:

[0151] (8).

[0152] (2) Principle of the equivalent sound speed profile method

[0153] In the post-processing of multi-beam data, due to the high accuracy of the ray tracing algorithm, it is widely used, but the model is relatively complex and the calculation amount is large. To address the above problems, Geng proposed the ESSP (Earth System Science Partnership) method. By keeping the beam incident angle and the surface sound speed unchanged, the multi-layer ray tracing problem is transformed into a single-layer constant gradient ray tracing problem. Only an equivalent sound speed profile with the same integral area as the actual complex sound speed profile is needed, which can greatly simplify the calculation. Therefore, the present invention performs an inversion of the equivalent sound speed.

[0154] As Figure 5 shown, in step S2, the optimal equivalent sound speed of each survey line is inverted, including:

[0155] The process of equivalent sound speed homing is:

[0156] (9)

[0157] In the formula, is a constant sound speed gradient, is the depth, is the sound speed, is the bottom sound speed, is the surface sound speed, is the bottom water depth value, is the incident surface sound speed, is a constant;

[0158] The expression for the propagation of sound waves in water conforming to Snell's law is:

[0159] (10)

[0160] In the formula, , are the upper and lower beam incident angles of a single layer respectively, is a constant obtained based on Snell;

[0161] The propagation time is obtained according to the principle of the equivalent sound speed profile method as follows:

[0162] (11)

[0163] Combining Equation (10) and Equation (11), the propagation time is derived as follows:

[0164] (12)

[0165] Since the propagation time is accurately measured and is a known quantity, Equation (12) is transformed to obtain as follows:

[0166] (13)

[0167] wherein are all one-way propagation times;

[0168] According to the theory of ray acoustics, the propagation trajectory of a sound ray in a constant gradient layer is a circular arc, and the radius of the circular arc is as follows:

[0169] (14)

[0170] Through derivation, the displacement of the sound ray in the depth direction and the displacement in the horizontal direction are respectively as follows:

[0171] (15).

[0172] It can be seen therefrom that to a certain extent, the equivalent sound speed profile method can be independent of the dependence on the original sound speed profile data. However, the calculation accuracy of the equivalent sound speed profile is closely related to the accuracy of the reference depth, and the reference depth in the actual situation is a position number that needs to be obtained through other means; from Figure 7 it can be seen that the equivalent sound speed gradient varies with the change of the reference depth, that is, the equivalent sound speed profile method holds for a specific depth and is not equivalent for other depths.

[0173] It can be seen from the above embodiments that the present invention proposes an automatic inversion method of equivalent sound speed based on multi-beam overlapping area constraints, which first uses the surface sound speed of multi-beam data to track the beam points of multi-beam raw data to obtain the overlapping area of ​​the strips; then the multi-beam strips are gridded, and the water depth deviation between the overlapping grids is minimized by using a genetic algorithm, and the optimal equivalent sound speed of each survey line is obtained by inversion; finally, the beam point tracking is performed again for each survey line using the corresponding equivalent sound speed obtained by inversion, thereby reducing the sound speed deviation existing in beam point tracking using the surface sound speed, so that the crying and laughing face-shaped curved terrain between the strips in the overlapping area is weakened, and a more continuous seabed terrain is obtained in the overlapping area of ​​adjacent strips, thereby reducing the dependence of beam point tracking on the measured sound speed profile. Therefore, the method proposed in the present invention effectively fills the gap in the current automatic correction of multi-beam sound speed errors.

[0174] Embodiment 3, as Figure 8 As shown, the equivalent sound velocity automatic inversion system based on multi-beam overlapping area constraint provided by the embodiment of the present invention includes:

[0175] The multi-beam strip overlapping area acquisition module 1 is used to track the beam points of the multi-beam raw data using the surface sound velocity of the multi-beam data to obtain the overlapping area of ​​the multi-beam strips;

[0176] The inversion equivalent sound velocity module 2 is used to grid the overlapping areas of the multi-beam strips, minimize the water depth deviation between the overlapping grids using a genetic algorithm, and invert the optimal equivalent sound velocity for each survey line;

[0177] The re-beam point tracking module 3 is used to re-track the beam points of each survey line using the corresponding equivalent sound speed obtained by inversion, reduce the sound speed deviation existing in beam point tracking using surface sound speed, weaken the curved terrain between strips in the overlapping area, and obtain a continuous seabed terrain in the overlapping area of ​​adjacent strips.

[0178] In order to verify the accuracy of the present invention in areas with drastic changes in sound speed, an experiment was conducted on a water area to verify the reliability of the present method. The experimental data contains four adjacent parallel survey lines. The experimental area data has terrain protrusions and relatively flat areas. This terrain undulation can better verify the accuracy of the inversion algorithm. The experimental area is as follows: Figure 9 shown.

[0179] The data were operated according to the process of this method, and the processed results were compared with the internal longitude. The depth deviation histogram of the overlapping area is as follows: Figure 10As shown, it can be seen that the water depth deviation in the overlapping area after being processed by this method is relatively small, no topographic distortion is generated, which meets the tolerance requirements of the hydrographic survey specifications, and it is possible to get rid of the dependence on the measured sound velocity profile when performing beam position tracking. As described above, only the relatively optimal specific implementation manner of the present invention is given, but the protection scope of the present invention is not limited thereto. Any modification, equivalent replacement, improvement, etc. made by any person skilled in the art within the technical scope disclosed by the present invention and within the spirit and principle of the present invention should be covered by the protection scope of the present invention.

Claims

1. An automatic inversion method for equivalent sound velocity based on multi-beam overlapping area constraint, characterized in that, The method includes the following steps: S1. Use the surface sound velocity of the multibeam data to perform beam point position tracking on the multibeam raw data, and obtain the overlapping area of the multibeam strips; S2. Grid the overlapping area of the multibeam strips, and use the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtain the optimal equivalent sound velocity of each survey line; S3. Use the corresponding equivalent sound velocity obtained by inversion to re-perform beam point position tracking for each survey line, weaken the sound velocity deviation existing in beam point position tracking using the surface sound velocity, weaken the topographic distortion between the strips in the overlapping area, and obtain the continuous real seabed topography of the adjacent strip overlapping area, so as to get rid of the dependence of beam point position tracking on the measured sound velocity information; Using the genetic algorithm to minimize the water depth deviation between the overlapping grids includes: Input the data of the overlapping area of the gridded multibeam strips into the genetic algorithm, and use the genetic algorithm to perform optimization; the genetic algorithm randomly generates a population within a range, and then performs crossover and mutation operations, and calculates the fitness function. The fitness function E(g) is: In the formula, g is the equivalent sound speed to be obtained, representing the sound speed gradient corresponding to N strips in the area; is the water depth of each grid node, is the inverse distance weighted average water depth of all measurement points in the grid, M is the total number of grids in the selected area, K m,n is the total number of measurement points of the nth survey line of the m grid in the area; Judge whether the average water depth deviation of the overlapping area calculated by g of the current fitness function is less than 1% of the average water depth. If it reaches, exit the loop and return the optimal equivalent sound velocity of each survey line. If it does not reach, perform iterative calculation again.

2. The equivalent sound velocity automatic inversion method based on multi-beam overlap area constraint according to claim 1, characterized in that In step S1, obtaining the overlapping area of the multibeam strips includes: S101. Parse the obtained multibeam raw data, and use the surface sound velocity information measured by the surface sound velocimeter in the data or assume a sound velocity profile to perform ray tracing on the data, and obtain the initial water depth and position of all measurement points; S102. Search for the overlapping area between each survey line along the direction perpendicular to the multibeam course.

3. The equivalent sound velocity automatic inversion method based on multi-beam overlapping area constraint according to claim 2, wherein, In step S101, ray tracing includes ray tracing based on constant sound velocity within a layer and ray tracing based on constant gradient within a layer.

4. The equivalent sound velocity automatic inversion method based on multi-beam overlap area constraint according to claim 3, characterized in that Ray tracing based on constant sound velocity within a layer includes: The sound speed propagates at a constant speed within the i-th layer. The depths of the upper and lower interfaces of layer i are z i and z i+1 , and the layer thickness is Δz i , θ i , C i are respectively the beam incident angle and the sound speed of the i-th layer. According to Snell's law, sinθ i / C i = p, then the horizontal displacement y i and the propagation time t i are respectively: The propagation horizontal distance y and propagation time t of the beam through the entire water column are: In the formula, N is the number of water layers; If the incident angle is 0, that is, in single-beam sounding, then there is: y＝0 5. The equivalent sound velocity automatic inversion method based on multi-beam overlap area constraint according to claim 4, wherein Ray tracing based on constant gradient within a layer includes: The speed of sound changes with a constant gradient g within layer i i The actual propagation trajectory of the beam within the layer is a continuous arc segment with a radius of curvature R i The expression is as follows: where C i+1 is the sound speed of the lower water layer; The horizontal displacement y experienced by the beam in this layer i and the arc length S i are as follows: where θ i+1 is the incident angle of the water depth in the next layer; According to Snell's law: Then the horizontal displacement y of the acoustic ray in the i-th layer i and the time t i are: In the formula, is the harmonic mean sound velocity in layer i; If the incident angle is 0, the propagation trajectory of the sound wave within each layer is a straight line, vertically downward, and the propagation time within each layer is:

6. The equivalent sound velocity automatic inversion method based on multi-beam overlapping area constraint according to claim 1, characterized in that The specific method for the genetic algorithm to randomly generate a population within a range, and then perform crossover and mutation operations, and calculate the fitness function is as follows: S201, randomly generate the initial value of the data volume according to the range (a, b) of the initial value, and call the generated initial value the initial population; ​ S202. Calculate the E(g) function values of each seed value of the generated initial population respectively, screen the initial population according to the size of each seed value, retain the seeds with small E(g) function values; then perform mutation and crossover on the seeds with deviations exceeding the threshold according to the mutation probability and crossover probability to form a new population; S203. Recalculate the respective E(g) values of the new population, genetically retain the excellent seed values, replace the inferior individuals, and generate new population individuals; S204. Perform mutation and crossover operations on the newly generated population individuals to generate a new population after mutation and crossover, and then calculate the E(g) of the new population; S205. Repeat steps S203 and S204 until the optimal value is found or the maximum number of iterations is reached, then exit and output the best equivalent sound velocity value.

7. The equivalent sound velocity automatic inversion method based on multi-beam overlap area constraint according to claim 1, characterized in that In step S2, the optimal equivalent sound velocity of each survey line is inversed, including: The process of equivalent sound velocity migration is as follows: In the formula, g is the sound speed gradient which is a constant, Z is the depth, C is the sound speed, C r is the bottom sound speed, C0 is the surface sound speed, Z r is the bottom water depth value, Z0 is the incident surface sound speed, and Constant is a constant value; The expression that the sound wave propagation in water conforms to Snell's law is: where θ0 and θ r are the incident angles of the upper and lower beams of a single layer respectively, and p is a constant obtained based on Snell's law; According to the principle of the equivalent sound velocity profile method, the propagation time t is obtained as: Combining equation (10) and equation (11), the propagation time t is deduced as: Since the propagation time is accurately measured and is a known quantity, Equation (12) is transformed to obtain θ r as follows: In the formula, t is the one-way propagation time; According to the theory of ray acoustics, the propagation trajectory of the sound ray in the constant gradient layer is a circular arc, and the radius R of the circular arc is: After derivation, the displacement Δz of the sound ray in the depth direction and the displacement Δx in the horizontal direction are respectively:

8. An automatic inversion system for equivalent sound velocity based on multi-beam overlap area constraint, characterized in that, This system implements the automatic inversion method of equivalent sound velocity based on multi-beam overlap area constraint as described in any one of claims 1-7. This system includes: An overlapping area acquisition module (1) of the multi-beam strip, which is used to track the beam point positions of the multi-beam raw data by using the surface sound velocity of the multi-beam data, and acquire the overlapping area of the multi-beam strip; An equivalent sound velocity inversion module (2), which is used to grid the overlapping area of the multi-beam strip, and use the genetic algorithm to minimize the water depth deviation between the overlapping grids, and inversely obtain the optimal equivalent sound velocity of each survey line; A re-beam point position tracking module (3), which is used to re-track the beam point positions of each survey line by using the corresponding equivalent sound velocity obtained by inversion, weaken the sound velocity deviation existing in the beam point position tracking by using the surface sound velocity, weaken the topographic distortion between the strips in the overlapping area, and obtain the continuous real seabed topography of the adjacent strip overlapping area splicing, so as to get rid of the dependence of the beam point position tracking on the measured sound velocity information.

9. The equivalent sound velocity automatic inversion system based on multi-beam overlapping area constraint according to claim 8, wherein The automatic inversion system of equivalent sound velocity based on multi-beam overlap area constraint is carried on a computer device. The computer device includes: at least one processor, a memory, and a computer program stored in the memory and running on the at least one processor. When the processor executes the computer program, the functions in the above-mentioned automatic inversion system of equivalent sound velocity based on multi-beam overlap area constraint are realized.

Images