A method for measuring underwater sound velocity profile based on time-of-flight drop-type

Abstract

The present invention discloses a method for measuring underwater sound velocity profile based on time flight, which relates to the field of underwater measurement technology and is used for drawing underwater sound velocity profile. The measuring method comprises the following steps: pre-processing the echo signal of sound velocity measurement by clock correction coefficient, then identifying the reference wave and obtaining underwater sound velocity data by variable threshold method; demodulating the underwater sound velocity data by differential phase shift keying modulation; eliminating abnormal data by local outlier factor algorithm, then filling the sound velocity value of the eliminated data position by extreme gradient boosting tree algorithm and partially predicting the sound velocity profile; drawing the sound velocity profile. The present invention can obtain high-precision underwater sound velocity by only a single measurement, and uses variable threshold method to identify the reference wave, avoiding the limitation of fixed threshold detection, effectively preventing the influence of noise on reference wave identification, and increasing the accuracy of time measurement.

Description

Technical Field

[0001] The invention relates to the technical field of underwater measurement, and in particular to a time-leap-based throwable underwater sound velocity profile measurement method. Background Art

[0002] The ocean sound speed measurement method can be divided into indirect measurement method and direct measurement method. The most common indirect measurement method is to calculate the ocean sound speed in the corresponding temperature, pressure and salinity range according to empirical formulas such as Digoso, Wilson, and Coppens after measuring the physical parameter information such as salinity, temperature and depth in seawater; the most common direct measurement method is the cross-correlation method, which uses the correlation between the ultrasonic transmission signal and the echo signal at different times on the same time axis to obtain the time difference between the ultrasonic packet transmission signal and the echo signal, effectively expressing the leading or lagging relationship between the two signals on the time axis. The time interval when the two signals are most correlated is the time difference between the transmission signal and the echo signal. Another direct measurement method is the pulse ringing method, which first drives the transducer to transmit an ultrasonic signal. The transducer receives the echo signal and again uses it as an excitation signal to transmit the next pulse signal. This process is repeated to form a "ringing". After reaching the preset number of times, the time for multiple round trips of the signal within a fixed distance is obtained.

[0003] In the indirect measurement method, the measurement of physical parameter information such as salinity, temperature and depth is affected by the accuracy of multiple sensors. Secondly, when using empirical formulas to calculate the ocean sound speed, each empirical formula has a different applicable range for temperature, pressure and salinity. Although the applicable ranges of some empirical formulas overlap, in the complex and changeable ocean environment, there are also differences in the accuracy of sound speed measurement, resulting in large errors when drawing sound speed profiles and the inability to transmit data in real time. The measurement accuracy of the cross-correlation method in the direct measurement method is related to the sampling rate of the system hardware. When the sampling rate of the system hardware is constant, the software algorithm can enhance the resolution of time measurement to a certain extent, but the time complexity of the cross-correlation algorithm is relatively high, and it is difficult to calculate the sound speed value quickly; the pulse ring method calculates the sound speed based on the time of multiple round trips within a fixed distance. This leads to additional delay time such as electro-acoustic conversion and circuit in the cyclic triggering process of the pulse ring method, and the distance of its sound wave transmission is generally short. After multiple reflections, the sound wave signal may be less attenuated and can still trigger the circuit to emit the sound wave signal, thereby causing errors. Moreover, due to limitations in manufacturing technology and measurement methods, this method is difficult to meet application scenarios where the sound speed changes rapidly. Therefore, a faster and more accurate sound speed measurement method is needed to achieve real-time transmission of sound speed data and draw a more accurate sound speed profile. Summary of the invention

[0004] In view of the problem that the current sound velocity measurement method is easily affected by noise, resulting in inaccurate time measurement and sound velocity data, the present invention provides a time-leap-based drop-type underwater sound velocity profile measurement method.

[0005] In order to solve the above technical problems, the present invention adopts the following technical solution: a method for measuring underwater sound velocity profile based on time leap, comprising the following steps:

[0006] S1: Preprocess the echo signal of the sound velocity measurement by the clock correction coefficient, and then use the variable threshold method to identify the reference wave and obtain the underwater sound velocity data;

[0007] S2: demodulate the underwater sound speed data using differential phase shift keying modulation;

[0008] S3: The local outlier factor algorithm is used to remove abnormal data, and then the extreme gradient boosting tree algorithm is used to fill in the sound velocity values ​​of the removed data positions and partially predict the sound velocity profile;

[0009] S4: Plot the sound velocity profile.

[0010] Further, in S1, two cycles of a 32.768kHz crystal oscillator are measured by a high-speed crystal oscillator to calculate a correction coefficient of the clock;

[0011] Get the maximum peak value of the echo signal , then we need to calculate the peak signal amplitude Ratio to maximum amplitude: , find the ratio of the two peaks with the largest ratio difference and The median value ;

[0012] Get the threshold function:

[0013] ;

[0014] After filtering, the value with amplitude greater than The first peak point after the search condition is used to find the required reference wave.

[0015] Furthermore, in S1, the corrected time interval is:

[0016] ;

