Genetic algorithms and smoothing filters in solving the geophysical inversion problem

Yugoslav Journal of Operations Research 12 (2002), Number 2, 215-226 GENETIC ALGORITHMS AND SMOOTHING FILTERS IN SOLVING THE GEOPHYSICAL INVERSION PROBLEM Vesna [E[UM Department of Mathematics, Faculty of Mechanical Engineering University of Belgrade, Belgrade, Yugoslavia Du{an TO[I] Faculty of Mathematics University of Belgrade, Belgrade, Yugoslavia Abstract: The combination of genetic algorithms, smoothing filters and geophysical tomography is used in solving the geophysical inversion problem. This hybrid technique is developed to improve the results obtained by using genetic algorithms only. The application of smoothing filters can improve the performance of GA implementation for solving the geophysical inversion problem. Some test-examples and the obtained comparative results are presented. Keywords: Genetic algorithms, smoothing filters, tomography. 1. INTRODUCTION Tomographic reconstruction techniques are used in different fields: geophysical exploration, medical imaging, astronomy etc. By using geophysical tomography, the characteristics of an underground region between two boreholes can be determined. The tomographic technique provides a means of estimating P-wave velocity in the region between two boreholes. In the time measuring procedure, wave transmitters are placed in one borehole and receivers in the other one in order to collect the first-arrival traveltime data. By using the first-arrival traveltime data, the distribution of wave velocities along a plan between two boreholes is calculated and displayed as a digital picture. Because of that, this problem is known as Image Reconstruction from Projections. The standard tomographic procedures for solving this problem are based on the decomposition of the cross-hole area into a number of cells and an assumption of straight raypaths. The details about the application of V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 216 geophysical tomography and basic notions related to this method can be found in: De Franco and Cavagna [4], Dines and Lytle [5], Ivanson [12], and Peterson et al. [17]. 1.1. Problem formulation The geophysical inversion problem is a problem of determining the characteristics of an underground region by using measurement data. The solution to this problem is often based on applying the tomographic reconstruction techniques. The examinated underground region, called the cross-hole region, is some rectangular area between two boreholes. The problem is to estimate the seismic velocities from first-arrival traveltimes in a cross-hole region using a straight-line ray model. To discretise the geophysical inversion problem, the usual tomographic procedure (Dines and Lytle [5], Peterson et al. [17]) is used. The cross-hole region is divided into a grid of cells (see Fig.1). Some constant, initial velocity is assumed in every cell. Traveltime and velocity structure are interrelated in the relation ( , )v x y ( , ) = ∫ k k R ds t v x y (1) where kt is the traveltime of the kth ray, is the differential raypath length of the k ds th ray, is the two-dimensional velocity function and ( , )v x y kR is the raypath trajectory of the kth ray (Dines and Lytle [5], Peterson et al. [17]). ∆ ijkS iju i kRay R j y x ( , )0 0 Figure 1: Discrete model of the cross-hole region V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 217 It is convenient to denote the reciprocal value 1 with u x , so we deal with slowness instead of velocity. Then the equation (1) becomes / ( , )v x y ( , )y ( , )= ∫ k k R t u x y ds ) (2) When the cross-hole region is divided into a grid of cells ( ×m n , the equation (2) can be approximated as = = = ∆∑ ∑ 1 1 m n k ij ijk i j t u s (3) where is the unknown slowness in the ( , -cell, iju )i j ∆ ijks is the length of a ray segment that intercepts the -cell, m is the number of the vertical cells and n is the number of the horizontal cells. It is understood that ( , )i j ∆ = 0ijks for those i and j values for which the associated cell is not intercepted by the ray kR . A detailed explanation about geophysical tomography can be seen in Dines and Lytle [5], Herman [10] and Ivanson [12]. 1.2. Genetic algorithms Genetic algorithms are a kind of a mathematical simulation for Darwin's theory of evolution. The basic notions and description of genetic algorithms may be found in Goldberg [8]. The starting point is a formation of the initial population either using some particular method or at random. The elements of the initial population, individuals, are points from the searching space for a giving problem. Every individual is uniquely determined by its genetic material. The adaptation of every individual has to be found according to the fact how good a solution that individual is. To that end, the appropriated value of the fitness function is assigned to every individual. Using selection and the values of the fitness function, "the best fitted" individuals are being chosen, while a new population is formed using crossover and mutation. The rules of genetics and evolution imply a greater probability that the new population have a better genetic material. Iterating this procedure, from generation to generation, the genetic material of the population becomes better and the process should converge to the optimal solution of a given problem. If the specified number of generation is reached or some criteria of convergence are fulfilled, the procedure stops. More information about GA can be found in: Beasley et al. [2], Beasley et al. [3], Goldberg [8], Ribeiro-Filho at al. [20], Srinivas and Patnaik [23]. 1.3. Smoothing filters The application of the appropriated smoothing filters can improve the performance of the used iterative reconstruction algorithm (Herman [10], Herman and Lent[11]) and form smooth areas in regions of the image. That implies better-defined boundaries in a reconstructed image (Balanis and Bentley [1]). The possibility of V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 218 achieving these improvements makes the combination of filters and the iterative reconstruction algorithm attractive (Herman [10], Herman and Lent[11]). Usually, the smoothing filter is implemented between the iterative steps of the reconstruction algorithm (Balanis and Bentley [1]). All cells in the image go throw filter. For every cell the window of the filter is formed consisting of that particular "center" cell and its neighboring cells. Some often-used shapes of window can be seen at Fig. 2. The value of the center cell is replaced with a new value and it is calculated using the values of cells within the window of that center cell. The new values become actual after all cells in the image have been passed through the filter. (More details are exposed in section 3.1, 3.2 and 3.3.) Figure 2: Square-shaped window and cross-shaped window 2. SOLVING THE GEOPHYSICAL INVERSION PROBLEM USING GENETIC ALGORITHMS This section contains a short review of solving the geophysical inversion problem with the use of genetic algorithms. Detailed explanation is presented in [e{um et al [24]. The encoding mechanism is a very important issue for the structure of a genetic algorithm. The genes representing coded numerical value of velocity in a cell of cross-hole region have the same fixed length. In the described implementation that length is 5 bits (this information has to be given in advance). The value of velocity is obtained by min= + ⋅v v dv tmp , (4) where is the velocity partition, dv (max min ) /= − 2nbitdv v v (5) min v is the minimum velocity, is the maximum velocity, tm is the value of the binary string that represents velocity over a cell and nbi is given in advance. max v p t V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 219 The objective function is ( ) δ = = − = ∑ ∑ 2 1 1 2 m n T C ij ij i j t t k (6) where is the measured traveltime, is the computed traveltime and k is the number of transmitters/receivers. The computed traveltimes are calculated from (3) and measured traveltimes are obtained after the measurement in a borehole. The objective function (6) arises as a measure of differences between the computed and measured traveltimes. The goal is to minimize this function. So, we have here a typical optimization problem and a genetic algorithm could be applied. T ijt C ijt The following genetic operators are used: rank-based selection, one-point crossover and mutation. The crossover rate is .= 0 85crossp and the mutation rate is in segment [0.001, 0.03], depending on the number of cells in the cross-hole region. The population size is 150 individuals and the maximum number of generations is 5000. The proposed method is tested on synthetic examples, where the dimensions of the cross-hole regions are 100m×100m and 80m×60m. This gave promising results and provided quite good performances. (More details about the application of genetic algorithms can be found in: Beasley at al. [2], Beasley at al. [3] and Goldberg [8].) 3. SOLVING THE GEOPHYSICAL INVERSION PROBLEM USING SMOOTHING FILTERS AND GENETIC ALGORITHMS This section describes a new hybrid technique - the combination of smoothing filters and genetic algorithms. The purpose of developing this technique was to improve the results obtained using genetic algorithm only. Experiences in researching filters and their application in geophysical tomography (Balanis and Bentley[1]) showed that the following filters gave relatively good results: − minimum variance partitioning-averaging, MVP-AVG − minimum variance partitioning-median filtering, MVP-MED − selective smoothing 3.1. MVP-AVG Instead of using all neighboring cells in the window (see Fig 2.), it is suggested to select a group of neighbors in the window (Balanis and Bentley [1]), in order to obtain better results. The selection for a group of neighbors can be done with MVP technique (Hall[9]). V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 220 For a given set of numbers, MVP technique finds the partition of set with minimum variance. The set of number is divided into several groups (usually two or three). To reduce the number of possible partitions, the numbers from a given set were previously sorted increasingly. For n numbers, there are −1n partitions and every partition has a fixed (in advance) number of groups. Each partition is processed in the following way: − for every partition's group, finding its mean , ,...,= 1im i number_of_ groups − calculating , ,..., - = ∈   = −   ∑ ∑ 2 1 1 j i numberofgroups l j i i u group S u m l = 1n (7) Finally, the minimum in minS ( ) −= 11nl lS is found and the partition that is joined with the value is the minimum variance partitioning for that set of numbers. This way provides an objective mechanism for determining the neighborhood for a given center cell. Then, a new value of the center cell can be determined in various ways. minS Example. Suppose that the following values are in the cells: 2,5,3,9,7 and the number of groups is 2. After the sorting, we have: 2,3,5,7,9. The partitions are: :{ }{ , , , }, . , . :{ , }{ , , }, . , . :{ , , }{ , }, . , . :{ , , , }{ }, . , . = = = = = = = = 1 1 2 1 3 1 4 1 2 3 5 7 9 2 00 6 00 2 3 5 7 9 2 50 7 00 2 3 5 7 9 3 33 8 00 2 3 5 7 9 4 25 9 00 S m m S m m S m m S m m 2 2 2 2 Applying the formula (7), we get: , . , . , .= = = =1 2 3 420 8 5 7 01 14 75S S S S :{ , , }{ , }3 2 3 5 7 9S . The minimum is S3 and the selected partition is . One approach is to replace the value of the center cell with the average value of group that contains the center cell. This technique is called MVP-AVG. More information about MVP technique, averaging technique and MVP-AVG can be found in: Balanis and Bentley [1], Rosenfeld [21], Radcliff et al. [19], Schowengerdt [22] and Hall [9]. 3.2. MVP-MED After determining the neighborhood for a given center cell using MVP technique, as described in 3.1, this technique replaces the value of the center cell with the median value of the group that contains the center cell (Balanis and Bentley [1], Ekstrom [7], Rosenfeld [21]). Example: Let us consider the square-shaped window from Fig 2. The center cell will take on the fifth highest value of the cells within the window. V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 221 3.3. Selective smoothing New value of the center cell is formed by that cell and those neighboring cells with values that satisfy the condition, implied by the threshold level. If the difference between the value u of the center cell and the value of the neighboring cell is less or equal to the threshold level, that neighboring cell will be included in calculating a new value of the center cell (Radcliff and Balanis [18], Herman [10], Balanis and Bentley [1]). 1 iu Example: Let us consider the square-shaped window from Fig. 2 with dimension 3×3 which is the neighborhood of the center cell . The new value of that center cell can be obtained from 1u = = = = + + = + + ∑ ∑ ∑ ∑ 5 9 1 1 2 3 2 6 1 5 9 1 2 3 2 6 i i i i I i i i i i i w u w f u w f u u w w f w f (8) where , ,..., ,  − ≤=  11 for 2 9 0 else i i u u t f =i 3 (9) , ,1 2w w w are smoothing weights (choosen in advance), t is threshold level. The successful application of this filter and the achievement satisfactory results are conditioned by the choice of the optimal values for w1,w2,w3 and t (Radcliff and Balanis [18]). 3.4. Genetic algorithms and filters The filters described in the previous section are implemented in the GA application for solving the geophysical inversion problem to improve the reconstruction of profile. Applying the appropriated filter between every iteration in the reconstruction algorithm (i.e. GA generation) can improve the convergency of that algorithm, if the system of equations (3) is consistent. Otherwise, this approach does not give good results (Balanis and Bentley [1]). In that case, the improvement can be obtained by using filter after a certain number of iterations, when the picture reaches some development (Balanis and Bentley [1]). A lot of numerical experiments are made with a different number of generations and different filters. In the following variants of the hybrid technique, the optimal number of generations (obtained by experiments) and optimal filters are used. a) 4000 GA generations were applied to the test-profile; then such partialy reconstructed profile is subjected to 100 passings through filter after every next GA generation; b) after every 100 GA generations, a partialy reconstructed profile is subjected to filtering for 10 times; V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 222 c) 3000 GA generations were applied to the test-profile; then such partialy reconstructed profile is subjected to 1 passing through filter after every next GA generation; d) after applying 5000 GA generations to the test-profile, the reconstructed profile (obtained by running of the pure GA) is subjected to filtering. Two-group MVP is used. Some conclusions about the parameters (presented below) in selective smoothing are accepted, as suggested in Radcliff and Balanis [18]. The investigated region may have one or more anomalies where values differ from background's values (see Fig 3.). For profiles that have high values in the anomaly's cells and low values in the background's cells, 1/6 of the difference between the highest value and the lowest value in the profile is assigned to the variable t. Reversely, for profiles that have low values in the anomaly's cells and high values in the background's cells, 2/5 of the difference between the highest value and the lowest value in the profile is assigned to the variable t. The values and are chosen to satisfy the relation w w for every profile. ,1 2w w 3w = =1 2 = 32w Figure 3: A hypothetical 'L' shaped anomaly in the background between two boreholes The implementation of MVP-AVG filtering technique, described in section 3.1. in the hybrid algorithm, can be described with Pascal pseudo-code: (* MVP-AVG variant(a) *) begin if (current_generation > 4000) then begin for filt := 1 to 100 do begin for i :=1 to m do for j :=1 to n do begin sort in asceding order cells in window that contains cell (i,j); apply MVP-AVG filter to the window of cell (i,j); new value of cell (i,j) put in mat(i,j); end; V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 223 replace old values of cells with the new ones from matrix mat; end end end. The other variants ((b), (c) and (d)) of MVP-AVG filter as well as MVP-MED filters and selective smoothings are implemented analoguously. 4. RESULTS 4.1. Test-examples Syntetic test-examples, where the examinated profile consists of part by part constant regions, are generated (see Fig. 4). The number of cells is 30 (6×5). In every borehole there are 10 transmitters/receivers. The distance between boreholes is 5.0m and the borehole's depth is 6.0m. The threshold level .= 0 1t , because in these test- examples high values are in the anomaly's cells and low values are in the background's cells. Figure4: Profiles of test-instances ,1 2P P and 3P where the values in anomaly's cells are 0.7 and the values in background's cells are 0.1 The testing has been done using PC compatible computer AMD 80486 at 120MHz with 8 MB of memory. Since genetic operators are nondeterministic, every test example was running 10 times and we have computed the average value. 4.2. Numerical results Table 1 contains results obtained by the following methods: • ART (Algebraic Reconstruction Technique)(see Dines and Lytle [4]) and different versions of ART i.e. M-ART (Modified-ART), DT-ART (Distance Transposed-ART), DTM-ART (Distance Transposed Modified-ART) and MDT-ART (Modified Distance Transposed-ART). These results are taken from Balanis and Bentley [1]. A detailed explanation about these techniques can be found in Balanis and Bentley [1]. V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 224 • GA ([e{um et al [24]). The measure of error δ2 given by relation: ( ) δ = = − = × ∑ ∑ 2 1 1 2 m n T C ij ij i j u u m n (10) is used, where is the true slowness, is the computed slowness and is the total number of cells in the profile. T iju C iju ×m n Table 1. 1P 2P 3P ( )δ2 GA 0.02901 0.03221 0.17492 ( )δ2 ART 0.13300 0.20100 0.22300 ( )δ −2 DT ART 0.08300 0.15600 0.81200 ( )δ −2 M ART 0.08100 0.16500 0.17600 ( )δ −2 MDT ART 0.08300 0.14800 0.17300 ( )δ −2 DTM ART 0.09200 0.14600 0.17900 The variants (a), (b), (c) and (d) of the implemented filters and the GA com- binations are described in section 3.4, while the obtained results are given in Table 2. Table 2. 1P 2P 3P ( ( ))δ −2 MVP AVG a 0.03009 0.01686 0.16008 ( ( ))δ −2 MVP AVG b 0.03053 0.22923 0.31102 ( ( ))δ −2 MVP AVG c 0.03356 0.23989 0.35823 ( ( ))δ −2 MVP MED a 0.05073 0.07503 0.41323 ( ( ))δ −2 MVP MED b 0.02984 0.05272 0.38797 ( ( ))δ −2 MVP MED c 0.03607 0.16241 0.24384 ( ( ))δ2 selective smooth. a 0.02640 0.04177 0.24936 ( ( ))δ2 selective smooth. b 0.03008 0.12567 0.37428 ( ( ))δ2 selective smooth. c 0.02778 0.09194 0.26610 ( / ( ))δ −2 GA MVP AVG d 0.00657 0.00933 0.14037 ( / ( ))δ −2 GA MVP MED d 0.00730 0.00000 0.14715 ( / ( ))δ2 selective smooth. GA d 0.02746 0.02992 0.16970 V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 225 Table 3 contains results of the application of the most successful ART modification, MDT-ART and filters (Balanis and Bentley [1]). Table 3. 1P 2P 3P ( / )δ − −2 MDT ART MVP AVG 0.050 0.187 0.195 ( / )δ − −2 MDT ART MVP MED 0.062 0.182 0.155 ( / )δ −2 selective smooth.MDT ART 0.076 0.142 0.170 From the comparison of the results obtained using GA with the results obtained using ART (see Table 1) we can conclude that GA give better results (in case of 1P and 2P ) and almost the same result (in case of 3P ). The comparison of the results obtained using the combination of GA and filters with the results obtained using the combination of ART and filters (from Table 2. and Table 3) shows that the combination of GA and filters gives better results for all three test-examples. (variant( ))d The values presented in Tables 1-3 confirm that the combination of GA and filters improves results obtained using GA only. 5. CONCLUSION Generally, GAs are very suitable for hybridizing with other methods. In this paper the new hybrid method (GA and filters) for solving the geophysical inversion problem is presented. The intention was to develop the technique that would improve the results obtained by using genetic algorithms only. The applying of smoothing filters can improve the performances of a main iterative reconstruction algorithm as well as the quality of reconstructed picture. The applying of new programming techniques, in this case, is not requested. The chosen GA structure provided relatively good results. We used the following filters: MVP-AVG, MVP-MED and selective smoothing. All of them showed good features in combination with GA. REFERENCES [1] Balanis, C.A., and Bentley, D., "Algorithm and filter selection in geophysical tomography", IEEE Transaction on Geoscience and Remote Sensing, 24(6) (1986) 983-996. [2] Beasley, D., Bull, D.R., and Martin, R.R., "An overview of genetic algorithms: Part 1, fundamentals", University Computing, 15 (2) (1993) 58-69. [3] Beasley, D., Bull, D.R., and Martin, R.R., "An overview of genetic algorithms: Part 1, research topics", University Computing, 15 (4) (1993) 170-181. [4] De Franco, R., and Cavagna, B., "Control of grouted tunnels through seismic topography", Gallerie e grandi pere otterrance, 47 (1995). [5] Dines, K.A., and Lytle, R.J., "Computerized geophysical tomography", Proc. of IEEE, 67 (7) V. [e{um, D. To{i} / Genetic Algorithms and Smoothing Filters 226 (1979) 1065-1073. [6] Dyke, V.A., and Young, P.R., "Physical characterization of rock masses using borehole methods", Geophysics, 50 (1985) 2530-2541. [7] Ekstrom, M.P., Digital Image Processing Technique, Academic, Orlando, 1984. [8] Goldberg, D.E., Genetic Algorithms in Search, Optimization and Machine Learning, Addison Wesley Pub. Comp., 1989. [9] Hall, E.L., Computer Image Processing and Recognition, Academic, New York, 1979. [10] Herman, G.T., Image Reconstruction From Projection, Academic, New York, 1980. [11] Herman, G.T., and Lent, A., "Iterative reconstruction algorithms", Comput. Biol. Med., 6 (1976) 273-294. [12] Ivanson, S., "Seismic borehole tomography - Theory and computational methods", Proc. of IEEE, 74 (2) (1986) 328-338. [13] Lytle, R.J., Lager, D.L., Laine, E.F., and Salisbury, J.D., "Monitoring fluid flow by using high- frequency electromagnetic probing", Lawrence Livermore Lab. Rep. UCRL-51979 (ERA 1:17178, available from National Technical Inform. Serv., Springfield, VA), 1979. [14] Lytle, R.J., Laine, E.F., Lager, D.L., and Okada, J.T., "Determination of the in situ high frequency electrical properties of permafrost rock", Radio Sci., 11 (4) (1976) 285-294. [15] Lytle, R.J., Dines, K.A., Laine, E.F., and Lager, D.L., "Electromagnetic cross-borehole survey of a site proposed for a future urban transit station", Lawrence Livermore Lab. Rep. UCRL- 52484, (available from National Technical Inform. Serv., Springfield, VA), 1978. [16] Mathisen, E.M., Vasiliou, A.A, Cunningam, P., Shaw, J., Justice, H.J., and Guinzy, J.N., "Time-lapse crosswell seismic tomogram interpretation: implications for heavy oil reservoir characterization, thermal recovery process monitoring and tomographic imagine technology", Geophysics, 60 (1995) 631-650. [17] Peterson, J.E., Paulsson, B.N.P., and McEvilly, T.V., "Application of algebraic reconstruction techniques to crosshole seismic data", Geophysics, 50 (10) (1985) 1566-1580. [18] Radcliff, R.D., and Balanis, C.A., "Electromagnetic geophysical tomography incorporating refraction and reflection", IEEE Trans. Antennas Propagat., 29 (2) (1981). [19] Radcliff, R.D., Balanis, C.A., and Hill, H.W., "A stable geotomography technique for refractive media", IEEE Trans. Geosci. Remote Sensing, 22 (6) (1984) 698-703. [20] Ribeiro-Filho, J.L., Treleaven, P.C., and Alippi, C., "Genetic algorithm programming environments", IEEE Computer, (1994) 28-43. [21] Rosenfeld, A., and Kak, A.C., Digital Picture Processing, Academic, Orlando, 1982. [22] Schowengerdt, R.A., Techniques for Image Processing and Classification in Remote Sensing, Academic, New York, 1984. [23] Srinivas, M., and Patnaik, L.M., "Genetic algorithms: A survey", IEEE Computer, (1994) 17- 26. [24] [e{um, V., Kratica, J., and To{i}, D., "Solving the geophysical inversion problem using genetic algorithms", YUJOR, 10 (2) (2000) 283-292. [25] Wong, J., Hurley, P., and West, G.F., "Crosshole seismology and seismic imaging in crystalline rocks", Geophys. Res. Lett., 10 (1983) 686-689. [26] Yamamoto, T., Nye, T., and Kuru, M., "Porosity, permeability, shear strength: crosswell tomography below an iron foundry", Geophysics, 59 (1994) 1530-1541.