Reconstructing High-Precision Coral Reef Geomorphology from Active Remote Sensing Datasets: A Robust Spatial Variability Modified Ordinary Kriging Method

Abstract

Active remote sensing technology represented by multi-beam and lidar provides an important approach for the effective acquisition of underwater coral reef geomorphological information. A spatially continuous surface model of coral reef geomorphology reconstructed from active remote sensing datasets can provide important geomorphological parameters for the research of coral reef geomorphological and ecological changes. However, the surface modeling methods commonly used in previous studies, such as ordinary kriging (OK) and natural neighborhood (NN), often represent a “smoothing effect”, which causes the strong spatial variability of coral reefs to be imprecisely reflected by the reconstructed surfaces, thus affecting the accurate calculation of subsequent geomorphological parameters. In this study, a spatial variability modified OK (OK-SVM) method is proposed to reduce the impact of the “smoothing effect” on the high-precision reconstruction of the complex geomorphology of coral reefs. The OK-SVM adopts a collaborative strategy of global parameter transformation, local residual correction, and extremum correction to modify the spatial variability of the reconstructed model, while maintaining high local accuracy. The experimental results show that the OK-SVM has strong robustness to spatial variability modification. This method was applied to the geomorphological reconstruction of the northern area of a coral atoll in the Nansha Islands, South China Sea, and the performance was compared with that of OK and NN. The results show that OK-SVM has higher numerical accuracy and attribute accuracy in detailed morphological fidelity, and is more adaptable in the geomorphological reconstruction of coral reefs with strong spatial variability. This method is relatively reliable for achieving high-precision reconstruction of complex geomorphology of coral reefs from active remote sensing datasets, and has potential to be extended to other geomorphological reconstruction applications.

1. Introduction

Global environmental changes have posed a strong threat to coral reef ecosystems, with increasing coral reef bleaching and death [1,2,3,4,5,6,7]. Especially in the process of exploitation and utilization of coral reef resources, human beings have caused continuous and strong disturbances to coral reef ecosystems, resulting in a significant trend of degradation of coral reef geomorphology and habitat [8,9,10,11]. Therefore, it is urgent to carry out investigation and exploration of coral reef geomorphology and continuously obtain information on coral reef geomorphology and its changes to provide data and a decision-making basis for coral reef research, development, utilization, and protection.

Remote sensing technology is regarded as the most promising way to carry out large-scale coral reef geomorphology survey and exploration [12,13,14,15,16], mainly including passive remote sensing and active remote sensing. Since there is no dedicated satellite for coral reef monitoring, most of the data used in passive remote sensing are mainly the optical images of airborne aerial platforms or land satellite platforms [17], and they are susceptible to the influence of the atmosphere and the water column [18]. Active remote sensing datasets of multi-beam and lidar are independent of solar radiation sources, which can make up for the limitations imposed by water transparency on passive remote sensing, and has been applied to the survey and evaluation of coral reef geomorphology [19,20,21,22]. Particularly, active remote sensing datasets provide the opportunity to conduct the geomorphological reconstruction of coral reefs. The parameters that characterize the spatial variability of coral reefs, such as surface roughness and morphological complexity, can then be derived, which is of great importance for the study of coral reef geomorphological and ecological changes [23,24,25,26]. Therefore, for coral reef geomorphology with strong spatial variability, a high-precision, spatially continuous geomorphological surface is the key to accurately obtain geomorphological parameters.

