Dynamic Optimization Method for Broadband ADCP Waveform with Environment Constraints

Broadband acoustic Doppler current profiler (ADCP) is widely used in agricultural water resource explorations, such as river discharge monitoring and flood warning. Improving the velocity estimation accuracy of broadband ADCP by adjusting the waveform parameters of a phase-encoded signal will reduce the velocity measurement range and water stratification accuracy, while the promotion of stratification accuracy will degrade the velocity estimation accuracy. In order to minimize the impact of these two problems on the measurement results, the ADCP waveform optimization problem that satisfies the environment constraints while keeping high velocity estimation accuracy or stratification accuracy is studied. Firstly, the relationship between velocity or distance estimation accuracy and signal waveform parameters is studied by using an ambiguity function. Secondly, the constraints of current velocity range, velocity distribution and other environmental characteristics on the waveform parameters are studied. For two common measurement applications, two dynamic configuration methods of waveform parameters with environmental adaptability and optimal velocity estimation accuracy or stratification accuracy are proposed based on the nonlinear programming principle. Experimental results show that compared with the existing methods, the velocity estimation accuracy of the proposed method is improved by more than 50%, and the stratification accuracy is improved by more than 22%.


Introduction
Acoustic Doppler current profiler (ADCP) is a popular Doppler sonar used for measuring water velocity and profile discharge in natural and human-made waterways. Due to the advantages of high estimation accuracy, multiple measurement parameters and low measurement cost, the ADCP has been widely used in the fields of hydrological survey, water resources exploration, underwater protection and navigation [1][2][3][4][5], etc. For example, 98% of United States Geological Survey (USGS) non-wading river discharge measurements were performed by ADCP in 2013 [6]. ADCP divides the river section into several vertical subsections. The average velocity and discharge of each subsection are obtained by measuring the velocity of multiple vertical layers in the subsection. Then, the discharge of all subsections is superimposed to get the discharge of the entire river section [7]. Therefore, high-precision velocity measurement is the cornerstone of accurate section discharge measurement. The following is a brief introduction to the concept and importance of stratification or layering. Stratification refers to dividing the entire echo into several segments during signal processing. Each segment is the echo generated by the scatterers in the corresponding depth water layer being illuminated by the transmitted sound wave. The velocity and discharge of each layer are estimated separately. Since the current velocity changes with depth, stratification can improve the accuracy of velocity and discharge measurement of entire river section and obtain the relationship between velocity and depth. The actual achievable layer thickness is generally determined by the transmitted pulse width of ADCP.
The thinner the layer, the more details about the current velocity change with depth can be obtained, which is very important for many hydrological survey applications.
The measurement methods used in ADCP include pulse incoherent method [8], pulse coherent method [9][10][11] and broadband method [12][13][14]. Since the broadband method combines the advantages of the incoherent method and the coherent method, it has become the most popular choice. The accuracy of velocity estimation directly determines the performance of velocity distribution estimation, profile discharge estimation, underwater navigation and other application requirements, while the stratification accuracy determines the observation scale of velocity variation with depth. Therefore, the velocity estimation accuracy and water stratification accuracy are the core performance parameters of ADCP. The broadband method usually uses a phase-encoded pulse as the transmitted signal, its velocity estimation accuracy and water stratification accuracy are determined by the pulse duration and the average time of phase difference, respectively.
The research on the performance and uncertainty of velocity estimation has always been the focus in the field of acoustic Doppler current measurement. In many studies, the measurement uncertainty guidelines [15], Monte Carlo simulation [16] and American Aerospace Society [17] methods were used to analyze the uncertainty of various types of ADCP measurements strictly. Aurélien Despax et al. studied the uncertainty of ADCP discharge measurement caused by cross-section selection and human operation by using repeated measures experiments and put forward some strategies to reduce the measurement uncertainty [18]. ADCP is usually installed on a moving carrier to complete the measurement. References [19,20] proposed a method to improve the velocity estimation accuracy by compensating the carrier motion. The bias error in moving-carrier ADCP current measurement can be separated into two classes: calibration error and application error; the major sources of calibration error and the influence of parameters on the discharge uncertainty were analyzed in Reference [21]. It is an efficient method to analyze and verify the estimation performance of ADCP through laboratory calibration system, which can simulate a variety of different measurement application environments in a short time [22]. At the same time, analyzing the estimation performance and error sources of ADCP through field comparison is also an important performance research method, especially for those new principle-based instruments [23]. Reference [24] uses fractional Fourier transform to separate the component of strong scatterers from ADCP echo, which improves the accuracy and stability of water velocity estimation. ADCP is widely used in current measurement, and there are many ways to improve the performance of velocity estimation, but there are few studies related to the optimization of transmitted waveform. Reference [25] points out that the measurement deviation of the phase-encoded signal is mainly determined by the energy of the autocorrelation function of a single baseband pulse (SBPAF) and gives the calculation method of the coding phase that makes the SBPAF energy maximum. Although this method can reduce the signal spectrum broadening caused by phase discontinuity, it also increases the complexity of waveform design, hardware implementation and signal processing. Reference [26] deduces the Cramer Rao bound of broadband Doppler sonar, points out that the factors affecting the accuracy of single beam velocity estimation are divided into SNR and echo randomness and gives the qualitative selection strategy of waveform parameters. Reference [27] proposes a waveform design method of dual band coherent phase-encoded transmitted signal by dividing the signal spectrum into two. This method enriches the form of the transmitted signal of ADCP, but its estimation performance is still limited by the lack of coherent method, and this method does not consider the adaptation to the application environment. To sum up, it is a relatively efficient strategy to adjust the velocity estimation accuracy and stratification accuracy by changing the code length and number of code repeats (hereinafter referred to as the number of repeats), and this strategy is easy to adapt to the application environment. However, improving the velocity estimation accuracy will reduce the velocity measurement range and water stratification accuracy, and improving the stratification accuracy will also degrade the velocity estimation accuracy. The application environment of acoustic Doppler Sensors 2021, 21, 3768 3 of 21 measurement is complex and diverse [28][29][30]. Due to the lack of specific strategies to adjust the waveform parameters according to the characteristics of the application environment, it is difficult to optimize the measurement performance of ADCP. In order to solve this problem, this paper introduces the constraints of sonar application environment, and researches the environment-adaptable waveform optimization methods of broadband ADCP based on the ambiguity function and nonlinear programming method.
The constraints of application environment on the waveform parameters of broadband ADCP are mainly reflected in the following aspects. In theory, the longer the code length and the more repeats, the higher the velocity estimation accuracy. In practice, it is necessary to consider the limitation of current velocity range, stratification accuracy requirement and other factors, and it is not advisable to blindly increase the code length and number of repeats. The velocity measurement range of the broadband method is inversely proportional to the code length and chip width (limited by the bandwidth of the transducer and receiver, the chip width is generally fixedly set to the reciprocal of the system bandwidth) [14]. In order to prevent the maximum current velocity from exceeding the velocity measurement range, there is an upper limit on code length. The layer thickness of the broadband method is proportional to the product of code length and number of repeats. The upper limit of the layer thickness requirement determines the upper limit of the product. The velocity estimation accuracy of the broadband method is inversely proportional to code length and number of repeats [13]. The lower limit of the velocity estimation accuracy requirement determines the lower limit of the product of code length and number of repeats. In order to ensure the output signal-to-noise ratio (SNR) of the matched filter of the receiver and the peak sidelobe level (PSL) of the autocorrelation function of the coded signal, there is a lower limit for code length. Therefore, the waveform design of the broadband ADCP must consider the above constraints of application environment.
In order to solve the problem that the existing waveform design methods lack the pertinence of application of environment characteristics and are difficult to realize the optimization of estimation performance, this paper introduces the constraints of the application environment, and uses an ambiguity function and nonlinear programming method to design two dynamic configuration methods of waveform parameters with environmental adaptability and optimal estimation performance. Section 2 introduces the basic operating principle of broadband ADCP. Section 3 employs the ambiguity function tool to broadband ADCP for the first time, studies the relationship between the velocity estimation accuracy of phase-encoded signal and its waveform parameters by using the ambiguity function, and gives the performance cost functions of velocity estimation and stratification. Then, aiming at the measurement application that prioritizes high velocity estimation accuracy or high stratification accuracy, a velocity estimation accuracy priority (VEAP) or a water stratification accuracy priority (WSAP) waveform optimization method are proposed. The constraints of the current velocity range, velocity estimation accuracy requirements and stratification accuracy requirements on the waveform parameters of the phase-encoded signal are studied in detail. Both methods model the design task of optimal waveform as a nonlinear programming problem with multiple inequality constraints and obtain waveform parameters that satisfy the constraints of the application environment and make the velocity estimation accuracy or stratification accuracy optimal. In Section 4, a principal prototype was built using an onshore ADCP calibration system. The experimental results demonstrate the correctness and practicability of the waveform design methods proposed by this paper. Finally, Section 5 highlights the conclusions of this study.

