The Improved Kriging Interpolation Algorithm for Local Underwater Terrain Based on Fractal Compensation

The interpolation-reconstruction of local underwater terrain using the underwater digital terrain map (UDTM) is an important step for building an underwater terrain matching unit and directly affects the accuracy of underwater terrain matching navigation. The Kriging method is often used in terrain interpolation, but, with this method, the local terrain features are often lost. Therefore, the accuracy cannot meet the requirements of practical application. Analysis of the geographical features is performed on the basis of the randomness and self-similarity of underwater terrain. We extract the fractal features of local underwater terrain with the fractal Brownian motion model, compensating for the possible errors of the Kriging method with fractal theory. We then put forward an improved Kriging interpolation method based on this fractal compensation. Interpolation-reconstruction tests show that the method can simulate the real underwater terrain features well and that it has good usability.


Introduction
The ever increasing capabilities of the autonomous underwater vehicle (AUV) allow for extended period long-range precision autonomous underwater navigation [1][2][3][4].Because of this, the underwater terrain matching navigation (UTMN) has been studied in depth.
In the study of UTMN, it is important to build the accurate underwater terrain matching unit by extracting the local underwater terrain data from the underwater digital terrain map (UDTM).The underwater terrain interpolationreconstruction is the most important step [5] and directly affects the accuracy of terrain matching [6].The underwater terrain real-time interpolation-reconstruction plays a key role in obtaining high-precision local terrain and in the analysis of the statistical features of terrain interpolation accuracy.The traditional terrain interpolation methods include Kriging [7], Gaussian-weighted average (GWA) [8], inverse distance weighting (IDW) [9], and bilinear interpolation [10].These methods do too much smoothing of the original data.The overall underwater terrain trends can be expressed in general, but too many details are lost and the interpolation accuracy is low.Therefore, these methods are difficult to apply in practice.Compared with these other methods, however, the Kriging based method has two distinct advantages [11,12]: (1) it considers the spatial correlation of a described object; (2) it shows the estimation error.Therefore, the Kriging based method gets in-depth research.However, the existing studies are limited to the application of Kriging in related fields, and the defects of Kriging are not corrected.
Numerous studies show that natural terrains have selfsimilarity features [13].Based on the fractal Brownian motion (FBM) model, the fractal compensation Kriging (FCK) method can be used for the interpolation-reconstruction of the underwater terrain.

