Using robust correlation matching to estimate sand-wave migration in Monterey Submarine Canyon, California

Marine Geology 376 (2016) 102–108

Contents lists available at ScienceDirect

Marine Geology journal homepage: www.elsevier.com/locate/margo

Using robust correlation matching to estimate sand-wave migration in Monterey Submarine Canyon, California Kai Zhang, Fanlin Yang ⁎, Chunxia Zhao, Chengkai Feng College of Geomatics, Shandong University of Science and Technology, Qingdao, China Key Laboratory of Surveying and Mapping Technology on Island and Reef, National Administration of Surveying, Mapping and Geoinformation, Qingdao, China Key Laboratory of Marine Surveying and Charting in Universities of Shandong, (Shandong University of Science and Technology), Qingdao, China

a r t i c l e

i n f o

Article history: Received 10 October 2015 Received in revised form 22 March 2016 Accepted 1 April 2016 Available online 6 April 2016 Keywords: Sandwave migration Robust correlation matching Outlier Adaptive trimming

a b s t r a c t Sandwave migration rate is of critical importance for the hydrodynamic research as well as the coastal engineering application. To acquire this information, the correlation matching algorithm is applied in this paper. By searching for similar seafloor areas between multi-temporal digital terrain models, the migration rate can be estimated. To account for the impact of the morphological distortion as well as outliers in the survey data, the correlation matching is adjusted based on the robust estimation theory. By executing the correlation matching based on the least median square criterion, a robust initial estimate is obtained. Afterwards, the potential outliers are detected based on the initial matching result, and the corresponding weights are adjusted to zero. Finally, correlation matching is executed with the updated weights. With the above procedure, the impact of outliers can be accounted for. The developed method is applied on the bathymetry time-series data collected in Monterey Canyon, California. Comparison shows that the presented algorithm significantly improves the credibility of the estimation. The generated migration vector field not only reveals the main tendency of the migration rate in the region, but also highlights the unique local migration patterns, which is helpful for understanding the mechanism of the sandwave migration. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The migration of sandwave plays a major part in underwater sediment transportation, and is of critical importance for both hydrodynamic researches (Best, 2005; Hulscher and Dohmen-Janssen, 2005; Van Landeghem et al., 2012) and coastal engineering applications (Besio et al., 2003, 2004). On the one hand, the knowledge of the sandwave migration rate is essential for understanding the behavior of the bedload surface current (Kubo and Nakajima, 2002). On the other hand, the sandwave migration is proved to be substantial for various hydrographic application, such as navigation (Knaapen et al., 2001) and submarine pipeline installation (Nemeth et al., 2003; Dorst, 2004; Van Dijk et al., 2008). The advent of multibeam echo-sounder system (MBES) has advanced the capability of bathymetric survey (Lurton, 2010). As a result, underwater digital terrain model (DTM) with higher accuracy can be generated in a more rapid and cost-effective way. With the collected data, the morphology of the seafloor is able to be described to a better degree, both in accuracy and resolution. In consequence, DTM timeseries data can provide valuable detailed information of the bedforms deformation, which leads to key advances in the analysis of the bedforms migration (Nemeth et al., 2002, 2007; Besio et al., 2003; Best, 2005). ⁎ Corresponding author at: No.579, Qianwangang Road, Qingdao, China. E-mail address: fl[email protected] (F. Yang).

http://dx.doi.org/10.1016/j.margeo.2016.04.002 0025-3227/© 2016 Elsevier B.V. All rights reserved.

When the migration of bedform is concerned, various approaches have been suggested over the past decades. Lindenbergh (2004) estimated the seabed dynamic information under a stable sandwave feature assumption. The approach was updated in Dorst et al. (2009) and (2011) by focusing on the uncertainty of the bathymetry data as well as the parameter estimates. Knaapen (2005) estimated the dune migration by subtracting crest lines and trough lines from the MBES data. Similar procedure was employed in Van Dijk and Kleinhans (2005), Cherlet et al. (2007) and Xu et al. (2008) to characterize the sandwave migration. Assisted by manual judgment, such straightforward method is practical and effective. However, when extracting information from multi-temporal DTM data, it was reported that the simple subtraction of 2 DTMs is sometimes not adequate (Duffy and Hughes-Clarke, 2005). Meanwhile, the manual judgment involved may introduce some level of subjectivity into the estimation result. For these two reasons, several works used correlation method to acquire the migration information of bedforms (Duffy and Hughes-Clarke, 2005; Buijsman and Ridderinkhof, 2008a, 2008b). The method involves comparing terrain nodes with equal window size from 2 different DTMs. The success of this kind of method is based on the premise that the bedform preserves its main morphological feature, which is often a reasonable assumption. The limitation of the method lies in that it addresses the disturbance of outliers by limiting the window size, so as to limit the number of outliers in the correlation process. However, the remained outliers in DTM data can still disturb the correlation. Moreover, irregular