Principle of Broadband ADCP
Broadband Doppler sonar usually uses a phase-encoded pulse as the transmitted signal for achieving higher range and velocity accuracy. The water profile is divided into multiple adjacent layers according to the length of the transmitted signal. First, the radial Doppler frequency shift of single beam in a single layer is estimated, and then the velocity profile in the geodetic coordinate system is obtained by combining the attitude sensing data Sensors 2021, 21, 3768 4 of 21 and the radial Doppler frequency shift of multiple beams. The common velocity estimation method of broadband Doppler sonar is complex covariance algorithm, which estimates the average phase change of two sub-echoes in a single layer to get the radial Doppler shift of the layer. The time delay of the two sub-echoes is called as correlation delay, and the duration of the single sub-echo is named as average time of phase difference. The current velocity estimation principle of broadband Doppler sonar is shown in Figure 1, where the code length and number of repeats are set as 7 and 3, respectively. The code selection criterion in the broadband method is that the main autocorrelation peak should be as sharp as possible, and the secondary autocorrelation peak should be as low as possible. There are some popular codes, such as Barker codes, minimum peak sidelobe (MPS) codes and pseudo-random noise sequences.

Principle of Broadband ADCP
Broadband Doppler sonar usually uses a phase-encoded pulse as the transmitted signal for achieving higher range and velocity accuracy. The water profile is divided into multiple adjacent layers according to the length of the transmitted signal. First, the radial Doppler frequency shift of single beam in a single layer is estimated, and then the velocity profile in the geodetic coordinate system is obtained by combining the attitude sensing data and the radial Doppler frequency shift of multiple beams. The common velocity estimation method of broadband Doppler sonar is complex covariance algorithm, which estimates the average phase change of two sub-echoes in a single layer to get the radial Doppler shift of the layer. The time delay of the two sub-echoes is called as correlation delay, and the duration of the single sub-echo is named as average time of phase difference. The current velocity estimation principle of broadband Doppler sonar is shown in Figure 1, where the code length and number of repeats are set as 7 and 3, respectively. The code selection criterion in the broadband method is that the main autocorrelation peak should be as sharp as possible, and the secondary autocorrelation peak should be as low as possible. There are some popular codes, such as Barker codes, minimum peak sidelobe (MPS) codes and pseudo-random noise sequences. In Figure 1, td represents the correlation delay, ta is the average time of the phase difference and the calculation formula of water layer thickness (oblique thickness along the beam direction) Δz is: where c is the underwater sound velocity. When the total chip number of the transmitted pulse is fixed, the more symbols ta contains, the smaller the standard deviation of velocity estimation [25]. Therefore, the ta is usually configured equal to the time width of N−1 repetitive coding and td is equal to the time width of single coding, namely: where N is the number of repeats, L is the code length and τ is the chip width. Substituting ta in Equation (2) into (1), we get the following result: In Figure 1, t d represents the correlation delay, t a is the average time of the phase difference and the calculation formula of water layer thickness (oblique thickness along the beam direction) ∆z is: where c is the underwater sound velocity. When the total chip number of the transmitted pulse is fixed, the more symbols t a contains, the smaller the standard deviation of velocity estimation [25]. Therefore, the t a is usually configured equal to the time width of N − 1 repetitive coding and t d is equal to the time width of single coding, namely: where N is the number of repeats, L is the code length and τ is the chip width. Substituting t a in Equation (2) into (1), we get the following result: The absolute value of velocity measurement range is called ambiguity velocity. In the broadband method, the ambiguity velocity is inversely proportional to the correlation delay. Assuming the phase change range of the two sub-echoes is [−π, π], the ambiguity velocity v max is calculated as follows [13]: where λ is the wavelength and f dmax is the maximum Doppler shift. According to the complex covariance algorithm [13], the Doppler angular shift can be calculated by the change rate of echo phase within the correlation delay t d , so the maximum Doppler angular frequency 2πf dmax is equal to the ratio of the maximum phase change π that may occur within t d to the time interval t d . Finally, f dmax = 1/(2t d ).

