Optimized Doppler Estimation and Symbol Synchronization for Mobile M-ary Spread Spectrum Underwater Acoustic Communication

In mobile underwater acoustic (UWA) communications, the Doppler effect causes severe signal distortion, which leads to carrier frequency shift and compresses/broadens the signal length. This situation has a more severe impact on communication performance in the case of low signal-to-noise ratio and variable-speed movement. This paper proposes a non-data-aided Doppler estimation method for M-ary spread spectrum UWA communication systems in mobile scenarios. The receiver uses the spread spectrum codes dedicated to transmitting signals with different frequency offsets as local reference signals. Correlation operations are performed symbol by symbol with the received signal. The decoding and Doppler estimation of the present symbol are achieved by searching the correlation maximum in the code domain and frequency domain. The length of the current symbol is corrected for the next symbol synchronization using the estimated Doppler coefficient. To optimize the process of Doppler estimation and symbol synchronization, a heuristic search method is used. By adjusting the Doppler factor search step size, setting the threshold value, and using the Doppler factor estimation of the previous symbol, the search range can be significantly reduced and the computational complexity decreased. The Fisher-Yates shuffle algorithm is used to traverse the search range to ensure reliability of the results. Simulation results show that enlarging the frequency-domain search step size in some degree does not affect the decoding accuracy. On 15 May 2021, a shallow-water mobile UWA spread spectrum communication experiment was conducted in Weihai, China. The horizontal distance between the transmitter and the receiver is 3.7–4.0 km, and the communication rate is 41.96 bits per second. The transmitting ship moves at a speed of 0–3 m/s, and the bit error rate (BER) is lower than 1e−3, which is better than that of the sliding correlation despreading method with average Doppler compensation.


