LSTM-Based Remote Sensing Inversion of Largescale Sand Wave Topography of the Taiwan Banks

: Shallow underwater topography has important practical applications in ﬁsheries, navigation, and pipeline laying. Traditional multibeam bathymetry is limited by the high cost of largescale topographic surveys in large, shallow sand wave areas. Remote sensing inversion methods to detect shallow sand wave topography in Taiwan rely heavily on measured water depth data. To address these problems, this study proposes a largescale remote sensing inversion model of sand wave topography based on long short-term memory network machine learning. Using multi-angle sun glitter remote sensing to obtain sea surface roughness (SSR) information and by learning and training SSR and its corresponding water depth information, the sand wave topography of a largescale shallow sea sand wave region is extracted. The accuracy of the model is validated through its application to a 774 km 2 area in the sand wave topography of the Taiwan Banks. The model obtains a root mean square error of 3.31–3.67 m, indicating that the method has good generalization capability and can achieve a largescale topographic understanding of shallow sand waves with some training on measured bathymetry data. Sand wave topography is widely present in tidal environments; our method has low requirements for ground data, with high application value. wave crest information more precisely, we used the Canny operator [34] to extract the sand ridge line data. The results are presented as a binary map.


Introduction
Extensive sand wave topography is distributed in tidal environments worldwide, such as the North Sea in Europe [1], the South Sea in Korea [2], San Francisco Bay in the Americas [3], and shallow shoals in Taiwan, China [4]. The study of sand wave shoal topography forms an important basis for coastal protection [5], navigation safety [6], submarine pipeline laying, and drilling platform construction [7]. The shallow sand wave topography of the Taiwan Banks, located offshore, is the largest sand wave shallow topography in the world [4]. The sand wave water depths are roughly distributed between 20 and 60 m [8], and water turbidity makes it difficult for sunlight to penetrate the seafloor, rendering seafloor reflectivity-based detection methods unsuitable [9]. Currently, there are three methods that can be used for sand wave shallow terrain detection: sonar multibeam bathymetry [10], synthetic aperture radar (SAR) detection, and sun glitter remote sensing. Among them, multibeam bathymetry uses vessels as the detection platform and has high efficiency, accuracy, and resolution in detecting trajectories. However, the scanning width is extremely limited, the coverage area is small, the measurement period is long, and the manpower and financial requirements are high [11].
Both SAR sounding and sun glitter remote sensing are ocean dynamics sounding techniques based on sea surface roughness (SSR) and are not constrained by seawater turbidity. They can be operated at large depths [12] and have a wide range of applications. A model for detecting shallow underwater topography (AH model) was developed using images were used to obtain high-resolution SSR information. A region with available measured water depth was selected, and the surface roughness and measured water depth were input into the LSTM model for training. Afterward, the trained model was applied to other areas to obtain the bathymetry of the inversion area. Finally, the model was validated using navigational survey lines to observe the migration capability of the model. The accuracy of the model and the factors affecting the accuracy were analyzed. The remainder of this paper is organized as follows: Section 2 introduces the study area and data sources and analyzes the mapping of water depth and roughness to time-series features. Section 3 introduces the LSTM model, data processing methods, and the training and inversion of the LSTM model. Section 4 analyzes the results of the experiments and the results applied to other regions. Section 5 discusses the results of the inversion. Finally, Section 6 concludes the paper.

Study Area
Marine areas with strong hydrodynamics and sufficient sand sources typically develop a rhythmic bottom bed morphology. The study area located in the southern Taiwan Strait is the world's largest submerged sand wave system [4]. As shown in Figure 1,the Taiwan Shoal extends from 117 • 14 to 119 • 26E and from 22 • 32.5 to 23 • 49N, spanning 228 km from east to west and 145 km from north to south. A total of 16,400 km 2 of the Taiwan Shoal, indicated by the red boundary in Figure 1, was observed using remote sensing imagery [8]. The sand wave directions (0 • due north) are distributed between 132 • and 264 • , with a statistically normal trend. The probability of wave directions at 180 • is the highest, and the sand ridge direction is primarily oriented east to west. The wavelengths of the sand waves in the Taiwan Banks range from 76 m to 2151 m, with a statistical mean wavelength of 721 m. Ridge lengths range from 254 to 19,992 m, while most are shorter than 3000 m, with a statistical mean ridge line length of 2453 m. The morphology of shallow sand waves is relatively diverse, with types including double-peaked sand waves, cosine sand waves [30], and pendulum sand waves [31]. The amplitude of sand waves in the Taiwan Banks is approximately 10 m, and the average water depth is approximately 20 m.
The remainder of this paper is organized as follows: Section 2 intro and data sources and analyzes the mapping of water depth and rou features. Section 3 introduces the LSTM model, data processing meth and inversion of the LSTM model. Section 4 analyzes the results of the results applied to other regions. Section 5 discusses the results of t Section 6 concludes the paper.

Study Area
Marine areas with strong hydrodynamics and sufficient sand velop a rhythmic bottom bed morphology. The study area located in Strait is the world's largest submerged sand wave system [4]. As s Taiwan Shoal extends from 117°14 to 119°26E and from 22°32.5 to 23°4 from east to west and 145 km from north to south. A total of 16,40 Shoal, indicated by the red boundary in Figure 1, was observed usin agery [8]. The sand wave directions (0° due north) are distributed be with a statistically normal trend. The probability of wave directions and the sand ridge direction is primarily oriented east to west. Th sand waves in the Taiwan Banks range from 76 m to 2151 m, with a s length of 721 m. Ridge lengths range from 254 to 19,992 m, while m 3000 m, with a statistical mean ridge line length of 2453 m. The mo sand waves is relatively diverse, with types including double-peake sand waves [30], and pendulum sand waves [31]. The amplitude o Taiwan Banks is approximately 10 m, and the average water depth m. ured bathymetry resolution was also reconstructed to 15 m according to the meth in [20], which was kept uniform with the remote sensing data. According to Zh [4], in the Taiwan Banks, large sand waves do not migrate and sediments only mo the wave crests, and we affirm that the remote sensing data and measured bat data correspond with this.