[0017] In the formula, For the measurement results, is the system clock period, is the number of base clock cycles, Indicates the time from the rising edge of the pulse start signal to the next rising edge of the system clock. Indicates the time from the rising edge of the pulse end signal to the next rising edge of the system clock. To measure time accurately, is a reference clock cycle, is two reference clock cycles, is the calibration clock period.

[0018] Furthermore, in S1, the resistance value of the temperature sensor is:

[0019] ;

[0020] In the formula, is the resistance of the reference resistor, is the discharge time of the capacitor through the reference resistor, is the discharge time of the capacitor after it is discharged through the reference resistor and then discharged through the resistor of the temperature sensor;

[0021] The gain error at each measurement point is calculated as:

[0022] ;

[0023] In the formula, The uncorrected temperature result measured at each temperature point, To calculate the ideal value based on the standard reference temperature and known physical model or reference table;

[0024] Taking the average of all gain errors yields the gain correction factor:

[0025] ;

[0026] In the formula, is the number of standard temperature points;

[0027] The gain-corrected temperature result is:

[0028] ;

[0029] In the formula, It is the uncorrected temperature result currently measured by the temperature sensor.

[0030] Furthermore, in S1, the resistance value of the pressure sensor is:

[0031] ;

[0032] In the formula, is the resistance of the reference resistor, To measure the discharge time of the pressure sensor resistor, The discharge time of the capacitor when measuring the reference resistance;

[0033] The uncorrected pressure result is obtained by the relationship between the resistance value of the pressure sensor used and the pressure, and then the corrected pressure result is obtained by using the gain correction coefficient.

[0034] Further, in S2, the original data code is differentially encoded to obtain a code element sequence, a square root raised cosine filter is used to obtain an output signal, and the output signal is multiplied by a carrier to obtain a differential phase shift keying modulated digital signal. ;

[0035] Differential Phase Shift Keying Modulated Digital Signal The time domain expression is:

[0036] ;

[0037] In the formula, is the carrier amplitude, is the carrier angular frequency, The value is +1 or -1. is the duration, For duration A gate function in It is equal to 1 during the time period and equal to 0 at other times;

[0038] In one symbol period, when hour, ,when hour, .

[0039] Furthermore, in S2, the differential phase-shift keying modulated digital signal is multiplied with the carrier recovered by the Costas loop, and coherent demodulation is realized after completing the first-level filtering. After coarse decoding, the lead-lag phase-locked loop method is used to achieve bit synchronization, and then the judgment result is differentially decoded under the bit synchronization clock to restore the original code data and obtain the underwater sound speed, temperature and pressure data.

[0040] Furthermore, in S3, the sound velocity value point As the center of the circle, Proximity distance A circle is formed with a radius, and the range of the circle is the k-distance neighborhood:

[0041] ;

[0042] In the formula, The sound velocity value point of Proximity distance, is the distance from the sound speed point Among the recent points, The nearest point and Point The distance between

[0043] Points by speed of sound is the center of the circle, the speed of sound value point To point No. The reachable distance is: ;

[0044] In the formula, The sound velocity value point of Proximity distance, is the distance from the sound speed point Among the recent points, The nearest point and Point The distance between

[0045] Sound speed value point The local reachable density is:

[0046] ;

[0047] In the formula, The sound velocity value point Centered A sound speed point in the distance neighborhood;

[0048] Get the sound speed value point The local outlier factor of is:

[0049] ;

[0050] if If it is less than 1, it means The density of is higher than that of its neighboring points. is a dense point; if If it is greater than 1, it means The density of is less than the density of its neighboring points. are outliers and should be removed.

[0051] Furthermore, in S3, the extreme gradient boosting tree algorithm is used to minimize the objective function:

[0052] ;

[0053] In the formula, For the The true value of the sound speed samples, For the The predicted value of the sound speed samples, is the loss function used to measure the gap between the predicted sound speed data and the actual sound speed data, A regularization term used to control the complexity of the model and prevent overfitting;

[0054] The regularization term is:

[0055] ;

[0056] In the formula, is the number of leaf nodes of the decision tree model, For the The weight of a leaf, , is the regularization parameter used to control the complexity of the model;

[0057] The extreme gradient boosting tree algorithm uses a forward distribution algorithm. Each iteration builds a new tree by optimizing the objective function. Given the prediction value of the current model , No. The goal of this iteration is to find a new function To minimize the objective function:

[0058] ;

[0059] In the formula, For the model The input features on the sound speed sample are used to generate the predicted value of the sound speed sample ;

[0060] The second-order Taylor expansion is used to approximate the objective function:

[0061] ;

[0062] In the formula, is the first-order derivative of the loss function with respect to the predicted value, is the second-order derivative of the loss function with respect to the predicted value;

[0063] The optimization problem is transformed into minimizing the objective function:

[0064] .

[0065] Furthermore, for a particular tree, the predicted value of its leaf node is:

[0066] ;

[0067] In the formula, is the mapping function from sample to leaf;

[0068] The objective function is expressed as:

[0069] ;

[0070] In the formula, is the weight of the leaf node, For the The gradient sum of leaf nodes, For the The sum of the diagonal elements of the second-order derivatives of the loss functions of leaf nodes;