Introduction
The underwater acoustic (UWA) channel is complex. When the source and receiving devices are in relative motion, the influence of Doppler and noise on the communication quality is severe [1][2][3]. The use of a spread spectrum communication system suppresses the multi-path effect, to a certain extent, and brings a spreading gain under the condition of a low signal-to-noise ratio (SNR) [4,5]. Nevertheless, for the UWA signal, the Doppler produces multiple effects of carrier frequency offset and compressed/expanded signal length [6,7]. Some researchers estimate Doppler based on the time-domain broadening/compression caused by the Doppler effect [8,9], which has a strong anti-noise ability. However, only the average Doppler of data block is obtained by the block estimation. For the autonomous undersea vehicle/unmanned undersea vehicle (AUV/UUV) that needs to realize reliable UWA communication while moving at a non-uniform velocity, the constant-velocity assumption is no longer applicable [10][11][12][13][14]. In this application scenario, to ensure the communication's robustness, Doppler estimation needs to be performed periodically [15]. Carrier frequency variation and synchronization error accumulation must be considered in the receiving process [16].
The research into mobile UWA communication mainly focuses on the Doppler estimation and equalization algorithm. Single carrier multiple phase shift keying (SC-MPSK) [17] and orthogonal frequency division multiplexing (OFDM) [18,19] are generally studied. Both have a weak anti-Doppler capability. Accurate symbol synchronization and channel equalization are the keys to achieving high-quality SC-MPSK communication; OFDM needs to conduct Doppler estimation and compensation with different accuracy step by step [20,21]. Considering the specific application limitations of AUV, we need to achieve reliable communication performance with the smallest hardware overhead possible. The spread spectrum does not have the problem of a high peak-to-average ratio in the OFDM and is not limited by the inter-symbol interference, as the SC-MPSK is. Without much demand on the communication rate, we believe that the spread spectrum technique is an appropriate choice. In [22], code modulation is used in the spread spectrum system so that the mapping bit does not reduce the Euclidean distance between transmission bits as it did in the traditional modulation, resulting in improved BER performance. The scheme proposed in this paper is based on M-ary spread spectrum (MSS) code modulation and combines symbol-by-symbol Doppler estimation and synchronization to meet the needs of mobile UWA communication. In addition, the MSS has a higher data rate and is more secure than the conventional direct sequence spread spectrum (DSSS), making it more practical [3,23].
To deal with the time-varying characteristics of the mobile UWA channel, the receiver needs to estimate the channel state information. The commonly employed Doppler estimation method is to use the thumb-tack ambiguity function of Doppler sensitive waveforms to calculate the cross-ambiguity function (CAF) on the two-dimensional grid of channel delay and Doppler compensation factor [24,25]. The maximum amplitude position of CAF correlation on Doppler grid is the estimated Doppler factor. However, this method has two obvious defects. One is a large amount of calculation. The other is the estimation accuracy limited by the time-bandwidth product of Doppler sensitive waveforms [26]. Inspired by [27], on the premise synchronization, we set up a group of discrete correlators. With different Doppler factors in the transmitted symbol as reference, we use these correlators to search in the expected Doppler range, obtaining the Doppler estimation. Simultaneously, it is assumed that the multi-path delay in the UWA channel is less than the duration of one symbol spreading code, ignoring the impact of inter-symbol interference on the propagated signal, and focusing on the signal distortion caused by Doppler [15,26]. In the spread spectrum system, the spreading sequences carrying information are Doppler sensitive waveforms [28]. Because the MSS method uses M distinctive sequences with the same length, it also needs to be set as a bank of discrete correlators. Therefore, the method presented needs to implement the correlation maximum search in the two-dimensional grid composed of spreading code numbers and different Doppler factors for Doppler factor estimation. The estimated Doppler is used to adjust the data length of the current symbol and synchronize the next.
When an idea similar to CAF is used, its disadvantages also become that of the proposed method, namely, a large amount of calculation and insufficient estimation accuracy. For estimation accuracy, the MSS technique is different from OFDM and SC-MPSK. It tolerates a certain degree of Doppler through the gain brought by spreading and despreading. Therefore, it does not need a high estimation accuracy if the residual Doppler is within its tolerance. The verification will be given in Section 3. This paper proposes a three-step optimization scheme to solve the problem of high computational complexity, which is also the innovation of this paper.
The first step is to adjust the frequency-domain search step size. Usually, a small search step size is required in the two-dimensional grid to get the maximum of all the correlation results in the traversal range. In this paper, we find that when searching in the frequency domain, within a certain range, this method is more tolerant to relatively large step sizes. Moreover, through this approach, we reduce the number of calculations to nearly 1/20 of the original.
The second step of the optimization scheme is to use a heuristic search method [29]. Observing the correlation results, we find that the maximum and the two adjacent values on the Doppler axis have obvious discrimination with other correlation results. This paper takes this as heuristic information and sets a threshold during the search. To reduce the amount of calculations, we present a random search method [30], which is different from the exhaustive method [31]. Some points are taken randomly from the search range in each iteration. We calculate the correlation results, obtain the temporary maximum, and compare it with the threshold. If the threshold is exceeded, only two times more calculations are needed neighboring the provisional maximum on the Doppler axis. The maximum among the results is the final decision. To prevent repeated search and ensure that the points are taken differently each time, we use the modern Fisher-Yates shuffle algorithm [32][33][34] to traverse the search results randomly. When the maximum number of iterations is reached, all points in the range are taken. The Fisher-Yates shuffle is an algorithm for generating random arrangements of finite linear arrays. It is unbiased, so the probability of each point being randomly selected is the same. The modern version is more computationally efficient and does not require additional storage space.
In the third step, we assume that the Doppler factor at the symbol level follows the Markov process [35]. Only the first symbol needs to be searched in the whole expected frequency domain. Then, the Doppler range of the following symbol is determined by the previous one, which dramatically reduces the search range of the two-dimensional grid. In this way, the calculation is further reduced.
The rest of this paper is organized as follows. In Section 2, we introduce the transmitter structure, UWA channel modeling, and the receiver structure of the MSS communication system. In Section 3, a detailed description of the receiving and processing optimization scheme is given, as are some simulation results. In Section 4, we provide a sea trial and data analysis. Section 5 summarizes this paper.