Dataset
In this study, the SSR inversion model based on multi-angle sun glitter, proposed by Zhang et al. [32], was used to invert the multi-angle sun glitter data obtained from the 3N and 3B channels of ASTER remote sensing, imaged on 16 July 2003. The inversion obtained the same SSR data with a resolution of 15 m as in [32]. ASTER is an advanced multispectral imager, and the ASTER visible NIR subsystem has three bands. The sky bottom telescope image (NVI) of channel 3N has the same bandwidth as the rear-view image (BVI) of channel 3B. The high spatial resolution (15 m), the tilt capability of the sensor, and the rear-view angle of channel 3B render the ASTER imager considerably advantageous for marine applications involving multiple angles. The location depicted in Figure 2 at the Taiwan Shoal is indicated by the black box in Figure 1.
In addition, this study also employed measured bathymetry data obtained using multibeam bathymetry methods in August 2017 and May 2012, as indicated by the blue and purple line segments, respectively, in Figure 2. In the May 2012 survey, 25 measured bathymetric lines occur at 500 m intervals in area A1 in Figure 2b. The small interval between the measured lines was filled using the joint interpolation method in [20] as the whole training area. They were interpolated from a dense set of measured points and fed into the model as real terrain for training. Areas A2, A3, and A4 have measured lines crossing them, as shown in the figure. These multibeam bathymetry data were acquired using the R2 Sonic 2024 multibeam bathymetry system. The system is a fifthgeneration physical multibeam system and can be used for the topographic mapping of the seabed at water depths of 1-500 m. The operating frequency was adjustable between 200 and 400 kHz, coverage width was 10 • to 160 • , and beam angle was 0.5 • × 1 • . In total, 256 effective detection beams could be formed, and the resolution was 1.25 cm. The measured bathymetry data with an average distance of 5 m were obtained after the sounding, and the measured bathymetry resolution was also reconstructed to 15 m according to the method of He in [20], which was kept uniform with the remote sensing data. According to Zhou et al. [4], in the Taiwan Banks, large sand waves do not migrate and sediments only move along the wave crests, and we affirm that the remote sensing data and measured bathymetry data correspond with this.

Sand Wave Terrain Feature Analysis
To explore both the general characteristics of the sand waves and roughness, we considered the profile line characteristics of the sand wave topography and SSR, respectively. The strongest signal variation was found perpendicularly to the sand ridge line. As shown in Figure 3a,b, a section of the profile lines was considered perpendicularly to the sand ridge line, and the variations in SSR and water depth were obtained as shown in Figure 3c.

Sand Wave Terrain Feature Analysis
To explore both the general characteristics of the sand waves and roughness, we con sidered the profile line characteristics of the sand wave topography and SSR, respectively The strongest signal variation was found perpendicularly to the sand ridge line. As show in Figure 3a,b, a section of the profile lines was considered perpendicularly to the san ridge line, and the variations in SSR and water depth were obtained as shown in Figure 3c  Figure 3c shows that along the direction perpendicular to the sand ridge line, bot SSR and bathymetry showed the same period of fluctuation; however, both showed di ferent characteristics of variation. The peaks and troughs were interspersed, and the var iations in bathymetry were always continuous, with the bathymetry at one point bein closely related to that of the surrounding areas. The surface roughness modulated by th sand wave topography varied with the same periodicity as the topography; however, it lightness and darkness did not correspond with the bathymetry data, and a compariso of spatial locations revealed a sudden change in surface roughness at the crest of the san wave.
The sequence characteristics displayed by the above-mentioned SSR and sand wav topography are similar to those of contextual sequences and periodic time sequences. Tex sequences, for example, indicate contextual relevance [33]; weather conditions on consec utive days have some connection with those over past days [34]. The average wavelengt of the sand waves in the Taiwan Banks is 721 m [8] and is characterized by long sequence In this study, we introduce an LSTM neural network, a time-series model for long se quences, to find the continuous and periodic variation characteristics of SSR and san wave topography in space, and to establish a prediction model to realize the remote sens ing inversion of largescale sand wave topography.

Methods Flow
SSR is not only modulated by underwater topography, but also by the wind field a the sea surface [32]. We extracted the SSR sequence modulated by underwater topograph as the input data corresponding to the respective bathymetric sequence. Next, we used reasonable filter sequence to process the wind streaks and enhance data abundance. W subjected the original roughness image to fast Fourier transform (FFT) filtering and cros filtering, extracted the sand ridge line data to calibrate the location of the sand wave cres  Figure 3c shows that along the direction perpendicular to the sand ridge line, both SSR and bathymetry showed the same period of fluctuation; however, both showed different characteristics of variation. The peaks and troughs were interspersed, and the variations in bathymetry were always continuous, with the bathymetry at one point being closely related to that of the surrounding areas. The surface roughness modulated by the sand wave topography varied with the same periodicity as the topography; however, its lightness and darkness did not correspond with the bathymetry data, and a comparison of spatial locations revealed a sudden change in surface roughness at the crest of the sand wave.
The sequence characteristics displayed by the above-mentioned SSR and sand wave topography are similar to those of contextual sequences and periodic time sequences. Text sequences, for example, indicate contextual relevance [33]; weather conditions on consecutive days have some connection with those over past days [34]. The average wavelength of the sand waves in the Taiwan Banks is 721 m [8] and is characterized by long sequences. In this study, we introduce an LSTM neural network, a time-series model for long sequences, to find the continuous and periodic variation characteristics of SSR and sand wave topography in space, and to establish a prediction model to realize the remote sensing inversion of largescale sand wave topography.