Relationship Between Velocity Estimation Accuracy and Waveform Parameters
The ambiguity function is a tool used in the radar field to analyze the effects of time delay and Doppler frequency shift on radar echoes. In view of the similarity between Doppler sonar and radar systems, it is also feasible to apply ambiguity functions to the field of sonar. By studying the ambiguity function of the transmitted waveform, the relationship between the waveform parameters and the resolution, ambiguity and estimation accuracy of the Doppler sonar can be determined. The definition of the complex ambiguity function is as follows: where t is the time delay relative to the expected matched filter peak output, f is the Doppler shift between the actual echo and the transmitted signal, x(t) is usually the complex envelope of the signal and the higher the value of χ(t, f ), the stronger the resolution of the system on this (t, f ). Ambiguity function is usually defined as the amplitude function of complex ambiguity function, namely |χ(t, f )|. The ambiguity function of the phase-encoded signal is based on the ambiguity function of the rectangular pulse signal. It is assumed that the complex envelope of the rectangular pulse is defined as follows: x r (t) = 1/ √ τ r , |t| ≤ τ r /2 0, |t| > τ r /2 (6) where τ r is the pulse width and the amplitude is normalized. The expression of x r (t) complex ambiguity function is as follows [31]: Calculate the amplitude of the complex ambiguity function and regard the result as the ambiguity function expression of the rectangular pulse, which is shown in Equation (8): The complex envelope of repeated binary phase-encoded signal can be expressed as: where θ n is equal to 0 or π, represents the polarity of biphase code, L is code length, τ is chip width and N is the number of repeats. By substituting Equation (9) into the definition formula of complex ambiguity function, the following results are obtained: where χ r (t, f ) is the complex ambiguity function of the rectangular pulse, the expression is as Equation (7) and τ r equals τ. Conducting double summation decomposition on Equation (10), the complex ambiguity function of N times repeated L-bit coded pulse is: Since the range of χ r (t, f ) on the t-axis is |t| ≤ τ, χ r (t, f ) with different delays in Equation (11) will not be aliased, and the ambiguity function |χ r (t, f )| is equal to the sum of the amplitudes of all the summation terms in Equation (11). Equation (11) shows that the ambiguity function of the phase-encoded signal is superimposed by rectangular pulse ambiguity functions with different delays after amplitude weighting, and the weighting coefficient is related to the autocorrelation characteristics of the encoding. The distance and speed resolution performance of |χ r (t, f )| is mainly determined by its main peak. Taking n = 0 in Equation (11) to obtain the expression of the ambiguity function in the main peak range is: The ambiguity function of the 13-bit Barker code phase modulation signal is shown in Figure 2, and the contour plot of Figure 2a is shown in Figure 2b. Figure 2 is the graph of the square of the modulus of Equation (12) when N = 1 and L = 13. The amplitude, frequency shift and delay in the figure are all normalized, and T is the transmitted pulse width (T equals NLτ). Based on the ambiguity function, the relationship between the velocity estimation accuracy and the waveform parameters such as code length and number of repeats can be concluded, and the cost function of velocity estimation performance can be estimated, which provides a theoretical basis for the optimal waveform design. Let t = 0 in Equation (12) to obtain the velocity ambiguity function of the phase-encoded signal: Based on the ambiguity function, the relationship between the velocity estimation accuracy and the waveform parameters such as code length and number of repeats can be concluded, and the cost function of velocity estimation performance can be estimated, which provides a theoretical basis for the optimal waveform design. Let t = 0 in Equation (12) to obtain the velocity ambiguity function of the phase-encoded signal: It can be seen from Equation (13) that the greater the product of code length and the number of repeats, the greater the PSL of the sin(πNLτf )/(πτf ) structure, the more obvious the "pushpin" shape of the velocity ambiguity graph, and the less the influence of the multi-valued ambiguity on the velocity estimation results.
The accuracy of frequency estimation is determined by the second derivative of the velocity ambiguity function |χ(0, f )| [32], and the relational expression is shown in Equation (14). Equation (14) means the sharpness of the frequency dimension of the ambiguity function at point (0,0). The smaller the value of σ 2 f , the sharper the peak, and the higher the accuracy of velocity estimation. The velocity estimation accuracy can be greatly improved through reasonable configuration of the waveform parameters.
where E is the energy of signal received by receiver, E n is the noise energy per unit bandwidth of receiver, E n = kT k F n , k is the Boltzmann constant (k = 1.38 × 10 −23 J/K), T k is the thermodynamic temperature corresponding to 17 • C and F n is the noise coefficient of receiver. Substituting Equation (13) into (14), the following results are obtained: The standard deviation of velocity estimation is obtained by substituting the relationship between Doppler shift and velocity v = λf /2 into Equation (15): Equation (16) shows that the standard deviation of velocity estimation is inversely proportional to N 3/2 (the 3/2 power of number of repeats), L 3/2 (the 3/2 power of code length) and τ (the chip width). That is to say, the wider the chip, the longer the code, and the more repeats, the higher the accuracy of velocity estimation. In practice, the design of the sonar waveform needs to consider the actual current characteristics of measurement applications. In the following, specific phase-encoded waveform design methods are proposed for two common applications with high velocity estimation accuracy requirements and high-water stratification accuracy requirements.