K. Zhang et al. / Marine Geology 376 (2016) 102–108

and inconsistent morphological distortion blurs the horizontal migration features, and thus makes the migration characterization by correlation matching challenging. The goal of this paper is to extend previous works by focusing on the interference of outliers in the process of the estimation. To this end, an algorithm termed the robust correlation matching (RCM) algorithm is developed to estimate the migration rate of the sandwave. By considering the local morphological distortion of bedform as outlier, and address it based on the robust theory, the main trend of migration can be estimated accurately with RCM. The remainder of this paper is organized as follows. Section 2 describes the methodology of the RCM method. Section 3 applies RCM to the multi-temporal bathymetric data in Monterey Canyon, so as to analyze the sandwave migration in the region. Finally, conclusions and discussions are given in Section 4. 2. Study area The time serial bathymetry data of Monterey Canyon are studied in this paper. Monterey Canyon is the largest submarine canyon in the west coast of United States. The region is supplied with a large amount of sediment every year from the river mouth at the northwest side of the canyon. This steady supply of sand and gravel moves gradually along the narrow pathway which is restricted to the canyon axis. As a result, extinct migration of sandwave organization can be observed in the area. In this paper, the study area is located in the head of the canyon (the left part of Fig. 1). The size of the area is around 900 m by 500 m, and the water depth of the area ranges from 30 m to 80 m. Benefiting from the large amount of annual sediment input, distinct sandwave organization has been formed in the region (the right part of Fig. 1). The seafloor mapping Laboratory of the California State University Monterey Bay conducted a series of bathymetric surveys at the head of the Monterey Submarine Canyon (Xu et al., 2008; Paull et al., 2010). In the survey, the bathymetric data were collected with a Reson 8101 multibeam sonar system. The horizontal position of the boat was measured with dual frequency differential GPS, which provided submeter horizontal positioning precision. Meanwhile, the tidal variation was also corrected with the same GPS systems so that the vertical bias was less than 0.2 m (Smith et al., 2005, 2007). After the manual editing, the bathymetic data were gridded into 3-m interval XYZ format. The multi-temporal DTMs used in this paper were collected in September and November of 2004 respectively. These two DTMs have the smallest

103

temporal gap among the whole DTM time series (around two months), which is ideal for sandwave migration estimation. 3. Methodology 3.1. Correlation matching for migration estimation Multi-temporal underwater DTM provides valuable detailed information about the change of bedforms. Based on these data, the spatial correlation technique can be applied to infer the migration of a particular patch of bedform during the different periods. To be more specific, for a particular patch of bedform P in the first DTM (termed foreground DTM), the most similar patch Q with the same size is searched in the second DTM (termed background DTM). Thus, by computing the horizontal displacement between P and Q, the migration information can be obtained. To find the most similar pair between the foreground and the background DTMs, the matching degree is used to quantify the similarity of 2 seafloor patches of different DTMs. To simplify the discussion, DTM is assumed to be in the grid format in this paper. In the foreground DTM, let P denote the terrain grids in a window with dimension (Wx, Wy). Likewise, Q is a patch with equal window size in the background DTM. Thus, to quantify the similarity degree between P and Q, the matching degree S between them is given as: S¼

2 N X ðdl −zl Þ 2 σ d;l þ σ 2z;l l¼1

ð1Þ

where {dl, l = 1 ~ n}and {zl, l = 1 ~ n} are the grid vectors of P and Q respectively, and σ2d , i and σ2z , i are the corresponding variances. The squared difference form given in Eq. (1) is a bit different from the criterion in (Duffy and Hughes-Clarke, 2005), but is used extensively in geostatistics and terrain navigation because of its analytically tractability (Nygren, 2005). To simplify the expression, a weight term is introduced in the equation S¼

n X

2

wl ðhl −zl Þ

ð2Þ

l¼1

where the weight wl ¼

1 . σ 2d;l þ σ 2z;l