Methods Flow
SSR is not only modulated by underwater topography, but also by the wind field at the sea surface [32]. We extracted the SSR sequence modulated by underwater topography as the input data corresponding to the respective bathymetric sequence. Next, we used a reasonable filter sequence to process the wind streaks and enhance data abundance. We subjected the original roughness image to fast Fourier transform (FFT) filtering and cross filtering, extracted the sand ridge line data to calibrate the location of the sand wave crest using the algorithm, and combined the three as input features to build a bathymetric inversion model through LSTM network training. In Figure 2, an iterative training model for A1 (255 km 2 ) area for accuracy verification is presented. The trained inversion model was used to calculate the water depths in the A2 (224 km 2 ), A3 (202 km 2 ), and A4 (348 km 2 ) areas, and the track survey line was used for verification. A flow chart detailing the entire experiment is shown in Figure 4.
x FOR PEER REVIEW 6 of 20 using the algorithm, and combined the three as input features to build a bathymetric inversion model through LSTM network training. In Figure 2, an iterative training model for A1 (255 km 2 ) area for accuracy verification is presented. The trained inversion model was used to calculate the water depths in the A2 (224 km 2 ), A3 (202 km 2 ), and A4 (348 km 2 ) areas, and the track survey line was used for verification. A flow chart detailing the entire experiment is shown in Figure 4.

Roughness Filtering
First, the effect of wind streaks on the accuracy of the terrain inversion by FFT filtering was excluded during preprocessing. Converting the roughness image from the spatial domain to the frequency domain revealed a significant difference in direction between the roughness modulated by the underwater terrain and that affected by the wind field. Data on the characteristic wavelengths and directions of the sand waves were selected in the frequency domain for FFT inversion transformation, and the results obtained are shown in the second image in Figure 5b. Thus, sand wave information was highlighted, and the effect of wind streaks was largely eliminated.

Roughness Filtering
First, the effect of wind streaks on the accuracy of the terrain inversion by FFT filtering was excluded during preprocessing. Converting the roughness image from the spatial domain to the frequency domain revealed a significant difference in direction between the roughness modulated by the underwater terrain and that affected by the wind field. Data on the characteristic wavelengths and directions of the sand waves were selected in the frequency domain for FFT inversion transformation, and the results obtained are shown in the second image in Figure 5b. Thus, sand wave information was highlighted, and the effect of wind streaks was largely eliminated.
The second step of preprocessing involved mean filtering with a series of cross windows to determine the relationship between the variations in the two sequences of roughness and water depth at different scales, and the data were then fed into the LSTM model as supplementary data for training the model. The shape of the cross-filtered window is shown in Figure 5a, which beneficially retained similar water depth in the horizontal direction and reduced the streaking in the inversion results. After experimental validation, we identified four sets of cross-filter templates, namely (3,3), (3,7), (3,15), and (3,25). Figure 5b shows the results of the original roughness, FFT filtering, and four sets of cross-mean filtering assessments of the blue line profile in Figure 3a.
x FOR PEER REVIEW 7 of 20 The second step of preprocessing involved mean filtering with a series of cross windows to determine the relationship between the variations in the two sequences of roughness and water depth at different scales, and the data were then fed into the LSTM model as supplementary data for training the model. The shape of the cross-filtered window is shown in Figure 5a, which beneficially retained similar water depth in the horizontal direction and reduced the streaking in the inversion results. After experimental validation, we identified four sets of cross-filter templates, namely (3,3), (3,7), (3,15), and (3,25). Figure  5b shows the results of the original roughness, FFT filtering, and four sets of cross-mean filtering assessments of the blue line profile in Figure 3a.

Sand Wave Crest Position Characteristics
The location of the sun glitter brightness and darkness reversal was the sand ridge line [19]. Precise extraction of sand ridges has been shown to facilitate the identification of wave crests. To retain sand wave crest information more precisely, we used the Canny operator [34] to extract the sand ridge line data. The results are presented as a binary map.

LSTM Networks
LSTM neural networks are a good solution to the problem of gradient disappearance and gradient explosion that can occur in recurrent neural networks when learning over long sequences [24,35]. The basic unit of an LSTM network implicit layer is called a storage block, and its structure is shown in Figure 6.

Sand Wave Crest Position Characteristics
The location of the sun glitter brightness and darkness reversal was the sand ridge line [19]. Precise extraction of sand ridges has been shown to facilitate the identification of wave crests. To retain sand wave crest information more precisely, we used the Canny operator [34] to extract the sand ridge line data. The results are presented as a binary map.

LSTM Networks
LSTM neural networks are a good solution to the problem of gradient disappearance and gradient explosion that can occur in recurrent neural networks when learning over long sequences [24,35]. The basic unit of an LSTM network implicit layer is called a storage block, and its structure is shown in Figure 6. The second step of preprocessing involved mean filtering with a series of cross wi dows to determine the relationship between the variations in the two sequences of roug ness and water depth at different scales, and the data were then fed into the LSTM mod as supplementary data for training the model. The shape of the cross-filtered window shown in Figure 5a, which beneficially retained similar water depth in the horizontal d rection and reduced the streaking in the inversion results. After experimental validatio we identified four sets of cross-filter templates, namely (3,3), (3,7), (3,15), and (3,25). Figu 5b shows the results of the original roughness, FFT filtering, and four sets of cross-mea filtering assessments of the blue line profile in Figure 3a.

Sand Wave Crest Position Characteristics
The location of the sun glitter brightness and darkness reversal was the sand ridg line [19]. Precise extraction of sand ridges has been shown to facilitate the identificatio of wave crests. To retain sand wave crest information more precisely, we used the Cann operator [34] to extract the sand ridge line data. The results are presented as a binary ma

LSTM Networks
LSTM neural networks are a good solution to the problem of gradient disappearan and gradient explosion that can occur in recurrent neural networks when learning ov long sequences [24,35]. The basic unit of an LSTM network implicit layer is called a stora block, and its structure is shown in Figure 6. The memory block contains the forget gate (f ), input gate (i), output gate (o), and memory cell (C). All three gating structures are identical and are implemented by the sigmoid neural layer dot product operation, and the elements of the sigmoid layer output take values in the range [0, 1], indicating the weight required to let the corresponding message through. The symbol ⊗ represents the dot product operation of two vectors; ⊕ represents the addition operation of two vectors; σ represents the sigmoid activation function; tanh is the hyperbolic tangent activation function; and x represents the input value. σ and tanh were calculated as follows.
The following characterized the flow of information in the LSTM storage block.
Output Gate:

5.
Cell State: 6. Hidden Gate: and W c represent the weight vectors of the forget gates, input gates, output gates, and memory cells, respectively; b f , b i , b o , and b c represent the bias vectors of the forget gates, input gates, output gates, and memory cells, respectively; and X t denotes the input to the t-node network.
The forget gate determined how historical information was retained; input gates determined how information from the input layer was passed to the memory unit; and output gates determined how information from the memory module was passed to the next moment of the storage block. Gate controllers described the proportion of messages that could pass through, and the σ function took values in the range [0, 1]. f t , i t , and o t are the outputs of the t-node σ function; C t is the output of the tanh function at node t, taking values in the range [−1, 1]. After the input sequence had been gated as described above, the t-node long-term memory (C t ) and short-term memory (h t ) were obtained and passed into the next memory block. In addition, h t was the output of node t (t = 1, 2, 3, . . . , n).

Training of LSTM Networks
The model was developed to minimize the distance between the predicted and actual water depths, and we used the root mean square error (RMSE), which is the distance between the two, as the indicator of the loss function. The updating of parameters, such as the weight vector (W) and the bias vector (b), during the training of the LSTM model depended not only on the current node but also on the previous node; this process is called back propagation in time (BPTT) [36]. In this process, the error derivatives are updated as time propagates, and matrices W and b are updated separately. In the experiment, the Remote Sens. 2021, 13, 3313 9 of 20 comparison of various gradient descent methods, such as Adam and SGD, revealed that the Adam algorithm had the highest accuracy. Thus, we used the Adam algorithm to update the parameters in our study.
The data obtained from preprocessing were used as input data and expressed as X = {X t }.t = 1,2,3,..., n with Y = {Y t }.t = 1,2,3,..., n, where t denotes the input data at position t and contains one set of FFT-filtered data, four sets of cross-filtered data on roughness, and one set of binary map data. The predicted water depth data obtained from X t are denoted as h t , and Y t denotes the actual water depth at position t. Given the window length L of the network, this parameter represents the prediction of the water depth at the end, based on an input sequence of length L. X t ,X t+1 . . . . . . X t+L−1 predicts the water depth at the end h t+L−1 . The topology of the neural network determined from L is shown in Figure 7, where LSTM t denotes an LSTM storage block at position t.
between the two, as the indicator of the loss function. The updating of parameters, such as the weight vector (W) and the bias vector (b), during the training of the LSTM model depended not only on the current node but also on the previous node; this process is called back propagation in time (BPTT) [36]. In this process, the error derivatives are updated as time propagates, and matrices W and b are updated separately. In the experiment, the comparison of various gradient descent methods, such as Adam and SGD, revealed that the Adam algorithm had the highest accuracy. Thus, we used the Adam algorithm to update the parameters in our study.
The data obtained from preprocessing were used as input data and expressed as X = {Xt}.t = 1,2,3,..., n with Y = {Yt}.t = 1,2,3,..., n, where t denotes the input data at position t and contains one set of FFT-filtered data, four sets of cross-filtered data on roughness, and one set of binary map data. The predicted water depth data obtained from Xt are denoted as ht, and Yt denotes the actual water depth at position t. Given the window length L of the network, this parameter represents the prediction of the water depth at the end, based on an input sequence of length L. Xt,Xt+1……Xt+L−1 predicts the water depth at the end ht+L−1. The topology of the neural network determined from L is shown in Figure 7, where LSTMt denotes an LSTM storage block at position t. Using this LSTM neural network structure, the input and output data were trained and predicted, respectively, using the following training process: (1). Two-dimensional images were transformed into one-dimensional images. To simulate continuously changing time streams, we connected the profile lines perpendicularly to the sand ridgeline according to the head and tail of the column. Correspondingly, we connected the topographic data head to tail along the profile lines. At this point, the two-dimensional image was converted into a one-dimensional continuous sequence of data. (2). Data normalization: The LSTM model learned the relative trends of two sequences, and, therefore, the Xt and yt data needed to be normalized. Their means were subtracted from both, and they were divided by the variance to obtain the final data set. The normalized depth value of the region was obtained by subtracting water depth data from the mean and dividing it by the variance.  Using this LSTM neural network structure, the input and output data were trained and predicted, respectively, using the following training process: (1). Two-dimensional images were transformed into one-dimensional images. To simulate continuously changing time streams, we connected the profile lines perpendicularly to the sand ridgeline according to the head and tail of the column. Correspondingly, we connected the topographic data head to tail along the profile lines. At this point, the two-dimensional image was converted into a one-dimensional continuous sequence of data. (2). Data normalization: The LSTM model learned the relative trends of two sequences, and, therefore, the X t and y t data needed to be normalized. Their means were subtracted from both, and they were divided by the variance to obtain the final data set. The normalized depth value of the region was obtained by subtracting water depth data from the mean and dividing it by the variance. (3). Network initialization: The weights (W) and bias vector (b) were set to 0 at initialization. hidden_size was set to indicate the number of hidden layer storage block dimensions. The mid layer was set to indicate the number of hidden layers contained in the network. L indicated the input window length, and Lr indicated the step length of each random gradient descent. (4). Data partitioning: The dataset was divided into a training set X tr = {X 1 ,X 2 , . . . X d } and a test set X te = {X d+1 , X d+2 , . . . , X n }. The training set was a subset organized according to the window length L. Each subset obtained was called a batch and was counted as {X tr1 ,X tr2 , . . . ,X trt , . . . ,X trd-L+1 }, where X trt = {X t ,X t+1 , . . . ,X t+L−1 } is the primary input to the LSTM network, and its corresponding output is {h t ,h t+1 , . . . ,h t+L−1 }, which is counted as an epoch according to the number of rounds of the network iteration. (5). The key steps in training the LSTM model are presented in Algorithm 1.