Velocity Estimation Accuracy Priority (VEAP) Waveform Optimization Method
For the current measurement application that prioritizes high velocity estimation accuracy, priority is given to achieve high velocity accuracy. In this section, a configuration method of code length and number of repeats is proposed, which can satisfy the constraints of stratification accuracy and velocity estimation range and optimize the velocity estimation accuracy. From the analysis in Section 3.1, we can see that the best velocity estimation accuracy can be obtained when we maximize the code length, the number of repeats and chip width of the transmitted signal. Limited by the bandwidth of the transducer and receiver, the chip width is generally fixed as the inverse of the system bandwidth. Extending the code length will reduce the ambiguity velocity, which may cause the current velocity exceeds the velocity estimation range of the sonar; increasing the number of repeats will increase the layer thickness, resulting in a decrease in the ranging resolution. Therefore, while improving the velocity estimation accuracy, it is necessary to take into account the constraints, such as stratification accuracy requirement and velocity estimation range, and not blindly increase the code length and the number of repeats. It can be seen from the above analysis that the waveform design in this application is a nonlinear programming problem constrained by multiple inequalities. The solution is to first transform the cost function from two-variable function to a one-variable function, and then transform inequality constrained optimization into equality constrained optimization by introducing slack variables, and then introduce Lagrangian multipliers to transform equality constrained optimization into unconstrained optimization. The flow chart of the VEAP method is shown in Figure 3.  Based on the relationship between layer thickness and waveform parameters, the performance cost function of velocity estimation is transformed into a one-variable function. The calculation formula of the layer thickness Δz in the broadband method is shown in Equation (3), assuming g(Δz) = 2Δz/(cτ), Equation (3) is equivalent to: By substituting Equation (17) into (16), the cost function σv is transformed into a single variable function: Equation (19) shows that the cost function and the number of repeats are non-linear. In addition, when N is fixed, the larger Δz, the higher the accuracy of velocity estimation, so Δz takes the maximum allowable layer thickness Δzm.
The following is a detailed analysis of the constraints on the waveform parameters, such as current velocity range and stratification accuracy requirement. The current velocity range constraint makes N have a lower limit. The correlation delay is inversely proportional to the ambiguity velocity, and the maximum correlation delay must ensure that no velocity estimation ambiguity occurs, that is, the ambiguity velocity should be greater than the to be measured maximum current velocity. In addition, Doppler tolerance of the Barker code is poor [31]. In order to control the Doppler mismatch loss within 1dB, the Based on the relationship between layer thickness and waveform parameters, the performance cost function of velocity estimation is transformed into a one-variable function. The calculation formula of the layer thickness ∆z in the broadband method is shown in Equation (3), assuming g(∆z) = 2∆z/(cτ), Equation (3) is equivalent to: By substituting Equation (17) into (16), the cost function σ v is transformed into a single variable function: Equation (19) shows that the cost function and the number of repeats are non-linear. In addition, when N is fixed, the larger ∆z, the higher the accuracy of velocity estimation, so ∆z takes the maximum allowable layer thickness ∆z m .
The following is a detailed analysis of the constraints on the waveform parameters, such as current velocity range and stratification accuracy requirement. The current velocity range constraint makes N have a lower limit. The correlation delay is inversely proportional to the ambiguity velocity, and the maximum correlation delay must ensure that no velocity Sensors 2021, 21, 3768 9 of 21 estimation ambiguity occurs, that is, the ambiguity velocity should be greater than the to be measured maximum current velocity. In addition, Doppler tolerance of the Barker code is poor [31]. In order to control the Doppler mismatch loss within 1dB, the phase change caused by the maximum current velocity generally does not exceed π/2, that is, the ambiguity velocity should be greater than 2 times the maximum current velocity to be measured. Thus, the inequality between maximum radial velocity v m and correlation delay t d is obtained: Substituting the t d expression in Equation (2) into (19), we get the following result: where L mid = λ/(8τv m ) represents the upper limit of the code length determined by the current velocity range. Combining Equations (17) and (20) to obtain a lower limit of N: The characteristics of the broadband method require that the number of repeats is not less than 2: The lower limit of the PSL of code autocorrelation function makes L have a lower limit. Combining Equation (17) to obtain the constraint of the minimum allowable code length L min on N is: Then, three slack variables r 2 1 , r 2 2 and r 2 3 are introduced to transform inequality constraints into equality constraints. According to Equations (21)−(23) the transformation result is: In order to transform the equality constrained optimization issue into the unconstrained optimization, Lagrangian multipliers µ 1 , µ 2 and µ 3 need to be introduced. According to Equation (24), the Lagrangian function f (N,r 2 1 , r 2 2 , r 2 3 , µ 1 , µ 2 , µ 3 ) is: After calculating the partial derivative of the Lagrangian function with respect to the respective variables, an equation group after simplification is obtained: where the expression of ∂σ v /∂N is: (27) According to the relationship between the expression on the right side of Equations (21) and (2), the solution of equation group (27) can be divided into two cases: v m ∆z m ≥ cλ/16 and v m ∆z m < cλ/16. When v m ∆z m ≥ cλ/16, the value range of N is shown in Figure 4a. The thin solid lines in the figure represent curve cluster obtained by Equation (17), all ∆z with different values satisfies ∆z ≤ ∆z m . The thick solid line represents the feasible region of N. The cost function σ v takes the global optimal value at the red circle in Figure 4a. At this point, the expressions of L and N are: where the brackets indicate rounding. Equation (28) shows that when v m ∆z m ≥ cλ/16, the standard deviation of velocity estimation can be minimized by first determining the code length by the velocity estimation range constraint and then determining the number of repeats by the stratification accuracy constraint. The minimum σ v obtained in this case is: Sensors 2021, 21, x FOR PEER REVIEW 11 of 22 where the expression of ∂σv/∂N is: According to the relationship between the expression on the right side of Equation (21) and 2, the solution of equation group (27) can be divided into two cases: vmΔzm ≥ cλ/16 and vmΔzm < cλ/16. When vmΔzm ≥ cλ/16, the value range of N is shown in Figure 4a. The thin solid lines in the figure represent curve cluster obtained by Equation (17) where the brackets indicate rounding. Equation (28) shows that when vmΔzm ≥ cλ/16, the standard deviation of velocity estimation can be minimized by first determining the code length by the velocity estimation range constraint and then determining the number of repeats by the stratification accuracy constraint. The minimum σv obtained in this case is:  When vmΔzm < cλ/16, the feasible region of N is shown in Figure 4b. The cost function σv takes the global optimal value at the black circle in Figure 4b. At this point, the expressions of L and N are: where the square brackets indicate rounding. Equation (30) shows that when vmΔzm < cλ/16, the standard deviation of velocity estimation can be minimized by first setting the number of repeats to 2 and then determining code length by the stratification accuracy constraint. The minimum σv obtained in this case is: When v m ∆z m < cλ/16, the feasible region of N is shown in Figure 4b. The cost function σ v takes the global optimal value at the black circle in Figure 4b. At this point, the expressions of L and N are: where the square brackets indicate rounding. Equation (30) shows that when v m ∆z m < cλ/16, the standard deviation of velocity estimation can be minimized by first setting the number of repeats to 2 and then determining code length by the stratification accuracy constraint. The minimum σ v obtained in this case is: In summary, the longest code length and the smallest number of repeats in the feasible region can minimize the velocity estimation standard deviation in the VEAP method.