Kriging Interpolation Method
The Kriging interpolation method is also called space autocovariance best interpolation method.It was proposed by South African geological engineer D. G. Krige and improved by French mathematicians G. Matheron.The core idea behind Kriging is that different sample points are weighted in their importance on the basis of their spatial location as well as their degree of correlation such that the estimation error is minimized.The Kriging method is the unbiased optimal estimation.This method is widely used in calculation of mineral reserves, remote sensing data processing, geology, hydrology, environmental science, agriculture, and forestry science.In these areas, the essence of the interpolation is the same.Using the structural features of the original regional variation data and variation function, do the best linear unbiased estimate for the regional variation values of interpolated points.The Kriging based underwater terrain interpolation method is similar to the other Kriging methods.Select the depth value of local electronic chart area for regionalized variables, determine the grid-node interpolation distance, and interpolate.
Suppose  = (,) is the depth value of regionalized variables, which is a second-order stationary random function.The mathematical expectation is , the covariance function is (ℎ), and the variation function is (ℎ): ( Assuming (ℎ) = (0) − (ℎ), suppose the estimated value  * (, ) of the true value (, ) is the linear combination of  known values (  ,   ) in estimated domain: ( The Kriging-based underwater terrain interpolation method involves finding the weight coefficients   ( = 1, 2, . . ., ), which makes  * (, ) the unbiased estimation of  * (, ).  ( = 1, 2, . . ., ) meet the following condition: In addition, the estimated variance is as follows: 2  =  ((, ) , (, )) In order to minimize  2  , apply the Lagrange multiplier principle: In Function (5),  is the Lagrange multiplier.Calculate /  and / and set them equal to 0: The Kriging functions can be obtained by Function (6): Based on the relationship between the covariance and variation functions, the Kriging functions are also expressed by variation function, as shown in Function (8): The Kriging-based underwater terrain interpolationreconstruction method is as follows.
(3) Find the mathematical expectation  where the covariance function is (ℎ) and the variation function is (ℎ).
The Kriging-based underwater terrain interpolationreconstruction method ignores the nonlinear terrain features so that the interpolated accuracy is relatively low.Therefore, it needs a nonlinear compensation.Studies have shown that natural terrain has self-similarity, or scale independence, and the feature can be expressed by a stationary random process.Fractal theory is often used for digital terrain modeling [14,15].Because of the self-similarity of terrain feature, fractal theory can approximate the terrain better and therefore has theoretical advantage.FBM is the important stochastic process in modern nonlinear time series analysis.FBM can effectively express the nonlinear phenomena in nature, so it is one of the most effective models used to describe the terrain features.Using analysis of fractal features of underwater terrain data based on the FBM principle, the terrain selfsimilarity features can be extracted and the Kriging based interpolation results can incorporate the nonlinear compensation.
The underwater terrain can be described by the random fractal model, but the actual results show that the degree of irregularity of terrain data in different areas has large difference.In other words, the underwater terrain has multifractal nature.Therefore, the fractal features calculated by Function (9) are local fractal features in approaching area.The "local" size can be defined by the UDTM grid-spacing size.Thus the size of  reflects the local FBF features.

Extraction of Local Fractal Features and Interpolation
Compensation.Fractal features can be represented using two parameters:  (Hurst exponent) and  (terrain variance) [19]. is the index of surface complexity or roughness, and  is the index of the adulatory features of local terrain.From the nature of the terrain parameters, we know that a larger Hurst exponent results in a simpler (smooth) terrain, and a smaller Hurst exponent results in a more complex (rough) terrain.We also know that a smaller  results in a smaller terrain undulation, and a larger  results in a larger terrain undulation.
The terrain generated by a fractional Brownian model has self-similarity in the statistical sense, and the fractal dimension can be obtained by  =  + 1 − ℎ, where  is the self-similarity parameter,  is the Euclidean dimension in three-dimensional space,  = 2, and  = 3 − .

The function for parameter extraction is
Letting  = 2/ √ 2, we obtain Taking the logarithm of both sides of Function (11) results in ln  ( ( + ) −  ()) =  ln ‖‖ + ln . ( It can be seen that there is a linear relationship between ln ((+)−()) and ln ‖Δ‖, so the regression function can be established as 2 can be obtained by minimizing  and .
The function is the interpolation in a one-dimensional curve, when the three-dimensional spatial interpolation is performed.Function ( 13) can be rewritten as In this function,  is the distance between a calculated grid node and nearby points, and  is the number of nearby points where the so-called nearby points are the points in a local fractal region.Using Function (14), the fractal feature number  and initial variance at each node of a UDTM can be obtained.
Using  and , the errors of traditional interpolation can be corrected using the fractal principle; using  for the scale factor of height and zone radius, the interpolation of any point can be expressed as In Function (15), (, ) is the base value obtained using Kriging based interpolation methods, Δ(, ) is the normal random increment that defers the distribution of (0,  2 0 (Δ 2 + Δ 2 )  ), 0 <  < 1, and Δ 2 + Δ 2 is the distance between the interpolation point and the nearest point within the local terrain region. is the average value obtained through the weighted average using the terrain features of nearby points.Function (15) is the underwater terrain interpolation function with the fractal correction term, which considers the effect of terrain self-similarity to interpolation points and takes into account correction for the interpolation error.

Expand the search radius
Num ≥ minimum interpolation number Num++ and store the eligible known point

The Underwater Terrain Interpolation
Method.Based on the above derivation, the steps of FCK based interpolationreconstruction method are as follows.
(3) Find the mathematical expectation  where the covariance function is (ℎ) and the variation function is (ℎ).
(5) Calculate the fractal feature number  and the initial terrain standard deviation  using Function ( 14).
The interpolation flow chart is as shown in Figure 1.

The Underwater Terrain Interpolation-Reconstruction Test
wherein  is the distance between the insertion point and search data points,  is the normalization factor to guarantee the sum of all weighs is 1, and  is the distance that reduced weight to 1/ of the maximum weight, usually taken to be 1/2 of the search radius.The distance weighted interpolation method (DW) uses the known points of nearby estimated points to determine the weights arithmetically.In order to obtain the unknown point values, we interpolate using Formula Note that the weighting is a function of distance: wherein ℎ  is the effective distance of grid node  around point .Ẑ is the depth of interpolated node ,   is the water depth near the known point,  is power of the distance, called weight parameters, and  is smoothing factor.DW is a precision interpolation method.When calculating the value of node, weights are assigned to each of the adjacent points.The weights sum up to 1.When there is a value on just one node, the weight is 1 and other weights are zero.This situation causes a "convex change" phenomenon.In calculating the node distance, we therefore introduce the smoothing factor to reduce "convex change" phenomena.
Typically, the value of the grid node is related to the distance, and sometimes it is related to the direction.When discrete points around grid point are evenly distributed, the direction is not considered.Parameters  are generally taken to be 2, and the weighting is calculated by inverse square of the distance.(The method is also called IDW.)The greater the weight parameters , the greater impact of the interpolated point by adjacent points.

The Results of Interpolation-Reconstruction Test.
Figures 2(a), (A), (B), (C), and (D) are the local terrain, FCK interpolation terrain, GMA interpolation terrain, and IDW interpolation terrain, respectively.As can be seen from the partial surface morphology, the FCK based method not only reflects the terrain overall trend but also shows the terrain local details.In comparison, the terrain obtained by GMA based method and IDW based method is relatively smooth.However, there are large differences in terrain local details.This is because GMA based method and IDW based method are approaching-point distance weighting.Both methods fail to consider terrain self-similarity in the local area and irregular features brought by nonlinear terrain; interpolation data are smooth mainly.FCK based method takes into account the terrain self-similarity and irregularities.The interpolated terrain by FCK method can show the terrain overall trend as the traditional methods.More importantly, the local details are similar with real terrain, and reflect the irregular features brought by non-linear terrain.So the FCK method has good usability.
In order to further test interpolation effects of FCK-based method under different terrain undulations as well as the differences with GMA and IDW, we use terrain standard deviation, selecting several different local terrains in UDTM and interpolating.In Figure 3, each figure only gives the contrast effect between FCK terrain and the real terrain.As can be seen from the visual effects, in a variety of terrain conditions, the FCK based method reflects terrain change trends better and has good ability to adapt to terrain features.Table 1 shows under four kinds of interpolation methods the statistical results of interpolation error mean and variance.As can be seen, with the terrain standard deviation increasing, the error mean and error covariance of three methods have certainly increased, and FCK based method is minimum.
We use the real terrain nodes to interpolate in a test so that the accuracy is consistent with the original node accuracy and the error testing standards ensure the consistency.The interpolation results showed that the FCK based interpolation method can reconstruct terrain features well and has good adaptability for different undulating extent terrain, and the average error is smaller than the GMA and IDW based methods.As can be seen from Figure 2, the FCK based method differs from the GMA based method and the IDW based method at a few points, but the average error is better.

Conclusion
For the defects of Kriging based methods when interpolating underwater terrain, the FCK based interpolation method has been proposed.Through extracting local terrain fractal features, this method can compensate for the deficiencies of the Kriging interpolation results.The following conclusions were reached.
(1) The proposed method can effectively improve the accuracy of the terrain interpolation, and this is evidenced in the underwater terrain interpolationreconstruction test.
(2) The target of the fractal compensation method is not limited to Kriging.Theoretically, it can be used for Set the interpolation grid spacing and format the interpolation areaTraverse each point to be interpolated iTraverse each known point jThe known point is in the search radius to the array raw dataEnd the loop jEnd the loop jInvoke the interpolation functions to interpolate and store the results

Figure 2 :
Figure 2: Comparison of four interpolation methods in one local terrain area: (a) graph of terrain surface and (b) histogram of mean error statistics.

Figure 3 :
Figure 3: Comparison between real terrain and interpolation terrain with different local terrain standard deviation.
[16][17][18]ch on Brownian motion, American mathematicians Mandelbrot and Van Ness proposed the FBM model in 1968[16][17][18], which is mainly used to describe irregular shapes in nature, such as mountains, clouds, terrain, and planet surfaces.The FBM based terrain modeling methods include the Midpoint displacement method, the Poisson step method, the Fourier filtering method, the successive random additions (SRA), 3.1.FBM of Underwater Terrain Data.and the Band-limited noise accumulation method.Among them, the SRA assumes the underwater terrain is generated by the FBF (, ).The original data (, , ) are the points on fractal Brown curved surface generated by FBF.The fluctuation degree of points (, ) can determine each part of the area and the  value of each point, so any point value can be decomposed into weight values of neighboring points and normal random fluctuation determined by terrain fractal features :  (, ) =  (, ) + Δ (, ) , [20]rder to illustrate the superiority of the proposed method, the GWA and IDW methods are compared with the FCK method.The GWA and IDW methods are conducted using Surfer 8[20]software.
4.1.The Interpolation-Reconstruction Test Methods.According to the conclusions of Section 3.3, the FCK based interpolation software can be built using C++ programming language, and different underwater terrain is used to execute underwater terrain interpolation-reconstruction tests.In these tests, the real terrain values are provided by 1 m × 1 m grid-spacing UDTM node values.After extracting UDTM fractal terrain features and taking into account the undulating terrain extent, sequentially select 10 m × 10 m size local terrain in eight different grid areas.Run the interpolation tests, respectively, interpolating for each 1 m × 1 m node in local terrain blocks.

Table 1 :
Statistical property of interpolation error with different method.