Evolutionary Programming for Systematic Evaluation of Aquifers: A Case Study from Dholera, Cambay Basin, Gujarat, India ()

Kriti Yadav^{*}, Anirbid Sircar^{}

School of Petroleum Technology, Pandit Deendayal Petroleum University, Gandhinagar, India.

**DOI: **10.4236/gep.2019.74009
PDF HTML XML
290
Downloads
439
Views
Citations

School of Petroleum Technology, Pandit Deendayal Petroleum University, Gandhinagar, India.

Joint inversion of different potentials improves subsurface model resolution. In this paper seismic refraction and magnetotelluric data are used to understand near subsurface features of Dholera, Gujarat, India. An extensive seismic and magnetotelluric survey was carried out in Dholera in order to delineate subsurface presence of aquifers. Ray Inversion for Near Surface Estimation (RINSE) is used for inversion of Dholera seismic data. The inversion output of seismic data is used as seed points for resistivity inversion of anomalies. Inversion of resistivity data is done using evolutionary programing method which is also a type of genetic algorithm. Here the optimization is done using four major steps, of evolutionary programing namely population generation, fitness function, crossover and mutation. This paper also compares the similarities between the natural and geophysical optimization. A Low Velocity Layer is identified up to a depth of 11 m from seismic refraction method. Three layers are identified after the interpretation of seismic and resistivity data. The average thicknesses of Layers one and two are calculated as 3.558 and 6.533 respectively.

Keywords

Evolutionary Programing, Seismic, Dholera, Resistivity, Magnetotelluric

Share and Cite:

Yadav, K. and Sircar, A. (2019) Evolutionary Programming for Systematic Evaluation of Aquifers: A Case Study from Dholera, Cambay Basin, Gujarat, India. *Journal of Geoscience and Environment Protection*, **7**, 139-155. doi: 10.4236/gep.2019.74009.

1. Introduction

Model resolution of subsurface features can be improved by joint inversion of different geophysical potential data (Ammon et al., 1990) . The study of joint inversion can be categorized into two groups; first is joint inversion of datasets which are sensitive to same physical parameters (Julia et al., 2000) and second one is joint inversion of datasets, which is essentially sensitive to different geophysical parameters (Gallardo & Meju, 2003, 2007). Magnetotelluric and seismic refraction techniques are the most effective and commercial methods for identification of aquifers. Inversion of these two potentials can be done by using various inversion methods like Gauss-Newton (GN) method, Quasi-Newton (QN) method, Genetic Algorithm, etc. Evolutionary programing or Genetic Algorithm method is stochastic method based on Darwin’s theory of “Natural selection and survival of the fittest” (Jamshidi & Mostafavi, 2013) . To achieve best solution Genetic Algorithm gives output based on responses obtained from environment and evolution operators. Genetic Algorithm generates near optimal solutions rapidly due to which it is good alternative for non-linear inversion of different geophysical potentials. The smaller size of population generated from premature convergence of nonlinear inversion problems can be avoided by either increasing the size of it or by re-scaling the parameters used (Gallagher et al., 1991; Gallagher & Sambridge, 1994) . It is observed that genetic algorithm approach for geophysical optimization problem is more efficient than other stochastic inversion techniques (Sambridge & Drijkoningen, 1992) . In this paper our study area is Dholera, Gujarat, India; many surveys have been done by co-researchers in this region. It has been found that there is compelling evidence of low enthalpy geothermal sources, which is identified by high gravity and magnetic anomalies in the region and manifestation of many hot water springs in the area (Shah et al., 2017) . These observations motivate our work, and our objective is to test these qualitative approaches by adding more potential data using a formal approach.

In this paper an attempt is made to apply joint inversion for magnetotelluric and seismic method by using Genetic Algorithm. The algorithm starts with the seismic refraction data inversion, where the head waves obtained from seismic refraction method are used to infer the subsurface structure. Inversion results of seismic data are used as a priori assumption for the resistivity. The interchangeability of potential data (velocity and resistivity) obtained from seismic and resistivity is used for evaluation of shallow depth aquifers. In this paper the algorithms for joint inversion of seismic and electromagnetic are coded in C++.

2. Study Area