Fig. 1. Left panel: the shaded relief image of the bathymetric data at the head of the Monterey Submarine Canyon, California. Right panel: the bathymetry of the study area.

104

K. Zhang et al. / Marine Geology 376 (2016) 102–108

Hence, for a particular patch P in the foreground DTM, the objective is to find a patch Q with equal window size from the background DTM based on the criterion Eq. (2). For this purpose, the position of the window is sliding over the background DTM, so as to find a horizontal position with the smallest S. When plotting the correlation sum Eq. (2) for every candidate position, the so-called correlation surface appears. Then, the minimum point on the correlation surface corresponds ^ that gives the best match: to the displacement D ^ ¼ arg min D D

n X

 2 wl hl −zl;x

ð3Þ

l¼1

where the term D denotes the position of Q in the background DTM. After the matching, the temporal displacement of a certain terrain patch can be deduced by subtracting the horizontal positions of the matching pair. Two implicit assumptions exist in the above procedure. Firstly, the success of estimation with correlation matching is based on the presence of morphological features in the different survey. Generally, the seafloor seldom preserves exact morphological expressions during different survey periods. However, general shapes of features are often preserved, which make the estimation possible. Secondly, the validity of measuring the similarity with Eq. (2) is based on the premise that there is little bias in depth dimension between multi-temporal DTMs. Nevertheless, this is seldom the case for real data. On the one side, the inaccurate tidal correction could lead to depth bias. On the other side, terrain changes happen during different survey periods. To deal with the depth bias, a modified algorithm, named the profile matching, is a natural choice. In the profile matching, the foreground grid vector {dl, l = 1 ~ n} is replaced by dcl ¼ dl −d

ð4Þ

where d denotes the mean of the foreground grid vector, given by 1 n dl ¼ ∑1 dl . Similarly, the background grid vector is replaced by n zcl ¼ zl −z:

ð5Þ

Therefore, by subtracting the mean depth from the grid vector, the absolute depth information is removed, and only the contour information is used for correlation matching. With the replacement, the influence of the bias in depth can be removed. It is noteworthy that the seafloor terrain is generally nonlinear. In consequence, the distribution of the matching result is seldom Gaussian (often multimodal), which prevents a straightforward interpretation of the matching result. However, as the number of grid points in the window increases, the distribution of the matching result will converge to Gaussian, owing to the central limit theorem (CLT) (Nygren, 2005). Moreover, variance R of the matching result can be approximated by the Cramer–Rao lower bound (CRLB): n o 2 R ¼ E ½^x−x 2 "  #2 XN ∂hi ⌢ x 6 ⌢ 6 l¼1 ∂x i 26 ¼ σe6  i  6 XN ∂hi ⌢ x ∂h ⌢ x 4 ⌢ ⌢ i¼1 ∂x ∂x j i

  i ⌢ 3−1 XN ∂hi ⌢ x ∂h x ⌢ ⌢ 7 l¼1 ∂x ∂x j 7 i 7 " i ⌢#2 7 7 XN ∂h x 5 ⌢ i¼1 ∂x j

ð6Þ

where x denotes the horizontal position of the window [xnorth, xeast], and h(.) the background DTM for searching. Therefore, when the number of the grid nodes in the window is large, the distribution of the matching error can be estimated with Eq. (6).

3.2. Robust correlation matching algorithm Obviously, the above procedure works efficiently under Gaussian noise assumption, because the matching converges CRLB as the number of matching points increases to infinity. However, outlier arises naturally in the terrain data, which violates the assumption. On the one side, outliers are introduced inevitably during the survey process, due to poor data collection and poor survey environment. On the other side, seafloor morphology features may be distorted between consecutive survey, owing to irregular erosion/deposit, sediment gravity flow, or even the sunken boats. All these factors will perform a disturbance on our investigation. Therefore, to estimate the bedform migration, the impact of outliers and terrain distortion has to be excluded. To this end, a robust procedure, termed robust correlation matching (RCM) algorithm is developed in this paper. The objective of the procedure is to estimate the displacement of the seabed in the presence of outliers, based on the robust estimation theory. Many robust estimators have been proposed during the last decades, and the most famous one is probably the M-estimator proposed by Huber (Huber, 2011). When applying robust theory to the correlation matching in the context discussed, two problems should be taken into consideration. Firstly, robust theory accesses the robustness of an estimator with the breakdown point (BP) criterion, which is defined as the largest fraction of outliers that the algorithm can tolerate. Since extensive variation of seafloor may occur between different survey periods, the developed matching algorithm should possess a high BP. Secondly, the estimator should have high efficiency at the same time. Based on the above considerations, the robust matching algorithm with both high BP and simple framework is introduced. To achieve these qualities, the method consists of 3 steps. Firstly, the least median square (LMS) matching (Rousseeuw, 1984), an analog to the correlation matching in Eq. (2), is executed to provide a robust initial matching result. Different from the least square form in Eq. (2), LMS method estimates the displacement based on LMS rule: ^ 0 ¼ arg min medianðdcl −zcl Þ2 D