At present, the surface modeling methods based on spatial interpolation are the main way to reconstruct the spatially continuous surface model of coral reef geomorphology from discrete lidar and multi-beam datasets [27,28]. The choice of spatial interpolation method often has a direct impact on the accuracy of surface modeling [29,30,31], which in turn will have impacts on the accuracy of various parameters derived from the surface model. In the existing studies, there are some direct empirical choices of kriging [32], natural neighborhood (NN) [33], and triangulation linear interpolation [21] to reconstruct the surface model of coral reef geomorphology, but these lack consideration for geomorphic adaptation. In contrast, other studies have conducted a comparison of the performance of spatial interpolation methods and given an analysis of adaptability for coral reef geomorphology. For example, Leon and Woodroffe [22] argue that the global interpolation method is clearly inappropriate in the geomorphological reconstruction of coral reefs with complex morphological features and that ordinary kriging (OK) is superior to radial basis function (RBF) for a given data distribution. In addition, Coleman et al. [28] compared the accuracy of NN, inverse distance weight (IDW), spline and OK in the reconstruction of coral reef geomorphology when there are gaps in the datasets. They believed that NN is unable to reproduce extreme surface morphologies such as peaks and valleys; the reconstructed surface model is smoother and more natural, with the best performance, followed by OK [28]. Similarly, Leon et al. [34] compared kriging and NN, and believed that NN preserves the range of the original values and generates a smooth and visually effective coral reef geomorphological surface, which is suitable for further geomorphological analysis.

Although NN and kriging have relatively excellent numerical accuracy performance in the studies of coral reef geomorphological reconstruction, it should be noted that coral reef geomorphology has strong spatial variability, especially for the reef flat, with diverse species composition and extremely complex geomorphology. Important parameters such as morphological complexity and roughness derived from the relatively “smooth” surface model generated by NN may not be able to accurately reflect the changes of the coral reef geomorphology or ecology. Similarly, kriging methods, especially OK, as a widely used and reliable interpolation method, lack sufficient adaptability in the studies of coral reef geomorphological reconstruction. Like other interpolation methods based on the least squares standard, kriging methods also have a “smoothing effect” [35]. The presence of the “smoothing effect” results in an underestimation of the spatial variability [35,36]. This “smoothing effect” will have a negative impact on the expression of geographical features with drastic spatial changes [37]. This will affect the accurate derivation of important geomorphological and ecological parameters such as the morphological complexity and roughness of coral reefs.

For the aforementioned problems, the specific objectives of this study are: (1) to develop a robust spatial variability modified OK method to reduce the influence of the “smoothing effect”; (2) to evaluate the robustness of the proposed method on spatial variability modification by experimental comparison with existing spatial interpolation methods; (3) to evaluate its adaptability in the application of active remote-sensing-data-based high-precision geomorphological reconstruction of coral reefs through case comparison with existing spatial interpolation methods.

2. Spatial Variability Modified Ordinary Kriging Method (OK-SVM)

Since OK presents the problem of the underestimation of the spatial variability, we developed an OK post-processing method. In this method, the spatial variability of OK estimation is modified by cooperating with global parameters transformation, local residual correction, and extremum correction, so as to achieve high-precision reconstruction of the surface model of geographic objects with strong spatial variability from discrete sampling point datasets. The method contains a total of four computational processes, as shown in Figure 1.

2.1. Normal Test and Normal Transformation

Before performing ordinary kriging interpolation and spatial variability modification, the normality test of the data distribution should first be performed. Because the theoretical basis of kriging interpolation is to require sampled data to obey a normal distribution, it is thus more sensitive to outliers in skewed data. Moreover, the skewed data distribution used for modification of interpolation estimation may cause larger errors, so the data should be normally distributed as much as possible [36]. Therefore, the significantly skewed original sample point (SP) data should be normalized according to its distribution characteristics using the following transformation (Equation (1)).

$$Z\left(x\right)=F\left(Z_o \left(x\right)\right)$$

(1)

where $Z\left(x\right)$ is the transformed sample point data, $Z_o \left(x\right)$ is the original sample point data, and $F$ is the normal cumulative distribution function transformation.

2.2. Ordinary Kriging Interpolation

The second process is to calculate the value of each estimated point (EP) based on the SPs using OK interpolation. The value of each EP in OK interpolation is the sum of the linear weight of the measured value of all SPs in the given neighborhood (Equation (2)). Its calculation focuses on the solution of weight coefficients, and its estimated value needs to be unbiased and locally optimal. Therefore, OK interpolation must satisfy the condition that the mean value of the EPs is equal to that of the SPs, the sum of the weights is 1 (Equation (3)), and the least square method must be used to minimize the variance of the EPs (Equation (4)). According to the Lagrange multiplier method, Equation (4) can be further derived to obtain the OK interpolation equations (Equation (5)). By solving the (n + 1)-order linear equations, obtaining the kriging weight coefficients, and substituting it into Equation (2), the value of the EPs can be obtained.

