Research on the Model Improvement of a DLA Fractal River Network

Computer simulations are introduced into the river dynamics field to provide a basis for studying fractal river networks. This paper improved on the defects of particle aggregation and growth methods in the standard diffusion-limited aggregation (DLA) model and proposed a fractal algorithm for simulating the evolution of the river network. To correct these defects, this model used watershed geomorphology as the entry point and analyzed the characteristics of runoff classification, aspect probability and contour lines. Thus, the unit-gradient growth calculation method, anisotropic Brownian motion and other related parameter optimization mechanisms (such as probability-position and step size) were established. In addition, the 2D fractal river network was simulated under different conditions. These simulations calculated the bifurcation ratio of the rendering based on Horton’s law, and its value was bounded between the values of natural river channels (3 ~ 5). The standard DLA algorithm was used for simulation, but the bifurcation ratio (±5.2) of the aggregate was higher than the value of the natural river network. Therefore, the improved model can reflect the structural characteristics of the river network more realistically based on quantitative indicators. The evolution mechanism of the natural watershed was evaluated by the cyber-physical systems theory, providing a model base to study the growth and evolution of river networks. The fractal dimension of the water system fitted with each simulated subgraph was based on the box counting method. The value boundary was between 1.457 and 1.747, providing a division standard for the development degree of the watershed geomorphology. Further, according to the relationship between the number of particles and the fractal dimension, the sensitivity test of the model was also analyzed and evaluated, and its sensitivity coefficient gradually increased from 2.75 to 8.38 based on a suitable data-set, which will make the predicted watershed water system closer to its development level. Finally, use of the improved DLA model was suggested for modeling river networks.


I. INTRODUCTION
Since their formation, river networks have been continuously eroded. Due to the acceleration of urbanization [1], the river network morphology has changed dramatically, and intense human activity has made people realize the importance of river networks [2]. For example, a quantitative description of fractal geometric simulations for river networks [3], large-scale dynamic simulation of river networks [4] and statistical models are used to simulate the geometric scale characteristics of self-similar river networks [5], but the simulation of nondirectional trajectories seems to be biased.
The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino .
Fractal self-similar coding and a distance calculation are used for natural river network simulation [6], however, for the prediction of runoff, it is necessary to use the corresponding algorithm to calculate topological distance of the links to determine the width function of the pattern based on the coding method. Combined with the geomorphic instantaneous unit hydrological model for estimation, the calculation method is more complicated. Developing a one-dimensional hydrodynamic model to simulate a modular water system [7], a perception model of a river corridor is established to simulate the dynamic expansion of the river network [8], and a predictive correction solver is proposed for the fast simulation of free water surface in dendritic and looped river networks [9]. However, certain hydrological factors of the above models are more difficult to obtain, such as river sediments, water loss, runoff, and water level values. Therefore, to provide a theoretical basis for the evolution of river networks, the combination of complex parameter characteristics and the physical mechanism of the river's evolution is the key direction for the future development of fractal chaos, water science, and computer science [10].
Diffusion-limited aggregation (DLA), proposed by Witten and Sander in 1978, is a stochastic model suitable for systems with diffusion as the main mode of motion [11]. Because the mathematical modeling of the morphological structure can help to understand the development of complex systems, the emergence of DLA models has attracted great attention from many scholars. For example, Balakirev et al. (2011) increased DLA seed particles and established a modified multicenter DLA model of dendrite growth in the ion beam synthesis geometry under external magnetic field conditions [12].  modified the starting point, trajectory, and stopping conditions of Brownian moving particles in the DLA model based on existing research methods for generating vascular curve structures [13]. In 2019, Sun et al. changed the step-size and angle of particle motion in the DLA model based on the membership function [14]. In 2020, Salcedo-Sanz et al. proposed a multifractal multiresolution structure based on the combination of DLA models and singular attractors of dynamic systems [15].
The improvement in the DLA model was also accompanied by some deficiencies. Failla et al. (2019) proposed a new dendrite growth model in batteries based on improving the anisotropy of DLA particles through theoretical and numerical analysis [16]. This model simulated the growth of dendrites under different electric field distributions, but only 4 directions of motion were defined for the particles. As a model of the fractal growth process, the DLA model provided an effective basis for two-dimensional simulation research in various computer graphics [17]. Guo et al. (2010) combined the DLA model with yarn and fabric simulation related technologies to perform a more realistic simulation of yarn hairiness [18], but the results lacked quantitative evaluation indicators. Simultaneously, the introduction of fractal theory also provided a basis for the 3D modeling of vegetation [19]. Zhao et al. (2012) improved this model based on the shortcomings of the DLA model and used it to simulate the water-borne growth of plant roots [20]. The results showed that the improved model can realistically and delicately simulate the hydrophobic growth morphology of the fibrous root system, further indicating the research direction for root system modeling. However, due to the many influencing factors in the growth process of plants, it is difficult to truly reflect the actual growth of the root system, resulting in a certain gap between the final root morphology and the true morphology. Based on improving the simulation effect for DLA, the optimization of its simulation time was also proposed. Kuijpers et al. (2014) proposed a new method for DLA simulation [21]. This method improved the time scale of the DLA simulation and reduced the number of calculations, but in some cases, the optimization still consumes considerable memory.
As stated, the DLA model can simulate large or local complex growth phenomena [22], but the modeling and simulation of fractal river networks have not been studied. Three-dimensional modeling of a natural scene is also an important application of graphics in the computer image field [23]. However, in this study, according to the formation characteristics of the river network itself, the standard DLA algorithm was used for a modeling simulation based on the growth mode of particle collision and aggregation. The simulation effect presented had the growth characteristics of fractal self-similarity, but for the construction method and the bifurcation ratio, the simulated characteristics were still far from being achieved. Hence, in order to overcome this shortcoming, this study needed to start with the particle aggregation method. Therefore, this paper analyzed the geomorphological characteristics of a river basin through the cyber-physical systems (CPS) modeling theory. The standard DLA model improved accordingly, and a fractal river network simulation algorithm combination of unit gradient growth, anisotropic motion and a parameter optimization mechanism was proposed. The fractal river network under different conditions was simulated separately by setting different parameter mechanisms for this improved model. This improved model provides a basis for future research on the growth and evolution of river networks in natural watersheds. Applying the idea of this modeling method to introduce the gravity gradient can be explored in the field of 3D modeling. This method will help to further understand the geometry and perform 3D simulations of the tributaries in a river network and study the erosional process of the river, which is of great significance to the theoretical research and the application of river dynamics.