[0071] All are assigned to leaf nodes The sample index set for is:

[0072] ;

[0073] Through Taking the derivative and setting the inverse to 0, we get the optimal solution for the leaf node weight:

[0074] ;

[0075] Will Substituting into the objective function, the minimized objective value is:

[0076] .

[0077] The beneficial effects of the present invention are as follows: the present invention can obtain high-precision underwater sound speed through only a single measurement, and adopts a variable threshold method to identify the reference wave, avoiding the limitation of fixed threshold detection, and can effectively prevent the influence of noise on the reference wave identification, increase the accuracy of time measurement, and be used to draw underwater sound speed profiles. The present invention adopts a square root raised cosine filter to limit the signal bandwidth and eliminate inter-symbol interference in modulation, and adopts Costas loop and lead-lag phase-locked method position synchronization technology in demodulation to complete the recovery of the original code element, so that the data transmission is more accurate. The present invention adopts a local outlier factor algorithm to remove abnormal values ​​from the obtained sound speed data, and then adopts an extreme gradient boosting tree algorithm to fill in the blanks of the removed sound speed data, and completes a certain prediction of the subsequent sound speed profile, which is conducive to the drawing of the sound speed profile. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] Figure 1 Shown is a schematic diagram of the principle of an embodiment of the present invention.

[0079] Figure 2 Shown is a schematic diagram of the system operation of the jettisonable underwater sound velocity profile measurement.

[0080] Figure 3 Shown is a block diagram of the working principle of the underwater data acquisition system.

[0081] Figure 4 Shown is a block diagram of the working principle of the long-distance data transmission communication system.

[0082] Figure 5 Shown is a block diagram of the working principle of the water data receiving and processing system. DETAILED DESCRIPTION

[0083] A method for measuring underwater sound velocity profile based on time-flight drop-type includes the following steps:

[0084] S1: Preprocess the echo signal of the sound velocity measurement by the clock correction coefficient, and then use the variable threshold method to identify the reference wave and obtain the underwater sound velocity data;

[0085] S2: demodulate the underwater sound speed data using differential phase shift keying modulation;

[0086] S3: The local outlier factor algorithm is used to remove abnormal data, and then the extreme gradient boosting tree algorithm is used to fill in the sound velocity values ​​of the removed data positions and partially predict the sound velocity profile;

[0087] S4: Plot the sound velocity profile.

[0088] Further, in S1, two cycles of a 32.768kHz crystal oscillator are measured by a high-speed crystal oscillator to calculate a correction coefficient of the clock;

[0089] Get the maximum peak value of the echo signal , then we need to calculate the peak signal amplitude Ratio to maximum amplitude: , find the ratio of the two peaks with the largest ratio difference and The median value ;

[0090] Get the threshold function:

[0091] ;

[0092] After filtering, the value with amplitude greater than The first peak point after the search condition is used to find the required reference wave.

[0093] Furthermore, in S1, the corrected time interval is:

[0094] ;

[0095] In the formula, For the measurement results, is the system clock period, is the number of base clock cycles, Indicates the time from the rising edge of the pulse start signal to the next rising edge of the system clock. Indicates the time from the rising edge of the pulse end signal to the next rising edge of the system clock. To measure time accurately, is a reference clock cycle, is two reference clock cycles, is the calibration clock period.

[0096] Furthermore, in S1, the resistance value of the temperature sensor is:

[0097] ;

[0098] In the formula, is the resistance of the reference resistor, is the discharge time of the capacitor through the reference resistor, is the discharge time of the capacitor after it is discharged through the reference resistor and then discharged through the resistor of the temperature sensor;

[0099] The gain error at each measurement point is calculated as:

[0100] ;

[0101] In the formula, The uncorrected temperature result measured at each temperature point, To calculate the ideal value based on the standard reference temperature and known physical model or reference table;

[0102] Taking the average of all gain errors yields the gain correction factor:

[0103] ;

[0104] In the formula, is the number of standard temperature points;

[0105] The gain-corrected temperature result is:

[0106] ;

[0107] In the formula, It is the uncorrected temperature result currently measured by the temperature sensor.

[0108] Furthermore, in S1, the resistance value of the pressure sensor is:

[0109] ;

[0110] In the formula, is the resistance of the reference resistor, To measure the discharge time of the pressure sensor resistor, The discharge time of the capacitor when measuring the reference resistance;

[0111] The uncorrected pressure result is obtained by the relationship between the resistance value of the pressure sensor used and the pressure, and then the corrected pressure result is obtained by using the gain correction coefficient.

[0112] Further, in S2, the original data code is differentially encoded to obtain a code element sequence, a square root raised cosine filter is used to obtain an output signal, and the output signal is multiplied by a carrier to obtain a differential phase shift keying modulated digital signal. ;

[0113] Differential Phase Shift Keying Modulated Digital Signal The time domain expression is:

[0114] ;

[0115] In the formula, is the carrier amplitude, is the carrier angular frequency, The value is +1 or -1. is the duration, For duration A gate function in It is equal to 1 during the time period and equal to 0 at other times;

[0116] In one symbol period, when hour, ,when hour, .