Algorithm 1: LSTM model iteration
Require: the initial value of W and b, number of hidden layers mid_layer, number of hidden layer storage block dimensions hidden_size, input window length L. for k = 1, 2, . . . Epoch do for t = 1,2,d + L − 1 do Calculate f t , i t , o t , C t , C t Get output h t end for Error calculation: E = (h t − y t ) 2 Calculate the Mean Square Error between the node output and the true value Node parameter update: Update W and b based on the error term E using the Adam gradient optimization algorithm. end for The optimal model was trained by adjusting the window length (L), the number of hidden layers (mid_layer), the number of nodes (hidden_size), the learning rate (Lr), and the number of epochs (epoch). Results were obtained when L was the column length of the training area, mid_layer was 4, hidden_size was 16, and Lr was 0.002. The accuracy of the model with iterations is shown in Figure 8. The horizontal coordinates indicate the number of iterations, and the vertical coordinates indicate the RMSE. The model remained largely stable.
for t = 1,2,d + L − 1 do Calculate , , , , Get output ℎ end for Error calculation: E = (ht − yt) 2 Calculate the Mean Square Error between the node output and th Node parameter update: Update W and b based on the error term gradient optimization algorithm. end for The optimal model was trained by adjusting the window le hidden layers (mid_layer), the number of nodes (hidden_size), th the number of epochs (epoch). Results were obtained when L w the training area, mid_layer was 4, hidden_size was 16, and Lr w of the model with iterations is shown in Figure 8. The horizontal number of iterations, and the vertical coordinates indicate the mained largely stable.

Validation and Application
After the training had been completed, the trained network predictions, using the following process: The first dataset (Xte1) of the test set was taken and merged of Xtrd−L+1 to form a new subset {Xd−L+1, Xd−L+2, …, Xd, Xte1}. This su network to obtain the first prediction of the validation set, hte1. Th to obtain the final prediction {hte1, hte2, …, hten}.
The data were normalized before being fed into the network of the mean divided by the variance from both the input and ou data were stitched into one dimension by column. After the norm had been outputted from the network, the variance was multiplie the predicted water depth, and the predicted water depth was stit

Validation and Application
After the training had been completed, the trained network was used to iterate the predictions, using the following process: The first dataset (X te1 ) of the test set was taken and merged with the last L−1 values of X trd−L+1 to form a new subset {X d−L+1 , X d−L+2 , . . . , X d , X te1 }. This subset was input into the network to obtain the first prediction of the validation set, h te1 . The process was repeated to obtain the final prediction {h te1 , h te2 , . . . , h ten }.
The data were normalized before being fed into the network by subtracting the value of the mean divided by the variance from both the input and output values. The image data were stitched into one dimension by column. After the normalized water depth data had been outputted from the network, the variance was multiplied by the mean to obtain the predicted water depth, and the predicted water depth was stitched into a topographic image by column.
The model was trained iteratively in the training area in region A1 shown in Figure 2, and accuracy validation was performed in the test area. The first 75% of the data were selected as the experimental set and the last 25% as the validation set, and good inversion results were obtained. The generalization ability of the model was tested by adjusting the training and test sets. Two parameters, the mean absolute value error (MAE) and the RMSE, were used to characterize the accuracy of the model.
The trained model was applied to a total area of 774 km 2 across the A2, A3, and A4 regions depicted in Figure 2 to obtain the predicted water depths for the three areas, and the accuracy was further verified by combining the presence of track lines in the area.

Model Evaluation
We achieved good inversion results within A1 in Figure 2 by selecting 75% of the data for the experimental set and 25% of the data for the validation set, called Model I. At this point, our training and test sets were placed adjacent to each other, and the training set was kept larger than the test set. Two new sets of experiments were designed to validate the generalization performance of the model. The first set took the first 50% of the data as the training set and the last 25% as the test set to obtain a model called Model II. At this point, the experimental and test regions were separated, and we could observe the effect on model accuracy when the data were adjacent to each other. The second group used 50% of the data as the training set and 50% of the data as the validation set to test the effect of changing the size and proportion of the sample set in the experiment; the resulting model was called Model III. The results of the above three sets of experiments are shown in Table 1. A comparison of the predicted water depth images with the real images is shown in Figure 9; scatter density distribution is shown in Figure 10; and the histogram of the absolute value difference frequency distribution is shown in Figure 11. The model was trained iteratively in the training area in region A1 shown in Figure  2, and accuracy validation was performed in the test area. The first 75% of the data were selected as the experimental set and the last 25% as the validation set, and good inversion results were obtained. The generalization ability of the model was tested by adjusting the training and test sets. Two parameters, the mean absolute value error (MAE) and the RMSE, were used to characterize the accuracy of the model.
The trained model was applied to a total area of 774 km 2 across the A2, A3, and A4 regions depicted in Figure 2 to obtain the predicted water depths for the three areas, and the accuracy was further verified by combining the presence of track lines in the area.

Model Evaluation
We achieved good inversion results within A1 in Figure 2 by selecting 75% of the data for the experimental set and 25% of the data for the validation set, called Model I. At this point, our training and test sets were placed adjacent to each other, and the training set was kept larger than the test set. Two new sets of experiments were designed to validate the generalization performance of the model. The first set took the first 50% of the data as the training set and the last 25% as the test set to obtain a model called Model II. At this point, the experimental and test regions were separated, and we could observe the effect on model accuracy when the data were adjacent to each other. The second group used 50% of the data as the training set and 50% of the data as the validation set to test the effect of changing the size and proportion of the sample set in the experiment; the resulting model was called Model III. The results of the above three sets of experiments are shown in Table 1. A comparison of the predicted water depth images with the real images is shown in Figure 9; scatter density distribution is shown in Figure 10; and the histogram of the absolute value difference frequency distribution is shown in Figure 11.   In Figure 9a, the blue line indicates the 50% area boundary, and the green line ind cates the 75% boundary. The images in Figure 9b-d indicate that all three sets of model essentially inverted the morphology of the sand wave topography and clearly identifie the location of the sand wave crests. Sand wave development was better below the imag than that above, and the inversion results were also better below than above. Figure 1 shows a scatter density plot where the horizontal coordinates are the predicted bathyme try values, the vertical coordinates are the measured bathymetry values, and the differ ence between the horizontal and vertical coordinates at each point represents the predic tion error. The scatter points for the three models were distributed roughly along the d agonal, with a higher density closer to the diagonal, indicating a greater number of point with lower deviations. Figure 11 shows a histogram of the different distributions of th three models. The errors of all three models were roughly normally distributed, with difference of zero, which verified the validity of the models. The accuracy of Model I wa slightly higher than that of Models II and III, but the difference was not significant. Th scatter distribution pattern of the three sets of experiments was similar, which indicates tha the models can be generalized in the area spaced from the training area or in a larger area.