Transmitter Structure
MSS is actually a kind of (N, K) encoding, which uses a group of pseudo-random codes with a quantity of M (M = 2 k ) and a length of N to represent the k-bit binary data information. MSS obtains processing gain K times higher than that of conventional DSSS with the same bandwidth. At the same time, MSS uses 2 K pseudo-random codes, and their corresponding relationship with K bits information is variable. Hence, the spectrum of the MSS signal is closer to white Gaussian noise and has stronger interception resistance and information confidentiality than the DSSS.
Among various pseudo-random sequences used in the DSSS system, the chaotic sequences have an excellent correlation performance. Compared with classical pseudonoise sequences, the length of the chaotic sequences can be adjusted unrestricted. For a fixed lengthsequence, there are many available chaotic sequences, which is especially suitable for the MSS system where there are more requirements for the number of codes [36].
After spread spectrum modulation, the signal is carrier modulated and transmitted into the UWA channel through a power amplifier and transducer. Assuming that the spreading sequences used by the MSS are c 1 , c 2 , . . . , c M , and every c i (1 ≤ i ≤ M) is an N × 1 vector, the chaotic sequence set is defined as where P is an N × M matrix. This denotes the sequences in P as 1, 2, . . . , M. The binary data, a[i], 1 ≤ i ≤ K are divided into a group every log 2 (M) bits. The bits of each group are converted into a decimal number, , by using a user-defined mapping relationship. We use the mapping method in [37]. For example, if M = 8, there are 3 bits in one group. If the data bits are [1 0 0], the corresponding decimal number is 4. Since the decimal number range is [0, 7], we need to change the range to [1,8]. We add one to the decimal number and get D[i]. Then, we choose the fifth spread sequence in P, i.e., c 5 . At the receiving end, we subtract one from the decoded number and convert it into the binary numbers to recover the original data.
Using the obtained D[i], we choose the chaotic sequence from P to construct the MSS baseband signal, which is expressed as where d is an N × 1 vector, α = [1, 0, . . . , 0] T M×1 , (·) T denotes transposition, and Q is the MSS mapping matrix defined as where I M−1 is the M − 1 order unit matrix and Q is an M × M matrix. The process of modulating data with M spreading sequences is achieved by left-multiplying P and Q. For example, if the decimal data converted by the mapping relationship is seven, the effect is that all the elements in P are shifted circularly from left to right seven times then multiplied by α, i.e., the sequence with the number of seven in vector α is chosen. Subsequently, the transmission signal of the nth symbol is expressed as where (n − 1)T s < t < nT, T s is the symbol duration, g(t) is the chip impulse shaping function, and f c is the carrier center frequency.
Because the pseudo-random sequence used in this paper is binary, the value of d[k] is −1 or 1. It depends on the current spreading sequence. d n is obtained by Equation (2) with a length of N and a chip duration of T c . {·} denotes the real part of the transmitted signal. With only one conversion rule, there are M! possibilities for the available mapping relationship, in theory. The security performance is obviously better than that of the DSSS. If the number of the transmitting symbols is L (L is an integer larger than zero), the transmitted signal is rewritten as

Channel Modeling
According to the description in [38], the influence of the UWA channel on communication is divided into two time-varying parts. One is expressed by slow varied channel impulse response, and the other is the effect of time-variant Doppler. We express the UWA channel by the following two equations where r(t) is the observation signal at the receiving end, h(t) is the impulse response of a slow time-varying channel, v(t) is the radial relative motion velocity of the transmitter and receiver, c is the sound propagation velocity, and n(t) is the additive white Gaussian noise. The UWA signals are affected by the multi-path effect. In general, different paths could have different Doppler factors. In [20], they use the assumption that all paths have the same Doppler scaling factor. They find that as long as the dominant Doppler shift is caused by the direct transmitter/receiver motion, the assumption seems to be justified. The same assumption is used in [15].
Due to the excellent autocorrelation performance of the pseudo-random sequences, it is easy to distinguish each multi-path signal. In addition, we assume that the time delay is less than the duration of one spreading symbol and ignore the impact of intersymbol interference. Hence, without losing generality, we focus on the Doppler's distortion on the propagated signals [15,26]. Define the Doppler factor as The channel model is simplified as The core of this paper is the received signal processing method under the change of motion over the duration of the transmitted signals. The duration of a one-symbol signal is relatively short, which is usually tens of milliseconds. Therefore, we assume that the relative motion velocity of transmitter/receiver remains changeless within one symbol, i.e., the Doppler factor of one symbol is a constant. During the whole transmission duration, the velocity of transmitter/receiver varies continuously, and the Doppler factors of different symbols are distinct. On this premise, the received signal model is modified, derived from Equations (4), (6), and (10) where ∆ n represents the Doppler factor on the nth symbol. It can be seen from Equation (11) that both the time-domain scale and carrier frequency of the signal are affected by Doppler effects. The time-domain scale variation on each symbol is −∆ n 1+∆ n T s , and the frequencydomain scale variation is ∆ n f c . In terrestrial wireless communications, the propagation speed of electromagnetic waves is quite large. According to Equation (9), the Doppler is small in general, so the time-domain amplitude variation of signal is completely ignored.
However, because the sound velocity in the water is only 1500 m/s, the Doppler factor underwater is five orders of magnitude higher than that on the land. Small relative motion causes a change of signal length that cannot be ignored. When each symbol is affected by different Doppler factors, the change of the symbol length is more complex. Coupled with the distortion of signal carrier frequency, it is not easy to recover data without corresponding processing at the receiving end.

Some Assumptions
Before receiving processing, some assumptions are given.

1.
Searching in the spreading code domain and Doppler domain, we only consider the correlation between the received signal and the passband reference signal with frequency offsets, omitting the signal length compression/expansion. This is because it has little effect on the correlation results for one single symbol, which is verified through simulations in Section 3. Moreover, although signal sampling conversion can simulate the influence of Doppler frequency offset and data length compression/expansion, the number of calculations is terrific, and it is easy to introduce additional errors.

2.
When intercepting the signal for processing, the Doppler factor of the current symbol is not known in advance. It is assumed that the length of the one-symbol-received signal is unchanged. After Doppler estimation, the length of the current symbol is adjusted according to the estimation results. Similarly, the motion speed of AUV is generally below 10 m/s. We believe that the data length error of the order is not enough to affect the despreading performance.