Water Stratification Accuracy Priority (WSAP) Waveform Optimization Method
For the current measurement application that prioritizes high stratification accuracy, priority is given to achieve small layer thickness. In this section, a configuration method of code length and number of repeats is proposed, which can satisfy the constraints of velocity estimation accuracy and velocity estimation range and optimize the stratification accuracy. It can be seen from Equation (3) that the thinnest layer thickness can be obtained by selecting the smallest number of repeats, code length and chip width. Limited by the bandwidth of the transducer and receiver, the chip width is generally fixed as the inverse of the system bandwidth. Shortening the code length will increase the ambiguity velocity, and the phase change caused by the same velocity will decrease, which will lead to the greater influence of phase noise; thus reducing the number of repeats will reduce the time bandwidth product of the echo segment from which the phase difference is estimated, which will lead to the decline of velocity estimation accuracy. Therefore, when increasing the hierarchical accuracy, it is necessary to consider constraints such as velocity estimation accuracy, velocity estimation range, etc., and cannot reduce the number of code length and repeats. From the above analysis, it can be seen that the waveform design in this case is also a nonlinear programming problem with multiple inequality constraints, but the cost function changes into the Equation (3). The solution is similar to that in Section 3.2 and will not be repeated here. The flow chart of the WSAP method is shown in Figure 5.
In summary, the longest code length and the smallest number of repeats in the feasible region can minimize the velocity estimation standard deviation in the VEAP method.