[0117] Furthermore, in S2, the differential phase-shift keying modulated digital signal is multiplied with the carrier recovered by the Costas loop, and coherent demodulation is realized after completing the first-level filtering. After coarse decoding, the lead-lag phase-locked loop method is used to achieve bit synchronization, and then the judgment result is differentially decoded under the bit synchronization clock to restore the original code data and obtain the underwater sound speed, temperature and pressure data.

[0118] Furthermore, in S3, the sound velocity value point As the center of the circle, Proximity distance A circle is formed with a radius, and the range of the circle is the k-distance neighborhood:

[0119] ;

[0120] In the formula, The sound velocity value point of Proximity distance, is the distance from the sound speed point Among the recent points, The nearest point and Point The distance between

[0121] Points by speed of sound is the center of the circle, the speed of sound value point To point No. The reachable distance is: ;

[0122] In the formula, The sound velocity value point of Proximity distance, is the distance from the sound speed point Among the recent points, The nearest point and Point The distance between

[0123] Sound speed value point The local reachable density is:

[0124] ;

[0125] In the formula, The sound velocity value point Centered A sound speed point in the distance neighborhood;

[0126] Get the sound speed value point The local outlier factor of is:

[0127] ;

[0128] if Less than 1, it means The density of is higher than that of its neighboring points. is a dense point; if If it is greater than 1, it means The density of is less than the density of its neighboring points. are outliers and should be removed.

[0129] Furthermore, in S3, the extreme gradient boosting tree algorithm is used to minimize the objective function:

[0130] ;

[0131] In the formula, For the The true value of the sound speed samples, For the The predicted value of the sound speed samples, is the loss function used to measure the gap between the predicted sound speed data and the actual sound speed data, A regularization term used to control the complexity of the model and prevent overfitting;

[0132] The regularization term is:

[0133] ;

[0134] In the formula, is the number of leaf nodes of the decision tree model, For the The weight of a leaf, , is the regularization parameter used to control the complexity of the model;

[0135] The extreme gradient boosting tree algorithm uses a forward distribution algorithm. Each iteration builds a new tree by optimizing the objective function. Given the prediction value of the current model , No. The goal of this iteration is to find a new function To minimize the objective function:

[0136] ;

[0137] In the formula, For the model The input features on the sound speed sample are used to generate the predicted value of the sound speed sample ;

[0138] The second-order Taylor expansion is used to approximate the objective function:

[0139] ;

[0140] In the formula, is the first-order derivative of the loss function with respect to the predicted value, is the second-order derivative of the loss function with respect to the predicted value;

[0141] The optimization problem is transformed into minimizing the objective function:

[0142] .

[0143] Furthermore, for a particular tree, the predicted value of its leaf node is:

[0144] ;

[0145] In the formula, is the mapping function from sample to leaf;

[0146] The objective function is expressed as:

[0147] ;

[0148] In the formula, is the weight of the leaf node, For the The gradient sum of leaf nodes, For the The sum of the diagonal elements of the second-order derivatives of the loss functions of leaf nodes;

[0149] All are assigned to leaf nodes The sample index set for is:

[0150] ;

[0151] Through Taking the derivative and setting the inverse to 0, we get the optimal solution for the leaf node weight:

[0152] ;

[0153] Will Substituting into the objective function, the minimized objective value is:

[0154] .

[0155] Combination Figure 1 and Figure 2 As shown, a throwable underwater sound velocity profile measurement system based on time flight includes an underwater data acquisition system, a data long-distance transmission communication system, and an above-water data receiving and processing system. The underwater data acquisition system includes a core control module, a time measurement module, a pressure measurement module, a temperature measurement module, and a system storage module. A throwable underwater sound velocity profile measurement method based on time flight includes the following specific steps:

[0156] like Figure 3 As shown, in the first step, the core control module controls the time measurement module through the SPI communication interface to complete its register configuration, and then sends instructions to calibrate the high-speed crystal oscillator. The high-speed crystal oscillator is used to measure two cycles of the 32.768kHz crystal oscillator, and the clock correction coefficient is calculated to correct the measured underwater sound speed data.

[0157] After obtaining the clock correction coefficient, the core control module sends the time measurement operation code. The time measurement module will transmit pulses of corresponding frequency and number according to the parameters in the register configuration, and transmit the ultrasonic drive signal from the Fire end into the transmitting module. The transmitting module matches the frequency of the ultrasonic drive signal with the transducer frequency through tuning and variable resistance, so that the mechanical wave signal converted and emitted by the transducer is large enough to increase the energy of the ultrasonic drive circuit to drive the transducer.

[0158] The time measurement module has two ports for processing ultrasonic signals, namely STOP1 and STOP2. STOP1 and STOP2 are connected to the same transceiver through capacitors. When the Fire end transmits the ultrasonic drive signal, the core control module controls the external analog switch to use the STOP1 port as the time measurement start signal. After the transducer sends the ultrasonic signal, the echo signal reflected by the reflector is received by the transducer. The echo signal is generally weak, usually tens of millivolts, and needs to be processed by parallel limiting, frequency selection amplification and signal amplification before it can be recognized. After the time measurement module recognizes the echo signal, the core control module opens the STOP2 port by controlling the external analog switch to receive the processed echo signal.