Dholera area falls under the Saurashtra Peninsula of Gujarat, India which is one of the three conspicuous physiographic divisions of the Gujarat state and lies between 20˚30'N to 22˚30'N latitude and 69˚00'E to 72˚30'E longitude. Dholera is situated in Gulf of Khambat, which lies 30 km south-west of Dhandhuka village in Ahmedabad and 60 km north of the city of Bhavnagar. It is surrounded by water from three sides, on north by Bavaliari creek, on south by Sonaria creek and on east by Gulf of Khambat (Aghil et al., 2014) . Geothermal springs of this area are located along the margins of Saurashtra Peninsula which falls under the vicinity of Western Marginal fault of Cambay Basin (Sharma 2013) . Terrain of Dholera is mainly covered by mudflats, while the basement is formed of Deccan Traps which is at a depth of 500 - 600 m. The area is also occupied by Quaternary soil deposits up to a thickness of about 100 m which is further followed by Tertiary sediment reposed over Deccan Traps.

3. Data Acquisition

This seismic survey has been deployed along different profiles in various orientations. The data is acquired along 4 profile lines, three east to west and one north to south. An array of 24 geophones (indicated by green triangles on the line in Figure 1) along a single profile is performed. A 10 kg gauge with a steel plate is used as source to generate P-waves. The geophone frequency used for refraction seismic is 28 Hz. 7 different shot points (indicated by red dots on the line) are selected at a particular profile. Group interval for seismic data acquisition was considered 2 m. MT survey was performed along 4 MT profiles and data were collected in the frequency range of 0.001 - 10,000 Hz. The orientation of the profiles was in WSW-ENE direction and one normal to three profiles. Measuring array was built for two orthogonal electric profiles namely: Ex, Ey electrical poles of the length 100 m. 3 magnetic sensors recorded two magnetic horizontal components Hx, Hy and one vertical Hz component.

Figure 1. Tectonic framework of Dholera and profile lines for data acquisition.

4. Ray Inversion Method

Inversion process allows us to map data from data space to model space. Estimation and appraisal are the two major steps of inversion, where the estimation stage involves the estimation of model parameters while the appraisal stage evaluate the accuracy of estimated model compared to the true model (Snieder, 1998) . Generally the inversion is categorised as direct and indirect, where the direct inversion outputs the model directly achieved after processing of data while in case of indirect inversion the model is obtained after optimization which minimizes the objective function that involves data set and model set (Weglein et al., 2009) . There are various methods for estimation of subsurface by seismic waves like Waveform Inversion, RINSE, etc. Theory of RINSE was proposed by Jones and Jovanovich (1985), which is used for the near surface estimation by ray inversion of seismic waves. The algorithm for RINSE is written in UNIX environment operating on “HP-UX 7.0”. This technique follows the method of interference lines plotted corresponding to arrivals from different layers. It helps us in understanding of number of layers at near subsurface. After determining the number of layer critical distance Xc is used to calculate the thickness of layer obtained from two way time analysis. Thickness value is determined by projecting the interference ray backward from surface point at a critical distance to the shot point. Intersection of rays exists on the refractor interface. Stripping travel time curves is done once all the depth points are obtained. The observation point for each refractor is assumed horizontal until maximum travel time curves are achieved. This phenomenon of stripping travel time curve is performed for each and every shot on same depth section. Number of trends in travel time curve represents number of layers. Thickness of layers can be calculated by using following equations. Thickness of layer can be calculated by intercept T1 which corresponds to time taken by wave to reach receiver/geophone, the equation can be expressed as:

${T}_{1}=\frac{2h\times \mathrm{cos}{\theta}_{c}}{{V}_{1}}$ (1)

where,

T_{1} = Time taken by wave to reach geophone;

V_{1} = Velocity of wave in first layer;

h = Thickness of layer.

Inversion of Magnetotelluric data by GA does not require forward solution to calculate derivatives of the fields with respect to model parameter changes. In this case Electric field (E) is calculated from Transverse Electric (TE) mode and Magnetic field (B) from Transverse Magnetic (TM) mode. This helps in estimating the apparent resistivity, phases and complex impedances for both TE and TM mode. According to Weaver’s method it has been found that the grids are generated automatically which were prepared on the information obtained from model and frequency (Taylor & Weaver, 1976; Poll, 1994) .

5. Genetic Algorithm