Water Stratification Accuracy Priority (WSAP) Waveform Optimization Method
For the current measurement application that prioritizes high stratification accuracy, priority is given to achieve small layer thickness. In this section, a configuration method of code length and number of repeats is proposed, which can satisfy the constraints of velocity estimation accuracy and velocity estimation range and optimize the stratification accuracy. It can be seen from Equation (3) that the thinnest layer thickness can be obtained by selecting the smallest number of repeats, code length and chip width. Limited by the bandwidth of the transducer and receiver, the chip width is generally fixed as the inverse of the system bandwidth. Shortening the code length will increase the ambiguity velocity, and the phase change caused by the same velocity will decrease, which will lead to the greater influence of phase noise; thus reducing the number of repeats will reduce the time bandwidth product of the echo segment from which the phase difference is estimated, which will lead to the decline of velocity estimation accuracy. Therefore, when increasing the hierarchical accuracy, it is necessary to consider constraints such as velocity estimation accuracy, velocity estimation range, etc., and cannot reduce the number of code length and repeats. From the above analysis, it can be seen that the waveform design in this case is also a nonlinear programming problem with multiple inequality constraints, but the cost function changes into the Equation (3). The solution is similar to that in Section 3.2 and will not be repeated here. The flow chart of the WSAP method is shown in Figure 5.  Based on the relationship between velocity estimation accuracy and waveform parameters, the cost function of stratification accuracy is transformed into a one-variable function. Assuming that the maximum allowable velocity standard deviation is σvm, the velocity estimation accuracy constraint can be expressed as: Based on the relationship between velocity estimation accuracy and waveform parameters, the cost function of stratification accuracy is transformed into a one-variable function.
Assuming that the maximum allowable velocity standard deviation is σ vm , the velocity estimation accuracy constraint can be expressed as: The Equation (32) is transformed into: where y(σ v ) is: By substituting Equation (33) into (3), the cost function ∆z is transformed into a single variable function: Equation (35) shows that the cost function and the number of repeats are non-linear. In addition, when N is fixed, the larger σ v , the higher the stratification accuracy, so σ v takes the maximum allowable standard deviation σ vm .
The following is a detailed analysis of the constraints on the waveform parameters, such as current velocity range and velocity estimation accuracy requirement. The restriction of current velocity range on N has been studied in detail in Section 3.2, and only the final expression is given: The characteristics of the broadband method require that the number of repeats is not less than 2: N ≥ 2 (37) The lower limit of the PSL of code autocorrelation function makes L have a lower limit. Combining Equation (32) to obtain the constraint of the minimum allowable code length L min on N is: The method of transforming inequality constrained optimization into unconstrained optimization is similar to that in Section 3.2. The equations obtained by partial derivation of Lagrange function are given directly: where the expression of ∂∆z/∂N is: According to the relationship between the expression on the right side of Equations (36) and (2), the solution of equation group (39) can be divided into two cases: y(σ vm ) ≥ 2L mid and y(σ vm ) < 2L mid . When y(σ vm ) ≥ 2L mid , the value range of N is shown in Figure 6a. The thin solid lines in the figure represent curve cluster obtained by Equation (33), all σ v with different values satisfies σ v ≤ σ vm . The thick solid line represents the feasible region of N. The cost function ∆z takes the global optimal value at the black circle in Figure 6a. At this point, the expressions of L and N are: where the brackets indicate rounding. Equation (41) shows that when y(σ vm ) ≥ 2L mid , the layer thickness can be minimized by first determining the code length by the velocity estimation range constraint and then determining the number of repeats by the velocity estimation accuracy constraint. The minimum layer thickness ∆z min obtained in this case is: According to the relationship between the expression on the right side of Equation (36) and 2, the solution of equation group (39) can be divided into two cases: y(σvm) ≥ 2Lmid and y(σvm) < 2Lmid. When y(σvm) ≥ 2Lmid, the value range of N is shown in Figure 6a. The thin solid lines in the figure represent curve cluster obtained by Equation (33), all σv with different values satisfies σv ≤ σvm. The thick solid line represents the feasible region of N. The cost function Δz takes the global optimal value at the black circle in Figure 6a. At this point, the expressions of L and N are: where the brackets indicate rounding. Equation (41) shows that when y(σvm) ≥ 2Lmid, the layer thickness can be minimized by first determining the code length by the velocity estimation range constraint and then determining the number of repeats by the velocity estimation accuracy constraint. The minimum layer thickness Δzmin obtained in this case is: When y(σvm) < 2Lmid, the feasible region of N is shown in Figure 6b. The cost function Δz takes the global optimal value at the black circle in Figure 6b. At this point, the expressions of L and N are: where the square brackets indicate rounding. Equation (43) shows that when y(σvm) < 2Lmid, the layer thickness can be minimized by first setting the number of repeats to 2 and then determining code length by the velocity estimation accuracy constraint. The minimum layer thickness Δzmin obtained in this case is: In summary, the longest code length and the smallest number of repeats in the feasible region can minimize the layer thickness in the WSAP method. When y(σ vm ) < 2L mid , the feasible region of N is shown in Figure 6b. The cost function ∆z takes the global optimal value at the black circle in Figure 6b. At this point, the expressions of L and N are: where the square brackets indicate rounding. Equation (43) shows that when y(σ vm ) < 2L mid , the layer thickness can be minimized by first setting the number of repeats to 2 and then determining code length by the velocity estimation accuracy constraint. The minimum layer thickness ∆z min obtained in this case is: In summary, the longest code length and the smallest number of repeats in the feasible region can minimize the layer thickness in the WSAP method.

Experiments and Verification
An experimental platform based on an onshore calibration system for ADCP is built to verify the effectiveness of the proposed waveform design methods.

Overview of Experimental Platform
The onshore calibration system can generate multi-channel simulated current measurement echo according to the set parameters and simulate underwater acoustic current measurement through seamless docking with the transducers of current measurement instruments in the laboratory. The onshore calibration system can flexibly set the transmitted waveform (including noncoherent pulse, coherent pulses train and phase-encoded pulse), water depth, scatterer distribution, current distribution, environmental noise etc., which provides convenience for the study of acoustic current measurement methods and the test of current measurement instruments. The shore calibration system consists of a computer (including acoustic Doppler echo simulation software), an arbitrary wave generator, a power amplifier, a transducer and a current measuring instrument (taking the 600 kHz ADCP as an example). The workflow chart of the onshore calibration system is shown in Figure 7a Table 1.

Experiments and Verification
An experimental platform based on an onshore calibration system for ADCP is built to verify the effectiveness of the proposed waveform design methods.

Overview of Experimental Platform
The onshore calibration system can generate multi-channel simulated current measurement echo according to the set parameters and simulate underwater acoustic current measurement through seamless docking with the transducers of current measurement instruments in the laboratory. The onshore calibration system can flexibly set the transmitted waveform (including noncoherent pulse, coherent pulses train and phase-encoded pulse), water depth, scatterer distribution, current distribution, environmental noise etc., which provides convenience for the study of acoustic current measurement methods and the test of current measurement instruments. The shore calibration system consists of a computer (including acoustic Doppler echo simulation software), an arbitrary wave generator, a power amplifier, a transducer and a current measuring instrument (taking the 600 kHz ADCP as an example). The workflow chart of the onshore calibration system is shown in Figure 7a. The echo simulation software running on the computer is responsible for generating the digital echoes according to the set parameters. The screenshot of the software interface is shown in the upper left corner of Figure 7b. The arbitrary waveform generator A and B convert the digital echoes into analog signals, and their actual photos are shown in the upper right corner of Figure 7b. Each arbitrary waveform generator can generate two channels of analog signals at the same time. The power amplifier is responsible for amplifying the analog signals to an appropriate power level and matching appropriate interface impedances for all the transducers. The photo of the power amplifier is shown in the lower right corner of Figure 7b. Finally, four transducers connected the output of the power amplifier are docked with transducers of an ADCP one by one. The photo of four groups of connected transducers is shown in the lower left corner of Figure  7b. The main parameters of the experimental platform are shown in Table 1.

VEAP Waveform Optimization Experiment Result
This experiment is aimed to compare the velocity estimation performance of the waveform obtained by the proposed VEAP method and other waveforms under the constraints of stratification accuracy requirement and current velocity range. This experiment contains two parts: a low velocity case and a high velocity case.