[0159] The first wave detection with a fixed threshold is to set a fixed amplitude threshold. After the echo signal reaches this value after filtering, the underwater data acquisition system automatically measures the half-cycle width of the first wave and the following waves until the half-cycle width is consistent, and the first wave with the same cycle width is considered as the first echo received by the underwater data acquisition system. If the signal is affected by noise or other external interference, the amplitude is attenuated, etc., the first wave detection will be difficult to perform accurately, and the identified reference wave is likely to be wrong. Therefore, in view of the limitations of fixed threshold detection, a variable threshold method is used to identify the reference wave. The variable threshold method can dynamically adjust the detection threshold according to the change of the actual signal amplitude, making it more flexible, thereby improving the accuracy of detection.

[0160] For the received ultrasonic echo signal, the signal amplitude will suddenly increase or decrease after being interfered by external noise, causing random changes in the amplitude of each periodic signal. The original maximum peak signal may be attenuated due to interference, and the amplitude of the signal close to it may surge due to interference. However, no matter how the position of the echo signal with the largest amplitude jumps, the ratio of the value of the maximum amplitude point of each signal when the envelope slope is a positive number to the amplitude of the echo signal with the largest amplitude remains basically unchanged. According to this rule, a threshold function can be set, and a certain echo signal can be accurately found based on it. The echo signal is set as the reference wave, and it can be accurately found every time among the other echo signals. The setting of the threshold function needs to be based on the actual signal amplitude.

[0161] First, we need to find the maximum peak value of the echo signal , then we need to calculate the peak signal amplitude Ratio to maximum amplitude , find the ratio of the two peaks with the largest ratio difference and The median value of:

[0162] ;

[0163] The threshold function is obtained as:

[0164] .

[0165] After filtering, the value with amplitude greater than The first peak point after the signal is used as the search condition to accurately find the required reference wave. The reference wave, as the reference point of the signal, provides a certain starting time. By identifying the reference wave, the system can eliminate the influence of noise and interference when receiving the signal and only measure the time of stable and reliable signals.

[0166] Assume that the start signal is Start, the end signal is Stop, and the entire measurement time interval is from the rising edge of the Start pulse to the rising edge of the Stop pulse. Indicates the number of system clock cycles that have passed. Indicates the time from the rising edge of the pulse Start to the next rising edge of the system clock. Indicates the time from the rising edge of the Stop pulse to the next rising edge of the system clock. The chip uses a coarse counter to calculate , calculated using an accurate counter and .

[0167] It can be concluded that the corrected time interval is:

[0168] ;

[0169] In the formula, For the measurement results, is the system clock period, is the number of base clock cycles, and The difference is the precise measurement time. is a reference clock cycle, is two reference clock cycles, and The difference is the calibration clock period. The measured data is stored in the result register. After the measurement is completed, an interrupt signal is sent to the core control module. The core control module reads the time data and calculates the sound speed.

[0170] The core control module sends the temperature measurement operation code, and the temperature measurement module starts the temperature measurement process. First, the capacitor passes through the reference resistor Discharge is performed and the discharge time is recorded as Then the capacitor is discharged through the resistance of the temperature sensor, and the discharge time is also measured and recorded as The core control module will read the measurement results stored in the result register. Since the time it takes for a capacitor to discharge through different resistors is proportional to the resistance value, the resistance value of the temperature sensor can be calculated by knowing the discharge time of the reference resistor:

[0171] .

[0172] When measuring temperature, the delay time of the Schmitt trigger will introduce a gain error. By querying the temperature table, you can obtain the uncorrected temperature result currently measured by the temperature sensor. ,In order to obtain more accurate temperature measurement results, gain correction is required.

[0173] Assume that standard temperature points, and the uncorrected temperature result measured at each temperature point is , standard temperature is , calculated from the standard reference temperature and known physical models or reference tables to obtain the ideal value , the gain error of each measurement point is calculated as:

[0174] ;

[0175] Taking the average of all gain errors yields the gain correction factor:

[0176] ;

[0177] The gain-corrected temperature result is:

[0178] .

[0179] The core control module sends the pressure measurement operation code, and the pressure measurement module starts the pressure measurement process. The process of measuring pressure data is similar to temperature measurement. Both convert the resistance change of the sensor into a high-precision time interval measurement. The piezoresistive pressure sensor is selected, and its resistance value changes with pressure. The discharge time when measuring the pressure sensor resistance , measure the reference resistance Capacitor discharge time , the resistance value of the pressure sensor is:

[0180] .

[0181] The uncorrected pressure result is obtained by the relationship between the resistance value of the pressure sensor used and the pressure, and the corrected pressure result is obtained by the gain correction coefficient.

[0182] like Figure 4As shown, in the second step, after obtaining the temperature, pressure and sound speed data signals of the underwater data acquisition system, due to the characteristics of the ocean channel, the original data of the unipolar non-return-to-zero code cannot be directly transmitted on the enameled wire, otherwise there will be serious waveform distortion and bit timing recovery cannot be performed, so differential phase-shift keying modulation and demodulation are used.