$$\widehat{Z\left(x\right)}={\sum }_{i=1}^n{\lambda }_{i}Z\left(x_i\right)$$

(2)

$${\sum }_{i=1}^n{\lambda }_i=1$$

(3)

$${\sigma }_E^2=\overline{cov}\left(x_0,x_0\right)+{\sum }_{i=1}^n{\sum }_{j=1}^n{\lambda }_i{\lambda }_j\overline{cov}\left(x_i,x_j\right)-2{\sum }_{i=1}^n\overline{cov}\left(x_0,x_i\right)$$

(4)

$$\{\begin{array}{c}{\sum }_{j=1}^n{\lambda }_j\overline{cov}\left(x_i,x_j\right)-\mu =\overline{cov}\left(x_0,x_i\right)\\ {\sum }_{j=1}^n{\lambda }_j=1 \end{array}$$

(5)

In the Equations (2)–(5), $\widehat{Z\left(x\right)}$ is the value of the EPs, $Z\left(x_i\right)$ is the value of the SPs in the given range $V\left(i=1,2,3\dots ,n\right)$, ${\lambda }_i$ is the weight coefficient, ${\sigma }_E^2$ is the variance of the EPs, and $\mu$ is the Lagrange multiplier.

Due to the least squares optimization method, the smaller values will be overestimated, and the larger values will be underestimated, which is the so-called “smoothing effect” [38,39]. The “smoothing effect” will also be intensified when the number of sample points in the estimation increases [39]. Moreover, this “smoothing effect” also shows non-uniformity in space, and the more distant the location from the SPs, the greater the “smoothing effect” will be [40].

2.3. Spatial Variability Modification

The reduction of the overall spatial variability after interpolation is essentially a deficit in global accuracy. The complex characteristics of geographic objects need to be reproduced by modifying the estimation results of spatial variability. However, the modification result often leads to a loss of local accuracy; that is, there is always a conflict between ensuring global spatial variability and local accuracy. This indicates that spatial variation modification is a nonlinear process. Therefore, a collaborative strategy of global parameters transformation, local residual correction, and extremum correction is proposed to achieve spatial variability modification without loss of local accuracy by coupling multiple linear processes.

2.3.1. Global Parameters Transformation (GPT)

The OK is an unbiased optimal estimate interpolator. However, in practice, cross-validation results show that the mean value of the EPs is almost equal to that of the SPs, and the variance of the EPs may be smaller than that of the SPs. This indicates that these two are two different sets with inconsistent global parameters (mean value and standard deviation). Therefore, in order to establish the modification relationship between them correctly, it is necessary to transform the two different sets into one set with uniform global parameters.

Z-score transformation (Equation (6)) is a commonly used data standardization method, which can transform two sets of data with different global parameters into a uniform standard score for comparing. Therefore, in the OK interpolation modification method proposed by Rezaee et al. [36] (OK-RM), the Z-score transformation is used for global parameters transformation of the EPs to improve its spatial variability, and the process is simple and easy to perform. This process is the same as OK-RM, which uses the Z-score transformation. The global parameters transformation based on Z-score consists of three steps. Firstly, the global parameters of the SPs and the EPs are calculated separately. Then, the Z-score transformation is performed on the EPs to obtain the standardized score (Equation (7)). Thirdly, the original global parameters of the EPs are replaced by those of the SPs, and then the inverse Z-score transformation is performed to obtain the new EPs (Equation (8)). After the aforementioned processing, the global parameters of the new EPs are consistent with those of the SPs. While ensuring unbiasedness, the spatial variability of the new EPs is significantly improved.

$$Z_{Score}=\frac{X-E\left[X\right]}{\sqrt{Var\left[X\right]}}$$

(6)

where $Z_{Score}$ is the standard score after Z-score transformation, $E\left[X\right]$ and $\sqrt{Var\left[X\right]}$ are the mean value and standard deviation of the variable, respectively.

$$\widehat{Z_s}=\frac{\widehat{Z\left(x\right)}-E\left[\widehat{Z\left(x\right)}\right]}{\sqrt{Var\left[\widehat{Z\left(x\right)}\right]}}$$

(7)

$$\widehat{\widehat{Z\left(x\right)}}=\widehat{Z_s}\times \sqrt{Var\left[Z\left(x\right)\right]}+E\left[Z\left(x\right)\right]$$

(8)

In Equations (7) and (8), $\widehat{Z_s}$ is the standard score of the EPs. $\widehat{\widehat{Z\left(x\right)}}$ is the value after GPT of the EPs. $E\left[\widehat{Z\left(x\right)}\right]$ and $E\left[Z\left(x\right)\right]$ are the mean values of the EPs and the SPs, respectively. $\sqrt{Var\left[\widehat{Z\left(x\right)}\right]}$ and $\sqrt{Var\left[Z\left(x\right)\right]}$ are the standard deviations of the EPs and the SPs, respectively.

2.3.2. Local Residual Correction (LRC)

Since the foregoing process is essentially a global linear transformation, there is still a certain residual error in the spatial variability between the EPs and the SPs in the local neighborhood, which shows some impact on the reproduction of the spatial structure characteristics of regional variables. Therefore, differing from OK-RM, we need to make further corrections to the local residuals.

Firstly, the residuals on the SPs can be obtained from the measured value and the value after GPT (Equation (9)). Then, OK interpolation is performed again on the residuals at the SPs to estimate the residuals at all the EPs. Finally, the local residual correction results can be obtained by removing the residual part from the EPs’ value (Equation (10)). It should be noted here that since increasing the number of sample points involved in the estimation can intensify the “smoothing effect”, the interpolation of the residuals should be performed with a minimum number of local neighbor SPs to avoid the “smoothing effect” [35] on the residual estimation. The four SPs nearest to each EP in the local neighborhood are used in the calculation to avoid the “smoothing effect”.

$$ReErr\left(x_s\right)=\widehat{\widehat{Z\left(x_s\right)}}-Z\left(x_s\right)$$

(9)

$$\stackrel{{\displaystyle ﹀}}{Z\left(x\right)}=\widehat{\widehat{Z\left(x\right)}}-\widehat{\widehat{ReErr\left(x\right)}}$$

(10)

In the Equations (9) and (10), $Z\left(x_s\right)$ and $\widehat{\widehat{Z\left(x_s\right)}}$ are the measured value and the value after GPT on the SPs, $ReErr\left(x_s\right)$ is the residual error on the SPs, $\widehat{\widehat{ReErr\left(x\right)}}$ is the estimated residual at the EPs, and $\stackrel{{\displaystyle ﹀}}{Z\left(x\right)}$ is the value after LRC of the EPS.

2.3.3. Extremum Test and Correction (ETC)

After the aforementioned LRC, it is possible that the values of the EPs fall outside the value range of their neighboring SPs. If so, the value after LRC of EPs also needs to be corrected. This is the other aspect that distinguishes it from OK-RM. It should be noted here that after the foregoing LRC, the global parameters of the EPs have changed. Therefore, differing from Yamamoto’s [35,41] extremum correction treatment, the above-mentioned transformation method needs to be used again to adjust the global parameters of the EPs in order to be able to compare with the SPs in its neighborhood. Then, by comparing with the maximum and minimum values of the SPs in the neighborhood (Equation (11)), the result after the ETC can be obtained (Equation (12)).

$$\Delta =\{\begin{array}{c}\stackrel{{\displaystyle ︾}}{Z\left(x\right)}-Z_{l\_max},\stackrel{{\displaystyle ︾}}{Z\left(x\right)}>Z_{l\_max}\\ \stackrel{{\displaystyle ︾}}{Z\left(x\right)}-Z_{l\_min}, \stackrel{{\displaystyle ︾}}{Z\left(x\right)}<Z_{l\_min}\end{array}$$

(11)

$$\widetilde{Z\left(x\right)}=\stackrel{{\displaystyle ︾}}{Z\left(x\right)}-\Delta$$

(12)

In the Equations (11) and (12), $\stackrel{{\displaystyle ︾}}{Z\left(x\right)}$ is the GPT result of the EPs after LRC, $Z_{l\_max}$ and $Z_{l\_min}$ are the extremum values of the SPs in the neighborhood, ∆ is the extremum difference between the EP and the SPs in the neighborhood, and $\widetilde{Z\left(x\right)}$is the value after ETC of the EPs.

Similar to the previous process, if the global parameters of the EPs deviate significantly from those of the SPs after the ETC, the GPT of the EPs needs to be performed again to obtain the final correction result (Equation (13)).

$$\widetilde{\tilde{Z\left(x\right)}}=\left(\frac{\widetilde{Z\left(x\right)}-E\left[\widetilde{Z\left(x\right)}\right]}{\sqrt{Var\left[\widetilde{Z\left(x\right)}\right]}}\right)\times \sqrt{Var\left[Z\left(x\right)\right]}+E\left[Z\left(x\right)\right]$$

(13)

where $\widetilde{\tilde{Z\left(x\right)}}$ is the final corrected value of the EPs and $E\left[\widetilde{Z\left(x\right)}\right]$ and $\sqrt{Var\left[\widetilde{Z\left(x\right)}\right]}$ are the mean value and standard deviation of the EPs after the ETC.

2.4. Inverse Normal Transformation

In most cases, the implementation objects are subjected to normal transformation to meet the application conditions of OK interpolation and Z-score transformation. Therefore, if the data are normalized before the modification, the result needs to be normalized inversely after the modification is completed to obtain the actual result (Equation (14)).

$$\overline{\overline{Z\left(x\right)}}=F^{-1}\left(\widetilde{\tilde{Z\left(x\right)}}\right)$$

(14)

where $\overline{\overline{Z\left(x\right)}}$ is the modification result of the EPs after the inverse normal transformation and $F^{-1}$ is the inverse normal transformation function.

3. Robustness Test for Spatial Variability Modification

To test the robustness of the proposed modification method in the reproduction of spatial variability, we performed an experimental comparison analysis with OK and a simple OK post-processing method (OK-RM) [36].

3.1. Experimental Data

Two regions with strong spatial variability (A and B) located on an island in China were selected as experimental objects. The raw data (RD) for them are both 52 × 63 raster format DEM (Figure 2) from the State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (http://www.lreis.ac.cn/ (accessed on 30 October 2021)). Sample data (50%) and test data (the remaining 50%) were generated by random sampling from the RD for modeling estimation and verification. The data in region A showed skewed distribution and those in region B showed normal distribution, as shown in Figure 3.

3.2. Test Methods

Firstly, OK, OK-RM, and OK-SVM were used to reconstruct the two experimental objects based on sample points, respectively. Then, the structure difference of spatial variability between the estimated value and the true value at the test points was verified by semi-variogram reproduction. In addition, the residual means square error (RMSE), the slope, and the R-squared (R²) of the regression line between the estimated value and the true value were used to reflect the loss of local accuracy.

3.3. Test Results

The three methods use the same model parameter settings for comparison. The search neighborhood type is circular, with four sectors and 45° offset, the number of samples in the neighborhood is 10 points, and the semi-variogram fitting models are all spherical models. The reconstruction results are consistent with the spatial resolution of the TVs, as shown in Figure 4. The visual comparison shows that OK has an obvious “smoothing effect”. The value range of OK-RM is closer to the RD, although with a large difference in local features, so the results of OK-SVM are more consistent with the RD.

The reproduction results of the semi-variogram on the test point are shown in Figure 5. The semi-variogram of the estimated value of OK differs greatly from that of the true value, which also indicates that the “smoothing effect” causes a large error between the estimated value of OK and the true value in the spatial variation structure. The semi-variograms of OK-RM and OK-SVM are relatively close, but the results of OK-SVM are more consistent with the true value. This indicates that OK-SVM shows better performance in reflecting the real spatial variation structure characteristics of regional variables through spatial variability modification.

The local accuracy evaluation is shown in Figure 6. At the test points with skewed distribution, the results of OK-SVM are very close to OK-RM and OK in terms of R², RMSE, and slope of the regression line. At the test points with normal distribution, the results of OK-SVM outperformed the other two methods. This shows that OK-SVM reduces the “smoothing effect” by modifying the spatial variability, while ensuring the local accuracy of the estimation results with strong robustness.

4. Case Study

4.1. Study Area and Datasets

The study area is located in a typical coral atoll (Figure 7b) in the Nansha Islands of the South China Sea (Figure 7a). An area covering about 15 km² in the north of the region was selected as the study area (Figure 7b). The study area contains most of the geomorphic types of reef flat and lagoon (Figure 7c), including bio-clusters zone, reef pit zone, lagoon slope, lagoon bottom, and patch reef, with strong spatial variability.

The sampling points of geomorphology in the study area were obtained from the State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (LREIS, IGSNRR, CAS). It was acquired by joint measurement of active remote sensing technology based on multi-beam and airborne lidar, with a total of 30,536,463 measurement points. The measured tidal data from the tide gauge stations in the study area were used for tide correction, and the outliers were checked and removed by manual editing. The corrected measurements were in the range of −25.65 m to 1.89 m.

To facilitate the subsequent evaluation of the accuracy and performance of the surface modeling methods, it is necessary to obtain a spatially continuous real geomorphological surface as the standard. However, both measurement sampling and modeling are generalizations of the real geomorphological surface, and it is impossible to obtain data that are completely consistent with the real geomorphological surface. Therefore, the real geomorphological surface is the logical object of expression and description, while the geomorphological generalization surface is the direct object of expression and description. The geomorphological generalization method is adopted here. First, the triangulated irregular network (TIN) surface model was constructed using all the measurement points. It was further converted into a grid DEM with a spatial resolution of 1 meter as the geomorphological generalization result (Figure 7d), which was used to replace the expression of the real geomorphological surface. Then, a random sampling method was used to extract 50% of the total measurement points for the reconstruction of the coral reef geomorphology. Finally, the reconstruction results and the generalization results were compared and analyzed to evaluate the geomorphological adaptability of different surface modeling methods.

4.2. Reconstruction of Coral Reef Geomorphology

The geomorphological surface of the coral reef in the study area was reconstructed based on the OK-SVM method, with a spatial resolution of 1 m. Meanwhile, the OK and NN with excellent performance in the published research were selected for comparison. Among them, the model parameters of OK and OK-SVM are the same, the search neighborhood is circular, with four sectors and 45° offset, the neighborhood search radius is 10 points, and the semi-variogram fitting model is the spherical model. The reconstruction results of the above three methods are shown in Figure 8. The visual comparison results show that the geomorphological surface reconstructed by OK is quite different from the geomorphological generalization surface, while the results reconstructed by NN and OK-SVM are very close to the geomorphological generalization surface, and the expression of local details is finer than that of OK. From the perspective of the bottom depth range, OK-SVM is the closest to the geomorphological generalization surface, followed by NN and OK.

4.3. Comparative Analysis of Accuracy

Geomorphological reconstruction is often influenced by morphological distortion, which is mainly reflected in spatial relationships, geomorphic features, and local morphology. To evaluate the approximation degree between the reconstructed geomorphological surface and the geomorphological generalization surface from the perspective of morphological fidelity, we improved the evaluation index system for the geomorphological accuracy of coral reefs proposed by Wang et al. [42]. From the perspective of morphological fidelity of “point, line, and area“, an index system that considers morphological numerical accuracy (NA) and attribute accuracy (AA) is constructed (

Table 1. Morphological fidelity evaluation index system (the calculation processes are provided in Appendix A).

Category	Indexes	Description and Morphological Meanings
NA	$RMSE\_LE=\sqrt{\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\left(Z^{*}-Z\right)}^2}$	$\mathrm{Where} Z^{*}$$, A^{*}$$, R^{*}$ are the local elevation (LE), local aspect (LA), and local relief (LR) of reconstructed surface at the test points, respectively, and Z, A, R$\mathrm{are} \mathrm{those} \mathrm{of} \mathrm{generalization} \mathrm{surface} \mathrm{at} \mathrm{the} \mathrm{test} \mathrm{points}, \mathrm{respectively}. \mathrm{The} \mathrm{smaller} RMSE\_LE$$, RMSE\_LA$$, \mathrm{and} RMSE\_LR$ are, the closer the reconstructed surface is to the generalization surface, and the higher the morphological numerical accuracy is.
	$RMSE\_LA=\sqrt{\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\left(A^{*}-A\right)}^2}$
	$RMSE\_LR=\sqrt{\frac{1}{n}{\displaystyle {\displaystyle \sum }_{i=1}^n}{\left(R^{*}-R\right)}^2}$
AA	$CR\_LP=\frac{N_{LP}}{N}$	Where N $\mathrm{is} \mathrm{the} \mathrm{number} \mathrm{of} \mathrm{the} \mathrm{test} \mathrm{points} \mathrm{and} N_{LP}$$, N_{LD}$$, N_{LS}$$\mathrm{are} \mathrm{the} \mathrm{changes} \mathrm{numbers} \mathrm{in} \mathrm{the} \mathrm{local} \mathrm{elevation} \mathrm{position} \left(\mathrm{LP}\right), \mathrm{local} \mathrm{direction} \left(\mathrm{LD}\right), \mathrm{and} \mathrm{local} \mathrm{shape} \left(\mathrm{LS}\right) \mathrm{of} \mathrm{the} \mathrm{reconstructed} \mathrm{surface} \mathrm{at} \mathrm{the} \mathrm{test} \mathrm{points}, \mathrm{respectively}. \mathrm{The} \mathrm{smaller} CR\_LP$$, CR\_LD$$, \mathrm{and} CR\_LS$ are, the closer the reconstructed surface is to the generalization surface, and the higher the morphological attribute accuracy is.
	$CR\_LD=\frac{N_{LD}}{N}$
	$CR\_LS=\frac{N_{LS}}{N}$

). The numerical accuracy indexes quantify the approximation degree of morphological features by numerical error, including the RMSE of local elevation (RMSE_LE), the RMSE of local aspect (RMSE_LA), and the RMSE of local relief (RMSE_LR). Attribute accuracy indexes quantify the change degree of morphological meaning by type attribute error, including the change rate of local elevation position (CR_LP), the change rate of local direction (CR_LD), and the change rate of local shape (CR_LS).

The numerical accuracy comparison (Figure 9) shows that the geomorphological surface reconstructed by OK-SVM is higher than OK in terms of local elevation, aspect, and convexity, and has a stronger correlation with the geomorphological generalization surface. Compared with the reconstruction results of NN, OK-SVM also has higher morphological numerical accuracy.

The comparison of attribute accuracy (Figure 10) shows that NN, OK, and OK-SVM have all changed 100% in terms of the attribute accuracy of the elevation location. This shows that no matter which method is used, it is difficult to reproduce the correct attributes of the elevation position. The shape attribute accuracy of OK-SVM is 39%, which is improved compared to NN (41%) and OK (47%). In terms of aspect attribute accuracy, OK-SVM (27%) is significantly higher than NN (75%) and OK (57%). This shows that the spatial variation modification has significantly improved the accuracy of the shape and aspect attributes of the reconstructed geomorphological surface. Comprehensive comparison can be drawn to show that OK-SVM has the highest morphological accuracy in the reconstruction of coral reef geomorphology, and its performance is better than NN and OK, and it can generate a geomorphological surface that is closer to reality.

5. Discussion

As with OK-RM, OK-SVM is also a post-processing method based on OK interpolation. Compared with OK and OK-RM, OK-SVM has an excellent performance in semi-variogram reproduction and local accuracy loss. This method shows high performance on both the reflection of the spatial variability structure characteristics and the local morphological characteristics, which provides a guarantee for the high-precision reconstruction of the coral reef geomorphology. It should be noted that the theoretical regression coefficient between the results after spatial variability modification and the true value should be equal to 1, while it is always close to 1 in practice. The main reason is that the method cannot completely modify its spatial variability to be consistent with the true one, and the results always have certain deviations [35]. The OK-RM proposed by Rezaee et al. [36] has made progress in terms of the simplicity of spatial variability modification compared with the classical post-processing method based on OK [35,41] due to the use of a global parameters transformation method based on Z-score. The OK-SVM proposed in this study also refers to the OK-RM in the spatial variability modification process, which can easily and conveniently improve the spatial variability of OK estimation. Compared with OK-RM, OK-SVM is more robust in spatial variability modification by coupling local residual correction and extremum correction processes, and the spatial variability is more consistent with reality without loss of local accuracy.

In cross-sectional comparison with OK and NN, which have excellent performance, the geomorphological surface reconstructed by OK-SVM is very close to the geomorphological generalization surface of coral reef, and it shows higher morphological numerical accuracy. A longitudinal comparison of the numerical accuracy indexes shows that the correlation between the reconstructed geomorphological surfaces and the geomorphological generalization surface in terms of LR is the lowest for NN, OK, or OK-SVM. The main reason is that the LR, the accuracy indicator of the “area” feature, is affected by the “smoothing effect” deeper than the other two numerical accuracy indicators. However, OK-SVM has relatively less accuracy loss in LR. In addition, OK-SVM has a relatively lower change ratio of morphological attributes. This is very critical, because morphological attributes are usually important parameters for the quantification and evolution of coral reef geomorphology [43,44]. Therefore, the collaborative spatial variability modification strategy of global parameters transformation, local residual correction, and extremum correction makes OK-SVM more adaptable in the geomorphological reconstruction of coral reefs with strong spatial variability. This will provide a guarantee for the subsequent accurate calculation of the multi-scale coral reef geomorphological parameters, which is beneficial to the research of active and passive remote sensing cooperative classification of coral reef geomorphology and analysis of coral reef ecological changes.

OK-SVM is a multi-process coupling method and contains two OK calculations, which will undoubtedly reduce the efficiency of the calculation for large data volume. Therefore, high-performance computing devices can be approached to solve this problem in large data volume computing scenarios. In addition, we need to conduct further research in the following aspects. In terms of the adaptability of high-precision coral reef geomorphological reconstruction, the OK and NN, which have been tested in the field of geomorphological reconstruction of coral reefs, were chosen for comparison. We will compare the adaptability of the proposed method with other methods for coral reef geomorphological reconstruction applications in our future research work. Because of the scale effect of interpolation on the representation of geographic phenomena [45,46], the reduction of spatial resolution will affect the accuracy of spatial variability representation. In this study, a high spatial resolution scale was selected for analysis, in order to lay the foundation for co-interpretation with high-resolution passive remote sensing datasets. Our next step will be to conduct related research in the selection of suitable scales. The lack of measurement data in some areas during coral reef geomorphology monitoring is a practical problem and one that cannot be ignored. Thus, the performance of the method in dealing with such cases needs to be tested in future work. Moreover, it is important to note that the OK-SVM proposed here has the potential to be applied to other surface modeling studies with strong spatial variability characteristics. Therefore, we will also test its performance in other potential modeling application areas in future research work.

6. Conclusions

Reconstruction of spatially continuous and high-precision coral reef geomorphological surfaces based on active remote sensing datasets is the basis for achieving accurate extraction of coral reef geomorphological parameters, which is crucial for coral reef geomorphological and ecological change studies. Surface modeling methods such as OK and NN, which have been frequently used in existing studies, have the problem of low estimation of spatial variability in the application of coral reef geomorphological reconstruction. This can have an impact on the accurate calculation of coral reef geomorphological parameters. To address the problem, an OK-SVM method is proposed in this study. This method adopts a collaborative strategy of global parameters transformation, local residual correction, and extremum correction to modify the underestimated spatial variability of the reconstructed surfaces, while avoiding the loss of local accuracy. The experimental results show that OK-SVM is highly robust in spatial variability modification. Meanwhile, the results of the application case study of coral reef geomorphological reconstruction show that, compared with OK and NN, the OK-SVM performs better in terms of detailed morphological fidelity, the reconstructed surface has higher morphological accuracy, and is more adaptable in the reconstruction of coral reef geomorphology with strong spatial variability.

This method can be used as a reliable tool for reconstructing high-precision coral reef complex geomorphology from active remote sensing data, compared with the commonly used methods in this field. Other potential modeling application areas of the proposed method will also be tested in the future.