Low Velocity Case
The onshore calibration system is set as follows: the current velocity of all layers is set to 0.37 m/s, so the maximum current velocity v m equals 0.37 m/s, and the maximum allowable layer thickness ∆z m is set to 0.8 m. The default code is Barker code, and the MPS code is used when there is no corresponding length Barker code. The four codes selected in this experiment are shown in Table 2 [33]. Three waveforms which are different from the waveform obtained by the VEAP method in this paper are selected, and the standard deviation of velocity estimation of the four waveforms under eight SNR (5-40 dB) is compared. Under each SNR condition, the sample number of velocity estimation is 500. The process of calculating the waveform parameter of VEAP method is described as follows: the chip width is 20 µs obtained from the system bandwidth, and the chip width is fixed. According to the value of ∆z m and v m , this measurement application belongs to the case of v m ∆z m ≥ cλ/16. According to Equation (28) The standard deviation of velocity estimation realized by the four waveforms is shown in Figure 9. The dotted line with star icon indicates the velocity standard deviation of the waveform (N = 2, L = 42) obtained by the VEAP method in this paper. The triangle icon indicates the deviation of (N = 3, L = 23) waveform. The diamond icon indicates the deviation of (N = 5, L = 13) waveform. The square icon indicates the deviation of (N = 8, L = 7) waveform. It can be seen from Figure 9 that the velocity standard deviation of four waveforms decreases with the increase of SNR. When SNR is greater than 25 dB, the velocity standard deviation tends to be stable, and the stable values ordered from largest to smallest are 0.22 cm/s, 0.14 cm/s, 0.1 cm/s and 0.05 cm/s. Under all SNR conditions, the velocity estimation accuracy of the waveform designed by the proposed VEAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method. The standard deviation of velocity estimation realized by the four waveforms is shown in Figure 9. The dotted line with star icon indicates the velocity standard deviation of the waveform (N = 2, L = 42) obtained by the VEAP method in this paper. The triangle icon indicates the deviation of (N = 3, L = 23) waveform. The diamond icon indicates the deviation of (N = 5, L = 13) waveform. The square icon indicates the deviation of (N = 8, L = 7) waveform. It can be seen from Figure 9 that the velocity standard deviation of four waveforms decreases with the increase of SNR. When SNR is greater than 25 dB, the velocity standard deviation tends to be stable, and the stable values ordered from largest to smallest are 0.22 cm/s, 0.14 cm/s, 0.1 cm/s and 0.05 cm/s. Under all SNR conditions, the velocity estimation accuracy of the waveform designed by the proposed VEAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method. The standard deviation of velocity estimation realized by the four waveforms is shown in Figure 9. The dotted line with star icon indicates the velocity standard deviation of the waveform (N = 2, L = 42) obtained by the VEAP method in this paper. The triangle icon indicates the deviation of (N = 3, L = 23) waveform. The diamond icon indicates the deviation of (N = 5, L = 13) waveform. The square icon indicates the deviation of (N = 8, L = 7) waveform. It can be seen from Figure 9 that the velocity standard deviation of four waveforms decreases with the increase of SNR. When SNR is greater than 25 dB, the velocity standard deviation tends to be stable, and the stable values ordered from largest to smallest are 0.22 cm/s, 0.14 cm/s, 0.1 cm/s and 0.05 cm/s. Under all SNR conditions, the velocity estimation accuracy of the waveform designed by the proposed VEAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method.

High Velocity Case
The onshore calibration system is set as follows: the current velocity of all layers is set to 2 m/s, so the maximum current velocity v m equals 2 m/s, and the maximum allowable layer thickness ∆z m is set to 0.8 m. According to Equation (28), the code length and number of repeats that minimize the velocity standard deviation are N = 8 and L = 7. In addition, other three waveforms are selected in the feasible region, and their waveform parameters are (N = 9, L = 6), (N = 10, L = 5) and (N = 12, L = 4).
The standard deviation of velocity estimation realized by the four waveforms is shown in Figure 10. The dotted line with star icon indicates the velocity standard deviation of the waveform (N = 8, L = 7) obtained by the proposed VEAP method. The triangle icon indicates the deviation of (N = 9, L = 6) waveform. The diamond icon indicates the deviation of (N = 10, L = 5) waveform. The square icon indicates the deviation of (N = 12, L = 4) waveform. Under all SNR conditions, the velocity estimation accuracy of the waveform designed by the proposed VEAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method.
ble layer thickness Δzm is set to 0.8 m. According to Equation (28), th number of repeats that minimize the velocity standard deviation are addition, other three waveforms are selected in the feasible region, an parameters are (N = 9, L = 6), (N = 10, L = 5) and (N = 12, L = 4).
The standard deviation of velocity estimation realized by the f shown in Figure 10. The dotted line with star icon indicates the veloci tion of the waveform (N = 8, L = 7) obtained by the proposed VEAP me icon indicates the deviation of (N = 9, L = 6) waveform. The diamond deviation of (N = 10, L = 5) waveform. The square icon indicates the de L = 4) waveform. Under all SNR conditions, the velocity estimation acc form designed by the proposed VEAP method is better than that of the which proves the effectiveness of the proposed method.

WSAP Waveform Optimization Experiment Result
This experiment aimed to compare the stratification performanc obtained by the proposed WSAP method and other waveforms under velocity estimation accuracy requirement and current velocity range. T tion system is set as follows: In order to simulate shear current appli velocity is set to be linearly distributed with depth. The velocity of top the velocity of bottom layer is 0.1 m/s, so the maximum current veloc is vm = 1 m/s. The maximum allowable velocity standard deviation is default code is Barker code, and the MPS code is used when there is length Barker code. The four codes selected in this experiment are show