A. THEORETICAL BACKGROUND 1) MATHEMATICAL PROPERTIES OF THE STANDARD DLA MODEL
The standard DLA (diffusion-limited aggregation) model is based on the Laplace [24] equation. Because the Laplace equation has scale invariance and growth probability consistency, it can reveal the relationship between fractals and growth. Scale invariance [25] refers to the symmetry of an object observed at different scales. Among them, the dynamic equation of the model is shown in formula (1) where is the particle number density in the standard DLA model, t is time, and r is a position vector.

2) HORTON'S LAW
Natural water systems follow certain laws during development, which are summarized into some river landform laws, such as Horton's law of river number and the Strahler river classification. The law of river number refers to the inverse geometric relationship between the number of rivers at any level and the VOLUME 8, 2020 level of the river system [26], [27]. The linear fitting formula of this law is The relationship between the regression coefficient k and the bifurcation ratio R B is R B = − log −1 k, where w is the level of the rivers and N w is the number of w rivers at the level. Related research shows [28], [29] that most water systems in the world conform to Horton's law. Water systems all have universal similarities. The bifurcation ratio is affected by geological structure, but most rivers in natural basins are less affected by complex geological structure, which makes the bifurcation ratio tend to be a constant. Hence, according to Horton's research, the bifurcation ratio of natural rivers is between 3 and 5. Consequently, this law is used in this paper to evaluate the simulation results of the model to test whether the simulation effect conforms to the structure of the river network in the natural watershed.

