# The Application of a Complex Composite Fractal Interpolation Algorithm in the Seabed Terrain Simulation

Authors: Chongyang Lv; Fei Yu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Seabed terrain modelling is one of the key technologies in the Subsea Environmental Information System, and this system is critical for underwater vehicle path planning. A composite fractal interpolation algorithm based on improved fractional Brownian motion (FBM) and an improved iterative function system (IFS) is proposed in this paper to increase the precision of the seabed terrain model for submarine topography and to account for the complexity and irregularity of fractal properties in each region. The MATLAB simulation experiment showed that fractal properties of the model built by the complex composite fractal interpolation algorithm were closer to real surface features. After calculation analysis, the model built by the complex composite fractal interpolation algorithm, when compared with the model built by the traditional interpolation algorithm or by the single fractal interpolation algorithm, had higher precision and was more suitable for path planning for underwater vehicles.


---

## Body

## 1. Introduction

Three-dimensional seabed terrain modelling is helpful to simulate the real seabed surface and to verify the static path planning technique for underwater vehicles under the sea [1]. It is very difficult to obtain the altitude data of the seabed, which requires substantial financial and material resources; thus, we use an interpolation method to establish a seabed terrain model based on the limited scattered data points [1].D.G. Krige, a South African scientist, first proposed the interpolation algorithm [2]. The Kriging interpolation algorithm has higher precision when compared with the traditional linear interpolation algorithm. This algorithm is still a smooth interpolation method with high geometric accuracy [3]. A surface modelled by the Kriging interpolation algorithm is smooth, whereas the real surface is disorganized and rugged. The model constructed with the Kriging interpolation algorithm hides the most important features of the real surface [4, 5]. Therefore, it is necessary to find an interpolation algorithm that can show real characteristics of the surface.In the 1970s, Mandelbrot, a famous mathematician, proposed fractal theory. The real surface can be approximated using fractal geometry; therefore, fractal interpolation is an effective tool [6]. Although the fractal interpolation algorithm based on fractional Brownian motion (FBM) has a single algorithm parameter and strong randomness, presently, it is most commonly used [7, 8]. M. F. Barnsley, an American mathematician, proposed the fractal interpolation algorithm based on IFS. The terrain model established by this fractal interpolation algorithm solves the previous problem well [9]. However, the fractal interpolation algorithm based on IFS cannot directly process the discrete data; it can only build a terrain model based on regular grid data [9].The original elevation data points of the seabed are scattered. This paper presents a composite fractal interpolation algorithm based on the improved FBM and IFS to establish a seabed terrain model that can reflect the real surface characteristics. The first step was to perform interpolation on the original elevation data of the seabed using the fractal interpolation algorithm of the improved FBM and then to obtain the preliminary seabed terrain model by improving the selection method of random quantity. The next step was to obtain regular grid data by interpolation processing to control the general direction and outline of the terrain. Finally, a local self-similar final seabed terrain model was established using regular grid data based on the improved IFS fractal interpolation.


## 2. Subsea Environmental Information System

The Subsea Environmental Information System is not only the application of GIS in the marine field but is also the result of the integration of science and earth science information [10]. This system contains various geological phenomena, geophysical field characteristics, and distribution of geographical features and attributes. This system is the embodiment of strategic thinking [11]. This system contains three core contents including the seabed information database, a GIS theoretical basis, and a technology supporting system [11] (Figure 1).Figure 1
Core content diagram.The structure of this system can be divided into three layers including the database layer, the application system layer, and the user terminal [12]. The database layer mainly refers to the standard seabed information database that is integrated. The application system layer is the centre of data scheduling, query, statistics, and access, as well as the establishment of the system model based on the seabed information database. The user terminal contains three levels including maintenance layer users, professional layer users, and browsing layer users [12].This paper focuses on the establishment of the system model in the application system layer. The marine environment is complex and dynamic, and the seabed terrain undulates. The difficulty of measuring ocean data has increased dramatically, and the available baseline data is limited and scattered [13]. Due to the above factors, it is very difficult to establish an underwater terrain model with greater certainty and to be able to show real seabed terrain features. Satellite imagery of the data area was used in this paper (Figure 2).Figure 2
Satellite imagery of the study area.Meanwhile, the contour map of this area, obtained from the 2009 Chinese electronic chart, can be seen in Figure3.Figure 3
Contour map.