WSAP Waveform Optimization Experiment Result
This experiment aimed to compare the stratification performance of the waveform obtained by the proposed WSAP method and other waveforms under the constraints of velocity estimation accuracy requirement and current velocity range. The onshore calibration system is set as follows: In order to simulate shear current application, the current velocity is set to be linearly distributed with depth. The velocity of top layer is 1 m/s, and the velocity of bottom layer is 0.1 m/s, so the maximum current velocity to be measured is v m = 1 m/s. The maximum allowable velocity standard deviation is σ vm = 0.05 m/s. The default code is Barker code, and the MPS code is used when there is no corresponding length Barker code. The four codes selected in this experiment are shown in Table 3 [33]. Three waveforms that are different from the waveform obtained by the WSAP method in this paper are selected, and the layer thickness of the four waveforms under 20 dB SNR is compared. The sample number of velocity estimation is 500. The process of calculating the waveform parameter of WSAP method is described as follows: The chip width is 20 µs obtained from the system bandwidth. y(σ vm ) = 16.8 is obtained from the maximum allowable velocity estimation error σ vm , and L mid = 15.6 is obtained from the maximum current velocity v m . Therefore, this measurement application belongs to the case of y(σ vm ) < 2L mid . According to Equation (43), the code length and number of repeats that minimize the layer thickness is N = 2 and L = 9. At the same time, the correlation delay is  Three waveforms that are different from the waveform obtained by the WSAP method in this paper are selected, and the layer thickness of the four waveforms under 20 dB SNR is compared. The sample number of velocity estimation is 500. The process of calculating the waveform parameter of WSAP method is described as follows: The chip width is 20 μs obtained from the system bandwidth. y(σvm) = 16.8 is obtained from the maximum allowable velocity estimation error σvm, and Lmid = 15.6 is obtained from the maximum current velocity vm. Therefore, this measurement application belongs to the case of y(σvm) < 2Lmid. According to Equation (43), the code length and number of repeats that minimize the layer thickness is N = 2 and L = 9. At the same time, the correlation delay is td = 0. 18  The layer thickness obtained by the four waveforms is shown in Figure 12. The black vertical bar indicates the layer thickness of the waveform (N = 2, L = 9) obtained by the WSAP method in this paper. The gray vertical bars indicate the layer thickness of the other waveforms (N = 3, L = 6), (N = 4, L = 5) and (N = 5, L = 4), respectively. It can be seen from Figure 12 that the layer thickness of the four waveforms decrease from left to right, and the values are 0.24 m, 0.23 m, 0.18 m and 0.14 m, respectively. The stratification accuracy of the waveform designed by the proposed WSAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method. The layer thickness obtained by the four waveforms is shown in Figure 12. The black vertical bar indicates the layer thickness of the waveform (N = 2, L = 9) obtained by the WSAP method in this paper. The gray vertical bars indicate the layer thickness of the other waveforms (N = 3, L = 6), (N = 4, L = 5) and (N = 5, L = 4), respectively. It can be seen from Figure 12 that the layer thickness of the four waveforms decrease from left to right, and the values are 0.24 m, 0.23 m, 0.18 m and 0.14 m, respectively. The stratification accuracy of the waveform designed by the proposed WSAP method is better than that of the other waveforms, which proves the effectiveness of the proposed method. Sensors 2021, 21, x FOR PEER REVIEW 20 of 22 Figure 12. Comparison of the layer thickness between the waveform obtained by the proposed WSAP method and the other waveforms. The waveform (N = 2, L = 9) obtained by the proposed WSAP method achieves the smallest layer thickness.

Conclusions
In order to overcome the shortcomings of the existing waveform design methods, such as lack of application environment pertinence and difficulties in optimizing the estimation performance, this paper designs two waveform optimization methods with environmental adaptability and optimal estimation performance. In these two methods, the environment characteristics such as current velocity range and distribution are taken as the constraints of the nonlinear optimization problems to realize the velocity estimation or stratification performance optimization.
The proposed methods improve the velocity estimation accuracy, stratification accuracy and environmental adaptability of broadband ADCP. Evaluation of algorithm performance in comparison with other existing methods indicate that the velocity estimation accuracy of the proposed VEAP method is improved by at least 50%, and the stratification accuracy of the proposed WSAP method is improved by at least 22%, which demonstrates the proposed methods are accurate and effective.
Obtaining high-precision velocity estimation and water stratification is the primary purpose of most ADCP measurement applications, so the proposed methods in this paper are valuable for many measurement activities. The waveform parameters can be dynamically configured according to the characteristics of the measurement environment by using the proposed methods, then higher accuracy and flexibility can be obtained than the fixed parameter measurement method.
However, there may be some limitations in practical applications. For example, it is sometimes not easy to obtain accurate characteristics of measurement environment. The environmental characteristics obtained by experience and rough estimation sometimes have certain errors. Therefore, research on obtaining accurate environmental characteristics will be interesting and valuable.

Conclusions
In order to overcome the shortcomings of the existing waveform design methods, such as lack of application environment pertinence and difficulties in optimizing the estimation performance, this paper designs two waveform optimization methods with environmental adaptability and optimal estimation performance. In these two methods, the environment characteristics such as current velocity range and distribution are taken as the constraints of the nonlinear optimization problems to realize the velocity estimation or stratification performance optimization.
The proposed methods improve the velocity estimation accuracy, stratification accuracy and environmental adaptability of broadband ADCP. Evaluation of algorithm performance in comparison with other existing methods indicate that the velocity estimation accuracy of the proposed VEAP method is improved by at least 50%, and the stratification accuracy of the proposed WSAP method is improved by at least 22%, which demonstrates the proposed methods are accurate and effective.
Obtaining high-precision velocity estimation and water stratification is the primary purpose of most ADCP measurement applications, so the proposed methods in this paper are valuable for many measurement activities. The waveform parameters can be dynamically configured according to the characteristics of the measurement environment by using the proposed methods, then higher accuracy and flexibility can be obtained than the fixed parameter measurement method.
However, there may be some limitations in practical applications. For example, it is sometimes not easy to obtain accurate characteristics of measurement environment. The environmental characteristics obtained by experience and rough estimation sometimes have certain errors. Therefore, research on obtaining accurate environmental characteristics will be interesting and valuable.