3) FRACTAL DIMENSION
The fractal dimension describes the complexity of the watershed. This paper attempts to test and analyze the sensitivity of the model through the fractal dimension to make the simulation effect more accurate. There are many ways to calculate the fractal dimension [30]. The appropriate calculation method is selected according to the specific situation of the physical scene. This paper applies the box counting [31] to the dimension calculation of two-dimensional fractal graphics.
A suitable square is selected to completely cover the measured figure. Given a positive number ε in the same proportion, ε is used as the side length to mesh the square above. The number of small squares in this area is recorded as N (ε). Finally, least-square fitting [32] is performed on ln ε and ln N (ε), as shown in formula (3) ln where D is the fractal dimension.
In a 2D plane, there is an N × N rectangle and is meshed. A particle is generated in the center of the grid as the seed of the aggregate. This seed particle is used as the center of the circle, and R is used as the radius to determine the circular domain. A particle on the boundary of the circular domain is randomly released so that this particle undergoes Brownian motion [33] in the circular domain. If this particle hits a seed particle, the two particles aggregate; if it crosses a circle with a fixed value from the center seed, it is absorbed by the outer boundary, which releases a new particle at the boundary of the circle. The same process is continuously repeated to obtain enough DLA clusters. R 1 is used as the initial circle domain boundary. R 2 = R 1 + r, r > 0, R max is the maximum boundary of the circular domain. The simulation process [34] is shown in Figure 1.

2) IMPROVING THE DLA MODEL
This paper combines the growth and evolution of fractal river networks [35] to improve the standard DLA model. The improvement mainly lies in 5 aspects: (1) Seed particle selection. River growth may occur at some point in the area.
(2) Selection of the particle area. Randomly generated particles, such as rain, can be generated throughout the area; the range of particle motion also varies according to the area where the river network grows in the natural watershed. (3) Particle aggregation mode. This paper proposes a growth method for particle aggregation, namely, unitgradient growth. (4) Particle motion step size. The step size of the particles varies with different terrain conditions in the area. (5) Brownian motion. Randomly generated particles have different directions and anisotropic probabilities of movement when they follow the aspect probability.
A flowchart was created based on the improved DLA model ( Figure 2).

III. MODEL FACTOR A. SEED POSITION SELECTION
The position of the seed particles is different, so the fractal map generated by the corresponding model is also different. However, the selection of the seed particle position is generally determined by the actual physical scene. The physical scene in this paper is a river network, so the location of the seed particles needs to be determined according to the characteristics of the river network.
The river network in a river basin includes the river sources, nodes, exits, external links, and internal links. The river source refers to the upstream position of the basin water system. The outlet refers to the downstream position of the basin water system. Because many rivers can converge into one river, a river system can have many river sources, but only one exit [36]. Because the formation of the river network gradually converges from the tributary to the main stream, the simulation of the improved DLA model is exactly the opposite, i.e., from the main stream to the tributary. Therefore, the location of the seed particles is the opposite (watershed outlet) of the river source in the basin.

B. PARTICLE AREA SELECTION
Two areas can be used to select the particle area: (1) the particle generation area. (2) the particle motion area. For the standard DLA model, the particle release area is on the boundary of the specified circle, and the motion area is within the predetermined circle. However, according to the factors of river network formation, particles randomly generated on the boundary of the circular domain cannot meet the simulation requirements. Since the formation of the river is simulated, the particles are generated in the entire range of the simulated area, not only on the boundary of the circle. The motion area of the particles can be adjusted according to the shape of the simulation area. Here, the motion area of the particles is set as a rectangular area.
Proceeding from the actual principle of river convergence, the probability that particles are generated in different areas is different according to the altitude. Among them, windward slopes are mostly distributed at high altitudes, while leeward slopes are mostly distributed at low altitudes. Because there is more precipitation on the windward slope and less precipitation on the leeward slope [37], the probability that particles are generated in the windward slope region is stipulated; otherwise, the probability of the leeward slope is small.

