3-D inversion of magnetic data from Sorkheh-Dizej region, Zanjan, using the nonlinear Marquardt-Levenberg algorithm

The aim of this study is 3-D inversion of magnetic data acquired from the Sorkheh-Dizej iron-bearing region in Zanjan province, Iran, in order to determine the geometrical distribution of the buried magnetic sources. For this purpose, our program discretizes the subsurface region into vertical right rectangular prisms and uses the nonlinear Marquardt-Levenberg algorithm to optimize the unknown parameters of the model, iteratively. To evaluate the applicability of the method, it has been first applied to modeling the synthetic data with added noise.
The nonlinear inversion of potential field data has long been used for determining the unknown parameters of the source bodies. One of the major applications of this method is to determine the topography of the basement relief (see, e.g., Bott, 1960; Pedersen, 1977; Bhattacharayya, 1980). In this study, we use a similar approach for modeling the magnetic sources with complicated geological shapes. To do this, the region is discretized into vertical right rectangular prisms. As a result, the unknown geometrical parameters will be the depth to the top and depth to the bottom of each prism. The only other unknown parameter is magnetization which is constant for each prism
A right rectangular prism is a widely used geometrical model for 3-D interpretation of magnetic source bodies. Bhattacharayya (1964) presented equations for computing the total field magnetic anomaly due to prismatic models. Kunaratnam (1981) simplified the logarithmic and arctangent terms in these equations using complex notation. The program used in this study uses these simplified expressions for forward modeling of magnetic data.
To obtain the body parameters that best describe the observed data, our program uses the Marquardt-Levenberg optimization algorithm (Marquardt, 1963) in an iterative way to minimize the difference between the observed and calculated anomalies. The objective function to be minimized is the sum of the square-roots of the errors at all the data points (equation 1). Marquardt-Levenberg algorithm can be regarded as an interpolation between the Gauss-Newton and Steepest Descent methods. At points distant from the solution, this algorithm acts like the Steepest Descent method (for being faster), but as it approaches the solution it acts more like a Gauss-Newton technique, for being more accurate. This behavior is controlled by a damping factor which is automatically adjusted at each iteration according to the success or failure of the iteration in reducing the objective function.
The program that is used in this study is the modified version of the FORTRAN-77 program presented by Rao and Babu (1993). For forming the Jacobian matrices, the program uses the expressions for derivatives from Rao and Babu (1991) and it uses the Cholesky decomposition technique for the factorization of the coefficient matrix (matrix D in equation 4). The input data to the program are the observed magnetic data and the magnetic properties of the region (i.e., magnetic declination and inclination and the regional constant). The program then generates an initial model of adjacent vertical prisms with similar upper and lower depths. There should be one data point above each prism on the surface of the earth. The goodness of fit between the observed data and calculated data is computed and then solving the inverse problem (equation 3) the increments that should be added to each body parameter to gain a smaller objective function in the next iteration are obtained. The iterations are continued until the desired objective function is reached or this residual is negligible. The geometrical configuration of the collection of the prisms at the end of the inversion shows the distribution of the magnetic source body.
In order to evaluate the applicability of the method, first we apply it to modeling the synthetic data. The synthetic model consists of three separate rectangular blocks which represent a complex geological configuration (Fig. 1). The blocks have the same magnetization intensity (10 A/m) and the declination and inclination of the magnetization are assumed to be zero and 45 degrees, respectively. In order to make the generated synthetic data resemble the realistic field data, random noise with zero mean and standard deviation of 5 percent of each datum magnitude has been added to the data set. The terrain is discretized into a grid of 10x10 prisms, each with a surface data centered above it. The magnetization is supposed to be known and constant throughout the body. The result of the inversion after 100 iterations has a good similarity to the original model although it can be seen that the accuracy reduces with increasing depth (Fig. 3).
We apply the inversion scheme to model the real field magnetic data from Sorkheh-Dizej region in Zanjan province, Iran. The shape of the anomaly map (the region isolated by the rectangle in Fig. 4) is typical of a large tabular body. As non-uniqueness is one of the major concerns in the inversion of potential field data, it is useful to limit the possible solutions by devising constraints on the variation of magnetization during the inversion procedure. The source body is completely buried with no outcrops, but the existence of several iron mines very close to the survey area, and the high amplitude of the anomaly intensity, suggest a similar genesis for the magnetic bodies of the region. Therefore, we performed magnetization intensity measurements on 15 core samples from the region which resulted in a range of 5-15 A/m for this parameter. We assume no remanent magnetization is present and the magnetization is only due to induction. The terrain is divided into a 17x15 grid of prisms and inversion performed for 100 iterations. Fig. 7 shows a view of the result body from East. The vertical extension of this dyke-like body is interpreted between -10 and -210 meters, and the average dip-angle of the body is 70 degrees towards North. The results have been compared with the results of the 3-D compact inversion of pseudo-gravity data. The good agreement between the two models shows the efficiency of the algorithm that is employed for the inversion of geomagnetic data in this study.