Receiving Processing
After receiving the signal, we use the chirp signal in front of the data frame for synchronization. Then, we intercept the signal symbol by symbol, as shown in Figure 1. The received signal most related to the local reference signal is searched in two dimensions of spread spectrum code domain and Doppler domain. The search range is using all of the spreading codes and the expected Doppler factor range. One symbol for a received signal affected by Doppler is expressed as where d n is the spreading sequence obtained by Equations (1)-(3). d n ∈ P. s ij r (t) represents the passband reference signal. In the additive white Gaussian noise channel, the optimal receiver is the one with the maximum amplitude in the correlator group (14) Then the cross-correlation result between the nth intercepted symbol and the passband reference signal is where R c (τ) is the cross-correlation between the one-symbol-received signal and a passband reference signal with spreading sequence, c i , and the Doppler factor, ∆ j . w(τ) represents the equivalent correlation noise term, and the superscript stands for conjugation. According to Equation (14), each symbol requires M 2 f d U correlation operations to find the maximum peak value and obtain the spreading sequence number, i, which is used by the reference signal when the correlation result reaches maximum value.
According to the mapping relationship between the spreading code number and the decimal number,D[i] is obtained. The data are recovered by the corresponding relationship between information bits and decimal numbers;â[i], 1 ≤ i ≤ K, and the decoding of the current symbol are implemented.
Using the Doppler estimation,∆ n , of the present symbol, we adjust the data length of the current symbol to accomplish the next symbol synchronization.
where V n is the starting position of the nth received signal, and len is the length of the single transmitted symbol. Then we intercept the next symbol according to the single transmitted symbol length, repeat the above process, and finish the decoding of all symbols.