C. PARTICLE AGGREGATION MODE
The particle aggregation mode in the standard DLA model is relatively simple: Randomly moving particles in a circular domain will aggregate when they encounter aggregates without touching the boundary. However, when rainfall and glacial meltwater form a river, part of it becomes surface runoff, part becomes underground runoff, and the rest evaporates into the atmosphere. Therefore, the standard aggregation mode cannot meet the needs of this study, and this paper proposes an aggregate growth method: unit-gradient growth. Figure 3(a) shows the unit-gradient growth mode of particles, which is described as follows: U is defined as an aggregate set. U contains u 1 , u 2 , · · · , and u m , which represent the adherent particles and the growing particles. m represents the total number of particles in U . Let T be a Brownian particle and P be a stationary particle that exceeds a limited number of steps. This translates to P after T rest. Use d = (P x − u x ) 2 + (P y − u y ) 2 to calculate the distance d 1 , d 2 , · · · , and d m of all particles u 1 , u 2 , · · · , and u m in P and U and take the smallest record as d min . Find the particle u min closest to P in U through d min so that the growing particles grow in a unit gradient manner. The particle calculation method is shown in Figure 3(b). Let the growth particle be u m+1 and add it to U . Ignore u m+1 for adhesion particles during this period.
In addition, the particle coordinate calculation formula for growing particles is shown in formula (4) Each unknown in formula (4) corresponds to Figure 3(b). R is the radius of the particle. 2R is 1 pixel.

D. PARTICLE MOTION STEP-SIZE
The particle motion steps are the same as those in standard DLA models. However, corresponding to the watershed in a complex environment, the surface runoff is different when it converges into the river. Therefore, the step-size problem is solved according to the contour lines for different river basins.
Contours [38] are the projections of the intersections between the horizontal plane and the actual ground at different altitudes on the horizontal plane. With the same contour distance, the denser the contours, the steeper the slope of the ground; conversely, the less dense the contours, the lower the slope of the ground. Therefore, this paper specifies the step size of particles according to the density of contours.

E. BROWNIAN MOTION 1) ISOTROPIC
Brownian motion is a typical feature in standard DLA models. The particles in this model undergo Brownian motion with an isotropic probability [39]. Let the probability in the up, right, down, and left directions be P u , P r , P d and P l , respectively. Among them, P u = P r = P d = P l = 0.25, P u + P r + P d + P l = 1.

2) ANISOTROPY
According to the geomorphological characteristics for the formation of a river network, the motion that specifies the isotropic probability for the particles is often inappropriate. In this paper, the direction and motion probability of particles are optimized based on the aspect of landform features.
The aspect [40] is the angle between the normal projection of any point on the ground surface perpendicular to the slope and the true north direction is on the horizontal plane. The angle ranges from 0-360 • . Therefore, according to the classification of the aspect, the particles are required to have 8 directions of movement based on the removal of the plane slope. The probability of each aspect is also different. The Brownian motion of particles with different anisotropy [41] is based on the number of aspects.
According to the analysis in section 3.1 ∼ 3.5, in this study, The basic parameters (Table 1) of the model factors of the standard and improved DLA models can be summarized in advance, with the parameters adjusted based on actual needs. Therefore, the quantitative parameters based on simulation and real river network maps are compared and evaluated to compare their performance so that the method can be applied on a larger practical scale.

A. RESULTS AND ANALYSIS OF STANDARD DLA MODELS
Analysis of the standard DLA model and its simulation effect is shown in Figure 4.
Because the seed particles of the standard DLA model are at the center of the circle (Figure 4(a)), the generated fractal graph can only be one type. Fractal self-similar [42] diffusion occurs from the center of the circle to the specified circle. However, to make reasonable improvements to the model, the initial optimization of the model in this paper is to change the seed particle coordinates according to the different water outlets of each river network basin without changing the original particle aggregation method in an attempt to simulate the starting point of the river source to simulate the river network. The simulation effect is shown in Figure 4(b). Changing the position coordinates of the seed particles in this model is not possible. Consequently, before analyzing other geomorphic factors to establish a related optimization mechanism, this article first needs to improve the aggregation mode of moving particles. According to the different coordinates of the seed particles, this improvement makes the simulated river network with water outlets of different  watersheds similar to the real river network in its construction method and related quantitative parameters. Therefore, this study combines the CPS theory to classify the runoff of the watershed to introduce unit gradient growth to change the particle aggregation method. Figure 4 is similar to the structural features of the self-similar multifractal river network [43], and river networks are often studied by this model as complex phenomena in the fractal growth process [17]. However, due to different environmental factors in the evolution of each fractal river network, the model will also optimize the other factors in the model in combination with environmental factors based on the introduction of a unit-gradient growth method. Subsequently, this model will also be based on the modification of limited conditions, which can simulate different shapes of river networks in different environments. Simultaneously, relevant research shows [44] that the bifurcation ratio A of the standard DLA model is 5.2 ± 0.2, which is not limited to the bifurcation ratio of natural river channels between 3 and 5. Therefore, based on the above factors, using this model to simulate fractal river networks in natural watersheds is still insufficient.