ð7Þ

By applying Eq. (7), LMS method possesses very high BP (up to 0.5) but poor efficiency, which makes it a suitable tool for outlier diagnose (Rousseeuw, 1984). During the computation, it is noteworthy that the mean values d and z involved in Eqs. (3) and (4) are also very sensitive to outliers. Hence, instead of computing d directly, we estimate its value by the trimmed mean, which is given as: dt ¼

U 1 X dðiÞ U−L i¼Lþ1

ð8Þ

where L = ⌊nδ⌋is the lower trimming bound, U = n − L the upper trimming bound, and δ the trimming ratio (20% for example). By trimming δ/2 percent data on both side, the trimmed mean can give a robust and efficient estimate of the mean value (Olive, 2008). Thus, by replacing d and z with dt and zt , the impact of outliers can be accounted for. Secondly, the contribution of each measurement pair is adjusted according to the initial estimate, so as to remove the impact of the ^ 0 , the residual values vi ¼ hi −zi ðD ^ 0 Þ are potential outlier. Based on D the main information for outlier identification. Under the no-outlier assumption, the distribution of residuals will converge to Gaussian distribution N(0, σv). Thus, when the absolute value of a certain residual is extremely large comparing to σv, it is reasonable to identify the corresponding grid pair as outlier. To handle the influence of outliers, the weight of each measurement pair is adjusted according to a threshold value M. Following the idea of (Gervini and Yohai, 2002), the empirical probability of each large residual is compared with its probability under the distribution N(0, σv).

K. Zhang et al. / Marine Geology 376 (2016) 102–108

To be more specific, for the value of vi which is larger than a certain value η (3σv for example), a measure of the proportion of outlier is defined as: dn ¼ sup fP ðvi Þ−F n ðvi Þgþ vi ≥ η

ð9Þ

where P(vi) is the probability of vi under N(0, σv), Fn(vi) is the empirical probability of vi , and {⋅}+ denotes the positive part. Note that if x (i) ≤ ··· ≤ x (n) are the order statistics of the residuals and i 0 = max {i : x(i) b η}, then     ði−1Þ þ dn ¼ max P vðiÞ − n i N i0

ð10Þ

At this stage, the ⌊ndn ⌋ observations with the largest residuals will be recognized as outliers and eliminated. As a result, the cutoff value is given as t þ ¼ minft : F n ðvi Þ≥1−dn g

ð11Þ

With the cutoff value t+, the potential outliers can be detected. After being identified, their weights are adjusted to zero. Instead of decreasing the weight to a decimal as in M-estimator, or setting it to zero directly when the residual exceeds a given threshold, the above procedure identifies the measurement pair as an outlier not only when its residual exceeds the threshold, but also when the residual is sufficiently larger than the corresponding order statistic. By adaptively calculating the cutoff threshold from the data, the method can achieve full efficiency at the Gaussian distribution, but retain a high BP at the same time (Gervini and Yohai, 2002). Obviously, the premise of the above procedure is the knowledge of σv. Nevertheless, this information is seldom available because the uncertainty information of the bathymetric data is rarely available. To deal with this problem, we estimate σv directly from the residual vectors. To account for the impact of outliers, instead of using sample standard deviation, the mean absolute deviation is used to estimate σv, which is given as: ^ v ¼ 1:4826  medianðabsðvi −medianðvi ÞÞ σ

ð12Þ