Basic principles of Genetic Algorithm/Evolutionary Programing are proposed by Holland (Holland, 1975) . In genetic algorithm method a class of adaptive algorithm are represented whose search methods are based on the simulation of natural genetics. It falls under the class of probabilistic algorithms. Evolutionary programing/Genetic Algorithm is a process of natural selection where the stronger individuals are the winners. It has been found that in GA the potential solution of problem is an individual which is represented by sets of parameters. The parameters are known as chromosomes which are structured as string in binary forms. A positive is represented as fitness value which reflects the degree of goodness of genes for solving the problems, this value represents the local minima for model set and data set. Good quality offspring is yielded by the fittest chromosome throughout its genetic evolution, which is a better solution to the problem.

In Figure 2 a multi-directional search is performed by genetic algorithm in order to maintain a population of potential solutions and encourage information formation along with exchange between these directions. A number of populated solutions are developed in this simulated evolution in which the relatively “good” solutions reproduce. The different solutions are distinguished on the basis of evaluation function which plays the role of an environment.

Figure 2. Flow diagram of Evolutionary Programing/Genetic Algorithm based inversion (modified after Sircar, 2000 ).

After the individuals are converted into population by determining the fitness function of individual, the simulation further proceeds to crossover selection and mutation. Probability of the simulation needs to be checked whether the termination conditions are satisfied or not. If it meets the desired output than the outputs are the optimization results otherwise it needs to be processed again from the fitness function determination point.

The parameters and steps which will be involved in GA inversion method of optimization are as follows.

5.1. Encoding

The optimization problem variables are represented by the encoding mechanism of GA. Each and every individual parameter in population consists its own genetic code. Particular genes of fixed length represents the velocities distributed over the cells in the cross-hole region of population. It is obtained from extraction of individuals genetic code. The length of genes are denoted by bits, the number of bits for representing velocity has to be given in advance. Velocity value over the cell can be obtained as:

$V=\mathrm{min}V+dv\ast tmp$ (2)

where,

minV = Minimum Velocity;

dv = Velocity partition;

tmp = Value of the binary string;

$dv=\left(\mathrm{max}V-\mathrm{min}V\right)/{2}^{nbit}$ (3)

where,

minV = Minimum Velocity;

nbit = Number of bit.

5.2. Population

The primary set of population is generated randomly and size of population is determined. The two important criteria for generation of population are population size and randomization of seed number (Haupt & Haupt, 2004) . In order to discover new clones (Rezaian et al., 2010) the initial population should be a large pool of different genes. Including different genes for initial population leads to an algorithm which has enough diversity in the population to get fast and good solutions.

5.3. Fitness Function

Fitness function is defined as the ratio of the assessment value of a particular clone to the average assessment of all the clones. The equation for fitness function probability selection (Chipperfield et al., 1994) can be expressed as:

${P}_{i}=\frac{{F}_{i}}{{\displaystyle {\sum}_{i=1}^{n}{F}_{j}}}$ (4)

where,

P_{i} = Fitness probability;

F_{i} = Individual parameters fitness.

In case of geophysical data optimization it is determined by using concept of Chi-square error from the observed and calculated apparent resistivity difference. The Chi-Square error is denoted as “ $\epsilon $ ” and can be expressed as:

$\epsilon =\frac{1}{N}{{\displaystyle \sum}}^{\text{}}\frac{{\rho}_{a}^{Cal}-{\rho}_{a}^{obs}}{{\rho}_{a}^{Cal}}$ (5)

where:

${\rho}_{a}^{Cal}$ = Calculated apperant resistivity;

${\rho}_{a}^{obv}$ = Observed apparent resistivity;

N = Number of iterations.

5.4. Crossover

Crossover is the process of generating better quality genes by exchanging the good information between the particular parents. The crossover probability can be calculated as the ratio of pairs of clones which will be selected for mating to the total number of pairs of clones.

5.5. Mutation

Mutation can generate new genes by flipping one or more gene values randomly in a clone. The mutation probability can be calculated as ratio of the bits to be flipped randomly to the total bits of clones (Thander & Sircar, 2014) . The simple mutation can be performed by using normal distribution, it leads to faster execution as only muted genes are processed. The number of mutations is performed on the basis of random pick in N(m, σ) distribution function.

$m=\left({n}_{pop}-{n}_{elite}\right)\ast {P}_{mut}$ (6)