Optimization and Performance Evaluation
Before introducing optimization, the assumptions made in Section 2 need to be verified. This paper compares the difference between the assumption-based correlation performance and the real situation. The simulation condition is a multi-path channel, and the channel impulse response is shown in Figure 2. Gaussian white noise is added to the signal, with −10 dB SNR. The maximum Doppler factor is 0.008. If the sound velocity in water is 1500 m/s, the relative motion velocity is 12 m/s. The signal bandwidth is 2000 Hz, the carrier center frequency is 3000 Hz, the spreading code length is 127, and the symbol duration is 63.5 ms. Eight different codes are used, with each carrying 3 bits of data. The pseudo-random sequence used is a binary quantized chaotic sequence generated by the second-order Chebyshev polynomial function with different initial values [39] x In the real case, the frequency offsets caused by the Doppler effect and the changes of data length in the time domain are both considered for the reference signal. Moreover, the received signal is intercepted with one symbol length affected by the actual present Doppler. Based on the assumptions, we only consider the influence of frequency offsets, and intercept the signal by one transmitted symbol length. The correlation results are obtained on the premise of accurate synchronization.
From Figure 3, it can be seen that the position and peak amplitude of the correlation results based on the hypothesis change slightly. As described in the hypothesis above, we only consider the frequency offsets, omitting the signal length compression/expansion. The length of the one-symbol-received signal is assumed to be unchanged. According to Figure 3, this deviation can be ignored. As for the synchronization error accumulation caused by the Doppler effect, it is corrected by the previous Doppler estimation.
We note that the passband spreading sequence has a large time-bandwidth product, which minimizes the error detection probability. The pseudo-random sequence also has a narrow ambiguity function on the Doppler axis [28]. Moreover, chaotic sequences with the same length have good autocorrelation and cross-correlation properties. Figure 4 shows the correlation search results of eight chaotic sequences under the condition of −10 Db SNR in the Doppler variation range of [−0.008, 0.008]. The search step is 0.3 Hz, and the real Doppler factor is 0.0069. Other simulation conditions are consistent with the above. It can be seen that the pseudo-random sequence has obvious discrimination in code and Doppler axis. This is the basis on which the characteristics of passband spreading signal are used for receiving and processing.  In the process of receiving, each symbol needs to obtain the maximum correlation through a large number of correlation operations. With the symbols increasing, the large number of calculations seriously affects the data processing speed. Hence, we consider establishing a grid based on the code number and the Doppler factor in the expected range to determine the search range. Taking the two-dimensional search results as the objective function, the search and the maximum decision process are modeled as an extreme value optimization, as shown in Equation (15). Algorithm optimization is performed next.
The optimization is carried out in three steps. In the first step, we reduce the number of computations by increasing the search step in the Doppler domain. Because the MSS system tolerates a small Doppler factor through despreading, we only need to make the search step within its tolerance range. The subsequent simulation results show that the search accuracy has little effect on the decoding performance within limits. Since there is an obvious difference between most search results and the optimal solution, in the second step, we set a threshold. The random search method is employed to take some points from the grid for correlation operation every cycle. When a result exceeds the threshold, we obtain the optimal solution only by searching near the result. Otherwise, the search scope is traversed. In the third step, we assume that the motion state of the underwater moving target conforms to the Markov process. The Doppler factor of the current symbol is calculated only in the adjacent area of the previous Doppler estimation, and the search range is reduced. A detailed description of all optimization steps is given next.
It can be seen from Figure 5 that if the search process is modeled as an extreme value optimization, the fitness function of correlations is a multimodal function. The simulation conditions are the same as before. When conducting correlation searching in the Doppler range of [−0.008, 0.008] in steps of 0.0001, we find the correlation result is a smooth and continuous three-dimensional surface. The peak position corresponds to the Doppler estimation and the estimation of the spread spectrum code number of the present symbol. The pseudo-random codes have an excellent correlation performance. They also have sharp peaks on the Doppler axis. These characteristics make it possible to estimate Doppler directly from the spreading signal without data assistance. Nevertheless, with such a Doppler step size and code numbers for only one symbol search, the correlation operations of 1e3 (8 × 161 = 1288) are required. The computations required for decoding are directly proportional to the symbols. Therefore, this paper first optimizes the search step size. We show that this optimization reduces the computations of correlation by an order of magnitude later.
With the help of Doppler estimation, we generally use the CAF method to realize accurate symbol synchronization and to suppress the influence of motion state on the UWA signal. Actually, the accuracy of Doppler estimation is not strict in the spread spectrum communication system. In a conventional DSSS communication system, the code synchronization is carried out by the sliding correlation method [40]. After coarse synchronization, they perform correlation operations on the one-symbol-length received signal and the reference signal, and find the position of the maximum correlation value. Then, taking this position as the starting position, they slide back half the symbol length, intercepting the signal with one symbol length again. The above steps should be repeated until the correlation operations are achieved for the whole frame signal. In incoherent UWA spread spectrum communication, the Doppler has an obvious effect on signal compression/expansion. The sliding correlation method can only suppress the cumulative effect of synchronization error caused by a small Doppler factor. Hence, the large Doppler factor is, firstly, modified to a smaller one by Doppler coarse estimation and compensation, then the residual Doppler is processed by sliding search.
The BER curves of the conventional MSS method under different motion scenarios are given in Figure 6. Assuming that the motion speed is constant within the duration of the whole frame signal, we simulate at the speeds of 0.6 m/s, 1.2 m/s, 1.8 m/s, and 2.4 m/s. The simulation results are shown as a blue dashed line, a purple solid line, an orange solid line, and an orange dashed line, respectively. We also give two comparisons of different conditions. One is that the relative motion speed changes continuously up and down at 1.5 m/s within the signal duration. The simulation result is shown as the black line at the top of Figure 6. The other is assuming that we perform accurate Doppler estimation and compensation processing at the receiver, when the uniform speed is 2.4 m/s. The result is shown as the black solid line at the lower left corner in Figure 6. Due to the narrow available bandwidth of underwater acoustic communication, the signal duration is generally in the order of a few seconds to tens of seconds, and the motion speed is difficult to maintain unchanged within the signal duration. This situation is an ideal state, in reality. We only regard it as a reference result. Through this comparison, we find the limitations of the conventional MSS method. Because the method proposed in this paper performs signal processing within the duration of symbol level, the method we use has a wider applicability.
The simulation parameters are consistent with the above. Two conclusions are drawn from the simulation results.

1.
If the relative motion speed changes slowly within the signal duration, with only spreading/despreading and sliding correlation symbol synchronization, the conventional MSS method tolerates the Doppler of 0.6-2.4 m/s (Doppler factor range, 1e−4-1.6e−3). If Doppler estimation and average Doppler compensation are added to the receiver, the MSS method can cope with a slightly larger motion velocity (2.4 m/s).