^ v , the Note that the median of the residual vector is 0. Thus, with σ outliers can be identified based on the adaptive trimming rule. Thirdly, the ordinary correlation matching (OCM) in Eq. (2) is executed again with the updated weights, and the improved matching result can be computed. By iterating weight updating based on previous estimated position and position updating based on updated weights until convergence, the final matching result is obtained. It is noteworthy that the matching result of the LMS matching in step 1 may be multimodal, because of the inherent nature of LMS estimation as well as correlation matching. To deal with this problem, it is necessary to track multiple local minimum values in LMS matching result simultaneously, so as to grantee the global minimum of the final result. For example, suppose that there are m minimum values in the LMS result. We execute the outlier detection based on every local minimum independently. Thus, after the correlation matching with the updated weights, we will have m matching results. To find the global minimum from these results, the residual vector {vi, i = 1 ~ n} is used as the indicator. Define the weighted sum of the residual square as: R¼

n n X w  ðvi Þ2 n−nt i¼1 i

105

Finally, the procedure of the robust correlation matching is summarized as: 1. Utilize LMS matching to calculate the initial value D0 based on Eq. (7); 2. Calculate the cutoff value based on the estimate D0, and adjust the weights accordingly; 3. Execute the ordinary correlation matching with the updated weights to obtain the updated estimate. 4. Iterate step 2 and step 3 until convergence. 5. Execute steps 2–4 on every local minimum in the LMS matching result simultaneously, and choose the matching result with the smallest R value in Eq. (13) as the final matching result. 4. Results and analysis In the computation, a first attempt is to estimate the migration of sandwaves with the ordinary weighted correlation matching algorithm based on Eqs. (1)–(3). During the matching, to make a compromise between matching accuracy and resolution, the window size is chosen as 20 grids by 20 grids (60 m by 60 m). Fig. 2 depicts the migration estimates of the correlation matching. As can be seen, for the most part of the area, the migration amplitude of the bedform is relatively small. While for some part of the area, the migration amplitude exceeds 10 m. In general, the main tendency of the bedform movement is in the down-canyon direction. While for a small portion of area, the migration shows an up-canyon tendency, which is an attracting feature when comparing to the neighboring area. To evaluate the quality of the matching result, the mean absolute value (MAV) of the residual vector is used. Clearly, when the matching result in the background is dissimilar to the terrain in the foreground window, MAV will increase. Fig. 3 reveals that MAV changes dramatically when the foreground window is sliding over the region. When the matching results are correlated with MAV, it turns out that high MAV (marked with red circles in Fig. 3) tempts to be related to the large migration amplitude (marked with red color in Fig. 2). Obviously, for a portion of large migration estimates, MAV indicates that their credibility is relatively low. Secondly, to evaluate the performance of RCM presented in Section 3.2, robust correlation matching is executed on the same foreground window as above. Fig. 4 shows the differences between the results of RCM and OCM. Generally, the main migration tendency is similar to the result of OCM. While for some local regions, the substantial differences (exceeds 10 m) show in the RCM results. This

ð13Þ

where {wi, i = 1 ~ n} are the weights derived by the outlier detection, and nt is the number of weights which are adjusted to 0 during the outlier detection. Thus, when multimodality occurs, the matching result with the smallest R value will be recognized as the final result.

Fig. 2. Migration estimates derived by OCM. The red ones denote the migration estimates with large MAV.

106

K. Zhang et al. / Marine Geology 376 (2016) 102–108

Fig. 3. MAV of the matching results derived by OCM. The red circles indicate MAV which are larger than 0.5.

indicates that in these local regions, the impact of outliers is not negligible. To compare the matching quality of RCM relative to OCM, MAV derived by two matching methods are compared (Fig. 5). It is obvious that MAVs of RCM are lower than OCM in general, which indicates the validity of RCM under the outlier environment. A detailed analysis of the relationship between MAV and RCM estimates shows a similar situation to OCM scenario. The bulk of large MAV correspond to large horizontal migration, which makes these large migration suspicious. Obviously, these suspicious migration estimates severally hinder the accurate characterization of the sandwave migration. To account for this problem, the error information of matching result has to be estimated, so as to provide the credibility degree of the matching result. Ideally, under the Gaussian noise environment, the theoretical error distribution of the matching can be computed based on Eq. (6). Thus, the theoretical distribution of the residual vector can be inferred by integrating this information with DTM. By comparing the actual residual vector with its theoretical distribution, the credibility can be assigned to the matching result. Nevertheless, the direct computation with Eq. (6) involves the uncertainty information of the bathymetry data, which are seldom available. Therefore,

Fig. 4. Difference in the migration estimates of RCM and OCM.

Fig. 5. Comparison of MAV derived by RCM and OCM.