${\sigma}^{2}=\left({n}_{pop}-{n}_{elite}\right)\ast {P}_{mut}\ast \left(1-{P}_{mut}\right)$ (7)

where,

m = Average number of mutations in the population;

σ = Standard Deviation;

${n}_{pop}$ = Population size;

${n}_{elite}$ = Number of elitist individuals.

Table 1 represents the mutation rate for different range of Chi-square error values. The rate of mutation used in Table 2 is a good trade-off for exploration and exploitation. After the selection of parameter from fitness function a random number is generated using the reference seed value. If the number is greater than list value of “ $\epsilon $ ” then the value of individual mutation is created by multiplying the individuals with probability Pm, otherwise it is estimated by dividing the individual with Pm.

Number of iterations are performed until modified “n” models are arranged

Table 1. Rate of mutation for different Chi-square error range value.

Table 2. Ranges of Velocity and Resistivity values calculated for study area.

from fitness values. The process is repeated until the population reaches high fitness value.

6. Integration of Seismic Refraction and Resistivity Inversion

Seismic refraction data is interpreted by using the first breaks in refraction surveys , the time section is prepared by using amplitudes and first-arrival travel times. Integration of seismic and resistivity data is sorted by low, high and targeted value of seismic refraction. Combination of different algorithms are classified into 6 categories which are described in Table 2. In inversion of seismic and resistivity the velocity is converted into resistivity for particular lithology from seismic database. The thickness of layers is obtained from seismic inversion which is perturbed at a fixed rate and provides upper and lower boundaries for resistivity inversion.

7. Results and Discussion

Seismic data is acquired along four profile lines. The group interval was taken to be 2 m with 4 shot points of shot interval 5 m. Data has been acquired by two direct and two reverse shots. Plot between first arrival time and offset is drawn for each seismogram obtained by varying source position along the profile. Plot is drawn for all four seismic profiles. Using each profile velocity of each individual layer is calculated with the help of inverse slope of ray path formed by joining first arrival time. The seismic data analysis suggest presence of three subsurface layers. Velocity and thickness of each layer is given in Table 3.

7.1. Profile 1

Low Velocity Zone is identified by seismic refraction survey, which need to be neglected in order to bring all the potential data to a common datum. By calculating velocity and thickness at each shot point the near subsurface model is constructed. The velocity calculated from profile 1 for Layer 1, Layer 2 and Layer 3 are 75.69 m/s, 107.15 m/s and 129.97 m/s respectively (Figure 3 & Figure 4). Thickness of Layer 1 and Layer 2 obtained from profile 1 is 4.182 m and 6.617 m respectively, while the thickness of third layer cannot be identified.

The nature of curve in Figure 5 is A-type, boundary condition for A-type curve is ρ1 < ρ2 < ρ3. The value of apparent resistivity for profile 1 is ρ1 = 60.27 Ωm in h1 = 0.36 m, ρ2 = 90.67 Ωm in h2 = 5.02 m, ρ3 = 120 Ωm in h1 = infinity. Chi-square computed for profile 1 is 0.5 at 1000 iteration (Figure 6).

Figure 3. Travel time curves―profile 1.

Figure 4. Subsurface geological layered cross-section obtained from Profile 1.

Figure 5. Apparent resistivity curve for profile 1.

Figure 6. Chi-Square error plot for profile 1.

Table 3. Velocity and thickness of different layers.

7.2. Profile 2

Second profile suggests presence of three layers in the identified subsurface. The velocity calculated from profile 2 (Figure 7 & Figure 8) for Layer 1, Layer 2 and Layer 3 are 82.85 m/s, 116 m/s and 138.18 m/s respectively. Thickness of Layer 1 and Layer 2 obtained from profile 2 is 4.156 m and 5.694 m respectively, while the thickness of third layer cannot be identified.

The nature of curve in Figure 9 is A-type. For profile 2 the apparent resistivity values are ρ1 = 73.27 Ωm and h1 = 0.46 m, ρ2 = 91.67 Ωm and h2 = 4.02 m, ρ3 = 110 Ωm and h1 = infinity. Chi-square computed for profile 1 is 0.005 at 1200 iterations (Figure 10).

Figure 7. Travel time graph―profile 2.

Figure 8. Subsurface geological layered cross-section obtained from profile 2.