2.
When the relative velocity changes irregularly in the whole signal duration, even if the average speed is only 1.5 m/s, the MSS method using sliding correlation for symbol synchronization cannot effectively deal with it.
The first conclusion is the reason why we can appropriately increase the frequencydomain search step. In Section 4, the second conclusion is confirmed.
The distinction between the coarse and fine Doppler estimation of the MSS method is shown in Figures 7 and 8. In the process of a two-dimensional correlation search, if the frequency shift of the reference signal is set to a small step, the correlation computations greatly increase. As shown in Figures 7 and 8, we note that even if the search step size is adjusted larger, it is enough to describe the signal variation in the frequency-domain. If setting the Doppler estimation step to 0.002, the difference between the real Doppler factor and the estimated value shall not exceed 0.001. The system should still be able to decode correctly according to the results in Figure 6.  Next, Monte Carlo simulations are conducted to verify the rationality of the hypothesis. As shown in Figure 9, the Doppler factor is roughly estimated (Doppler step size is 0.002) and accurately estimated (Doppler step size is 0.0001) in a fixed frequency domain. The SNR ranges from −15 dB to−8 dB. Each method is calculated 3000 times, and other simulation parameters remain unchanged. The BER performance of coarse Doppler estimation is the same as that of fine Doppler estimation. For each code number of one symbol, [−0.008, 0.008] is taken as the Doppler factor range for the search. A total of 161 correlation operations are required for fine Doppler estimation and only nine for coarse Doppler estimation. In other words, simply changing the Doppler search step size reduces the correlation computations to 9/161. It is also found from Figure 10 that the objective function value of the optimal solution is highly differentiated from the general solution, which contains the local optimum solution, i.e., the peak with a small amplitude. Using this characteristic, this paper proposes the second optimization, the heuristic search method. Unlike the general extreme value optimization process, the optimization proposed in this paper only accepts the best solution in the limited two-dimensional search range. If the code number is incorrect, it directly affects decoding. If the Doppler factor estimation is wrong, it causes the accumulation of symbol synchronization errors. This affects the subsequent symbol correlation results and may indirectly lead to decoding errors.
As is seen from Figure 10, on the Doppler axis the difference between the maximum and its two adjacent values is small. The discrimination between these three points and all the other general solutions is obvious. Using this feature, the paper sets a threshold. When optimizing, we first determine whether the maximum exceeds the threshold. If so, the search cycle is stopped. Only two more points, which are next to the maximum on the Doppler axis, are calculated. The max value between the three is taken as the final decision result. If the threshold is not exceeded, we take out a fixed number of points randomly again, and their correlations with the received signal are calculated. They are still compared with the threshold until all points in the range are taken, and the point with the max correlation result is chosen. The search process is shown in Figure 11.  It should be explained here that this paper uses the Fisher-Yates shuffle algorithm [32][33][34] to design the search process. In this way, the search is not only random, but also has a memory function. Random points are taken to find the solution exceeding the threshold with a certain probability, but while using the minimum number of searches so as to reduce the correlation computations. The memorability of the algorithm is reflected in excluding the previous points in every next cycle, so as to avoid repeated searches. It ensures that the maximum number of the random search does not exceed the exhaustive method. To avoid frequent comparison between the temporary optimal solution and the threshold due to too many cycles, the number of cycles is set to three. According to the simulation parameters, in the MSS system, there are nine Doppler factors and eight spreading code numbers to be searched for each symbol. According to the exhaustive method, there are 72 times correlation calculations needed for every single symbol. Then the maximum is obtained.
It is assumed that the true value position of each symbol in the search area follows a uniform distribution. It is not considered that the search cannot exceed the threshold temporarily (when the threshold is set reasonably, the optimal solution will be found with a high probability). Since the search cycle is three times, the correlations of 1/3 points in the search range are calculated each time. We can stop the cycle as long one of the three points that exceed the threshold is found.
Therefore, from the perspective of probability calculation, let the number of all combinations of m (m ≤ n) elements be taken out from n different elements, expressed by C m n , and the number of all possible permutations of m elements be taken out from n elements, represented by A m n , and the equation is where x! denotes factorial of x.
Taking the above parameters as an example, there are 72 elements, and 24 elements are taken out each cycle. Since the points are taken randomly, there are C 24 72 possibilities in total. As long as one, two, or three elements are taken out from the three points with high correlation, and the remaining 23, 22, and 21 elements are selected from other general solutions, the optimal solution can be found in the first cycle. The probability of finding the temporary optimal solution exceeding the threshold in the first cycle is calculated as In the second search, because we use the Fisher-Yates shuffle method to exclude the points searched in the first iteration, we only need to extract 24 points from the remaining 48 points to calculate the correlations. Then the probability that the optimal solution is not found in the primary cycle and instead found in the second cycle is as follows The probability of finding the optimal solution in the first cycle is 71%, and the calculations are reduced to about 1/3 of the original. The probability of finding the peak in two cycles is about 96%. That means that through this optimization step, the number of calculations is probably reduced to at least 2/3 of the original. If the threshold decision fails, the computation will be higher than the above results, but it is still similar to that of the exhaustive method. The calculation results are only an example of specific parameter conditions. In fact, as long as the distribution of search results conforms to the heuristic information, this optimization method can be used for different search ranges to reduce the calculations.
Based on the motion feature of the UUVs, this paper puts forward a further optimization scheme. Due to the considerable resistance of water, the speed and acceleration of UUVs are limited to a certain extent, since these factors do not change as fast as for surface ships or aircrafts. In this paper, the motion state of an underwater moving target is modeled as a Markov process [35]. It is assumed that the Doppler factor of each symbol is affected by the previous Doppler estimation. Coupled with the premise that the velocity of the underwater moving target does not change much in a short time, we only need to search the Doppler region neighboring the prior estimated Doppler.
For example, if the estimated Doppler factor of one symbol is 0.004, the search range of the next is all the code numbers with the Doppler factors of 0.002, 0.004, and 0.006. Hence, compared with the whole search area, [−0.008, 0.008], the search area is reduced to the original 1/3. If the Doppler factor estimation of one symbol is 0.008, the next search area is only 0.006 and 0.008. In this way, the number of calculations is reduced to about 1/3 of the original.
Using the three-step optimization method, according to the search range parameters utilized in the simulation, the correlation times of the method with no optimization are reduced to The original two-dimensional search for each symbol required 1288 correlation operations; now about 16 correlation operations are needed. The MSS method using sliding correlation only needs M correlation operations for each symbol. It takes the synchronization signal estimation as the average Doppler of the whole frame signal. Therefore, the conventional method does not search on the Doppler axis and cannot adapt to the Doppler change within the signal duration, as shown in Figure 6. Compared with this method, the optimized method needs 2M correlation operations for each symbol. In other words, the computation volume of the optimized method is only about twice that of the conventional method, and the computing time is still at the same order of magnitude. By increasing a small number of calculations, the optimization gains a reliable communication ability under the condition of high moving speed and acceleration.