instead of estimating the theoretical distribution of the residuals, the empirical distribution of the actual residual vectors is analyzed. For this purpose, the mean value of the residual vector (MRV) is used to access the estimates empirically. Fig. 7 shows the empirical probability distribution of MRV series derived by RCM. As can be seen, the major part of MRV distribution shows a Gaussian shape. While for the tail section, the heavy tail phenomena is obvious, which indicates the existence of the dubious matching results. Based on the empirical distribution information, the adaptive trimming rule described in Section 3.2 is used to identify the suspicious matching results. To be more specific, the mean absolute deviation given in Eq. (12) is used to estimate the stand deviation of MRV. Afterwards, MRV with large absolute values is tested according to Eqs. (9)–(11). For those MRV which fail the test, the corresponding matching results will be regarded as false. The spatial distribution of such migration estimates is depicted in Fig. 6 (marked in red). Although some of these marked results seem reasonable when comparing with the nearby estimates, the corresponding residual vectors are too large to be trusted. Meanwhile, it is noteworthy that some irregular migration estimates pass the test (marked with the blue circle), which indicates the unique behavior of the local sandwave.

Fig. 6. Migration estimates derived by RCM and tested by the adaptive trimming rule. The estimates which don't pass the test are marked in red.

K. Zhang et al. / Marine Geology 376 (2016) 102–108

Fig. 7. Empirical probability distribution of MRV (the probability distribution curve is obtained by the kernel density estimation).

5. Discussion A method for characterizing sandwave migration is presented. The method extends existed methods by addressing the impact of outliers based on the robust estimation theory. Based on the experimental results in Section 3, the following statements can be made. (1) The migration tendency is an inherent property of the sandwave (Besio et al., 2004). During the migration, the shape of the crest changes inevitably. However, the main morphological feature of sandwave is often preserved. For this reason, it is possible to estimate the migration rate of sandwave by comparing the multi-temporal DTMs of the same area. Based on this idea, the correlation matching algorithm is proved to be effective. Nevertheless, the ordinary correlation matching is based on the least square criterion, which is the most efficient one under the Gaussian noise assumption, but extremely sensitive to outliers. Unfortunately, in the problem of estimating the sandwave migration rate, the violation of Gaussian assumption is inevitable. On the one hand, outliers occur naturally in DTM data, due to the hostile environment when collecting bathymetry data. On the other hand, during the migration process, a portion of sandwave may change irregularly, which leads to the distortion of the sandwave features. Therefore, when these effects are significant, the result of the ordinary correlation matching could be degraded seriously. (2) From a statistical perspective, morphological distortion of the sandwave can be regarded as outliers. Thus, to account for the impact of outliers, RCM algorithm is developed. In the matching process, the correlation matching based on LMS criterion is applied to provide a robust initial estimate. LMS matching can tolerate nearly 50% of outliers in the measurements, but has a low efficiency. Hence, the adaptive trimming rule is applied according to the initial LMS estimate, and OCM is executed afterwards to provide a better estimate. Noteworthily, like correlation matching, LMS criterion is accompanied with multimodality. Therefore, it is necessary to trace all of the local minimum values in LMS matching result simultaneously. Afterwards, for each local minimum, the normalized residual vector is computed, so as to find the global minimum. In contrast to some other popular robust estimation routines such as M estimators which adjust the weight of the outlier to be decimal, the adaptive trimming rule directly set it to zero. Thus, the following normalized step of residual values is straightforward, which make the adaptive trimming very suitable for the problem herein discussed.

107

(3) When estimating the sandwave migration with the correlation matching algorithm, outliers do have non-negligible influence due to the irregular morphological change of the sandwave migration. The analysis of the residual vectors reveals that such irregular change is especially distinct in the region where migration amplitude is large. Meanwhile, a careful analysis of the matching result shows that even with RCM, such morphological distortion can't be fully addressed. To account for this problem, the quality controlling of RCM results is necessary. Since the uncertainty knowledge of bathymetry data is rarely available, the empirical analysis is employed for the purpose. In the process, MRV is chosen as an indicator to assess the matching quality. This procedure is somewhat arbitrary as the distributions of matching results of different positions are seldom identical to each other. However, under the similar survey environment, it is reasonable that large MRV indicates a suspicious matching result. Hence, we believe that the procedure still provides a practical way for error analysis, at least in a conservative sense. (4) The primary driving force of sandwave migration is the nearbottom current induced by tide or tumidity difference (Hulscher, 1996; Wynn et al., 2002; Besio et al., 2004). As far as Monterey Canyon is concerned, the predominated part of the strong current (50 cm or higher) is the down-canyon current, as reported in (Xu et al., 2008). Accordingly, the down-canyon direction also dominates the sandwave migration. Migration estimates obtained by RCM also coincident with such expectation. In addition to this main tendency, the upstream migration is also highlighted in the matching result. The similar phenomenon is also discussed in (Besio et al., 2004; Duffy and Hughes-Clarke, 2005; Paarlberg et al., 2009). For this unique migration pattern, more detailed DTM time-serial data (higher survey frequency) will be helpful to understand its behavior and to validate the existed model.