B. IMPROVED DLA SIMULATION RESULTS AND ANALYSIS
Through the comprehensive set of (1) seed particle selection; (2) particle area selection; (3) particle condensation mode; (4) particle motion step size; and (5) Brownian motion, fractal river networks of various shapes are generated. Therefore, setting up the improved model according to the difference in different natural conditions can better simulate the river network in the actual watershed.
According to the analysis, based on the subsequent improvement in the standard DLA model, and the reasonable setting of the parameters according to different conditions, the simulated river network diagrams obtained are different ( Figure 5). The simulation subgraph meets the criteria of the river network hierarchical classification and fractal self-similarity. The fractal river network constructed by this improved model will be more reasonable because the cyber-physical systems modeling theory is used to analyze  Figure 5 feature parameter statistics. and combine the geomorphological characteristics of the real watershed, which will achieve a realistic response to the simulation effect. The quantitative description is based on morphological analysis. The bifurcation ratio R B of each simulation subgraph is bounded between 3 and 5 ( Table 2). This R B falls in the range of that of the real river. It is similar in structure to natural rivers, strictly following Horton's tributary order law. The number of different orders of branches tends to obey the inverse geometric series ( Figure 6).
The bifurcation ratio quantitatively reflects the degree of bifurcation of the river system and the degree of landscape development in the basin. The larger the bifurcation ratio, the larger the number of downstream tributaries and the more complicated the bifurcation conditions [45]. A bifurcation ratio between 3 and 5 indicates that although the basin has geological structure, it has little effect on the water system. A bifurcation ratio is less than 3, it indicates that the terrain of the basin is almost flat without any structural influence [46]. A bifurcation ratio greater than 5 indicates that the terrain of the basin is strongly controlled by the local complex geological structure [47]. Comparing the bifurcation ratios of rivers in related studies, it is found that the slope influences the size of the bifurcation ratio [48]. Complicated terrain makes Horton's law deviate [49], and human activities also affect Horton's law [50]. Therefore, when using this model  for simulation, it is important that the bifurcation ratio of the rendering is close to the bifurcation ratio of the natural watershed.
Therefore, according to the different actual conditions, it is reasonable to apply this model to simulate the generation of real fractal river networks with limited modification.