## 3. Simulation Method and Process

### 3.1. Fractal Interpolation and Improvement Based on Fractional Brownian Motion (FBM)

FBM can be generated in many ways, such as the midpoint displacement; this is the simplest and most classic method [14]. The important point of random midpoint displacement method is the determination of the random displacement quantity. In terrain modelling, the size of this random displacement is related to the statistical characteristics of the elevation data, such as how to describe the self-similar parameter of rough terrain or the mean and standard deviation of the elevation data. For terrain changes, the terrain model flattened out from rough when the index increased. The least squares method was used to calculate the magnitude of the random displacement. The calculation formula is as follows: (1)Δi=diH×δ×1-22H-2×gausswhere i is the iterative times of random midpoint displacement and di is the distance between two points after the i iteration.Although the rhombus was the basic unit interpolated, there were “crease problems” in all these measures. “Crease problems” represent a constant common edge between two adjacent units, and the midpoint of this common edge has two different elevation values. A cross-section finally appears on the terrain model [15]. The simplest way to solve this problem is to use a line segment to take the midpoint elevation value. The key point is the connection and selection of the line segments.The original data contained many scattered points, and the best way to manage these data was to connect them using a triangular irregular network (TIN) that can describe terrain surfaces at different levels of resolution [16]. Compared with the grid model, the TIN model can show a complex terrain surface more effectively. Previous research indicates that the Delaunay triangular network is the most prominent in terrain fitting among all TIN [16]. The performance of random midpoint displacement interpolation based on the Delaunay triangular network can achieve better results. The improved random midpoint displacement method not only maintains the advantages of the Delaunay triangular network for fitting topography but also solves the “crease problem.”The selection of the scale-free zones and theH index is relatively unitary. This paper intends to use statistical analysis of the terrain data to determine the range of the scale-free zone. The entire data interval was divided into several communities based on the scale-free zone; then, the H index and δ parameter were determined separately.The method to select the minizone is shown below. For A(x, y), near the centre of the entire interval plane, the elevation value is Z. The influence of point(xi,yi,zi) on point A is related to d, the plane distance between two points. σ is the relevance between two points. When the distance d increased, σ exponentially decreased. The relevance σ=1 when two points are close together and when the elevation of these two points is the same.(2)σ=e-d2The relationship between relevanceσ and distance d is reflected well in this function. For a given threshold ω, we assumed that there is no interaction between two points when the relevance σ<ω.Among all the points(xi,yi) that related to point (x,y), the farthest point was defined as dmax, the midpoint as A(x,y), and the length of the side as 2dmax. A rectangle was drawn to find index H and parameter δ of the data point in this area, and the random displacement of this region was calculated.The sea area coverage of the original data points in this paper was the north latitude 10.11°-29.6° and the east longitude 120.11°-132.8°. The interpolation terrain model in this region was constructed by the random midpoint displacement method, which improved with three iterations (Figure4).Figure 4
Three iterations of the terrain model.The random displacement measured by this method is better than the single random displacement of the whole region. It takes too long to run three or more iterations due to the large quantity of data, the time complexity, and the space complexity. Meanwhile, insufficient memory in the computing platform results in a large error. Too few iterations, however, will result in low resolution of the terrain model and will enlarge the random amplification in this model. Therefore, it is necessary to conduct more detailed modelling based on this model.


### 3.2. Acquisition of Regular Grid Data