6. Conclusion In this paper, a correlation matching algorithm is developed, so as to provide a statistical framework for the estimation of sandwave migration. By concentrating on the morphological distortion as well as outliers in the data, a robust migration estimation can be obtained. With the empirical quality control, each migration estimate is assigned with a confidence level. Thus, a credible migration vector field of the sandwave migration can be reproduced. The developed algorithm is applied on the serial bathymetry data of Monterey Canyon. The computation result demonstrates that the presented method can provide a better estimation solution when dealing with the impact of above described factors. Consequently, the proposed method not only reveals the detailed migration rate information, but also highlights the unique migration behavior effectively. As suggested by the experimental results, the presented method improves the migration estimation performance, and provides a better characterization of the bedform deformation. Acknowledgment Data used in this study were acquired, processed, archived, and distributed by the Seafloor Mapping Lab of California State University Monterey Bay. This work was supported by the project of the National Natural Science Foundation of China (41506111, 41376108), and the Public science and technology research funds projects of surveying, mapping and geoinformation (201512034). References Besio, G., Blondeaux, P., Brocchini, M., Vittori, G., 2003. Migrating sandwaves. Ocean Dyn. 53, 232–238. http://dx.doi.org/10.1007/s10236-10003-10043-x. Besio, G., Blondeaux, P., Brocchini, M., Vittori, G., 2004. On the modeling of sand wave migration. J. Geophys. Res. 109, C04018. http://dx.doi.org/10.1029/2002JC001622.

108

K. Zhang et al. / Marine Geology 376 (2016) 102–108

Best, J., 2005. The fluid dynamics of river dunes: a review and some future research directions. J. Geophys. Res. 110, F04S02. http://dx.doi.org/10.1029/2004JF000218. Buijsman, M., Ridderinkhof, H., 2008a. Long-term evolution of sand waves in the Marsdiep inlet. I: high-resolution observations. Cont. Shelf Res. 28, 1190–1201. Buijsman, M., Ridderinkhof, H., 2008b. Long-term evolution of sand waves in the Marsdiep inlet. II: relation to hydrodynamics. Cont. Shelf Res. 28 (9), 1202–1215. Cherlet, J., Besio, G., Blondeaux, P., Van Lancker, V., Verfaillie, E., Vittori, G., 2007. Modeling sand wave characteristics on the Belgian Continental Shelf and in the Calais-Dover Strait. J. Geophys. Res. 112. Dorst, L.L., 2004. Survey plan improvement by detecting seafloor dynamics in archived echo sounder survey. Int. Hydrogr. Rev. 5, 49–63. Dorst, L.L., Roos, P.C., Hulscher, S., Lindenbergh, R.C., 2009. The estimation of sea floor dynamics from bathymetric surveys of a sand wave area. J. Appl. Geophys. 3 (3), 97–120. Dorst, L.L., Roos, P.C., Hulscher, S., 2011. Spatial differences in sand wave dynamics between the Amsterdam and the Rotterdam region in the Southern North Sea[J]. Cont. Shelf Res. 31 (10), 1096–1105. Duffy, G.P., Hughes-Clarke, J.E., 2005. Application of spatial cross correlation to detection of migration of submarine sand dunes. J. Geophys. Res. 110 (F4). http://dx.doi.org/10. 1029/2004JF000192. Gervini, D., Yohai, V.J., 2002. A class of robust and fully efficient regression estimators. Ann. Stat. 583-616. Huber, P.J., 2011. Robust Statistics. Springer, Berlin Heidelberg. Hulscher, S.J., 1996. Tidal-induced large-scale regular bedform patterns in a threedimensional shallow water model. J. Geophys. Res. 101 (20), 727–20,744. Hulscher, S.J., Dohmen-Janssen, C.M., 2005. Introduction to special section on marine sand wave and river dune dynamics. J. Geophys. Res. 110. Knaapen, M.A., 2005. Sandwave migration predictor based on shape information. J. Geophys. Res. 110 http://dx.doi.org/10.1029/2004JF000195. Knaapen, M.A., Hulscher, S., Vriend, H.J., Stolk, A., 2001. A new type of sea bed waves. Geophys. Res. Lett. 28 (7), 1323–1326. Kubo, Y., Nakajima, T., 2002. Laboratory experiments and numerical simulation of sediment-wave formation by turbidity currents. Mar. Geol. 192 (1–3), 105–121. Lindenbergh, R.C., 2004. Parameter Estimation and Deformation Analysis of Sand Waves and Mega Ripples, Paper Presented at Second International Conference on Marine Sandwave and River Dune Dynamics (MARID2004). Univ. Of Twente, Enschede, Netherlands. Lurton, X., 2010. An Introduction to Underwater Acoustics. Springer, London.

