Three-dimensional Magnetotelluric Inversion Based On Limited Memory Quasi-Newton Method

The three-dimensional magnetotelluric inversion technology has been greatly improved and proven to be a reliable data interpretation method due to continual struggling researches by professors worldwide in recent years. Both modeling and synthesized data inversion shows 3-D inversion can give or recover the true resistivity model beneath the surface to let the interpreters get more confidential on the results they have. But the technology has still not been applied widely on real data inversion because of two aspects: first, many 3D inversion program need huge amount of memory and runs slowly which can only be installed and launched on high level workstation or pc-clusters, another reason is most of the algorithms been published were still on the synthesized data inversion stage, which means they do not have a convenient mechanism to introduce a prior information acquired by geologists and geophysicists before real data inversion. Aim on that two shortages, we have improved the 3-D staggered grid finite difference method to do data modeling, which means the memory requirement has been greatly shorten while speed of solving the finite difference linear system was heighten a lot. We using the limit memory Quasi-Newton method with simple upper and lower bound constraints to make inversion objective functional minimum, the converge rate is much better than other method based on gradient direction.The frequency domain Maxwell equation has been expanded under regular mesh condition first, and then we use the scaling transform to get the expanded form of equations on ill-regular mesh; the modeling linear systems’ matrix can be calculated at last. Because this matrix cannot satisfy the band width minimum, we decide to solve it by traditional iterative method instead of multi-front direct solving method which is most popular in recent years. We use the standard CSR format to stored only the non-zero elements of matrix’s upper tri-angle part, all matrix and vector computing function has been reprogrammed such as incomplete Cholesky decompose of the symmetric hermitain matrix and its’ fast forward & backward substitution method which only need the upper triangle non-zero elements. We also deeply researched a problem of how frequently to do static divergence correction is the best to get the most efficient scheme of iterative solver of finite difference linear system. After electric field of two polarity modes been solved, a set of linear interpolate vectors is been multiplied to get every components of both electric and magnetic field right on the observing site, then we using 3-D tensor impedance and tipper formula to calculate magnetotelluric response. We use a 1-D multi-layered models analytic solution, a 2-D vertical fault system’s pseudo-analytic solution and a 3-D international standard test-model synthesized data from other professors to test our 3-D code, the results shows our data modeling program have good accuracy to fulfill the requirement of 3-D data inversion.Based on the data modeling code mentioned before, different from other researchers who chose typical OCCAM or NLCG method, after compare analyze of traditional inversion theory, we chose the BFGS formula of Quasi-Newton method to apply our 3-D magnetotelluric inversion. Consider of 3D inversion is a big scale minimize problem, the limited memory BFGS strategy has been introduced to lower the memory requirement as to NLCG. Then, to avoid the shortage of transforming the minimize problem to unconstrained problem by log transform, we developed a direct simple upper and lower bound constrained LBFGS algorithm by imply a projection sub-space, this procedure first use unconstrained minimum direction to get a set of parameters which could be extend on upper or lower bounds, then the hessian matrix has been appropriate in the sub-space of parameters left. We designed the inversion kernel as a callback function to maintain loose coupling between the minimizing method and the code to calculate objective functional and its’ gradient. This design makes our inversion kernel can easily been modified to apply on other geophysical inversion problems. We developed a more general regularized operator to suit on geological interpretation which is calculated base on pre-set cells’ properties, and a frequency based parallel strategy under Open MP protocol was implied to efficiently use every core of CPU. All parallel compute related functions have been designed to have independent memory management to convenient the modification adapt to other parallel computation platform. During the research, we also did a lot of primal test on line search method under Wolfe condition, self-adapt regularizing factor, and precondition of decent direction. Finally, our 3-D inversion code have two more constraints, less memory requirement, less mount of line search times, and also more efficient converge rate.After the development of 3-D magnetotelluric modeling & inversion code, we use a synthesized data from simple low resistivity cube model which used by other researchers to test and analyze our inversion algorithm, the dependence of initial guess model and rationality of regularizing factor self-adjusting strategy have been discussed. Then we chose a standard model above which the observing sites can not cover the surface completely to test ability of our inversion code to adopt the real geological exploration mission. Our test results shows that after trial inversion computation to obtain simple constraint, the 3-D inversion still can recover the real resistivity model even the observed data cannot full cover the whole model. By the way, through analyzing of trail inversions, because of the gradient of objective functional calculation is deeply related with surface resistivity, we can make a conclusion that 3-D inversion’s initial model should keep the surface resistivity as accurate as possible, meanwhile, we found after we obtain some rough space interval of partial resistivity changing, the inversion results can be remarkable improved.Finally, applying our 3-D modeling & inversion code, we analyze a real 3D AMT data set from one Kimberley exploration project of southwest China. First, we provide a simple correction technology of impedance phase invariant based on 3-D modeling, which can recover the low resistivity’s phase information overwhelmed by complex topography. After that, the corrected phase invariant can be used to qualitatively analyze the traditional 1-D & 2-D inversion results of this area to avoid the pseudo images of these inversion methods in real 3-D model’s data. Then we use different initial model to make trial 3-D inversion, all results shows a low resistance pipe shape abnormality, based on the information we get from the trial, the last constrained 3-D inversion completely foreclose pseudo images of 1-D & 2-D inversion. The 3-D results can perfectly coherent with qualitative analyzing result by corrected impedance phase invariant.