The fractal interpolation based on IFS is an important method for further modelling, but it can only handle regular grid data. The processing of data is necessary.The seabed elevation data is obtained after iteration calculation from the preliminary model. These data are interpolated to get the required regular grid data. Triangle interpolation was used in this research and the algorithm steps were as follows:(1)
The scattered data points were found in the grid after dividing the whole data area grid and the grid interval wast.(2)
The grid node is defined as the centre and the Matts, whose side is 2t, as the unit; then, a search was applied for the data point n in the unit.(3)
If n≥3, three points were taken to build a triangle that contains the centre node in Matts and the triangle with the smallest area was chosen; then, the elevation value of the plane made by the interpolation node in this triangle was calculated. The result was the interpolation elevation value of the node.(4)
If n<3, a grid interval was added to the scope of the search unit, the search for data points continued until n≥3, and the interpolation elevation value of the node was calculated.The interpolation process for the triangle is shown in Figure5.Figure 5
Grid interpolation.The blue grid node was within the triangle formed by the three data points. Another interval was added if the data points within Matts were less than 3, until three data points within the triangle were found.Compared with the triangular mesh interpolation algorithm, the search of the grid nodes and the maker in the triangle mesh were omitted in this interpolation algorithm. By finding the triangle through grid nodes, the algorithm was simplified and the programming was easier to implement. Rule grid data was obtained based on the preliminary terrain model through the interpolation algorithm. The range of rule grid data was 22-28 degrees north latitude and 124-130 degrees east longitude. The topographic map by grid data is presented in Figure6.Figure 6
Regular grid diagram.


### 3.3. Fractal Interpolation and Improvement Based on the Iterative Function System (IFS)

The IFS fractal interpolation is based on the self-similarity principle of fractal [17].Regular grid data was previously obtained, andI is the data point set. Compressional transformation was adopted on x, y, and z directions, (xn,ym,zn,m) and is defined as the point set: I;n=0,1,⋯,N;m=0,1,⋯,M and a≤xn≤b,c≤ym≤d. Δx and Δy are step-size and divided into grids, as follows: (3)a=x0<x1<⋯<xN=b(4)c=y0<y1<⋯<yM=dAfter calculation and derivation, the compression transformation of each direction is as follows:(5)ϕnx=xn-1+xn-xn-1x-x0xN-x0φmy=ym-1+ym-ym-1y-y0yM-y0Fn,mx,y,z=en,mx+fn,my+gn,mxy+αn,mz+kn,mAnd(6)gn,m=zn-1,m-1-zn-1,m-zn,m-1+zn,mx0y0-xNy0-x0yM+xNyM-αn,mz0,0-zN,0-z0,M+zN,Mx0y0-xNy0-x0yM+xNyMen,m=zn-1,m-1-zn,m-1-αn,mz0,0-zN,0x0-xN-gn,mx0y0-xNy0x0-xNfn,m=zn-1,m-1-zn-1,m-αn,mz0,0-z0,My0-yM-gn,mx0y0-x0yMy0-yMkn,m=zn,m-en,mxN-fn,myM-αn,mzN,M-gn,mxNyMThe compression transformation ofx and y was determined. The vertical scale factor was defined by adding a free variable α to the compression transformation of z(0≤α<1). The role of the vertical scale factor was to control the surface roughness of the terrain. The more α grows, the rougher the terrain is and the more serrated the surface is [17].Iterative operation was used on these three compression transformations. Based on the contracting mapping principle, we can get a fixed point known as the attractor of the IFS, and the image of this attractor represents the three-dimensional seabed topographic map.When the overall data area is larger, the changes in seabed topography are larger (gentle or rugged). It is not reasonable to use theα value to reflect the whole seabed topography. The model is closer to the real terrain if the whole area is divided into small pieces and then the α value for each small piece is calculated [17].The vertical scale factorα represents the compression ratio for that feature of the terrain change, mapped from the whole area to a smaller range. The coefficient of variation was used to calculate the vertical scaling factor ω for the local plots in this research. The algorithm of ω is shown below.As for(xi,yi,zi), the grid data on a local interval, the least square method was used to fit the linear equation of trend surface.(7)z^=b0+b1x+b2ywhere z^i(xi,yi) was the trend elevation value for each point in the region,which was used to define(8)z^-=∑i=1nz^ixi,yinwherez^- is recorded as the overall average trend in elevation.As is known,zi(xi,yi) is the real elevation value of each grid node and(9)ei=zi-z^iwhich was used to define (10)δ=∑i=1nei2nTake δ as the similar standard deviation of the value of grid node elevation in the region.The vertical scale factorω was considered as(11)ω=δz^-Compared with the selection method that uses deviation alone, the selection method of vertical scale factor relates to the similar standard deviation and mean values for data point elevation values in the region. It not only accounts for characteristics of the grid data but also the geomorphic change trend in the region. The selection method of the vertical scale factor can reflect the characteristics of topography changes in the region.