Nemeth, A.A., Hulscher, S., de Vriend, H.J., 2002. Modelling sand wave migration in shallow shelf seas. Cont. Shelf Res. 22 (18–19), 2795–2806. Nemeth, A.A., Hulscher, S., de Vriend, H.J., 2003. Offshore sand wave dynamics, engineering problems and future solutions. Pipeline Gas J. 230 (4), 67–69. Nemeth, A.A., Hulscher, S., Van Damme, R.M.J., 2007. Modelling offshore sand wave evolution. Cont. Shelf Res. 27 (5), 713–728. Nygren, I., 2005. Terrain navigation for underwater vehicles. Doctoral dissertation. School of Electrical Engineering, Royal Institute of Technology (KTH), Sweden. Olive, D.J., 2008. Applied robust statistics. Preprint M-02-006. http://www.math.siu.edu. Paarlberg, A.J., Dohmen-Janssen, C.M., Hulscher, S.J., Termes, P., 2009. Modeling river dune evolution using a parameterization of flow separation. J. Geophys. Res. 114. http://dx. doi.org/10.1029/2007JF000910. Paull, C.K., Ussler, W., Caress, D.W., Lundsten, E., Covault, J.A., Maier, K.L., Xu, J., Augenstein, S., 2010. Origins of large crescent-shaped bedforms within the axial channel of Monterey Canyon, offshore California. Geosphere 6, 755–774. http://dx. doi.org/10.1130/GES00527.1. Rousseeuw, P.J., 1984. Least median of squares regression. J. Am. Stat. Assoc. 79 (388), 871–880. Smith, D.P., Ruiz, G., Kvitek, R., Iampietro, P.J., 2005. Semiannual patterns of erosion and deposition in upper Monterey Canyon from serial multibeam bathymetry. Geol. Soc. Am. Bull. 117 (9–10), 1123–1133. Smith, D.P., Kvitek, R., Iampietro, P.J., Wong, K., 2007. Twenty-nine months of geomorphic change in upper monterey canyon (2002–2005). Mar. Geol. 236 (236), 79–94. Van Dijk, T.A., Kleinhans, M.G., 2005. Processes controlling the dynamics of compound sandwaves in the North Sea, Netherlands. J. Geophys. Res. 110, F04S10. http://dx. doi.org/10.1029/2004JF000173. Van Dijk, T.A., Lindenbergh, R.C., Egberts, P.J., 2008. Separating bathymetric data representing multiscale rhythmic bed forms: a geostatistical and spectral method compared. J. Geophys. Res. 113 http://dx.doi.org/10.1029/2007JF000950. Van Landeghem, K.J., Baas, J.H., Mitchell, N.C., Wilcockson, D., Wheeler, A.J., 2012. Reversed sediment wave migration in the Irish Sea, NW Europe: a reappraisal of the validity of geometry-based predictive modelling and assumptions. Mar. Geol. 295, 95–112. Wynn, R.B., Piper, D.J.W., Gee, M.J.R., 2002. Generation and migration of coarse-grained sediment waves in turbidity current channels and channel-lobe transition zones. Mar. Geol. 192 (1–3), 59–78. Xu, J.P., Wong, F.L., Kvitek, R., Smith, D.P., Paull, C.K., 2008. Sandwave migration in Monterey submarine canyon, central California. Mar. Geol. 248 (3), 193–212.