Figure 9. Apparent resistivity curve for profile 2.

Figure 10. Chi-square error plot for profile 2.

7.3. Profile 3

Three layer subsurface model is identified by third profile of seismic refraction survey. The velocity calculated from profile 3 (Figure 11 & Figure 12) for Layer 1, Layer 2 and Layer 3 are 83.85 m/s, 116.83 m/s and 139.75 m/s respectively. Thickness of Layer 1 and Layer 2 obtained from profile 3 are 4.04 m and 7.29 m respectively, while the thickness of third layer cannot be identified.

The nature of curve in Figure 13 is A-type. For profile 3 the apparent resistivity values are ρ1 = 68.27 Ωm and h1 = 0.5 m, ρ2 = 89.67 Ωm and h2 = 8.02 m, ρ3 = 100 Ωm and h3 = infinity. Chi-square computed for profile 3 is 0.15 at 1800 iteration (Figure 14).

7.4. Profile 4

As the length of the geophone array taken for profile 4 was half of the profile 1, 2 & 3, hence only two layers were identified in case of profile 4 (Figure 15 & Figure 16). The velocity calculated from profile 4 for layer 1 is 77.69 m/s and layer 2 is 106.04 m/s. Thickness of layer 1 is 3.558 m.

Average thickness of layer one and two are identified to be 3.984 m and 6.533 m respectively. The maximum depth of Low velocity zone which needs to be neglected in order to set datum level for all potential data is about 11m. Average thickness of layer 1, 2 and 3 are calculated 80.02 m/sec, 111.505 m/sec and 135.97 m/sec respectively. By analysing the velocity range of this region it can be calculated that the subsurface is unsaturated sand.

The apparent resistivity value for profile 4 (Figure 17) are ρ1 = 77.27 Ωm and h1 = 0.49 m, ρ2 = 96.67 Ωm and h2 = 10.02 m, ρ3 = 115 Ωm and h1 = infinity. After the analysis of data it has been found that the apparent resistivity value in Dholera for first layer ranges from 60.27 Ωm to 77.27 Ωm, for second layer 90.67 Ωm to 96.67 Ωm and for third layer 110 Ωm to 120 Ωm. Chi-square computed for profile 4 is 0.1 at 1200 iteration (Figure 18).

Figure 11. Travel time graph―profile 3.

Figure 12. Subsurface geological layered cross-section obtained from profile 3.

Figure 13. Apparent resistivity curve for profile 3.

Figure 14. Chi-Square error plot for profile 3.

Figure 15. Travel time graph―profile 4.

Figure 16. Subsurface geological layered cross-section obtained from profile 4.

Figure 17. Apparent Resistivity curve for profile 4.

Figure 18. Chi-Square error plot for profile 4.

8. Conclusion

An extensive seismic and magnetotelluric survey is carried out in Dholera in order to understand the subsurface model. Evolutionary programing method which is also known as Genetic Algorithm method is used for joint inversion of these potential data. Evolutionary programings (genetic algorithm) are class of search algorithms which are used mainly for optimization of problems. In this paper a basic evolutionary programing method is used for cross correlation as fitness function where multiple point crossover method is used. The values of initial population size, probabilities of crossover, mutation and up-gradation of data were used for developing the subsurface model. In this paper joint inversion of seismic refraction is done and the seed point obtained from this is used as primary guess for subsurface resistivity calculation. From seismic refraction method the Low Velocity Layer is obtained up to a depth of 11 m which is neglected in order to bring all the potential data to common datum. Three layers are identified from seismic refraction survey of average thickness of layer one as 3.984 m and layer two as 6.533 m respectively with thickness range of 3.558 m to 7.29 m. The resistivity curves of Dholera region are mainly A-type which means that the resistivity of layers increases with depths. Average resistivity of first layer ranges from 60.27 Ωm to 77.27 Ωm, for second layer 90.67 Ωm to 96.67 Ωm and for third layer 110 Ωm to 120 Ωm.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

[1] | Aghil, T. B., Mohna, K., Srinivas, Y., Rahul, P., Paul, J., Alby, E. R., Nair, N. C., & Chandrasekar, N. (2014). Delineation of Electrical Resistivity Structure Using Magnetotellurics: A Case Study from Dholera Coastal Region, Gujarat, India. Journal of Coastal Sciences, 1, 41-46. |