## 3.1. Fractal Interpolation and Improvement Based on Fractional Brownian Motion (FBM)

FBM can be generated in many ways, such as the midpoint displacement; this is the simplest and most classic method [14]. The important point of random midpoint displacement method is the determination of the random displacement quantity. In terrain modelling, the size of this random displacement is related to the statistical characteristics of the elevation data, such as how to describe the self-similar parameter of rough terrain or the mean and standard deviation of the elevation data. For terrain changes, the terrain model flattened out from rough when the index increased. The least squares method was used to calculate the magnitude of the random displacement. The calculation formula is as follows: (1)Δi=diH×δ×1-22H-2×gausswhere i is the iterative times of random midpoint displacement and di is the distance between two points after the i iteration.Although the rhombus was the basic unit interpolated, there were “crease problems” in all these measures. “Crease problems” represent a constant common edge between two adjacent units, and the midpoint of this common edge has two different elevation values. A cross-section finally appears on the terrain model [15]. The simplest way to solve this problem is to use a line segment to take the midpoint elevation value. The key point is the connection and selection of the line segments.The original data contained many scattered points, and the best way to manage these data was to connect them using a triangular irregular network (TIN) that can describe terrain surfaces at different levels of resolution [16]. Compared with the grid model, the TIN model can show a complex terrain surface more effectively. Previous research indicates that the Delaunay triangular network is the most prominent in terrain fitting among all TIN [16]. The performance of random midpoint displacement interpolation based on the Delaunay triangular network can achieve better results. The improved random midpoint displacement method not only maintains the advantages of the Delaunay triangular network for fitting topography but also solves the “crease problem.”The selection of the scale-free zones and theH index is relatively unitary. This paper intends to use statistical analysis of the terrain data to determine the range of the scale-free zone. The entire data interval was divided into several communities based on the scale-free zone; then, the H index and δ parameter were determined separately.The method to select the minizone is shown below. For A(x, y), near the centre of the entire interval plane, the elevation value is Z. The influence of point(xi,yi,zi) on point A is related to d, the plane distance between two points. σ is the relevance between two points. When the distance d increased, σ exponentially decreased. The relevance σ=1 when two points are close together and when the elevation of these two points is the same.(2)σ=e-d2The relationship between relevanceσ and distance d is reflected well in this function. For a given threshold ω, we assumed that there is no interaction between two points when the relevance σ<ω.Among all the points(xi,yi) that related to point (x,y), the farthest point was defined as dmax, the midpoint as A(x,y), and the length of the side as 2dmax. A rectangle was drawn to find index H and parameter δ of the data point in this area, and the random displacement of this region was calculated.The sea area coverage of the original data points in this paper was the north latitude 10.11°-29.6° and the east longitude 120.11°-132.8°. The interpolation terrain model in this region was constructed by the random midpoint displacement method, which improved with three iterations (Figure4).Figure 4
Three iterations of the terrain model.The random displacement measured by this method is better than the single random displacement of the whole region. It takes too long to run three or more iterations due to the large quantity of data, the time complexity, and the space complexity. Meanwhile, insufficient memory in the computing platform results in a large error. Too few iterations, however, will result in low resolution of the terrain model and will enlarge the random amplification in this model. Therefore, it is necessary to conduct more detailed modelling based on this model.


## 3.2. Acquisition of Regular Grid Data