C. FRACTAL DIMENSION RESULTS AND ANALYSIS
The fractal dimension of the river network plays an important role in the watershed and is a key to choosing an appropriate method to calculate the fractal dimension of its water system. First, the box counting method can be used to fit the corresponding river network map (Figure 7).
Second, according to the fitted curve in Figure 7, the fractal dimension of the water system for each simulated river network diagram can be obtained. The fitting formula is shown in formula (3). The parameters in the fitting formula are shown in Table 3. Table 3, shows the corresponding fractal dimension of each river network diagram. The determination coefficient R 2 is an indicator for measuring the fitting effect of the model. If R 2 of each linear fitting formula is greater than 0.99, all are close to 1, and correlation is high, which shows that it is feasible to calculate the fractal dimension of the river network using the box counting method. Simultaneously, the fractal dimension of the water system can also reflect the development degree of the water system [51]. When D ≤ 1.6, the watershed is judged to be in its infancy; when 1.6 < D ≤ 1.9, the watershed is judged to be in its prime; when 1.9 < D ≤ 2.0, the watershed is judged to be in old age. Zakharov et al. (2019) reconstructed the water system based on the digital elevation model, and used its calculated locally changed fractal dimensions (1∼1.53) to explain the formation of riverbeds and valleys [52]. Tunas et al. (2016) verified the stability of the fractal dimension based on the fractal characteristics of the watershed [53], and the fractal dimension (1∼2) of the water system is basically the same, providing variables for the integrated unit hydrological model. Hence, when the model is simulated, the fractal dimension of the simulation rendering is very important. Therefore, the model parameters of the algorithm need to be controlled so that the fractal dimension of the effect map is closer to the fractal dimension of the natural water system evolution.
Finally, based on the above analysis, to make the simulation effect of the model more accurate, this paper analyzes the relationship between the fractal dimension and the number of particles. According to the fractal dimension and particle number of each simulated river network diagram, the functional relationship can be obtained as formula (5). Y = 6.161 × 10 5 X 3 − 2.912 × 10 6 X 2 + 4.589 × 10 6 X − 2.406 × 10 6 (5) In formula (5), X represents the fractal dimension and Y represents the number of particles.
Select the particle fitting method according to the size of the determination coefficient R 2 of the curve. The curve is a cubic fit. From the perspective of quantitative analysis, all other things being equal, set the different number of particles. This process will obtain the correlation coefficient R of the fitted curve, where R is 0. 99943 and is greater than 0. This correlation coefficient indicates that the number of particles is positively related to the fractal dimension. Therefore, sensitivity analysis [54] can be performed on the model based on the one-variable-at-a-time approach (OAT) [55], [56]. The basic principle of OAT is to calculate the change rate of the model output caused by small changes (such as 10% increase or decrease) in each parameter near its best estimate. The absolute value of the change rate represents the sensitivity of the parameter. In this paper, X (1.457, 1.557, 1.697, 1.723, 1.747) in the fitting curve is increased by 10%, and the sensitivity coefficients are 2.75, 3.05, 8.12, 8.25, 8.38, respectively. Consequently, within the appropriate range of the data-set, as the fractal dimension increases, the number of particles increases and the sensitivity of the model also increases. Therefore, the fractal dimensions of the real and simulated river networks are compared using polynomial fitting for the number of particles in the simulated river network and the fractal dimension of the real river network. The number of particles predicted by the river network is calculated according to the fitting formula, which will bring the predicted river system closer to its developmental level.

V. CONCLUSION
This paper used the fractal features of the river network to improve the standard DLA model to simulate the construction of a natural river network and analyze the characteristics of the river network in natural watersheds from the perspective of fusion system modeling. This model added a unit-gradient growth calculation method to optimize Brownian motion and other parameters. Different model factors were set to generate different fractal river networks affected by natural factors. The simulation results between standard DLA models and the improved models were compared. Based on the quantitative evaluation, the improved model simulated the growth structure of the natural river network. This structure truly reflected the fractal features of the river network that were irregular, randomly changing, and self-similar. In addition, the fractal dimensions of river systems in different simulated river network diagrams were calculated using the box counting method, which provided a quantitative criterion for dividing the geomorphological characteristics of the river basin. Fractals can be applied to various fields of science. The improvement of the DLA model in this paper is still in the early stage of analysis. This model needs to be further combined with CPS modeling to extract watershed geomorphological features and improve prediction applications. The research work provides a great opportunity for the development of knowledge such as watershed hydrology, which is beneficial to the future exploration of the hydrodynamic dispersion equation and the diffusion of Brownian motion molecules in Fick's law. This knowledge will enable the understanding of runoff movement and solute migration in groundwater and provide some theoretical basis for the natural system to predict the evolution trend of fractal river networks and early warning of river basin disasters.
HAO JI was born in Shangluo, Shanxi, in 1997. He is currently pursuing the master's degree with the College of Computer Science and Engineering, Northwest Normal University. His main research directions are modeling and simulation of complex systems.
YULIN ZHAN was born in Linxia, Gansu, in 1996. He is currently pursuing the master's degree with the College of Computer Science and Engineering, Northwest Normal University. His main research directions are pattern recognition and image processing.
HONGHONG LI was born in Lvliang, Shanxi, in 1995. She is currently pursuing the master's degree with the College of Computer Science and Engineering, Northwest Normal University.
Her main research directions are information management and technology.
PING LI was born in Jinan, Shandong, in 1996. She is currently pursuing the master's degree with the College of Computer Science and Engineering, Northwest Normal University. Her main research direction is the geographic information systems. VOLUME 8, 2020