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

Authors: Pengyun Chen; Ye Li; Yumin Su; Xiaolong Chen
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

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.


---

## Body

## 1. Introduction

The ever increasing capabilities of the autonomous underwater vehicle (AUV) allow for extended period long-range precision autonomous underwater navigation [1–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 interpolation-reconstruction 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 self-similarity 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.


## 2. 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.Supposez=fx,y is the depth value of regionalized variables, which is a second-order stationary random function. The mathematical expectation is m, the covariance function is ch, and the variation function is γh:
(1)Efx,y=m,ch=Efx,yfx+h,y+h-m2,γh=12Efx,y-fx+h,y+h2.Assumingγh=c0-ch, suppose the estimated value f*x,y of the true value fx,y is the linear combination of n known values fxi,yi in estimated domain:
(2)f*x,y=∑i=1nλifxi,yi.The Kriging-based underwater terrain interpolation method involves finding the weight coefficientsλii=1,2,…,n, which makes f*x,y the unbiased estimation of f*x,y. λii=1,2,…,n meet the following condition:
(3)∑i=1nλi=1.In addition, the estimated variance is as follows:(4)σE2=cx,y,x,y+∑i=1n∑j=1nλiλjcxi,yi,xj,yj-2∑i=1nλicxi,yi,x,y.In order to minimizeσE2, apply the Lagrange multiplier principle:
(5)F=σE2-2μ∑i=1nλi-1.In Function (5), μ is the Lagrange multiplier. Calculate ∂F/∂λi and ∂F/∂μ and set them equal to 0:
(6)∂F∂λi=2∑j=1nλjcxi,yi,xj,yj-2cxi,yi,x,y-2μ=0,∂F∂μ=-2∑i=1nλi-1=0.The Kriging functions can be obtained by Function (6):
(7)∑j=1nλjcxi,yi,xj,yj-μ=cxi,yi,x,y,∑i=1nλi=1.Based on the relationship between the covariance and variation functions, the Kriging functions are also expressed by variation function, as shown in Function (8):
(8)∑j=1nλjγxi,yi,xj,yj+μ=γxi,yi,x,y,∑i=1nλi=1.The Kriging-based underwater terrain interpolation-reconstruction method is as follows.(1)
Extract the local water depth data and corresponding grid numbersxi,yi,zi, i=1,2,…,n, from UDTM and read the grid spacing.(2)
Set the interpolated grid spacing.(3)
Find the mathematical expectationm where the covariance function is ch and the variation function is γh.(4)
Calculate the weight coefficientsλii=1,2,…,n based on Functions (7) and (8).(5)
Calculate the interpolated depth value using Function (2).The Kriging-based underwater terrain interpolation-reconstruction 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 self-similarity features can be extracted and the Kriging based interpolation results can incorporate the nonlinear compensation.


## 3. The Underwater Terrain Interpolation Method Based on Fractal Compensation

### 3.1. FBM of Underwater Terrain Data

Using research on Brownian motion, American mathematicians Mandelbrot and Van Ness proposed the FBM model in 1968 [16–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), and the Band-limited noise accumulation method. Among them, the SRA assumes the underwater terrain is generated by the FBF Z(x,y). The original data (x,y,z) are the points on fractal Brown curved surface generated by FBF. The fluctuation degree of points z(x,y) can determine each part of the area and the H 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 H:
(9)zx,y=fx,y+Δξx,y,
wherein Δξ(x,y)~N(0,(dx2+dy2)H), 0<H<1.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 H reflects the local FBF features.


### 3.2. Extraction of Local Fractal Features and Interpolation Compensation

Fractal features can be represented using two parameters:H (Hurst exponent) and σ (terrain variance) [19]. H 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 byD=n+1-h, where H is the self-similarity parameter, n is the Euclidean dimension in three-dimensional space, n=2, and D=3-H.The function for parameter extraction is(10)Ey=Efx+δx-fxδxH,Ey=2∫0∞y12πσ·exp⁡y22σ2dy=2σ2π.LettingA=2σ/2π, we obtain
(11)Ef(x+δx)-f(x)=AδxH.Taking the logarithm of both sides of Function (11) results in
(12)ln⁡Efx+δx-fx=Hln⁡δx+ln⁡A.It can be seen that there is a linear relationship betweenln⁡Ef(x+δx)-f(x) and ln⁡Δx, so the regression function can be established as
(13)e2=∑δx=δxmin⁡δxmax⁡ln⁡Efx+δx-fx-Hln⁡δx-ln⁡A2.e2 can be obtained by minimizing H 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
(14)e2=∑r=r1rnln⁡Efx+δx,y+δy-fx,y-Hln⁡δx2+δy2-ln⁡A2.In this function,r is the distance between a calculated grid node and nearby points, and n 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 H and initial variance at each node of a UDTM can be obtained.UsingH 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
(15)zx,y=fx,y+Δξx,y,fx,y=∑i=1nλifxi,yi.In Function (15), f(x,y) is the base value obtained using Kriging based interpolation methods, Δξx,y is the normal random increment that defers the distribution of N0,σ02Δx2+Δy2H, 0<H<1, and Δx2+Δy2 is the distance between the interpolation point and the nearest point within the local terrain region. H 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.


### 3.3. The Underwater Terrain Interpolation Method

Based on the above derivation, the steps of FCK based interpolation-reconstruction method are as follows.(1)
Extract the local water depth data and corresponding grid numbersxi,yi,zi, i=1,2,…,n from UDTM and read the grid spacing.(2)
Set the interpolated grid spacing.(3)
Find the mathematical expectationm where the covariance function is ch and the variation function is γh.(4)
Calculate the weight coefficientsλii=1,2,…,n based on Functions (7) and (8).(5)
Calculate the fractal feature numberH and the initial terrain standard deviation σ using Function (14).(6)
Calculate the interpolated depth value using Function (15).The interpolation flow chart is as shown in Figure1.Figure 1
The interpolation flow chart.


## 3.1. FBM of Underwater Terrain Data

Using research on Brownian motion, American mathematicians Mandelbrot and Van Ness proposed the FBM model in 1968 [16–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), and the Band-limited noise accumulation method. Among them, the SRA assumes the underwater terrain is generated by the FBF Z(x,y). The original data (x,y,z) are the points on fractal Brown curved surface generated by FBF. The fluctuation degree of points z(x,y) can determine each part of the area and the H 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 H:
(9)zx,y=fx,y+Δξx,y,
wherein Δξ(x,y)~N(0,(dx2+dy2)H), 0<H<1.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 H reflects the local FBF features.


## 3.2. Extraction of Local Fractal Features and Interpolation Compensation

Fractal features can be represented using two parameters:H (Hurst exponent) and σ (terrain variance) [19]. H 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 byD=n+1-h, where H is the self-similarity parameter, n is the Euclidean dimension in three-dimensional space, n=2, and D=3-H.The function for parameter extraction is(10)Ey=Efx+δx-fxδxH,Ey=2∫0∞y12πσ·exp⁡y22σ2dy=2σ2π.LettingA=2σ/2π, we obtain
(11)Ef(x+δx)-f(x)=AδxH.Taking the logarithm of both sides of Function (11) results in
(12)ln⁡Efx+δx-fx=Hln⁡δx+ln⁡A.It can be seen that there is a linear relationship betweenln⁡Ef(x+δx)-f(x) and ln⁡Δx, so the regression function can be established as
(13)e2=∑δx=δxmin⁡δxmax⁡ln⁡Efx+δx-fx-Hln⁡δx-ln⁡A2.e2 can be obtained by minimizing H 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
(14)e2=∑r=r1rnln⁡Efx+δx,y+δy-fx,y-Hln⁡δx2+δy2-ln⁡A2.In this function,r is the distance between a calculated grid node and nearby points, and n 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 H and initial variance at each node of a UDTM can be obtained.UsingH 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
(15)zx,y=fx,y+Δξx,y,fx,y=∑i=1nλifxi,yi.In Function (15), f(x,y) is the base value obtained using Kriging based interpolation methods, Δξx,y is the normal random increment that defers the distribution of N0,σ02Δx2+Δy2H, 0<H<1, and Δx2+Δy2 is the distance between the interpolation point and the nearest point within the local terrain region. H 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.


## 3.3. The Underwater Terrain Interpolation Method

Based on the above derivation, the steps of FCK based interpolation-reconstruction method are as follows.(1)
Extract the local water depth data and corresponding grid numbersxi,yi,zi, i=1,2,…,n from UDTM and read the grid spacing.(2)
Set the interpolated grid spacing.(3)
Find the mathematical expectationm where the covariance function is ch and the variation function is γh.(4)
Calculate the weight coefficientsλii=1,2,…,n based on Functions (7) and (8).(5)
Calculate the fractal feature numberH and the initial terrain standard deviation σ using Function (14).(6)
Calculate the interpolated depth value using Function (15).The interpolation flow chart is as shown in Figure1.Figure 1
The interpolation flow chart.


## 4. The Underwater Terrain Interpolation-Reconstruction Test

### 4.1. The Interpolation-Reconstruction Test Methods

According to the conclusions of Section3.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. In order 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.1. GWA Method

The GWA is an interpolation method typically used for batch processing multibeam data. It uses a Gaussian function as the interpolation weight function. The form is similar to that used in the IDW method. The weighting function of GWA is(16)ωr=Fexp⁡-rd2,
wherein r is the distance between the insertion point and search data points, F is the normalization factor to guarantee the sum of all weighs is 1, and d is the distance that reduced weight to 1/e of the maximum weight, usually taken to be 1/2 of the search radius.


#### 4.1.2. IDW Method

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 (17). Note that the weighting is a function of distance:
(17)Z^j=∑i=1nzi/hijα∑i=1n1/hijα,hij=dij2+δ2,
wherein hij is the effective distance of grid node j around point i. Z^j is the depth of interpolated node j, Zi 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.


### 4.2. The Results of Interpolation-Reconstruction Test

Figures2(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.Comparison of four interpolation methods in one local terrain area: (a) graph of terrain surface and (b) histogram of mean error statistics.
(a)
(b)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 Figure3, 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.Table 1
Statistical property of interpolation error with different method.
Terrain standard deviation  (m)
0.0936
0.3056
0.615
1.072
1.316
1.535
1.828
2.219
The average error
GWA
0.0334
0.0698
0.1622
0.1623
0.152
0.191
0.1901
0.1889
IDW
0.0478
0.0793
0.2314
0.2758
0.2854
0.3953
0.4557
0.4946
FCK
0.0276
0.0676
0.1332
0.1579
0.1489
0.1849
0.1746
0.1658
The error covariance
GWA
0.0024
0.0038
0.048
0.0407
0.042
0.0854
0.0916
0.0698
IDW
0.0043
0.0106
0.0965
0.1068
0.1135
0.2064
0.3059
0.3629
FCK
0.0016
0.0039
0.0393
0.034
0.0325
0.0796
0.1042
0.0742Comparison between real terrain and interpolation terrain with different local terrain standard deviation.
(a)
σ
=
0.0936
(b)
σ
=
0.3056
(c)
σ
=
0.615
(d)
σ
=
1.072
(e)
σ
=
1.316
(f)
σ
=
1.5348
(g)
σ
=
1.828
(h)
σ
=
2.2188We 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 Figure2, 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.


## 4.1. The Interpolation-Reconstruction Test Methods

According to the conclusions of Section3.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. In order 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.1. GWA Method

The GWA is an interpolation method typically used for batch processing multibeam data. It uses a Gaussian function as the interpolation weight function. The form is similar to that used in the IDW method. The weighting function of GWA is(16)ωr=Fexp⁡-rd2,
wherein r is the distance between the insertion point and search data points, F is the normalization factor to guarantee the sum of all weighs is 1, and d is the distance that reduced weight to 1/e of the maximum weight, usually taken to be 1/2 of the search radius.


### 4.1.2. IDW Method

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 (17). Note that the weighting is a function of distance:
(17)Z^j=∑i=1nzi/hijα∑i=1n1/hijα,hij=dij2+δ2,
wherein hij is the effective distance of grid node j around point i. Z^j is the depth of interpolated node j, Zi 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.


## 4.1.1. GWA Method

The GWA is an interpolation method typically used for batch processing multibeam data. It uses a Gaussian function as the interpolation weight function. The form is similar to that used in the IDW method. The weighting function of GWA is(16)ωr=Fexp⁡-rd2,
wherein r is the distance between the insertion point and search data points, F is the normalization factor to guarantee the sum of all weighs is 1, and d is the distance that reduced weight to 1/e of the maximum weight, usually taken to be 1/2 of the search radius.


## 4.1.2. IDW Method

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 (17). Note that the weighting is a function of distance:
(17)Z^j=∑i=1nzi/hijα∑i=1n1/hijα,hij=dij2+δ2,
wherein hij is the effective distance of grid node j around point i. Z^j is the depth of interpolated node j, Zi 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.


## 4.2. The Results of Interpolation-Reconstruction Test

Figures2(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.Comparison of four interpolation methods in one local terrain area: (a) graph of terrain surface and (b) histogram of mean error statistics.
(a)
(b)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 Figure3, 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.Table 1
Statistical property of interpolation error with different method.
Terrain standard deviation  (m)
0.0936
0.3056
0.615
1.072
1.316
1.535
1.828
2.219
The average error
GWA
0.0334
0.0698
0.1622
0.1623
0.152
0.191
0.1901
0.1889
IDW
0.0478
0.0793
0.2314
0.2758
0.2854
0.3953
0.4557
0.4946
FCK
0.0276
0.0676
0.1332
0.1579
0.1489
0.1849
0.1746
0.1658
The error covariance
GWA
0.0024
0.0038
0.048
0.0407
0.042
0.0854
0.0916
0.0698
IDW
0.0043
0.0106
0.0965
0.1068
0.1135
0.2064
0.3059
0.3629
FCK
0.0016
0.0039
0.0393
0.034
0.0325
0.0796
0.1042
0.0742Comparison between real terrain and interpolation terrain with different local terrain standard deviation.
(a)
σ
=
0.0936
(b)
σ
=
0.3056
(c)
σ
=
0.615
(d)
σ
=
1.072
(e)
σ
=
1.316
(f)
σ
=
1.5348
(g)
σ
=
1.828
(h)
σ
=
2.2188We 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 Figure2, 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.


## 5. 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 interpolation-reconstruction test.(2)
The target of the fractal compensation method is not limited to Kriging. Theoretically, it can be used for the error of other traditional interpolation methods as well.(3)
How to further improve the terrain interpolation accuracy and interpolation efficiency is still an area worthy of study.



---

*Source: 289521-2014-12-18.xml*