Model Application
The predicted data from the network were multiplied by the variance plus the mea to obtain the final predicted bathymetry values. The measured bathymetry data in are A1 were obtained using a survey vessel, which facilitated the acquisition of the varianc and mean values. Using Model I, predictions were made within areas A2, A3, and A4 i Figure 2. Accuracy validation was performed using the measured line data within th respective regions. We outputted the normalized water depth values using the model.
After the normalized bathymetry had been obtained by inversion, calculating th variance in the bathymetry across the whole area and the regional mean bathymetry wa  In Figure 9a, the blue line indicates the 50% area boundary, and the green line indicates the 75% boundary. The images in Figure 9b-d indicate that all three sets of models essentially inverted the morphology of the sand wave topography and clearly identified the location of the sand wave crests. Sand wave development was better below the image than that above, and the inversion results were also better below than above. Figure 10 shows a scatter density plot where the horizontal coordinates are the predicted bathymetry values, the vertical coordinates are the measured bathymetry values, and the difference between the horizontal and vertical coordinates at each point represents the prediction error. The scatter points for the three models were distributed roughly along the diagonal, with a higher density closer to the diagonal, indicating a greater number of points with lower deviations. Figure 11 shows a histogram of the different distributions of the three models. The errors of all three models were roughly normally distributed, with a difference of zero, which verified the validity of the models. The accuracy of Model I was slightly higher than that of Models II and III, but the difference was not significant. The scatter distribution pattern of the three sets of experiments was similar, which indicates that the models can be generalized in the area spaced from the training area or in a larger area.

Model Application
The predicted data from the network were multiplied by the variance plus the mean to obtain the final predicted bathymetry values. The measured bathymetry data in area A1 were obtained using a survey vessel, which facilitated the acquisition of the variance and mean values. Using Model I, predictions were made within areas A2, A3, and A4 in Figure 2. Accuracy validation was performed using the measured line data within the respective regions. We outputted the normalized water depth values using the model. After the normalized bathymetry had been obtained by inversion, calculating the variance in the bathymetry across the whole area and the regional mean bathymetry was necessary to obtain the specific sand wave bathymetry. Although the sand waves in the In Figure 9a, the blue line indicates the 50% area boundary, and the green line indicates the 75% boundary. The images in Figure 9b-d indicate that all three sets of models essentially inverted the morphology of the sand wave topography and clearly identified the location of the sand wave crests. Sand wave development was better below the image than that above, and the inversion results were also better below than above. Figure 10 shows a scatter density plot where the horizontal coordinates are the predicted bathymetry values, the vertical coordinates are the measured bathymetry values, and the difference between the horizontal and vertical coordinates at each point represents the prediction error. The scatter points for the three models were distributed roughly along the diagonal, with a higher density closer to the diagonal, indicating a greater number of points with lower deviations. Figure 11 shows a histogram of the different distributions of the three models. The errors of all three models were roughly normally distributed, with a difference of zero, which verified the validity of the models. The accuracy of Model I was slightly higher than that of Models II and III, but the difference was not significant. The scatter distribution pattern of the three sets of experiments was similar, which indicates that the models can be generalized in the area spaced from the training area or in a larger area.