[2] |
Ammon, C. J., Randall, G. E., & Zandt, G. (1990). On the Nonuniqueness of Receiver Function Inversions. Journal of Geophysical Research, 95, 15, 303-15, 318. https://doi.org/10.1029/JB095iB10p15303 |

[3] | Chipperfield, A., Fleming, P., Pohlheim, H., & Fonseca, C. (1994). Genetic Algorithm Toolbox (pp. 1-105). Sheffield: University of Sheffield. |

[4] |
Gallagher, K., & Sambridge, M. (1994). Genetic Algorithms: A Powerful Tool for Large-Scale Nonlinear Optimization Problems. Computers & Geosciences, 20, 1229-1236. https://doi.org/10.1016/0098-3004(94)90072-8 |

[5] |
Gallagher, K., Sambridge, M., & Drijkoningen, G. (1991). Genetic Algorithms: An Evolution from Monte Carlo Methods for Strongly Non-Linear Geophysical Optimization Problems. Geophysical Research Letters, 18, 2177-2180. https://doi.org/10.1029/91GL02368 |

[6] | Haupt, R. L., & Haupt, S. E. (2004). Practical Genetic Algorithms (2nd ed.). Hoboken: Wiley Interscience Publication. |

[7] | Holland, J. H. (1975). Adaption in Natural and Artificial Systems. Cambridge, MA: MIT Press. |

[8] | Jamshidi, E., & Mostafavi, H. (2013). Soft-Computation Application to Optimize Drilling Bit Selection Utilizing Virtual Intelligence and Genetic Algorithms. In International Petroleum Technology Conference. |

[9] | Poll, H. E. (1994). Automatic forward Modeling of Two-Dimensional Problems in Electromagnetic Induction. PhD Thesis, Victoria: University of Victoria. |

[10] |
Rezaian, A., Alipanah, M., Pour, H. N., & Kazemi, H. (2010). A Genetic Algorithm Approach to Determination of Optimum Diameter of Gas Transmission Pipes. In 34th Annual SPE International Conference and Exhibition (p. 5). Society of Petroleum Engineers. https://doi.org/10.2118/140671-MS |

[11] |
Sambridge, M. S., & Drijkoningen, G. G. (1992). Genetic Algorithms in Seismic Waveform Inversion. Geophysical Journal International, 109, 323-342. https://doi.org/10.1111/j.1365-246X.1992.tb00100.x |

[12] | Shah, M., Sircar, A., Sahajpal, S., Sarkar, P., Sharma, D., Garg, S., Mishra, T., & Shukla, Y. (2017). Geochemical Analysis for Understanding Prospectivity of Low Enthalpy Geothermal Reservoirs of Dholera. In 42nd Workshop on Geothermal Reservoir Engineering (pp. 1-16). Stanford, CA: Stanford University. |

[13] | Sharma, N. (2013). Physico-Chemical Characterization of the Thermal Spring Waters Occurring in Gujarat, India (pp. 1-16). CEGE Report. |

[14] | Sircar, A. (2000). Integrated Seismic Refraction and Geoelectric Sounding Inversion Algorithm for Systematic Evaluation of Aquifers. Acta Geophysica Polonica, 48, 359-375. |

[15] |
Snieder, R. (1998). The Role of Nonlinearity in Inverse Problems. Inverse Problems, 14, 387-404. https://doi.org/10.1088/0266-5611/14/3/003 |

[16] |
Taylor, B. C. R., & Weaver, J. (1976). On the Finite Difference Solution of Two-Dimensional Induction Problems. Geophysical Journal of the Royal Astronomical Society, 47, 375-396. https://doi.org/10.1111/j.1365-246X.1976.tb01280.x |

[17] | Thander, B., & Sircar, A. (2014). Hydrocarbon Resource Estimation: A Stochastic Approach. Journal of Petroleum Exploration and Production Technology, 5, 445-452. |

[18] |
Weglein, A. B., Zhang, H., Ramırez, A. C., Liu, F., & Lira, J. E. M. (2009). Clarifying the Underlying and Fundamental Meaning of the Approximate Linear Inversion of Seismic Data. Geophysics, 74, WCD1-WCD13. https://doi.org/10.1190/1.3256286 |

Journals Menu

Copyright © 2021 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.