[0183] First, the original data code Perform differential encoding to obtain the codeword sequence In the communication system, inter-symbol interference will cause distortion of the received signal, so the square root raised cosine filter is used to limit the signal bandwidth and eliminate the inter-symbol interference. The phase selection module multiplies the signal output by the square root raised cosine filter with the carrier of a specific frequency, and the phase selection module outputs the differential phase shift keying modulated digital signal. .

[0184] Phase shift keying modulated digital signal The time domain expression is:

[0185] ;

[0186] In the formula, is the carrier amplitude, is the carrier angular frequency, The value is +1 or -1. is the duration, For duration A gate function, exist It is equal to 1 during the time period and equal to 0 at other times. The probability of taking 1 is , The probability of taking -1 is Then the DPSK signal in one symbol period is:

[0187] when hour, ;when hour, .

[0188] The digital-to-analog converter converts the above signal into a digital-to-analog converter, and then connects the signal to the drive amplifier circuit to improve the signal driving capability, and transmits the data to the water receiving circuit via the enameled wire in a differential transmission mode. After the water receiving circuit receives the data transmitted by the enameled wire, it first passes the signal through a high-pass filter circuit to filter out the high-voltage DC component used to start the underwater system and the low-frequency interference mixed in the channel, and then sends the signal to the differential input end of the instrument amplifier for amplification. After the signal is amplified by the instrument amplifier, it is low-pass filtered to filter out the high-frequency interference. The high-performance analog-to-digital converter converts the received differential phase-shift keying analog electrical signal into a digital signal and inputs it into the demodulation system for demodulation processing.

[0189] The demodulation system first filters out the noise outside the effective frequency band through a bandpass filter to improve the signal-to-noise ratio of the received signal. The differential phase-shift keying modulated signal is multiplied with the carrier recovered by the Costas loop, and coherent demodulation is completed after completing a first-level filtering. However, the signal after coherent demodulation is not the original code element sequence. Next, coarse decoding is performed to convert the signal into 0 and 1 signals. The coarse decoding method is to calculate the average value of hundreds of adjacent sampling points. When the sampling point is greater than the average value, it is judged as 1, otherwise it is judged as 0. At this time, the coarse decoding result at a certain sampling rate will have many consecutive 1s or 0s. According to the sampling clock, it is impossible to determine which sampling points belong to a symbol. At this time, in order to realize the symbol bit sampling decision, bit synchronization technology is required to obtain an accurate sampling decision clock to determine the decision sampling point.

[0190] The bit synchronization technology obtains the symbol bit decision clock, and the lead-lag phase-locked loop method is used to achieve bit synchronization. The falling edge of the bit synchronization clock is the middle moment of a complete code element cycle. The sampling at this point is used as the final decision output. After the bit synchronization clock is inverted, the sampling decision can be completed at the rising and falling edges of the bit synchronization clock. The decision result can be differentially decoded under the bit synchronization clock to restore the original code element data. After demodulation, the sound velocity, temperature and pressure data obtained by the underwater acquisition system can be obtained.

[0191] like Figure 5 As shown in the figure, in the third step, after the water data receiving and processing system receives the data transmitted by the long-distance transmission communication system, the upper computer and the lower computer communicate through the serial port, and the communication protocol is implemented in the upper computer software and the lower computer, so that the upper computer controls the lower computer and reads and processes the data. In the measurement of sound velocity profile data, the sound velocity in each small area changes little, but it will also be affected by various marine environmental noises, causing a few sound velocity data to be abnormal.

[0192] The local outlier factor algorithm is used to remove abnormal data. Among the recent points, The nearest point and Point The distance between them is:

[0193] ;

[0194] This distance is called the sound speed value point of Proximity distance.

[0195] Points by speed of sound As the center of the circle, Proximity distance A circle is formed with a radius of Distance Neighborhood:

[0196] ;

[0197] In the formula, is a certain sound speed value point within the circle.

[0198] Points by speed of sound As the center of the circle, a certain sound speed value point To point No. The reachable distance is:

[0199] ;

[0200] The local reachable density at this sound speed value is:

[0201] ;

[0202] In the formula, The sound velocity value point Centered A sound speed point in the distance neighborhood;

[0203] From this, we can get the local outlier factor of the sound speed value point as:

[0204] .

[0205] If the local outlier factor is closer to 1, it means The density of the neighborhood points is similar; if the local outlier factor is less than 1, it means The density of is higher than that of its neighboring points. is a dense point; if the local outlier factor is greater than 1, it means The density of is less than the density of its neighboring points. It is an outlier and should be removed.

[0206] The seawater sound velocity gradient has certain stratified characteristics, mainly including the surface layer, the transition layer and the deep sea isotherm layer. The change law of the sound velocity is different in different layers, but the change trend of the sound velocity in the same layer is usually consistent. In order to improve the accuracy of the measurement data, the local outlier factor algorithm can be used to identify and remove outliers of the sound velocity, thereby enhancing the reliability of the data.

[0207] After removing abnormal values, the water data receiving and processing system uses the extreme gradient boosting tree algorithm to fill in the sound speed values ​​of the removed data positions and partially predict the sound speed profile. The goal of the extreme gradient boosting tree algorithm is to minimize the objective function:

[0208] ;

[0209] In the formula, For the The true value of the sound speed samples, For the The predicted value of the sound speed samples, is the loss function used to measure the gap between the predicted sound speed data and the actual sound speed data, A regularization term used to control the complexity of the model and prevent overfitting.

[0210] The regularization term is in the form of:

[0211] ;

[0212] In the formula, is the number of leaf nodes of the decision tree model, For the The weight of a leaf, and is the regularization parameter used to control the complexity of the model.

[0213] The extreme gradient boosting tree algorithm uses a forward distribution algorithm. Each iteration builds a new tree by optimizing the objective function. Given the prediction value of the current model , No. The goal of this iteration is to find a new function To minimize the objective function:

[0214] ;

[0215] The second-order Taylor expansion is used to approximate the objective function:

[0216] ;

[0217] In the formula, is the first-order derivative of the loss function with respect to the predicted value, is the second-order derivative of the loss function with respect to the predicted value;

[0218] Then the optimization problem can be transformed into minimizing the objective function:

[0219] ;

[0220] For a particular tree, the predicted value of its leaf node is:

[0221] ;

[0222] In the formula, is the mapping function from sample to leaf, and the weight of leaf node is , so the objective function can be expressed as:

[0223] ;

[0224] In the formula, For the The gradient sum of leaf nodes, For the The sum of the diagonal elements of the second-order derivatives of the loss functions of leaf nodes.

[0225] All are assigned to leaf nodes The sample index set for is:

[0226] ;

[0227] Through Taking the derivative and setting the inverse to 0 gives the optimal solution for the leaf node weight:

[0228] ;

[0229] Will Substituting into the objective function, the minimized objective value is:

[0230] ;

[0231] The extreme gradient boosting tree algorithm is used to predict the sound speed value at the removed data to improve the accuracy of the data.

[0232] The fourth step is to use the sound velocity data at a known depth to predict the sound velocity data at a certain subsequent depth and draw a sound velocity profile using the extreme gradient boosting tree algorithm.

[0233] Embodiment: First, obtain the system clock correction coefficient in the underwater data acquisition system. The core control module is responsible for controlling the measurement unit to complete the acquisition of sound velocity, temperature and pressure data. The collected data is transmitted to the data receiving and processing system on the water by means of a remote data transmission system, and the data is finally processed to generate a sound velocity profile. Table 1 shows the ultrasonic propagation time and sound velocity data measured at a fixed distance between the ultrasonic wave and the baffle under the condition that the temperature of the constant temperature water tank is 13.868 degrees Celsius. In the experiment, a certain length of enameled wire is used to simulate remote transmission, and the entire device is placed in a constant temperature water tank. After collecting the data, the abnormal values ​​in the data are eliminated with the help of the upper computer software, and the average value of the data in the interval time is calculated. In the data measured in the constant temperature water tank, the data fluctuation is small, and the sound velocity data value is accurate and stable.

[0234] Table 1 Ten sets of sound velocity data values ​​measured at a water temperature of 13.868 degrees

[0235]

[0236] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for measuring underwater sound velocity profile based on time-flight, characterized in that: It includes an underwater data acquisition system, a data long-distance transmission communication system and an above-water data receiving and processing system. The underwater data acquisition system includes a core control module, a time measurement module, a pressure measurement module, a temperature measurement module and a system storage module. Specifically, it includes the following steps: S1, pre-processing the echo signal of the sound velocity measurement by the clock correction coefficient, and then using the variable threshold method to identify the reference wave and obtain the underwater sound velocity data; S2, demodulates the underwater sound speed data using differential phase shift keying modulation; S3, using the local outlier factor algorithm to remove abnormal data, and then using the extreme gradient boosting tree algorithm to fill in the sound velocity values ​​of the removed data positions and partially predict the sound velocity profile; S4, plot the sound velocity profile.

2. The method for measuring underwater sound velocity profile based on time leap according to claim 1 is characterized in that: In S1, two cycles of a 32.768kHz crystal oscillator are measured by a high-speed crystal oscillator, and the correction coefficient of the clock is calculated; Get the maximum peak value V of the echo signal max , then we need to calculate the peak signal amplitude V i Ratio to maximum amplitude: Find the ratio R of the two peaks with the largest ratio difference diff1 and R diff2 The median value Get the threshold function: V threshold =KV max ; After filtering, the amplitude is greater than V thrsshold The first peak point after the search condition is used to find the required reference wave.

3. The method for measuring underwater sound velocity profile based on time leap according to claim 2 is characterized in that: In S1, the corrected time interval is: Where Time is the measurement result, T ref is the system clock period, C c is the number of reference clock cycles, F1 represents the time from the rising edge of the pulse start signal to the next system clock rising edge, F2 represents the time from the rising edge of the pulse end signal to the next system clock rising edge, F1-F2 is the precise measurement time, C1 is one reference clock cycle, C2 is two reference clock cycles, and C2-C1 is the calibration clock cycle.

4. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 3 is characterized in that: In S1, the temperature measurement module measures the temperature through a temperature sensor, and the resistance value of the temperature sensor is: In the formula, R ref is the resistance of the reference resistor, t Rref is the discharge time of the capacitor through the reference resistor, t temp is the discharge time of the capacitor after it is discharged through the reference resistor and then discharged through the resistor of the temperature sensor; The gain error at each measurement point is calculated as: In the formula, The uncorrected temperature result measured at each temperature point, To calculate the ideal value based on the standard reference temperature and known physical model or reference table; Taking the average of all gain errors yields the gain correction factor: Where N is the number of standard temperature points; The gain-corrected temperature result is: Where, T uncorr It is the uncorrected temperature result currently measured by the temperature sensor.

5. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 4 is characterized in that: In S1, the pressure measurement module measures pressure through a pressure sensor, and the resistance value of the pressure sensor is: In the formula, R preref is the resistance of the reference resistor, t pressure To measure the discharge time of the pressure sensor resistor, t preref The discharge time of the capacitor when measuring the reference resistance; The uncorrected pressure result is obtained by the relationship between the resistance value of the pressure sensor used and the pressure, and then the corrected pressure result is obtained by using the gain correction coefficient.

6. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 5, characterized in that: In S2, the original data code is differentially encoded to obtain a code element sequence, a square root raised cosine filter is used to obtain an output signal, and the output signal is multiplied by the carrier to obtain a differential phase shift keying modulated digital signal e DPSK (t); Differential Phase Shift Keying modulated digital signal DPSK The time domain expression of (t) is: e DPSK (t)=[∑ n a n g(t-nT s )]Acos(ω c t); Where A is the carrier amplitude, ω c t is the carrier angular frequency, a n The value is +1 or -1, T s is the duration, g(t) is the duration T s A gate function, at 0-T s It is equal to 1 during the time period and equal to 0 at other times; In one symbol period, when a n =1, e DPSK (t) = Acos(ω c t), when a n = -1, e DPSK (t) = -Acos(ω c t).

7. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 6, characterized in that: In S2, the differential phase-shift keying modulated digital signal is multiplied with the carrier recovered by the Costas loop, and coherent demodulation is realized after completing the first-level filtering. After coarse decoding, the lead-lag phase-locked loop method is used to achieve bit synchronization, and then the judgment result is differentially decoded under the bit synchronization clock to restore the original code data and obtain the underwater sound speed, temperature and pressure data.

8. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 7, characterized in that: In S3, with the sound velocity point O as the center, the distance d k (O) forms a circle with radius, and the range of the circle is the k-distance neighborhood: R k (O)={d(O,P′)≤d k (P)}; Where, d k (P) is the k-neighboring distance of the sound speed value point P, and d(O, P′) is the distance between the kth nearest point P′ and point O among the points closest to the sound speed value point O; With the sound speed value point O as the center, the kth reachable distance from the sound speed value point P to point O is: Distance k (O, P) = max{d k (O), d(O, P)}; Where, d k (O) is the k-neighboring distance of the sound speed value point O, and d(O, P) is the distance between the kth nearest point P and point O among the points closest to the sound speed value point O; The local reachable density of the sound speed value point O is: Where X is a sound speed point in the k-distance neighborhood of the sound speed point O as the center of the circle; The local outlier factor of the sound speed value point O is obtained as: If Factor k (O) is less than 1, which means that the density of X is higher than the density of its neighboring points, and X is a dense point. k If (O) is greater than 1, it means that the density of X is less than the density of its neighborhood points. X is an outlier and should be removed.

9. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 8, characterized in that: In S3, the extreme gradient boosting tree algorithm is used to minimize the objective function: Where yi is the true value of the i-th sound speed sample, is the predicted value of the i-th sound speed sample, is the loss function used to measure the gap between the predicted sound speed data and the actual sound speed data, Ω(f k ) is a regularization term used to control the complexity of the model and prevent overfitting; The regularization term is: Where T is the number of leaf nodes in the decision tree model, ω j is the weight of the jth leaf, γ and λ are regularization parameters used to control the complexity of the model; The extreme gradient boosting tree algorithm uses a forward distribution algorithm. Each iteration builds a new tree by optimizing the objective function. Given the prediction value of the current model The goal of the t+1th iteration is to find a new function f t+1 To minimize the objective function: In the formula, x i is the input feature of the model on the i-th sound speed sample, which is used to generate the predicted value of the sound speed sample The second-order Taylor expansion is used to approximate the objective function: In the formula, is the first-order derivative of the loss function with respect to the predicted value, is the second-order derivative of the loss function with respect to the predicted value; The optimization problem is transformed into minimizing the objective function:

10. The method for measuring underwater sound velocity profile based on time-flight drop-type according to claim 9, characterized in that: For a particular tree, the predicted value of its leaf node is: Where q is the mapping function from sample to leaf; The objective function is expressed as: Where ω is the weight of the leaf node, is the gradient sum of the jth leaf node, is the sum of the diagonal elements of the second-order derivative of the loss function of the j-th leaf node; The set of all sample indexes assigned to leaf node j is: I j ={i|q(x i )=j}; Through the ω j Taking the derivative and setting the inverse to 0, we get the optimal solution for the leaf node weight: Will Substituting into the objective function, the minimized objective value is:

Images