Experimental Data Analysis and Processing
Because the proposed receiving and processing method is based on some assumptions, it is necessary to carry out experiments in a real marine environment to verify the effectiveness of the scheme. Because of this, a mobile UWA spread spectrum communication trial was conducted in Weihai, China in May of 2021. In the experiment, the length of one data packet is 7.15 s. It includes a 0.2 s linear frequency modulation synchronization signal, 0.2 s protection interval, 6.35 s data frame, and 0.4 s protection interval at the end of the data. The frequency bandwidth is 2 kHz (2-4 kHz). The carrier center frequency is 3 kHz. The sampling rate is 18 kHz, and each spreading symbol carries 3 bits of data. (Kaddoum has confirmed in [41] that using an MSS symbol to carry 3 bits of data is the best compromise between communication rate and system hardware complexity. If more bits are carried, the complexity of this method will exceed that of the spread spectrum M-ary PSK algorithm.) The chip length is 127 for each symbol. The duration of each symbol is 63.5 ms, with 300 bits of data per packet, and a total of 40 packets of data are transmitted. The trial scenario is shown in Figure 12; the receiving end (point R in Figure 12) is fixed, the transmitting ship (point T in Figure 12) moves with the ocean current at a speed of 0-3 m/s, and the distance is 3.7-4.0 km. In order to reflect the BER performance of the proposed method, channel coding is not used. The measurement of coherence time is of great significance for the design of the UWA communication system. For the UWA channel which changes slowly (the coherence time is tens of seconds or even more), the length of the data frame can be increased to improve the effectiveness of the system without frequent use of a detection signal and protection interval. For the time-varied channel (the coherence time is hundreds of milliseconds or even less), the symbol-by-symbol data processing method we proposed is a better choice.
In Figures 13 and 14, we find that the channel structure changes fast, resulting in a short channel coherence time. Within 0.667 s, the temporal coherence decreases below 0.707, so the coherence time of the UWA channel is about 0.667 s. In the experiment, the duration of one data frame is 7.15 s, which is far longer than the coherence time. It is inappropriate to treat the whole data frame as affected by a constant Doppler factor. However, the one symbol length is only 0.0635 s. Therefore, within a symbol duration, we consider that the channel is time invariant (including the multipath structure and the Doppler factor). Utilizing symbol-by-symbol Doppler estimation and synchronization, we can decode correctly. The experimental results prove this.  The time-varying channel impulse response is shown in Figure 13. We do not perform Doppler compensation processing at the receiving end. It can be seen that the moving speed is relatively consistent in one direction. The time delay is about 10 ms, and the energy of the main path is intense, which is in line with the assumptions of the channel modeling. The transducer is placed 8 m undersea, and the hydrophone is about 10 m underwater. The sound speed profile diagram is shown in Figure 15. Due to the downward bending sound line, the communication distance cannot be too far. We process the collected 40 packets of data and compare the proposed method with the sliding correlation method in a conventional incoherent MSS system. Moreover, the average Doppler compensation is added to the MSS method using the sliding correlation. In Figure 16, the green blocks show the BER performance of the proposed method. Our method does not need the Doppler compensation technique. The blue blocks display the BER of the MSS method using the sliding correlation, and there is no Doppler compensation processing in the receiver. We abbreviate this method as "sliding correlation without compensation". The purple blocks exhibit the BER of the MSS method with the sliding correlation. The difference is, before de-spreading, we perform Doppler coarse estimation and compensation for this method. We abbreviate it as "sliding correlation with compensation". We try to explore whether the absence of the Doppler compensation affects the performance of the proposed method. As shown in Figure 16, the proposed method has 33 packets of data decoded correctly, five packets with the BER of 1e−3, and two packets with the BER about 1-2%. The optimized method obviously outperforms the two sliding correlation methods without the help of the Doppler compensation technique.
The BER performance of the sliding correlation method with Doppler compensation is slightly better than that without Doppler compensation. In most cases, Doppler compensation before sliding correlation despreading improves the system performance. However, when the local instantaneous mutation of the Doppler factor is serious, the Doppler compensation affects the overall decoding effect of the signal. It is explained in the next diagram.
Taking the 14th packet data with error in the proposed method as an example, in Figure 17, there is one error bit in the 94th symbol and two in the 100th symbol. The difference between the Doppler factor estimation of the 94th and 100th symbols and their previous symbols are 0.003 and 0.005, respectively. Because the Doppler estimation range of the proposed method is not set that large, it causes bit error. In fact, the error can be corrected by expanding the search range. For this sudden change of velocity in variablespeed motion, the BER of the sliding correlation despreading method without and with Doppler compensations is 2% and 1.67%, respectively, which are both higher than that of the proposed method. The number of error bits of the sliding correlation method without Doppler compensation is six. Using Doppler compensation in the sliding correlation method, there are five error bits.
For the ninth packet, we can see more abrupt changes in Doppler occurring for this packet in Figure 18. By using the proposed optimization scheme, this method is well adapted to this mutation and complete error-free decoding. However, the BER of the two conventional methods is 3.33% and 5%, respectively. Without using Doppler compensation, the sliding correlation method obtains ten error bits. While using Doppler compensation, there are 15 error bits, which is even worse. When the underwater acoustic channel changes slowly, the Doppler compensation based on an average Doppler factor generally has a positive impact on the sliding correlation method. However, when the Doppler varies within the duration of one signal frame, this Doppler compensation method is counterproductive because the Doppler compensation based on an average Doppler factor often causes a more serious Doppler effect in some parts of the whole signal frame. In Figure 18, the Doppler varies more violently in the 6th and 79th symbols by using Doppler compensation, and there are five error bits in the 6th and 79th symbols using the sliding correlation method with Doppler compensation, while there is only one error bit in the 79th symbol without Doppler compensation.  Moreover, we can find that the sliding correlation method has no good solution for the case where the Doppler change exceeds 1e−3. As shown in the figure, the wrong symbols are behind the sudden change of the Doppler. In other words, although Doppler compensation plays a role in the mutation position, it intensifies the Doppler effect of subsequent symbols, thus affecting the normal decoding.
The BER of 40 data packets in the proposed method is 0.0917%. Compared with 1.91% (sliding correlation with compensation) and 2.42% (sliding correlation without compensation), the BER of this scheme is one order of magnitude lower.