The fractal interpolation based on IFS is an important method for further modelling, but it can only handle regular grid data. The processing of data is necessary.The seabed elevation data is obtained after iteration calculation from the preliminary model. These data are interpolated to get the required regular grid data. Triangle interpolation was used in this research and the algorithm steps were as follows:(1)
The scattered data points were found in the grid after dividing the whole data area grid and the grid interval wast.(2)
The grid node is defined as the centre and the Matts, whose side is 2t, as the unit; then, a search was applied for the data point n in the unit.(3)
If n≥3, three points were taken to build a triangle that contains the centre node in Matts and the triangle with the smallest area was chosen; then, the elevation value of the plane made by the interpolation node in this triangle was calculated. The result was the interpolation elevation value of the node.(4)
If n<3, a grid interval was added to the scope of the search unit, the search for data points continued until n≥3, and the interpolation elevation value of the node was calculated.The interpolation process for the triangle is shown in Figure5.Figure 5
Grid interpolation.The blue grid node was within the triangle formed by the three data points. Another interval was added if the data points within Matts were less than 3, until three data points within the triangle were found.Compared with the triangular mesh interpolation algorithm, the search of the grid nodes and the maker in the triangle mesh were omitted in this interpolation algorithm. By finding the triangle through grid nodes, the algorithm was simplified and the programming was easier to implement. Rule grid data was obtained based on the preliminary terrain model through the interpolation algorithm. The range of rule grid data was 22-28 degrees north latitude and 124-130 degrees east longitude. The topographic map by grid data is presented in Figure6.Figure 6
Regular grid diagram.


## 3.3. Fractal Interpolation and Improvement Based on the Iterative Function System (IFS)

The IFS fractal interpolation is based on the self-similarity principle of fractal [17].Regular grid data was previously obtained, andI is the data point set. Compressional transformation was adopted on x, y, and z directions, (xn,ym,zn,m) and is defined as the point set: I;n=0,1,⋯,N;m=0,1,⋯,M and a≤xn≤b,c≤ym≤d. Δx and Δy are step-size and divided into grids, as follows: (3)a=x0<x1<⋯<xN=b(4)c=y0<y1<⋯<yM=dAfter calculation and derivation, the compression transformation of each direction is as follows:(5)ϕnx=xn-1+xn-xn-1x-x0xN-x0φmy=ym-1+ym-ym-1y-y0yM-y0Fn,mx,y,z=en,mx+fn,my+gn,mxy+αn,mz+kn,mAnd(6)gn,m=zn-1,m-1-zn-1,m-zn,m-1+zn,mx0y0-xNy0-x0yM+xNyM-αn,mz0,0-zN,0-z0,M+zN,Mx0y0-xNy0-x0yM+xNyMen,m=zn-1,m-1-zn,m-1-αn,mz0,0-zN,0x0-xN-gn,mx0y0-xNy0x0-xNfn,m=zn-1,m-1-zn-1,m-αn,mz0,0-z0,My0-yM-gn,mx0y0-x0yMy0-yMkn,m=zn,m-en,mxN-fn,myM-αn,mzN,M-gn,mxNyMThe compression transformation ofx and y was determined. The vertical scale factor was defined by adding a free variable α to the compression transformation of z(0≤α<1). The role of the vertical scale factor was to control the surface roughness of the terrain. The more α grows, the rougher the terrain is and the more serrated the surface is [17].Iterative operation was used on these three compression transformations. Based on the contracting mapping principle, we can get a fixed point known as the attractor of the IFS, and the image of this attractor represents the three-dimensional seabed topographic map.When the overall data area is larger, the changes in seabed topography are larger (gentle or rugged). It is not reasonable to use theα value to reflect the whole seabed topography. The model is closer to the real terrain if the whole area is divided into small pieces and then the α value for each small piece is calculated [17].The vertical scale factorα represents the compression ratio for that feature of the terrain change, mapped from the whole area to a smaller range. The coefficient of variation was used to calculate the vertical scaling factor ω for the local plots in this research. The algorithm of ω is shown below.As for(xi,yi,zi), the grid data on a local interval, the least square method was used to fit the linear equation of trend surface.(7)z^=b0+b1x+b2ywhere z^i(xi,yi) was the trend elevation value for each point in the region,which was used to define(8)z^-=∑i=1nz^ixi,yinwherez^- is recorded as the overall average trend in elevation.As is known,zi(xi,yi) is the real elevation value of each grid node and(9)ei=zi-z^iwhich was used to define (10)δ=∑i=1nei2nTake δ as the similar standard deviation of the value of grid node elevation in the region.The vertical scale factorω was considered as(11)ω=δz^-Compared with the selection method that uses deviation alone, the selection method of vertical scale factor relates to the similar standard deviation and mean values for data point elevation values in the region. It not only accounts for characteristics of the grid data but also the geomorphic change trend in the region. The selection method of the vertical scale factor can reflect the characteristics of topography changes in the region.