Model Application
The predicted data from the network were multiplied by the variance plus the mean to obtain the final predicted bathymetry values. The measured bathymetry data in area A1 were obtained using a survey vessel, which facilitated the acquisition of the variance and mean values. Using Model I, predictions were made within areas A2, A3, and A4 in Figure 2. Accuracy validation was performed using the measured line data within the respective regions. We outputted the normalized water depth values using the model.
After the normalized bathymetry had been obtained by inversion, calculating the variance in the bathymetry across the whole area and the regional mean bathymetry was Remote Sens. 2021, 13, 3313 13 of 20 necessary to obtain the specific sand wave bathymetry. Although the sand waves in the Taiwan Banks formed under complex environmental conditions, the sand wave formation conditions in adjacent areas were similar. We used the roughness of bathymetry data in A1 as an approximation for the three validation regions.
There were 30,000, 60,000, and 130,000 measured bathymetry points on the track lines in areas A2, A3, and A4, respectively, as shown in Figure 2. Although these data points did not cover the entire map, they were representative. The average bathymetry data from the three areas were taken as the average bathymetry of the area. The normalized bathymetry was multiplied by the variance of the bathymetry in area A1, and the mean bathymetry of the survey line was added to obtain the predicted bathymetry. The predicted bathymetry was reduced to a topographic map, as shown in Figure 12. A section of the profile was taken from the measured line in each of the three areas and compared to the predicted bathymetry, as shown in Figure 13. The position of the profile is indicated by the black line in Figure 12. The scatter density distributions of the predicted and measured bathymetry are shown in Figure 14, and the difference distribution is shown in Figure 15. The prediction errors for the three areas are listed in Table 2.
Taiwan Banks formed under complex environmental conditions, the sand wave formation conditions in adjacent areas were similar. We used the roughness of bathymetry data in A1 as an approximation for the three validation regions.
There were 30,000, 60,000, and 130,000 measured bathymetry points on the track lines in areas A2, A3, and A4, respectively, as shown in Figure 2. Although these data points did not cover the entire map, they were representative. The average bathymetry data from the three areas were taken as the average bathymetry of the area. The normalized bathymetry was multiplied by the variance of the bathymetry in area A1, and the mean bathymetry of the survey line was added to obtain the predicted bathymetry. The predicted bathymetry was reduced to a topographic map, as shown in Figure 12. A section of the profile was taken from the measured line in each of the three areas and compared to the predicted bathymetry, as shown in Figure 13. The position of the profile is indicated by the black line in Figure 12. The scatter density distributions of the predicted and measured bathymetry are shown in Figure 14, and the difference distribution is shown in Figure 15. The prediction errors for the three areas are listed in Table 2.       According to Figure 12, on comparing the three sets of roughness images with the inverse bathymetric images, the inversion results reproduced the sand wave pattern. The position of the sand wave crests and of the length and direction of the sand waves were accurately determined. The presence of a trough was due to an impact in the upper right of the A3 area, which is represented by the corresponding position of the inverse bathymetry. Based on the roughness image of area A4, the sand waves were sparse in the upper half of the area and dense in the lower half, with the same effect being seen in the bathymetry image. Though all distinctive sand waves were represented in the inverse bathymetry data, the small and poorly characterized sand waves present in Figure 12c are not present in Figure 12d. The is because the sand waves in this region were shaped too indistinctly to be captured by the model at this time. From Figure 13, the measured and predicted water depths on the three sets of profile lines basically showed the same trend variation. The results of the inversions accurately identified the locations of significant sand wave crests. Although there was a small number of unmatched or incorrectly matched peaks and troughs in all three sets of profiles, and there were some differences between the inverse and real bathymetry values, the overall trend in bathymetry from peak to trough could be reflected.
As indicated by the results in Table 2, the MAE of the inversion results for the three regions ranged from 2.56 to 2.90 m, and the MSE ranged from 3.31 to 3.67 m. As indicated in Figure 14, the predicted and measured values for the three regions were roughly distributed along the diagonal, with an increasing density closer to the diagonal. The difference in values in Figure 15 shows a roughly normal distribution, indicating the validity of the inversion model.
The errors in the inversion results come from three sources: first, the error in the accuracy of the model itself; second, the error caused by the difference between the variance of the water depth in area A1 and in other areas; and third, the error in the mean value of the water depth in the sampled survey line and that in the mean value of the water depth in the whole area. To reduce these errors, the accuracy error of the model itself can be reduced using a higher-accuracy roughness image, or by enhancing the roughness image.  According to Figure 12, on comparing the three sets of roughness images with the inverse bathymetric images, the inversion results reproduced the sand wave pattern. The position of the sand wave crests and of the length and direction of the sand waves were accurately determined. The presence of a trough was due to an impact in the upper right of the A3 area, which is represented by the corresponding position of the inverse bathymetry. Based on the roughness image of area A4, the sand waves were sparse in the upper half of the area and dense in the lower half, with the same effect being seen in the bathymetry image. Though all distinctive sand waves were represented in the inverse bathymetry data, the small and poorly characterized sand waves present in Figure 12c are not present in Figure 12d. The is because the sand waves in this region were shaped too indistinctly to be captured by the model at this time. From Figure 13, the measured and predicted water depths on the three sets of profile lines basically showed the same trend variation. The results of the inversions accurately identified the locations of significant sand wave crests. Although there was a small number of unmatched or incorrectly matched peaks and troughs in all three sets of profiles, and there were some differences between the inverse and real bathymetry values, the overall trend in bathymetry from peak to trough could be reflected.
As indicated by the results in Table 2, the MAE of the inversion results for the three regions ranged from 2.56 to 2.90 m, and the MSE ranged from 3.31 to 3.67 m. As indicated in Figure 14, the predicted and measured values for the three regions were roughly distributed along the diagonal, with an increasing density closer to the diagonal. The difference in values in Figure 15 shows a roughly normal distribution, indicating the validity of the inversion model.
The errors in the inversion results come from three sources: first, the error in the accuracy of the model itself; second, the error caused by the difference between the variance of the water depth in area A1 and in other areas; and third, the error in the mean value of the water depth in the sampled survey line and that in the mean value of the water depth in the whole area. To reduce these errors, the accuracy error of the model itself can be reduced using a higher-accuracy roughness image, or by enhancing the roughness image.  According to Figure 12, on comparing the three sets of roughness images with the inverse bathymetric images, the inversion results reproduced the sand wave pattern. The position of the sand wave crests and of the length and direction of the sand waves were accurately determined. The presence of a trough was due to an impact in the upper right of the A3 area, which is represented by the corresponding position of the inverse bathymetry. Based on the roughness image of area A4, the sand waves were sparse in the upper half of the area and dense in the lower half, with the same effect being seen in the bathymetry image. Though all distinctive sand waves were represented in the inverse bathymetry data, the small and poorly characterized sand waves present in Figure 12c are not present in Figure 12d. The is because the sand waves in this region were shaped too indistinctly to be captured by the model at this time. From Figure 13, the measured and predicted water depths on the three sets of profile lines basically showed the same trend variation. The results of the inversions accurately identified the locations of significant sand wave crests. Although there was a small number of unmatched or incorrectly matched peaks and troughs in all three sets of profiles, and there were some differences between the inverse and real bathymetry values, the overall trend in bathymetry from peak to trough could be reflected.
As indicated by the results in Table 2, the MAE of the inversion results for the three regions ranged from 2.56 to 2.90 m, and the MSE ranged from 3.31 to 3.67 m. As indicated in Figure 14, the predicted and measured values for the three regions were roughly distributed along the diagonal, with an increasing density closer to the diagonal. The difference in values in Figure 15 shows a roughly normal distribution, indicating the validity of the inversion model.
The errors in the inversion results come from three sources: first, the error in the accuracy of the model itself; second, the error caused by the difference between the variance of the water depth in area A1 and in other areas; and third, the error in the mean value of the water depth in the sampled survey line and that in the mean value of the water depth in the whole area. To reduce these errors, the accuracy error of the model itself can be reduced using a higher-accuracy roughness image, or by enhancing the roughness image. It is also possible to improve the representativeness of the sampling to obtain a more accurate mean water depth and to reduce the error introduced by sampling.