Conclusions and Perspectives
In this paper, a mobile MSS UWA communication scheme is proposed. The quasi orthogonality of a pseudo-random sequence and its thumb-tack peak on the Doppler axis is utilized. The received signal is matched symbol by symbol with the reference signals, which have different frequency shifts and spreading code numbers, to find the maximum correlation value. Then we complete the decoding and Doppler estimation of the current symbol. The conventional Doppler compensation is not used in receiving and processing. By utilizing only the large time-bandwidth product of the spreading signals, the symbollevel synchronization is simply carried out after the Doppler estimation. The sea trial results show a satisfying decoding effect. More importantly, on the premise of ensuring the robustness of communication, we optimize the receiving process. A large Doppler search step and some heuristic information are used to reduce the correlation operations. On the premise of ensuring the real-time implementation of hardware, this method obtains good robustness in a variable speed motion scenario.
Many researchers are using CAF for channel information estimation, and two-dimensional search for Doppler and time delay in a certain range. In fact, since the Doppler sensitive waveform on time delay and Doppler axis has a relatively sharp peak, it is not necessary to traverse all correlation results in the search range. Using some heuristic conditions to narrow the search range expands the application scope of this method. If the threshold is reasonably set, it is easier to find the optimal solution accurately in a shorter time. Although the search domain proposed in this paper is different from the conventional CAF method, the idea is actually universal.
In some studies, for better Doppler compensation, more effective methods are used to improve the accuracy of Doppler estimation. The pursuit of extreme Doppler estimation accuracy is not necessarily the only solution in mobile UWA communication. Considering the hardware overhead of underwater equipment research and development, this simple and stable system design needs to be paid attention to.