- Email: [email protected]
A wavelet adaptive-homotopy method for inverse problem in the fluid-saturated porous media
Applied Mathematics and Computation 208 (2009) 189–196
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A wavelet adaptive-homotopy method for inverse problem in the fluid-saturated porous media q Ying He *, Bo Han Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e
i n f o
Keywords: Wavelet Adaptive-homotopy method Inversion Fluid-saturated porous media
a b s t r a c t In this paper, we consider a parameter identification problem for the elastic wave equations in the fluid-saturated porous media. A wavelet adaptive-homotopy method is proposed for the recovery of porosity. This algorithm combines the wavelet multiscale inversion idea with the adaptive-homotopy method. And the adaptive-homotopy method provides a simple way to adapt computational refinement to the choice of the homotopy parameter. Numerical results illustrate that the wavelet adaptive-homotopy method is globally convergent and highly effective for solving the inverse problem. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction The elastic wave propagation in porous media has received considerable attention in resent years because of increasing applications in science and engineering. Applications of porous media include Mechanics, Engineering, Geosciences, Material science. In many engineering problems, such as those found in petroleum and geotechnical engineering, rocks have to be modeled as porous materials saturated with fluid. In 1956, Biot [1,2] firstly considered the propagation of elastic wave in the fluid-saturated porous media. He summarized a general and systematic theory to describe a porous solid containing a viscous fluid. Since then, various numerical techniques have appeared in the literatures of the inverse problem in porous media. The inverse problem in the fluid-saturated porous media can be viewed as a parametric data-fitting problem. It is possible to formalize such a problem in the framework of optimization where a functional defined in terms of discrepancy between observed and computed data is minimized over a model space. In general, such problem is very difficult to solve, since it is nonlinear and ill-posed. The nonlinear property causes the presence of numerous local minimum, while the illposedness often affects the efficiency of the numerical scheme significantly. For the nonlinear and ill-posed problem, traditional linearized inversion algorithms, such as conjugate gradient method, time convolution iteration method [3] and perturbation method [4], are locally convergent. The widely popular recent approaches such as simulated annealing algorithm, genetic algorithm, trust region method, and neural networks [5], are globally convergent. However, along with the searching space decreasing, the efficiency of these are worse than before. Therefore, the shortages of the above inversion methods motive us to design a highly efficient, numerically stable and globally convergent algorithm. Wavelet multiscale method is a newly developmental inversion strategy which accelerate convergent rate, enhance stability, and overcome disturbance of local minimum. Because of these features, the applications of the wavelet multiscale inversion approach are giving rise to much attention. Liu [6] proposed wavelet multiresolution method for distributed parameter estimation of the elliptic problem. Ebrahimi and Sahimi [7] has described scale up of geological models of field-scale porous media using a method based on the wavelet transformation. Zhang et al. [8] firstly recovered the porosity q
Project supported by the National Natural Science Foundation of China (No. 40774056). * Corresponding author. E-mail address: [email protected] (Y. He).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.11.033
190
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
in the fluid-saturated porous media by using wavelet method associated with the Gauss–Newton iterative method, and the results showed that the Gauss–Newton method performed much better when employed with wavelet multiscale method. Homotopy method is a powerful tool for solving inverse problems due to its globally convergent property. It is a very interesting research to apply homotopy method to solve the inverse problem. Various formulations for the homotopy method have been used in the literature [9–13]. Although the homotopy method is applied in many fields, the literatures about the availability of the homotopy method in the fluid-saturated porous media have rarely been found except [14], where Han et al. recovered the parameter of 1D wave equation. All these works testified that the homotopy method is effective for the inverse problem. As mentioned above, in this paper, we will design an inversion strategy called wavelet adaptive-homotopy method, for parameter identification problem in the fluid-saturated porous media. This new algorithm combines the property of good stability and fast convergent rate of wavelet multiscale method with the advantage of global minimum of homotopy method. In order to search the global minimum efficiently, we provide a simple way to adapt computational refinements to the choice of the homotopy parameter. The structure of the paper is organized as follows: we start with presenting the elastic wave equations and correlative parameters in the fluid-saturated porous media based on Biot theory in Section 2. In Section 3, an adaptive-homotopy method is proposed and the adaptive property of this algorithm provides a simple way to adapt a computational refinements to the choice of the homotopy parameter. In Section 4, a wavelet adaptive-homotopy is proposed for solving the inverse problem in the fluid-saturated porous media. Then we describe how this algorithm can be carried out by combining the wavelet multiscale method with the adaptive-homotopy method. Finally, in order to illustrate the performance of this algorithm, we carry out some numerical simulations for the inversion of porosity in fluid-saturated porous media. The results obtained show that the wavelet adaptive-homotopy method is effective and feasible to the inverse problem. 2. Governing equations The propagation of elastic wave in the fluid-saturated porous media is different from that in the single-phase media. Based on continuum mechanics and macroscopic constitutive relationship, Biot developed a theory of wave motion in an elastic solid saturated with a viscous compressible fluid. We ignore viscous effects and consider the 2D elastic wave propagation equations in fluid-saturated porous media as follows:
2
o ox
o ox
oux o ou ou o oux ouz o ow ow o2 u o2 w l xþl z þ aM x þ aM z ¼ q 2x þ qf 2 x f1 ; k þ þk þ oz ox ox ox oz ox ox oz ox oz ot ot
ð2:1Þ
oux ouz o ou o oux ouz o ow ow o2 u o2 w l z þ aM x þ aM z ¼ q 2z þ qf 2 z f2 ; k þl þ2 þk þ oz oz oz oz ox oz ox oz ox oz ot ot
ð2:2Þ
l
l
o ou ou ow ow o2 u o2 w aM x þ aM z þ M x þ M z ¼ qf 2x þ m 2 x f1 ; ox ox oz ox oz ot ot
ð2:3Þ
o ou ou owx owz o2 uz o2 wz aM x þ aM z þ M þM ¼ qf 2 þ m 2 f2 : oz ox oz partialx oz ot ot
ð2:4Þ
where u ¼ ðux ; uz ÞT denotes the solid-frame displacement, w ¼ ðwx ; wz ÞT is the fluid displacement relative to solid-frame, k is coefficient, the quantity j is the Darcy permeability coefficient, g is the viscosity of the pore fluid, qf is the density the Lame of the pore fluid, qs is the density of the solid grain, q is the bulk density defined by q ¼ bqf þ ð1 bÞqs , and b is the porosity. The relationship of the tortuosity a, added mass density m, and coupling constant M can be given by
a¼1
Ks ; Kr
M¼
K 2r ; Dr K s
Kr Dr ¼ K r 1 þ b 1 ; Kf
m¼
qf b
;
where K r is the bulk modulus of the grain, K f is the bulk modulus of the pore fluid, K s is the bulk modulus of the skeletal frame. The boundary conditions for the problem considered here can be given by
oux ðx; z; tÞ oux ðx; z; tÞ ¼ ¼ 0; ox ox x¼L x¼0 oux ðx; z; tÞ ¼ 0; oz z¼H ouz ðx; z; tÞ ouz ðx; z; tÞ jx¼0 ¼ ¼ 0; ox ox x¼L
ð2:5Þ ð2:6Þ ð2:7Þ
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
ouz ðx; z; tÞ ¼ 0; oz z¼H owx ðx; z; tÞ owx ðx; z; tÞ ¼ ¼ 0; ox ox x¼0 x¼L owx ðx; z; tÞ ¼ 0; oz z¼H owz ðx; z; tÞ owz ðx; z; tÞ ¼ ¼ 0; ox ox x¼0 x¼L owz ðx; z; tÞ ¼ 0; oz z¼H
191
ð2:8Þ ð2:9Þ ð2:10Þ ð2:11Þ ð2:12Þ
where ðx; zÞ 2 ½0; L ½0; H; t 2 ½0; T. The initial conditions are given as follows:
oux ðx; z; 0Þ ¼ 0; ot ouz ðx; z; 0Þ uz ðx; z; 0Þ ¼ 0; ¼ 0; ot owx ðx; z; 0Þ ¼ 0; wx ðx; z; 0Þ ¼ 0; ot owz ðx; z; 0Þ ¼ 0; wz ðx; z; 0Þ ¼ 0; ot ux ðx; z; 0Þ ¼ 0;
ð2:13Þ ð2:14Þ ð2:15Þ ð2:16Þ
and the additional condition is
ux ðx; 0; tÞ ¼ ux ðx; tÞ;
ð2:17Þ
where ux is the observed data. If given a value of the porosity b, Eqs. (2.1)–(2.6), (2.8), (2.7), (2.9)–(2.16) enable to form the direct problem for wave fields simulations, i.e. we can solve the wave field functions ux ; uz ; wx ; wz using a numerical method. Consequently, the additional condition ux ðx; 0; tÞ can be obtained. Note that since the porosity b is unknown in practice, governing equations (2.1)–(2.4) with the boundary conditions (2.5), (2.6), (2.8), (2.7), (2.9)–(2.12), the initial conditions (2.13)–(2.16) and the additional condition (2.17) completely define the inverse problem for the elastic wave equations in the fluid-saturated porous media. 3. An adaptive-homotopy method In general, the inverse problem in the fluid-saturated porous media can be viewed as a parametric data-fitting problem. It is possible to formalize such a problem in the framework of optimization where a function defined in terms of discrepancy between observed and computed data is minimized over a model space. Now, we discretize Eqs. (2.1)–(2.6), (2.8), (2.7), (2.9)–(2.17) by the second order finite difference method as follows:
ux ði; j; k þ 1Þ ¼ wx ði; j; k þ 1Þ ¼
mU 1 qf U 2 mqði; jÞ q2f
;
uz ði; j; k þ 1Þ ¼
mU 3 qf U 4 mqði; jÞ q2f
;
qði; jÞU 2 qf U 1 qði; jÞU 4 qf U 3 ; wz ði; j; k þ 1Þ ¼ ; 2 mqði; jÞ qf mqði; jÞ q2f
ux ð0; j; kÞ ¼ ux ð1; j; kÞ;
ux ðp 1; j; kÞ ¼ ux ðp; j; kÞ;
ux ði; n 1; kÞ ¼ ux ði; n; kÞ;
uz ð0; j; kÞ ¼ uz ð1; j; kÞ;
uz ðp 1; j; kÞ ¼ uz ðp; j; kÞ;
uz ði; n 1; kÞ ¼ uz ði; n; kÞ;
wx ð0; j; kÞ ¼ wx ð1; j; kÞ;
wx ðp 1; j; kÞ ¼ wx ðp; j; kÞ;
wx ði; n 1; kÞ ¼ wx ði; n; kÞ;
wz ð0; j; kÞ ¼ wz ð1; j; kÞ;
wz ðp 1; j; kÞ ¼ wz ðp; j; kÞ;
wz ði; n 1; kÞ ¼ wz ði; n; kÞ;
ux ði; j; 0Þ ¼ ux ði; j; 1Þ ¼ 0; wx ði; j; 0Þ ¼ wx ði; j; 1Þ ¼ 0; ux ði; 0; kÞ ¼
uz ði; j; 0Þ ¼ uz ði; j; 1Þ ¼ 0; wz ði; j; 0Þ ¼ wz ði; j; 1Þ ¼ 0;
ux ði; kÞ;
where b is the model to be recovered, ux is the observed data, ux ði; j; kÞ ¼ ux ði hx ; j hz ; k sÞ, hx ; hz are the step sizes of the rectangle grid in the x- and z-direction, respectively, s is the time step size, p ¼ hLx ; n ¼ hHz ; l ¼ Ts. The above difference equations define an operator equation
F : b ! U;
i:e: FðbÞ ¼ U;
192
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
where b and U denote vectors, which are composed of bði; jÞ and ux ði; kÞ in the following sequence:
bT ¼ ðbð1; 1Þ; . . . ; bð1; nÞ; bð2; 1Þ; . . . ; bð2; nÞ; . . . ; bðp; nÞÞ; U T ¼ ðux ð1; 1Þ; . . . ; ux ðp; 1Þ; ux ð2; 1Þ; . . . ; ux ðp; 2Þ; . . . ; ux ðp; lÞÞ; where p; n; l denote the number of mesh nodes for the x; z; t components, respectively. In order to identify the porosity b from the observed data U containing noise, we define a nonlinear operator equation
TðbÞ , FðbÞ U ¼ 0:
ð3:1Þ
Eq. (3.1) is usually highly nonlinear and ill-posed, therefore we consider a globally convergent algorithm, namely an adaptive-homotopy method, to solve the nonlinear operator equation (3.1). Let us construct a homotopy regularization method
Hðb; aÞ ¼ ð1 aÞkFðbÞ Uk2 þ akb b0 k2 ;
ð3:2Þ
0
where b is the initial guess, the parameter a 2 ½0; 1 is not only the homotopy parameter bur also the regularization parameter. When a ¼ 1, the solution of Hðb; 1Þ ¼ 0 is b0 , which is known. When a ¼ 0, the solution of Hðb; 0Þ ¼ 0 is the desired solution of the nonlinear operator equation (3.1). Above all, we divide the interval [0, 1] into 1 ¼ a0 > a1 > > aN ¼ 0. For some ak , we use some iterative method to solve Hðb; ak Þ ¼ 0 sequently. As the solution b0 of Hðb; 1Þ ¼ 0 is known, it can be used as the initial approximation of the next equation. Assume that the approximation solution bk of the Hðb; ak Þ ¼ 0 has already been found. For computing the solution bkþ1 of Hðb; akþ1 Þ ¼ 0, we use the damped Gauss–Newton method. In practice, we need not compute the accurate solution of each equation. So we only iterate one step when using the damped Gauss–Newton method to solve the equation Hðb; ak Þ ¼ 0. The iterative formulation is given as follows:
bkþ1 ¼ bk ½ð1 ak ÞF 0 ðbk Þ F 0 ðbk Þ þ ak I1 ½ð1 ak ÞF 0 ðbk Þ ðFðbk Þ UÞ þ ak bk ;
k ¼ 0; 1; . . . ; N 1:
ð3:3Þ
From formulation (3.3), the approximation solution bN of the Hðb; 1Þ ¼ 0 can be solved, which can be taken as the initial guess to the problem (3.1). Then iterate repeatedly by using the damped Gauss–Newton algorithm
bkþ1 ¼ bk ½F 0 ðbk Þ F 0 ðbk Þ1 ½F 0 ðbk Þ ðFðbk Þ UÞ;
k ¼ N; N þ 1; . . .
ð3:4Þ
The formulations (3.3) and (3.4) constitute a globally convergent method for the inverse problem. As we all know, the homotopy parameters ak ðk ¼ 0; 1; . . .Þ plays an important role. But the homotopy parameters always can be given by ak ¼ Nk in the existing literatures. Such choice may loss some convergent approximation solutions of the Eq. (3.1). Therefore, in this paper, we propose a modified homotopy method named by the adaptive-homotopy method and apply the new approach to find the global minimum of the inverse problem. This algorithm can adaptively choose the homotopy parameters ak . Algorithm 1. Choose an initial model b0 , small threshold e and initial homotopy step size h. Set a1 ¼ a0 h, a2 ¼ a1 h. Calculate the solution b1 of Hðb; a1 Þ ¼ 0 (using damped Gauss–Newton method). For k ¼ 2; . . ., do:
Calculate the solution bk of Hðb; ak Þ ¼ 0 (using damped Gauss–Newton method). kbk bk1 k . Calculate j ¼ Vertb k1 bk2 k Check if j < e, then set h ¼ 2h, else if j P e, then set h ¼ 2h. Check the convergence, if not converged: set k ¼ k þ 1 ak ¼ ak1 h: and go to step 1.
4. Wavelet adaptive-homotopy method The wavelet multiscale method is a technique that improves the performance of inversion method by decomposing the problem into subproblem on different scales. At coarser scale there are fewer local minima and those that remain are further apart from each other, and then inversion method can get the solution closer to the neighborhood of the global minimum. As the scale increases, the complexity of the original problem is reduced and the size of the corresponding domain of convergence is decreasing. In this section, we try to use wavelet multiscale method to decompose the original operator equation (3.1). Once the problem has been decomposed into different subproblems, it is natural to solve first the subproblem on Table 1 The decomposition of the original problem by wavelet. Scale space
Source
Wave field
Solved problem
Solution of the problem
Initial guess
V0 V1 V2
f f1 f2
ux u1 x u2 x
FðbÞ ¼ U F 1 ðbÞ ¼ U 1 F 2 ðbÞ ¼ U 2
b b1 b2
b1 b2 b0
193
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0
20
40
60
80
100
120
140
160
180
200
Time axis (0.001s) Fig. 1. The source function.
Table 2 The results of inversion by wavelet adaptive-homotopy method. Number
Initial value
Real value
Iterative value
Iterative times
Relative error (%)
1 2 3 4
0.045 0.237 0.384 0.659
0.266000 0.266000 0.266000 0.266000
0.26472 0.26579 0.26678 0.26713
27 14 17 30
0.483 0.079 0.294 0.425
Table 3 The results of inversion by homotopy method. Number
Initial value
Real value
Iterative value
Iterative times
Relative error (%)
1 2 3 4
0.045 0.237 0.384 0.659
0.266000 0.266000 0.266000 0.266000
0.26464 0.26554 0.26714 0.26758
34 21 26 39
0.512 0.174 0.427 0.593
0
0
5
10
Surface axis (×10m) 15
20
25
0.4 0.38
5
Depth axis (x10m)
30
0.36 0.34
10
0.32 15
0.3 0.28
20
0.26 25
0.24 0.22
30
0.2 Fig. 2. The true value of porosity model.
194
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
the coarsest scale by using adaptive-homotopy method, and to use the resulting solution as initial guess for the subproblem on the finer scale. Without loss of generality, we describe the wavelet adaptive-homotopy method by decomposing the original problem into two scales. Assume that the original problem is on the scale space V 0 and decomposed on the coarser scale space V 1 as well as the coarsest scale space V 2 . We utilize the compactly supported orthogonal wavelet ‘DB4’ to decompose the wave field functions and source function. The algorithm can be described as follows:
Step 3. Step 4.
Surface axis (x10m) 0
0
5
10
15
20
25
30 0.55
5
Depth axis (×10m)
Step 2.
Given an initial model b0 , the true porosity model b , the source function f. Calculate the forward problem to obtain the model data U. Decompose f and U at each scale, such as V 1 and V 2 by the compactly supported orthogonal wavelet ‘DB4’ (using Mallat algorithm [15]), i.e. the original problem (3.1) can be decomposed as F 1 ðbÞ ¼ U 1 and F 2 ðbÞ ¼ U 2 , which is shown in the following Table 1. At the coarsest scale space V 2 , calculate the solution b2 of F 2 ðbÞ ¼ U 2 with the initial guess b0 (using the adaptivehomotopy method). At the finer scale space V 1 , we regard the solution b2 as the initial guess, and then calculate the solution b1 of F 1 ðbÞ ¼ U 1 (using the adaptive-homotopy method).
0.5 0.45
10
0.4 0.35
15
0.3 0.25
20
0.2 0.15
25
0.1 0.05
30 Fig. 3. The inversion result at scale 3.
Surface axis (×10m) 00
5
10
15
20
25
30 0.45
5 0.4
Depth axis (×10m)
Step 1.
10 0.35 15
0.3
20
0.25
25
0.2 0.15
30 Fig. 4. The inversion result at scale 2.
195
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
Step 5.
At the initial space V 0 , we regard the solution b1 as the initial guess, and then calculate the solution b of FðbÞ ¼ U (using the adaptive-homotopy method).The obtained solution b is the solution of original problem (3.1).
5. Numerical simulations In this section, we deal with numerical experiments using the wavelet adaptive-homotopy method. We choose to recover the porosity for the elastic wave equations in the fluid-saturated porous media. In the first example, let us consider a homogenous model. The source function we selected is Ricker function with the amplitude 0.8 m and the frequency 40 Hz, which is shown in Fig. 1. The parameters describing the physical properties of the medium are given by
kb ¼ 3:3568 106 ;
l ¼ 2:32 106 ; kf ¼ 1:25 106 Pa; ks ¼ 6:296106 Pa;qf ¼ 1:000 kg=m3 ; qs ¼ 2:400 kg=m3 : In this numerical test, we choose the spatial step size hx ¼ hz ¼ 10 m, the temporal step size s ¼ 0:001 s and 32 grid points in the horizonal and vertical directions, respectively. We decompose the original problem into five scales using the wavelet
Surface axis (×10m) 0
0
5
10
15
20
25
30
Depth axis (×10m)
5
0.4
10
0.35
15
0.3
20
0.25
25 0.2 30 Fig. 5. The inversion result at scale 1.
Surface axis (×10m) 0
0
5
10
15
20
25
30 0.4 0.38
Depth axis (×10m)
5
0.36 10
0.34 0.32
15
0.3 0.28
20
0.26 0.24
25
0.22 0.2
30 Fig. 6. The inversion result at the original scale.
196
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
‘DB4’, and then apply the wavelet adaptive-homotopy method proposed in Section 4 to recover the porosity in the fluid-saturated porous media. The inversion results are shown in Table 2. In order to perform the advantages of the wavelet adaptivehomotopy method, we also apply the homotopy method for recovering the porosity and the results are shown in Table 3. Comparison of Table 2 with Table 3 on the aspects of relative errors and iterative times, we find that the wavelet adaptive-homotopy method is faster than the homotopy method when searching the global minimum. In the second example, a complicatedly exact porosity model is considered in Fig. 2. The parameters are the same as the first example. Note that the digitization process would have added 5% random noisy to the observed data, and then we recovery the porosity from noise data. We decompose the original problem into three scales and the inversion results at different scales are shown in Figs. 3–6. These results illustrate the fact that the method has the ability of noise suppression. The above numerical results show that the wavelet adaptive-homotopy method is globally convergent and it is feasible to use this method to recover the porosity in the fluid-saturated porous media. 6. Conclusions In order to accelerate convergence and overcome the problem of local minimum, we have designed the wavelet adaptivehomotopy method consisting of the wavelet multiscale technique and the adaptive-homotopy method. And then we carry out numerical simulations by using the proposed method. The present results indicate that the method has clear advantages in comparison with other well established numerical algorithms. Numerical results demonstrate the ability of the method to recovery the parameter for the elastic wave equation in the fluid-saturated porous media. References [1] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid: low-frequency range, J. Acoust. Soc. Am. 28 (1956) 168–178. [2] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid: higher-frequency range, J. Acoust. Soc. Am. 28 (1956) 179–191. [3] K.A. Liu, H.W. Liu, B.Q. Guo, Time convolution regularization iteration method for wave equation porosity inversion in 2-D two-phase medium, Oil Geophys. Prospect. 31 (1996) 410–414 (in Chinese). [4] K.A. Liu, H.W. Liu, Method of perturbation for inversion of 2-D elastic wave equation in a two-phase medium, J. Harbin Univ. Architect. Eng. 29 (1996) 80–84 (in Chinese). [5] P.J. Wei, Z.M. Zhang, H. Han, Inversion of fluid-saturated porous media based on neural networks, Acta Mech. Solid Sin. 15 (2002) 342–349. [6] J. Liu, A multiresolution method for distributed parameter estimation, SIAM J. Sci. Comput. 14 (1993) 389–405. [7] F. Ebrahimi, M. Sahimi, Multiresolution wavelet scale up of unstable miscible displacements in flow through heterogeneous porous media, Transport Porous Med. 57 (2004) 75–102. [8] X.M. Zhang, K.A. Liu, J.Q. Liu, The wavelet multiscale method for inversion of porosity in the fluid-saturated porous media, Appl. Math. Comput. 180 (2006) 419–427. [9] S. Li, B. Wang, J.H. Hu, Homotopy solution of the inverse generalized eigenvalue problems in structural dynamics, Appl. Math. Mech.-Engl. 25 (2004) 580–586. [10] B. Han, H.S. Fu, Z. Li, A homotopy method for the inversion of a two-dimension wave equation, Inverse Probl. Sci. Eng. 13 (2005) 411–431. [11] B. Richter, Homotopy algebras and the inverse of the normalization functor, J. Pure Appl. Algebra 206 (2006) 277–321. [12] H.J. Su, J.M. McCarthy, A polynomical homotopy formulation of the inverse static analysis of planar compliant mechanisms, J. Mech. Des. 128 (2006) 776–786. [13] B. Han, H.S. Fu, H. Liu, A homotopy method for well-log constraint waveform inversion, Geophysics 72 (2007) R1–R7. [14] H. Han, Z.M. Zhang, Y.S. Wang, Homotopy method for inversing the porosity of 1D wave equation in porous media, Acta Mech. Sin. 35 (2003) 235–239. [15] S.G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Mach. Intell. 11 (1989) 674–693.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
A wavelet adaptive-homotopy method for inverse problem in the fluid-saturated porous media q Ying He *, Bo Han Department of Mathematics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e
i n f o
Keywords: Wavelet Adaptive-homotopy method Inversion Fluid-saturated porous media
a b s t r a c t In this paper, we consider a parameter identification problem for the elastic wave equations in the fluid-saturated porous media. A wavelet adaptive-homotopy method is proposed for the recovery of porosity. This algorithm combines the wavelet multiscale inversion idea with the adaptive-homotopy method. And the adaptive-homotopy method provides a simple way to adapt computational refinement to the choice of the homotopy parameter. Numerical results illustrate that the wavelet adaptive-homotopy method is globally convergent and highly effective for solving the inverse problem. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction The elastic wave propagation in porous media has received considerable attention in resent years because of increasing applications in science and engineering. Applications of porous media include Mechanics, Engineering, Geosciences, Material science. In many engineering problems, such as those found in petroleum and geotechnical engineering, rocks have to be modeled as porous materials saturated with fluid. In 1956, Biot [1,2] firstly considered the propagation of elastic wave in the fluid-saturated porous media. He summarized a general and systematic theory to describe a porous solid containing a viscous fluid. Since then, various numerical techniques have appeared in the literatures of the inverse problem in porous media. The inverse problem in the fluid-saturated porous media can be viewed as a parametric data-fitting problem. It is possible to formalize such a problem in the framework of optimization where a functional defined in terms of discrepancy between observed and computed data is minimized over a model space. In general, such problem is very difficult to solve, since it is nonlinear and ill-posed. The nonlinear property causes the presence of numerous local minimum, while the illposedness often affects the efficiency of the numerical scheme significantly. For the nonlinear and ill-posed problem, traditional linearized inversion algorithms, such as conjugate gradient method, time convolution iteration method [3] and perturbation method [4], are locally convergent. The widely popular recent approaches such as simulated annealing algorithm, genetic algorithm, trust region method, and neural networks [5], are globally convergent. However, along with the searching space decreasing, the efficiency of these are worse than before. Therefore, the shortages of the above inversion methods motive us to design a highly efficient, numerically stable and globally convergent algorithm. Wavelet multiscale method is a newly developmental inversion strategy which accelerate convergent rate, enhance stability, and overcome disturbance of local minimum. Because of these features, the applications of the wavelet multiscale inversion approach are giving rise to much attention. Liu [6] proposed wavelet multiresolution method for distributed parameter estimation of the elliptic problem. Ebrahimi and Sahimi [7] has described scale up of geological models of field-scale porous media using a method based on the wavelet transformation. Zhang et al. [8] firstly recovered the porosity q
Project supported by the National Natural Science Foundation of China (No. 40774056). * Corresponding author. E-mail address: [email protected] (Y. He).
0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.11.033
190
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
in the fluid-saturated porous media by using wavelet method associated with the Gauss–Newton iterative method, and the results showed that the Gauss–Newton method performed much better when employed with wavelet multiscale method. Homotopy method is a powerful tool for solving inverse problems due to its globally convergent property. It is a very interesting research to apply homotopy method to solve the inverse problem. Various formulations for the homotopy method have been used in the literature [9–13]. Although the homotopy method is applied in many fields, the literatures about the availability of the homotopy method in the fluid-saturated porous media have rarely been found except [14], where Han et al. recovered the parameter of 1D wave equation. All these works testified that the homotopy method is effective for the inverse problem. As mentioned above, in this paper, we will design an inversion strategy called wavelet adaptive-homotopy method, for parameter identification problem in the fluid-saturated porous media. This new algorithm combines the property of good stability and fast convergent rate of wavelet multiscale method with the advantage of global minimum of homotopy method. In order to search the global minimum efficiently, we provide a simple way to adapt computational refinements to the choice of the homotopy parameter. The structure of the paper is organized as follows: we start with presenting the elastic wave equations and correlative parameters in the fluid-saturated porous media based on Biot theory in Section 2. In Section 3, an adaptive-homotopy method is proposed and the adaptive property of this algorithm provides a simple way to adapt a computational refinements to the choice of the homotopy parameter. In Section 4, a wavelet adaptive-homotopy is proposed for solving the inverse problem in the fluid-saturated porous media. Then we describe how this algorithm can be carried out by combining the wavelet multiscale method with the adaptive-homotopy method. Finally, in order to illustrate the performance of this algorithm, we carry out some numerical simulations for the inversion of porosity in fluid-saturated porous media. The results obtained show that the wavelet adaptive-homotopy method is effective and feasible to the inverse problem. 2. Governing equations The propagation of elastic wave in the fluid-saturated porous media is different from that in the single-phase media. Based on continuum mechanics and macroscopic constitutive relationship, Biot developed a theory of wave motion in an elastic solid saturated with a viscous compressible fluid. We ignore viscous effects and consider the 2D elastic wave propagation equations in fluid-saturated porous media as follows:
2
o ox
o ox
oux o ou ou o oux ouz o ow ow o2 u o2 w l xþl z þ aM x þ aM z ¼ q 2x þ qf 2 x f1 ; k þ þk þ oz ox ox ox oz ox ox oz ox oz ot ot
ð2:1Þ
oux ouz o ou o oux ouz o ow ow o2 u o2 w l z þ aM x þ aM z ¼ q 2z þ qf 2 z f2 ; k þl þ2 þk þ oz oz oz oz ox oz ox oz ox oz ot ot
ð2:2Þ
l
l
o ou ou ow ow o2 u o2 w aM x þ aM z þ M x þ M z ¼ qf 2x þ m 2 x f1 ; ox ox oz ox oz ot ot
ð2:3Þ
o ou ou owx owz o2 uz o2 wz aM x þ aM z þ M þM ¼ qf 2 þ m 2 f2 : oz ox oz partialx oz ot ot
ð2:4Þ
where u ¼ ðux ; uz ÞT denotes the solid-frame displacement, w ¼ ðwx ; wz ÞT is the fluid displacement relative to solid-frame, k is coefficient, the quantity j is the Darcy permeability coefficient, g is the viscosity of the pore fluid, qf is the density the Lame of the pore fluid, qs is the density of the solid grain, q is the bulk density defined by q ¼ bqf þ ð1 bÞqs , and b is the porosity. The relationship of the tortuosity a, added mass density m, and coupling constant M can be given by
a¼1
Ks ; Kr
M¼
K 2r ; Dr K s
Kr Dr ¼ K r 1 þ b 1 ; Kf
m¼
qf b
;
where K r is the bulk modulus of the grain, K f is the bulk modulus of the pore fluid, K s is the bulk modulus of the skeletal frame. The boundary conditions for the problem considered here can be given by
oux ðx; z; tÞ oux ðx; z; tÞ ¼ ¼ 0; ox ox x¼L x¼0 oux ðx; z; tÞ ¼ 0; oz z¼H ouz ðx; z; tÞ ouz ðx; z; tÞ jx¼0 ¼ ¼ 0; ox ox x¼L
ð2:5Þ ð2:6Þ ð2:7Þ
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
ouz ðx; z; tÞ ¼ 0; oz z¼H owx ðx; z; tÞ owx ðx; z; tÞ ¼ ¼ 0; ox ox x¼0 x¼L owx ðx; z; tÞ ¼ 0; oz z¼H owz ðx; z; tÞ owz ðx; z; tÞ ¼ ¼ 0; ox ox x¼0 x¼L owz ðx; z; tÞ ¼ 0; oz z¼H
191
ð2:8Þ ð2:9Þ ð2:10Þ ð2:11Þ ð2:12Þ
where ðx; zÞ 2 ½0; L ½0; H; t 2 ½0; T. The initial conditions are given as follows:
oux ðx; z; 0Þ ¼ 0; ot ouz ðx; z; 0Þ uz ðx; z; 0Þ ¼ 0; ¼ 0; ot owx ðx; z; 0Þ ¼ 0; wx ðx; z; 0Þ ¼ 0; ot owz ðx; z; 0Þ ¼ 0; wz ðx; z; 0Þ ¼ 0; ot ux ðx; z; 0Þ ¼ 0;
ð2:13Þ ð2:14Þ ð2:15Þ ð2:16Þ
and the additional condition is
ux ðx; 0; tÞ ¼ ux ðx; tÞ;
ð2:17Þ
where ux is the observed data. If given a value of the porosity b, Eqs. (2.1)–(2.6), (2.8), (2.7), (2.9)–(2.16) enable to form the direct problem for wave fields simulations, i.e. we can solve the wave field functions ux ; uz ; wx ; wz using a numerical method. Consequently, the additional condition ux ðx; 0; tÞ can be obtained. Note that since the porosity b is unknown in practice, governing equations (2.1)–(2.4) with the boundary conditions (2.5), (2.6), (2.8), (2.7), (2.9)–(2.12), the initial conditions (2.13)–(2.16) and the additional condition (2.17) completely define the inverse problem for the elastic wave equations in the fluid-saturated porous media. 3. An adaptive-homotopy method In general, the inverse problem in the fluid-saturated porous media can be viewed as a parametric data-fitting problem. It is possible to formalize such a problem in the framework of optimization where a function defined in terms of discrepancy between observed and computed data is minimized over a model space. Now, we discretize Eqs. (2.1)–(2.6), (2.8), (2.7), (2.9)–(2.17) by the second order finite difference method as follows:
ux ði; j; k þ 1Þ ¼ wx ði; j; k þ 1Þ ¼
mU 1 qf U 2 mqði; jÞ q2f
;
uz ði; j; k þ 1Þ ¼
mU 3 qf U 4 mqði; jÞ q2f
;
qði; jÞU 2 qf U 1 qði; jÞU 4 qf U 3 ; wz ði; j; k þ 1Þ ¼ ; 2 mqði; jÞ qf mqði; jÞ q2f
ux ð0; j; kÞ ¼ ux ð1; j; kÞ;
ux ðp 1; j; kÞ ¼ ux ðp; j; kÞ;
ux ði; n 1; kÞ ¼ ux ði; n; kÞ;
uz ð0; j; kÞ ¼ uz ð1; j; kÞ;
uz ðp 1; j; kÞ ¼ uz ðp; j; kÞ;
uz ði; n 1; kÞ ¼ uz ði; n; kÞ;
wx ð0; j; kÞ ¼ wx ð1; j; kÞ;
wx ðp 1; j; kÞ ¼ wx ðp; j; kÞ;
wx ði; n 1; kÞ ¼ wx ði; n; kÞ;
wz ð0; j; kÞ ¼ wz ð1; j; kÞ;
wz ðp 1; j; kÞ ¼ wz ðp; j; kÞ;
wz ði; n 1; kÞ ¼ wz ði; n; kÞ;
ux ði; j; 0Þ ¼ ux ði; j; 1Þ ¼ 0; wx ði; j; 0Þ ¼ wx ði; j; 1Þ ¼ 0; ux ði; 0; kÞ ¼
uz ði; j; 0Þ ¼ uz ði; j; 1Þ ¼ 0; wz ði; j; 0Þ ¼ wz ði; j; 1Þ ¼ 0;
ux ði; kÞ;
where b is the model to be recovered, ux is the observed data, ux ði; j; kÞ ¼ ux ði hx ; j hz ; k sÞ, hx ; hz are the step sizes of the rectangle grid in the x- and z-direction, respectively, s is the time step size, p ¼ hLx ; n ¼ hHz ; l ¼ Ts. The above difference equations define an operator equation
F : b ! U;
i:e: FðbÞ ¼ U;
192
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
where b and U denote vectors, which are composed of bði; jÞ and ux ði; kÞ in the following sequence:
bT ¼ ðbð1; 1Þ; . . . ; bð1; nÞ; bð2; 1Þ; . . . ; bð2; nÞ; . . . ; bðp; nÞÞ; U T ¼ ðux ð1; 1Þ; . . . ; ux ðp; 1Þ; ux ð2; 1Þ; . . . ; ux ðp; 2Þ; . . . ; ux ðp; lÞÞ; where p; n; l denote the number of mesh nodes for the x; z; t components, respectively. In order to identify the porosity b from the observed data U containing noise, we define a nonlinear operator equation
TðbÞ , FðbÞ U ¼ 0:
ð3:1Þ
Eq. (3.1) is usually highly nonlinear and ill-posed, therefore we consider a globally convergent algorithm, namely an adaptive-homotopy method, to solve the nonlinear operator equation (3.1). Let us construct a homotopy regularization method
Hðb; aÞ ¼ ð1 aÞkFðbÞ Uk2 þ akb b0 k2 ;
ð3:2Þ
0
where b is the initial guess, the parameter a 2 ½0; 1 is not only the homotopy parameter bur also the regularization parameter. When a ¼ 1, the solution of Hðb; 1Þ ¼ 0 is b0 , which is known. When a ¼ 0, the solution of Hðb; 0Þ ¼ 0 is the desired solution of the nonlinear operator equation (3.1). Above all, we divide the interval [0, 1] into 1 ¼ a0 > a1 > > aN ¼ 0. For some ak , we use some iterative method to solve Hðb; ak Þ ¼ 0 sequently. As the solution b0 of Hðb; 1Þ ¼ 0 is known, it can be used as the initial approximation of the next equation. Assume that the approximation solution bk of the Hðb; ak Þ ¼ 0 has already been found. For computing the solution bkþ1 of Hðb; akþ1 Þ ¼ 0, we use the damped Gauss–Newton method. In practice, we need not compute the accurate solution of each equation. So we only iterate one step when using the damped Gauss–Newton method to solve the equation Hðb; ak Þ ¼ 0. The iterative formulation is given as follows:
bkþ1 ¼ bk ½ð1 ak ÞF 0 ðbk Þ F 0 ðbk Þ þ ak I1 ½ð1 ak ÞF 0 ðbk Þ ðFðbk Þ UÞ þ ak bk ;
k ¼ 0; 1; . . . ; N 1:
ð3:3Þ
From formulation (3.3), the approximation solution bN of the Hðb; 1Þ ¼ 0 can be solved, which can be taken as the initial guess to the problem (3.1). Then iterate repeatedly by using the damped Gauss–Newton algorithm
bkþ1 ¼ bk ½F 0 ðbk Þ F 0 ðbk Þ1 ½F 0 ðbk Þ ðFðbk Þ UÞ;
k ¼ N; N þ 1; . . .
ð3:4Þ
The formulations (3.3) and (3.4) constitute a globally convergent method for the inverse problem. As we all know, the homotopy parameters ak ðk ¼ 0; 1; . . .Þ plays an important role. But the homotopy parameters always can be given by ak ¼ Nk in the existing literatures. Such choice may loss some convergent approximation solutions of the Eq. (3.1). Therefore, in this paper, we propose a modified homotopy method named by the adaptive-homotopy method and apply the new approach to find the global minimum of the inverse problem. This algorithm can adaptively choose the homotopy parameters ak . Algorithm 1. Choose an initial model b0 , small threshold e and initial homotopy step size h. Set a1 ¼ a0 h, a2 ¼ a1 h. Calculate the solution b1 of Hðb; a1 Þ ¼ 0 (using damped Gauss–Newton method). For k ¼ 2; . . ., do:
Calculate the solution bk of Hðb; ak Þ ¼ 0 (using damped Gauss–Newton method). kbk bk1 k . Calculate j ¼ Vertb k1 bk2 k Check if j < e, then set h ¼ 2h, else if j P e, then set h ¼ 2h. Check the convergence, if not converged: set k ¼ k þ 1 ak ¼ ak1 h: and go to step 1.
4. Wavelet adaptive-homotopy method The wavelet multiscale method is a technique that improves the performance of inversion method by decomposing the problem into subproblem on different scales. At coarser scale there are fewer local minima and those that remain are further apart from each other, and then inversion method can get the solution closer to the neighborhood of the global minimum. As the scale increases, the complexity of the original problem is reduced and the size of the corresponding domain of convergence is decreasing. In this section, we try to use wavelet multiscale method to decompose the original operator equation (3.1). Once the problem has been decomposed into different subproblems, it is natural to solve first the subproblem on Table 1 The decomposition of the original problem by wavelet. Scale space
Source
Wave field
Solved problem
Solution of the problem
Initial guess
V0 V1 V2
f f1 f2
ux u1 x u2 x
FðbÞ ¼ U F 1 ðbÞ ¼ U 1 F 2 ðbÞ ¼ U 2
b b1 b2
b1 b2 b0
193
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0
20
40
60
80
100
120
140
160
180
200
Time axis (0.001s) Fig. 1. The source function.
Table 2 The results of inversion by wavelet adaptive-homotopy method. Number
Initial value
Real value
Iterative value
Iterative times
Relative error (%)
1 2 3 4
0.045 0.237 0.384 0.659
0.266000 0.266000 0.266000 0.266000
0.26472 0.26579 0.26678 0.26713
27 14 17 30
0.483 0.079 0.294 0.425
Table 3 The results of inversion by homotopy method. Number
Initial value
Real value
Iterative value
Iterative times
Relative error (%)
1 2 3 4
0.045 0.237 0.384 0.659
0.266000 0.266000 0.266000 0.266000
0.26464 0.26554 0.26714 0.26758
34 21 26 39
0.512 0.174 0.427 0.593
0
0
5
10
Surface axis (×10m) 15
20
25
0.4 0.38
5
Depth axis (x10m)
30
0.36 0.34
10
0.32 15
0.3 0.28
20
0.26 25
0.24 0.22
30
0.2 Fig. 2. The true value of porosity model.
194
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
the coarsest scale by using adaptive-homotopy method, and to use the resulting solution as initial guess for the subproblem on the finer scale. Without loss of generality, we describe the wavelet adaptive-homotopy method by decomposing the original problem into two scales. Assume that the original problem is on the scale space V 0 and decomposed on the coarser scale space V 1 as well as the coarsest scale space V 2 . We utilize the compactly supported orthogonal wavelet ‘DB4’ to decompose the wave field functions and source function. The algorithm can be described as follows:
Step 3. Step 4.
Surface axis (x10m) 0
0
5
10
15
20
25
30 0.55
5
Depth axis (×10m)
Step 2.
Given an initial model b0 , the true porosity model b , the source function f. Calculate the forward problem to obtain the model data U. Decompose f and U at each scale, such as V 1 and V 2 by the compactly supported orthogonal wavelet ‘DB4’ (using Mallat algorithm [15]), i.e. the original problem (3.1) can be decomposed as F 1 ðbÞ ¼ U 1 and F 2 ðbÞ ¼ U 2 , which is shown in the following Table 1. At the coarsest scale space V 2 , calculate the solution b2 of F 2 ðbÞ ¼ U 2 with the initial guess b0 (using the adaptivehomotopy method). At the finer scale space V 1 , we regard the solution b2 as the initial guess, and then calculate the solution b1 of F 1 ðbÞ ¼ U 1 (using the adaptive-homotopy method).
0.5 0.45
10
0.4 0.35
15
0.3 0.25
20
0.2 0.15
25
0.1 0.05
30 Fig. 3. The inversion result at scale 3.
Surface axis (×10m) 00
5
10
15
20
25
30 0.45
5 0.4
Depth axis (×10m)
Step 1.
10 0.35 15
0.3
20
0.25
25
0.2 0.15
30 Fig. 4. The inversion result at scale 2.
195
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
Step 5.
At the initial space V 0 , we regard the solution b1 as the initial guess, and then calculate the solution b of FðbÞ ¼ U (using the adaptive-homotopy method).The obtained solution b is the solution of original problem (3.1).
5. Numerical simulations In this section, we deal with numerical experiments using the wavelet adaptive-homotopy method. We choose to recover the porosity for the elastic wave equations in the fluid-saturated porous media. In the first example, let us consider a homogenous model. The source function we selected is Ricker function with the amplitude 0.8 m and the frequency 40 Hz, which is shown in Fig. 1. The parameters describing the physical properties of the medium are given by
kb ¼ 3:3568 106 ;
l ¼ 2:32 106 ; kf ¼ 1:25 106 Pa; ks ¼ 6:296106 Pa;qf ¼ 1:000 kg=m3 ; qs ¼ 2:400 kg=m3 : In this numerical test, we choose the spatial step size hx ¼ hz ¼ 10 m, the temporal step size s ¼ 0:001 s and 32 grid points in the horizonal and vertical directions, respectively. We decompose the original problem into five scales using the wavelet
Surface axis (×10m) 0
0
5
10
15
20
25
30
Depth axis (×10m)
5
0.4
10
0.35
15
0.3
20
0.25
25 0.2 30 Fig. 5. The inversion result at scale 1.
Surface axis (×10m) 0
0
5
10
15
20
25
30 0.4 0.38
Depth axis (×10m)
5
0.36 10
0.34 0.32
15
0.3 0.28
20
0.26 0.24
25
0.22 0.2
30 Fig. 6. The inversion result at the original scale.
196
Y. He, B. Han / Applied Mathematics and Computation 208 (2009) 189–196
‘DB4’, and then apply the wavelet adaptive-homotopy method proposed in Section 4 to recover the porosity in the fluid-saturated porous media. The inversion results are shown in Table 2. In order to perform the advantages of the wavelet adaptivehomotopy method, we also apply the homotopy method for recovering the porosity and the results are shown in Table 3. Comparison of Table 2 with Table 3 on the aspects of relative errors and iterative times, we find that the wavelet adaptive-homotopy method is faster than the homotopy method when searching the global minimum. In the second example, a complicatedly exact porosity model is considered in Fig. 2. The parameters are the same as the first example. Note that the digitization process would have added 5% random noisy to the observed data, and then we recovery the porosity from noise data. We decompose the original problem into three scales and the inversion results at different scales are shown in Figs. 3–6. These results illustrate the fact that the method has the ability of noise suppression. The above numerical results show that the wavelet adaptive-homotopy method is globally convergent and it is feasible to use this method to recover the porosity in the fluid-saturated porous media. 6. Conclusions In order to accelerate convergence and overcome the problem of local minimum, we have designed the wavelet adaptivehomotopy method consisting of the wavelet multiscale technique and the adaptive-homotopy method. And then we carry out numerical simulations by using the proposed method. The present results indicate that the method has clear advantages in comparison with other well established numerical algorithms. Numerical results demonstrate the ability of the method to recovery the parameter for the elastic wave equation in the fluid-saturated porous media. References [1] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid: low-frequency range, J. Acoust. Soc. Am. 28 (1956) 168–178. [2] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid: higher-frequency range, J. Acoust. Soc. Am. 28 (1956) 179–191. [3] K.A. Liu, H.W. Liu, B.Q. Guo, Time convolution regularization iteration method for wave equation porosity inversion in 2-D two-phase medium, Oil Geophys. Prospect. 31 (1996) 410–414 (in Chinese). [4] K.A. Liu, H.W. Liu, Method of perturbation for inversion of 2-D elastic wave equation in a two-phase medium, J. Harbin Univ. Architect. Eng. 29 (1996) 80–84 (in Chinese). [5] P.J. Wei, Z.M. Zhang, H. Han, Inversion of fluid-saturated porous media based on neural networks, Acta Mech. Solid Sin. 15 (2002) 342–349. [6] J. Liu, A multiresolution method for distributed parameter estimation, SIAM J. Sci. Comput. 14 (1993) 389–405. [7] F. Ebrahimi, M. Sahimi, Multiresolution wavelet scale up of unstable miscible displacements in flow through heterogeneous porous media, Transport Porous Med. 57 (2004) 75–102. [8] X.M. Zhang, K.A. Liu, J.Q. Liu, The wavelet multiscale method for inversion of porosity in the fluid-saturated porous media, Appl. Math. Comput. 180 (2006) 419–427. [9] S. Li, B. Wang, J.H. Hu, Homotopy solution of the inverse generalized eigenvalue problems in structural dynamics, Appl. Math. Mech.-Engl. 25 (2004) 580–586. [10] B. Han, H.S. Fu, Z. Li, A homotopy method for the inversion of a two-dimension wave equation, Inverse Probl. Sci. Eng. 13 (2005) 411–431. [11] B. Richter, Homotopy algebras and the inverse of the normalization functor, J. Pure Appl. Algebra 206 (2006) 277–321. [12] H.J. Su, J.M. McCarthy, A polynomical homotopy formulation of the inverse static analysis of planar compliant mechanisms, J. Mech. Des. 128 (2006) 776–786. [13] B. Han, H.S. Fu, H. Liu, A homotopy method for well-log constraint waveform inversion, Geophysics 72 (2007) R1–R7. [14] H. Han, Z.M. Zhang, Y.S. Wang, Homotopy method for inversing the porosity of 1D wave equation in porous media, Acta Mech. Sin. 35 (2003) 235–239. [15] S.G. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Trans. Pattern Anal. Mach. Intell. 11 (1989) 674–693.