## 4. Results and Discussion

### 4.1. Simulation Results

As previously shown, after preliminary modelling, the original scattered seabed altitude data in this area was derived from rule grid data; finally, the required altitude data of rule grid was obtained. The final seabed terrain model was obtained from an improved IFS interpolation algorithm (Figure7).Figure 7
Terrain model of composite fractal interpolation.The seabed terrain model was established by traditional linear interpolation and satellite images of this area, as shown in Figure8.Figure 8
Traditional linear interpolation model.When compared with the traditional linear interpolation model or FBM fractal interpolation model, the composite fractal interpolation model, with higher resolution, had more realistic terrain features.


### 4.2. Results Analysis

Contour reconstruction method was used in this paper for the quantitative analysis comparison between the composite fractal interpolation model and the traditional linear interpolation model. Bathymetric charts of the composite fractal interpolation model and the traditional linear interpolation model are presented in Figures9 and 10, respectively.Figure 9
Bathymetric chart of the composite fractal interpolation model.Figure 10
Bathymetric chart of the traditional linear interpolation model.The contour map of this region obtained from the 2009 Chinese electronic chart represents the “truth value” contour line. The error value was calculated by comparing the contour lines for the first two interpolation models at contour lines of -500 m, -1000 m, and -1500 m. The maximum error, mean error, and standard deviation of each group were then calculated (Figure11).Figure 11
Error analysis.Figure11 shows that the maximum error for the two models was not significantly different. Whereas the average error for the composite fractal interpolation model was better than that of the traditional linear interpolation model, the standard deviation for the composite fractal interpolation model was larger than that of the traditional linear interpolation model at -500 m and -1000 m. The main reason for this situation may be that the topography of this region undulates. Partial fractal dimension is used in the composite fractal interpolation model. The topographics are very close to the real surface of the earth, but they resulted in a large error for terrain altitude.


## 4.1. Simulation Results

As previously shown, after preliminary modelling, the original scattered seabed altitude data in this area was derived from rule grid data; finally, the required altitude data of rule grid was obtained. The final seabed terrain model was obtained from an improved IFS interpolation algorithm (Figure7).Figure 7
Terrain model of composite fractal interpolation.The seabed terrain model was established by traditional linear interpolation and satellite images of this area, as shown in Figure8.Figure 8
Traditional linear interpolation model.When compared with the traditional linear interpolation model or FBM fractal interpolation model, the composite fractal interpolation model, with higher resolution, had more realistic terrain features.


## 4.2. Results Analysis

Contour reconstruction method was used in this paper for the quantitative analysis comparison between the composite fractal interpolation model and the traditional linear interpolation model. Bathymetric charts of the composite fractal interpolation model and the traditional linear interpolation model are presented in Figures9 and 10, respectively.Figure 9
Bathymetric chart of the composite fractal interpolation model.Figure 10
Bathymetric chart of the traditional linear interpolation model.The contour map of this region obtained from the 2009 Chinese electronic chart represents the “truth value” contour line. The error value was calculated by comparing the contour lines for the first two interpolation models at contour lines of -500 m, -1000 m, and -1500 m. The maximum error, mean error, and standard deviation of each group were then calculated (Figure11).Figure 11
Error analysis.Figure11 shows that the maximum error for the two models was not significantly different. Whereas the average error for the composite fractal interpolation model was better than that of the traditional linear interpolation model, the standard deviation for the composite fractal interpolation model was larger than that of the traditional linear interpolation model at -500 m and -1000 m. The main reason for this situation may be that the topography of this region undulates. Partial fractal dimension is used in the composite fractal interpolation model. The topographics are very close to the real surface of the earth, but they resulted in a large error for terrain altitude.


## 5. Conclusions

We conclude that good visualization effects and better precision were obtained from the seabed terrain model built by the composite fractal interpolation algorithm. Compared with the traditional linear interpolation algorithm, topographic features of the model built by the composite fractal interpolation algorithm were much closer to the real surface of the earth. This model better reflects the complexity and irregularity of the seabed features and is more suitable for planning paths for underwater vehicles.



---

*Source: 8641471-2018-07-09.xml*