Errors at Different Depths
As shown in Figure 13, most of the bathymetric profile was a gentle trough. Wave crests were located in areas where bathymetry varied considerably and covered less of the overall location. To investigate the inversion accuracy of the model for different areas, the first 10% of the shallowest water depths were classified as the crest region (<27 m), the deepest 50% (>33 m) were classified as the trough region, and the middle part was termed the transition region (27-33 m). We statistically calculated the errors at different depths within regions A2, A3, and A4, and the results are shown in Table 3. In all three inversion regions, the error in the wave crest values was significantly higher than that in the other values. There are two reasons for this: first, variance migration has a much greater effect on the peaks and troughs than on the middle part. Second, there was a small number of peaks or troughs that were not matched, and the direct errors in the unmatched peaks and troughs can be greater than 10 m, which significantly increased the errors across the interval. Among the three zones, only the accuracy of the trough data in the A2 zone was higher than the mean value, and the deviations in the positions of the troughs in the other two zones, A3 and A4, were also greater than those in the excess zone. In Figure 12b, several downward peaks from the predicted bathymetry values are evident, and these peaks increased the error in the trough positions. However, because of the wide distribution at trough locations, the sensitivity was lower than that at the crest locations.

Errors at Different Wavelengths
According to Zhang et al. [8], the wavelength of the sand waves in the Taiwan Banks ranged from 76 to 2151 m, with an average wavelength of 721 m, and their peak occurred at 500 m. We considered sand waves less than 200 m as small sand waves, those between 200 m and 500 m as medium sand waves, and those longer than 500 m as large sand waves. The accuracies at different wavelengths on the oblique survey lines in regions A2, A3, and A4 are listed in Table 4. The error increased in all regions as the wavelength increased. The reasons for this are two-fold. First, the wavelength of the sand waves is highly correlated with the height of the sand waves. According to the conclusions, higher sand waves resulted in greater errors. Second, larger wavelengths entail more complex structures, such as bimodal sand waves or pendulum sand waves, and the characteristics of the immediately adjacent bimodal peaks make inversion difficult.

Sensitivity to the Size of the Training Dataset
Neural networks require sufficient training samples. To test the sensitivity of the model to the samples, we used 100%, 75%, 50%, and 25% of region A1 to train the model, and the four groups of trained models were applied to regions A2, A3, and A4. The results are listed in Table 5. As the size of the training set increased, there was an upward trend in accuracy in the validation region. Using 75% of the training samples, the accuracy was 0.2-0.3 m lower than when using the whole region for training. Each subsequent 25% reduction in training data resulted in a greater error. Therefore, the experimental area needed to be as large as possible; we used the whole A1 area to train the model, and made predictions in the A2, A3, and A4 areas. The experimental dataset was limited to 255 km 2 , and an increase in the experimental area did not reach the bottleneck of model accuracy. Adding new bathymetric analysis areas to the training set would improve the inversion accuracy.

Strengths and Weaknesses of the Method
The proposed method is mainly based on SSR for inversion, and remote sensing images with high resolution and high accuracy SSR are available as long as they can be obtained. If an RMSE of 3.31-3.67 m is acceptable for areas without measured water depth data, we can conclude that the proposed LSTM-based inversion method can significantly extend the area of the inversion region. Zhang et al. obtained good results for several sections [8]; however, in their study, the presence of more than two control points located at the crest and trough of the wave on the section was required, and detachment from the control points could not be achieved. He et al. [20] used a joint reconstruction method to achieve the inversion of the water depth; however, accuracy could only be achieved at a maximum of 1 km from the region of measured water depth. Though these data are often used to fill the gaps in multibeam bathymetry data, inversion is not possible in regions without multibeam bathymetry data.
In this study, we proposed the model based on the serial permutation relationship between the results of bathymetry and surface roughness assessment perpendicularly to the sand ridge line. After training on a small amount of measurement data, it can be migrated to other regions for application. For example, in this experiment, bathymetry data of a total sand wave area of 774 km 2 were obtained across three areas, A2, A3 and A4. The core feature of this model is the ability to invert large areas without measured bathymetry data with a guaranteed degree of accuracy and the ability to migrate.
The shallow sand wave topography is widely distributed in the tidal environment, and the model proposed herein can be transferred for use in other sand wave areas. However, when transferring the model for other shallow sand wave areas, such as the North Sea in Europe and San Francisco Bay in the USA, local measured water depth data are still required to iterate the model and obtain a locally applicable sand wave sounding model.

Conclusions
In this study, LSTM models were introduced to predict the water depth of sand waves, and the following conclusions were obtained.
In this study, a time-series LSTM model was introduced to exploit the continuous, periodic relationship between SSR and bathymetry sequences oriented perpendicularly to the sand waves. The relationship in the spatial domain was mapped to the temporal domain to establish a bathymetric inversion model. The model was designed based on the characteristic relationship between water depth and roughness and has a low dependence on measured water depth data. This addresses the concerns of high cost and small coverage associated with current sun glitter bathymetric inversion models, which rely heavily on the measured water depth data.
The effect of the sea-surface wind field was excluded by experimentally selecting suitable features. The iterative model was trained in the experimental area, and the generalization of the model to non-adjacent areas was verified by changing the experimental area. The model was applied in three non-contiguous areas with a total area of 774 km 2 , and the accuracy was verified using existing track lines, for which the RMSE was between 3.31 and 3.67 m, and the observed inverse bathymetric images essentially restored the sand wave pattern, enabling the acquisition of low-cost sand wave topography.
The model is based on the relationship between the continuity and periodicity of SSR and sand wave topography and is generalizable in shallow sand wave regions. It can be extended to other areas of sand wave topography, such as the North Sea in Europe, San Francisco Bay in the USA, the South Sea in Korea, and the Bohai Sea and northern South China Sea in China. The model requires a small number of measured bathymetric points to obtain the mean bathymetry of the region, and the reliance on measured bathymetry can be further reduced using other methods to roughly obtain the background bathymetry information for the region.  Data Availability Statement: ASTER remote sensing data from publicly available datasets.