A Rapidly Direct Algorithm of Explicit Fréchet Derivatives for MCSEM in a 3-D TI Formation Using Finite Volume of Coupled Potentials

This paper presents a rapidly direct algorithm (RDA) to calculate explicit Fréchet derivatives (EFD) for marine controlled-source electromagnetic measurement (MCSEM) in a three-dimensional (1-D) transversely isotropic (TI) formation. By discretizing the Helmholtz equations about coupled potentials on Yee’s staggered grids by the 3-D finite volume method (FVM), we obtain a complex linear system about the unknown potentials excited by a mass of moving electric current sources. To efficiently determine electromagnetic (EM) fields, we introduce an interpolation operator and projection operator per receiver by using the direct solver PARDISO and 3-D Newtown interpolation. Based on this, the perturbation in goal conductivity is expressed as a piece-wise constant function according to block or pixel model. The spatial scattered electric currents will be decomposed into a series of electric current elements distributed on Yee’s grids due to the conductivity perturbation. We then discretize the scattered electric currents by 3-D FVM and obtain the new right-hand terms about the unknown scattered potentials. This allows for fast production of the linear relationship between the changes in the EM fields and the relative conductivity perturbation per block or pixel, ultimately resulting in the EFD of MCSEM responses. Numerical results demonstrate the efficiency and accuracy of this method. The 3-D pixel sensitivities are presented to further investigate the response characteristics of MCSEM in several cases.


I. INTRODUCTION
The marine controlled-source electromagnetic measurement (MCSEM) is mainly applied in seabed structure exploration, characteristic identification, and quantitative evaluation of offshore oil and gas reservoirs to help reduce offshore blind drilling rate and exploration costs.The typical measurement method is to use the low-frequency (0.1-10 Hz) horizontal The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan .electric dipole transmitting antenna, which is towed by a ship, to continually excite alternating EM fields.The receivers, which are laid on the seafloor, are used to measure the induction EM fields with multi-frequency and multi-offset for determination of the spatial distribution of underground conductivity [1], [2].Due to the complex topographic structure of the seabed and the inhomogeneous conductivity distribution of the sediments, a lot of numerical simulation, Fréchet derivatives (FD) calculation, and inversion are required in the optimization design of signal acquisition and the interpretation of MCSEM data.Therefore, in recent decades, the corresponding theories and methods of MCSEM have been developed rapidly [3], [4].
In the forward modeling of MCSEM, analytical and numerical methods both have been widely studied and applied.For example, analytic methods, such as the transmission line method (TLM) [5], the propagation matrix method (PMM) [6], and the numerical mode matching method (NMM) [7], can be used to simulate the MCSEM responses with high accuracy in 1-D layered models or some axisymmetric models.For the arbitrary 3-D models, various 3-D numerical methods are usually required.The numerical algorithms chiefly include the 3-D finite element method (FEM) [8], the 3-D finite volume method (FVM) [9], [10], and the integral equation method (IEM) [11].In terms of the inversion of CSEM data, the various iterative methods (such as conjugate gradient or Gauss-Newton techniques) are often applied to realize the best fit between the theoretical synthetic data and the actual input data.Up to now, almost all inversion methods have been systematically studied in MCSEM and a lot of corresponding literatures have been published.For example, several representative papers about 1-D inversion include [12], [13], [14], and those about 2-D and 3-D inversions do [15], [16], [17], [18], [19], [20].During conductivity inversion imaging, the two different inversion models: pixels and blocks are usually used to describe the spatial distribution of formation conductivity [20], [21], [22].Furthermore, it must be pointed out that determining the descent direction of the objective function has been the key technique in EM inversion.Except that the FD in 1-D stratified formations and 2-D axisymmetric conductivity formations can be effectively and accurately calculated by the analytical method [13], [14], [23] or semi-analytical methods [24], the FD in the 2-D and 3-D EM inversion are usually determined by only some approximate methods, such as nonlinear conjugate gradient (NLCG) or quasi-Newton techniques [4], [20], [25], [26], [27] because the accurate algorithm of FD of EM response is usually much time-consuming.
In this paper, we will utilize the 3-D FVM of coupled potentials and a direct solver to establish a set of accurate and efficient algorithms of explicit Fréchet derivatives (EFD) for the 3-D pixel or block model of MCSEM in a 3-D transversely isotropic (TI) formation.We firstly discretize the Helmholtz equations about coupled potentials on Yee's staggered grids by the 3-D FVM and obtain a complex linear system about the unknown coupled potentials excited by a mass of moving electric current sources.To efficiently determine EM fields, we will introduce an interpolation operator and a projection operator per receiver by the combination of the direct solver PARDISO with the 3-D Newtown interpolation.Based on this, the perturbation in goal conductivity is expressed as a piece-wise constant function for both block and pixel models.The scattered electric current density can be decomposed into a series of electric current elements distributed on Yee's staggered grids.Each scattered current element is represented by the product of the relative perturbation in conductivity with the electric current elements on the grid.To establish a rapidly direct algorithm (RDA) of EFD, we then discretize the equations about scattered EM fields and obtain the new right-hand terms about the unknown scattered potentials.By using a projection operator per receiver, we achieve the linear relationship between the changes in the EM fields and the relative perturbation in conductivity on each block or pixel, which ultimately results in the EFD for MCSEM responses.We apply numerical results to validate the efficiency and accuracy of this new method.The 3-D pixel EFD are presented to further investigate the characteristics of MCSEM in several different cases.

II. THEORY
In this section, we first study the discretization and direct solution of coupled potential equations by the 3-D FVM of MCSEM in 3-D TI formation.Especially, we will introduce an interpolation operator and a projection operator through 3-D Newtown interpolation to enhance the efficiency of 3-D forward modeling.We then give the expressions of scattered electric currents according to block and pixel models and discretize the equations about scattered EM fields.By projection operator per receiver, we establish a rapidly direct approach of block and pixel EFD of MCSEM responses.

A. 3-D FINITE VOLUME FORWARD MODELING AND PROJECTION OPERATOR
In order to overcome the low induction number problem, we introduce the coupled potentials (A, φ) and obtain the following Helmholtz equations (the e iωt is assumed) [10]: and Here, J S = ILê x , IL is the electric dipole moment of the horizontal transmitting antenna, êx is a unit vector in the is the Laplace operator, and σ * (r) = diag(σ H (r) + iωε, σ H (r) + iωε, σ V (r) + iωε) is the complex conductivity tensor.σ H and σ V are respectively the horizontal and vertical conductivity.ε is the dielectric constant.µ 0 is the vacuum permeability.ω = 2πf is the angular frequency.r S is the transmitter position.
We use the symbol to denote a sufficiently large computation region.On its outer boundary ∂ , n×A(r, r S ) | ∂ = 0 and φ(r, r S ) | ∂ = 0 are assumed.Here, n is the unit normal vector on ∂ .The domain is divided into a series of staggered grids V i+1 / 2,j,k , V i,j+1 / 2,k , V i,j,k+1 / 2 , and V i,j,k .We then define the three mutually orthogonal components A x i+1/2,j,k , A y i,j+1/2,k , A z i,j,k+1/2 of the vector potential A(r) at the centers of the grids V i+1 / 2,j,k , V i,j+1 / 2,k , and V i,j,k+1 / 2 , and the scalar potential φ i,j,k at the center of the grid V i,j,k .We further calculate the equivalent conductivities σ * i+1 / 2,j,k , σ * i,j+1 / 2,k , σ * i,j,k+1 / 2 and σ * i,j,k of σ * (r) in all cells by the volume average [28].
By using the 3-D FVM, we can discretize (1) and (2) into the following linear algebraic equation about unknown coupled potentials [10], [28]: Here, F is a N × N complex sparse asymmetric coefficient matrix, x(r S ) is the N -dimensional unknown column vector composed of all discrete potentials A x i+1/2,j,k (r S ), A y i,j+1/2,k (r S ), A z i,j,k+1/2 (r S ), and φ i,j,k (r S ).N = m 1 + m 2 + m 3 + m 4 is the total number of unknowns where m 1 = (N x + 1)N y N z , m 2 = N x (N y + 1)N z , m 3 = N x N y (N z + 1), and m 4 = N x N y N z .b(r S ) represents the discrete vector of the right-hand terms of ( 1) and (2).
Because the MCSEM usually adopts moving sources, J S δ(r − r S,k ), (k = 1, 2, • • • , K ), the right-hand terms in (3) can be combined to form a N ×K -order sparse banded matrix B = b(r S,1 ), b(r S,2 ), • • • , b(r S,K ) .Similarly, all unknown vectors consisting of the left side of (3) are combined into a N × K -order matrix X = x(r S,1 ), x(r S,2 ), • • • , x(r S,K ) .Since the matrix F is fixed in the domain , we can solve the coupled potentials of all the transmitters simultaneously by PARDISO [29]: Here, the inverse matrix F−1 is obtained by PARDISO.
In addition, we assume that r R,l , l = 1, 2, • • • , L are all receiver positions in and L ≪ K .According to the central difference formula, we can discretize the following equations: at r R,l .We then obtain the interpolation operator Q E x (r R,l ) and Q H y (r R,l ), (l = 1, 2, • • • , L), to compute the values of E x (r R,l ) and H y (r R,l ) from the coupled potentials X.
Therefore, the values of E x (r R,l ) and H y (r R,l ) at r R,l by all transmitters are given by: Substituting ( 4) into (6) yields: Here, are respectively the electric and magnetic projection operators and as well the N -dimensional row vectors.
Because above interpolation operator in (6) and projection operator in (8) are independent of the transmitter position, we can calculate them in advance and apply them repeatedly.Especially, under the condition of L ≪ K , (7) can largely enhance the simulation efficiency.

B. THE SOLUTION OF THE PERTURBATION EQUATION AND THE EFD
When the conductivity happens to small perturbation δ σ (r) = diag( δσ H (r) δσ H (r) δσ V (r) ), we use the perturbation principle to obtain the following Helmholtz equations about the scattered potentials δA(r, r S ) and δφ(r, r S ): and Here, δJ(r, r S ) = δ σ (r)E(r, r S ) is the scattered electric currents generated by δ σ (r).It is easy to see that the left-hand terms of ( 9) and ( 10) are the same as that of ( 1) and ( 2).We try to only discretize the right-hand sides of ( 9) and ( 10) by 3-D FVM on the same Yee's grids because the discretized matrix of the left-hand term of ( 9) and ( 10) is also equal to F in (3).[20], [21].In the pixel model, we suppose the integer grid V i,j,k on Yee's staggered grid as a basic pixel element.Its conductivity σ i,j,k per pixel is also mutually independent.In the block model, the known conductivities σ Block,γ (r) is distributed in different 3-D regions V γ , (γ = 1, 2, • • • , ), and each region V γ may include multiple different basic pixels with the same conductivity.

1) EFD IN THE PIXEL MODEL
For the pixel model in Fig. 1(b), we assume that in the anomalous bodies or target domains there are . The small perturbation δ σ i α 1 ,j α 2 ,k α 3 in the conductivity of V i α 1 ,j α 2 ,k α 3 will lead to a change in the conductivity function: Here, (r, , which are respectively the relative variables in horizontal and vertical conductivities, the scattered currents in ( 9) and ( 10) can be described as follows, where, are the volumes of the corresponding staggered grids, respectively, while are the current densities at the center of half-integer grids.Substituting (12), as shown at the bottom of the next page, into (9) and integrating on V i+1/2,j,k , V i,j+1/2,k and V i,j,k+1/2 respectively give the following discretized right-hand terms of (9): Here, δ i α ,i = 1, i α = i 0, i α ̸ = i .Similarly, integral on V i,j,k yields: Because δν H,i α 1 ,j α 2 ,k α 3 and δν V,i α 1 ,j α 2 ,k α 3 in ( 13) and ( 14), as shown at the bottom of the next page, are any variables enough small, we can properly sort out the above discrete results and obtain a compact form. Here, is a 2M -dimensional column vector, and composed of relative perturbations in horizontal and vertical conductivities in M pixels, V(r S ) is a N × 2M -order matrix derived from the discrete coefficients in( 13) and ( 14).We then apply a N ×2M order of unknown matrix δx(r S ) to represent the discrete scattered potentials of the transmitter at r S due to the relative perturbation δν.Finally, we obtain the following system of ( 9) and (10): Insertion of ( 8) into ( 16) yields the following linear relationship of both δE x (r R,l , r S,k ) and δH y (r R,l , r S,k ) with δν: Here, are the 2M -dimensional row vectors and respectively represent the pixel EFD of δE x (r R,l , r S,k ) and δH y (r R,l , r S,k )

2) EFD IN THE BLOCK MODEL
For the block model in Fig. 1(c), the conductivities in block abnormal bodies γ , The relative perturbations in horizontal and vertical conductivities are denoted by ν H,Block,γ = δσ H,Block,γ /σ H,Block,γ and ν V,Block,γ = δσ V,Block,γ /σ V,Block,γ , respectively.Then the perturbations in conductivities of all anomalous bodies can also be expressed as the piece-wise constant function: Here, δ νBlock,γ = diag( ν H,Block,γ ν H,Block,γ ν V ,Block,γ ), (r, γ ) is the characteristic function of γ .In addition, it is assumed that in the anomalous body γ there are the integer grids 143288 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Then the scattered currents on the right-hand side of ( 9) can be approximately expressed as: Similarly, plugging of ( 20) into ( 9) and ( 10) and discretization by FVM give: Here, δx Block (r S ) is an unknown matrix of order N × , whose column vector corresponds to the discrete scattered potentials due to the perturbation in conductivity of γ .
T is the -order column vector composed of relative perturbation in conductivities of all block bodies, and VBlock (r S ) is the N × 2 -order matrix formed by the combination of discrete coefficients in the right-hand side of ( 9) and (10).
By using ( 7) and ( 8), we obtain the block EFD of δE x (r R,l , r S,k ) and δH y (r R,l , r S,k ) with respect to δν Block : The linearization formula of δE x (r R,l , r S,k ) and δH y (r R,l , r S,k ) with respect to δν Block can be expressed as It must be pointed out that the difference between pixel EFD and block EFD lies in that the former describes the influence of relative perturbation in conductivities of the single integer grid ) and H y (r R,l , r S,k ), while the latter does the influence of relative perturbation in conductivities of blocks including multiple integer grids on E x (r R,l , r S,k ) and H y (r R,l , r S,k ).

3) EFD OF AMPLITUDE AND PHASE
In MCSEM, the amplitude and phase of EM fields are often used as output results.For example, the horizontal electric intensity is expressed by: where A E x (r R,l , r S,k ) and θ E x (r R,l , r S,k ) are the amplitude and phase of E x (r R,l , r S,k ), respectively.Differentiating the following equation: with respect to δν yields the EFD of A Ex (r R,l , r S,k ) and θ Ex (r R,l , r S,k ) Inserting ( 17) into (26) results in the following pixel EFD S Pixel,A Ex and S Pixel,θ Ex of the amplitude and phase of E x (r R,l , r S,k ) with respect to δν in the pixel model S Pixel,θ Ex (r R,l , r S,k ) Similarly, block EFD S Block,A Ex and S Block,θ Ex of the amplitude and phase of E x (r R,l , r S,k ) with respect to δν Block in the block model can be expressed as: The EFD of the amplitude and phase of H y (r R,l , r S,k ) can be obtained by a similar method.For brevity, their corresponding contents will be ignored.

III. NUMERICAL EXAMPLES
In this section, we first give the numerical results of block EFD of the amplitude and phase of E x and H y in the block model by the previous RDA and finite difference algorithm (FDA), respectively, to validate RDA and investigate the variable characteristics and computational efficiency of the EFD.
In addition, the validity of the FVM employed in the forward simulation has been rigorously assessed through comparisons with both TLM and NMM [7].Then, we present several numerical results of pixel EFD of the amplitude and phase of E x with respect to the horizontal and vertical conductivity under the seafloor at different frequencies and offsets to further investigate the spatial distribution and variable characteristics of EFD of E x .Unless otherwise stated, the transmitting antenna is always placed 50 m above the seafloor and moves along the y = 0 main line.The electric dipole moment IL is assumed to be equal to 1 A•m.All the numerical results are calculated on the workstation with the CPU model of Intel Xeon Platinum 8269CY.A.
To validate the previous RDA about EFD, we apply the model shown in Fig. 2. The model consists of a three-layer background and a single abnormal body.The background includes air (R Air = 10 6 •m), seawater (R Seawater = 0.3 •m and 1km thick), and sediment (R H, Sedm = 1 • m with several different anisotropic coefficients λ Sedm ) from up to down.The single abnormal body is regarded as a simplified hydrocarbon reservoir.Its lengths in x, y, and z directions are 3 km, 3 km, and 0.3 km, respectively, the central position is placed at (0,0,1) km, and the resistivity is R H = 100 • m and R V = 200 • m.The computation domain in x, y, and z direction is 100km × 100km × 100km, and divided with (N x , N y , N z ) = (134, 34, 70) grid nodes.The distribution and spacing of grids are given in Table 1.Uniform grids with a spacing of 250 m in x-direction range from -10 km to 21 km, uniform grids with a spacing of 250 m in y-direction range from -3 km to 3 km, and grids with spacings of 100 m and 25 m in z-direction range from -1.1 km to 2 km.Near the seawater surface, the seafloor, and the abnormal body, the mesh spacing is reduced to 25 m in the z-direction for enhancement of discretization precision of ( 1) and ( 2).The outermost several gradient grids are selected by Lebedev grids [28].The total of discrete vector and scalar potentials is N = 1291996.The positions of the receiver and transmitter are shown in Fig. 2. The fixed receiver can measure horizontal EM components E x and H y at (-4,0,0) km.The transmitter moves from -6 km to 22 km along the x-axis with an interval of 250 m and operates at three frequencies: 0.25 Hz, 0.5 Hz, and 1 Hz.
In Fig. 3, we compare the block EFD of amplitude and phase of ∂ν H, Block , and ∂θ E x /∂ν V, Block computed by the RDA with that by FDA in the model shown in Fig. 2. Here, the results by RDA are identified with lines, and those by FDA are labeled by scatter symbols.The results by RDA and FDA are well in agreement.From the results, we can see that the EFD of E x with respect to the horizontal conductivity of the abnormal body are about three orders of magnitude less than that of the vertical conductivity.Therefore, it can be inferred that E x is much more sensitive to the change in vertical conductivity of the abnormal body.Furthermore, by comparison of FD from different frequencies, we can see that operating frequency has an obvious influence on the detection range of E x .The higher frequency, the smaller the detection range.In Fig. 4, the relative errors of the FD of amplitude and phase of E x by RDA and that by FDA are almost less than 1% except for the case of the transmitter position larger than 16 km, where the FD have become very small.Fig. 5 shows the comparison of the block EFD of amplitude and phase of H y : ∂A H y /(A H y ∂ν H, Block ), ∂A H y /(A H y ∂ν V, Block ), ∂θ H y /∂ν H, Block , and ∂θ H y /∂ν V, Block by the RDA with that by FDA with respect to horizontal and vertical conductivity of a block abnormal body in the same model.Here, the results by RDA are identified with lines, and those by FDA are labeled by scatter symbols.
The results by RDA and FDA are well in agreement.Similarly, we can observe that the EFD of H y with respect to horizontal conductivity are largely less than that of vertical conductivity.Therefore, H y is also much more sensitive to the change in vertical conductivity.The operating frequency has an obvious influence on the detection range of H y .In Fig. 6, the relative errors of the four FD of H y between RDA and FDA are almost less than 1% except for the case of the transmitter position larger than 16 km, where the FD have become very small as well.Fig. 4 and Fig. 6 demonstrate that the relative errors of FD for the vertical conductivity of the abnormal body are smaller and more stable compared to those for the horizontal conductivity, particularly in the section larger than 12 km.This is due to a difference of three orders of magnitude in value.Regarding the relative error of FD for the horizontal conductivity of the abnormal body, it can be observed from Fig. 3 (a) (b) and Fig. 5 (a) (b), specifically in the range between 12km-16km, that FD gradually diminishes to nearly zero across all three frequencies, corresponding to relative errors shown in Fig. 4 (a) (b) and Fig. 6 (a) (b), respectively.We can see that the relative error of 1.0Hz is the largest, but the absolute difference of the part with the relative error greater than 10% is very small, close to 0, thereby the reference significance of the relative error is reduced.models: (δv H, Block , δv V, Block ) = 0, (δv H, Block , δv V, Block ) = (0.005, 0), and (δv H, Block , δv V, Block ) = (0, 0.005).The results prove that the computational efficiency of the RDA is enhanced more than triples.Obviously, when the dimension of δν Block or δν Pixel increases a large value, FDA will become impractical.

B. BLOCK EFD IN SEDIMENT WITH VARIABLE ANISOTROPIC COEFFICIENTS
The above results prove that the block EFD of amplitude and phase of E x and H y are chiefly influenced by the vertical conductivity in the abnormal body.In this section, we will investigate the influence of different anisotropic coefficients of the sediment on the block EFD of E x and H y .We still apply the same mode in Fig. 2. The horizontal resistivity of the sediment is fixed at R H, Sedm = 1 • m while its anisotropic coefficients λ Sedm will be equal to three different values of 1, 2, and 3.The operating frequency is equal to 0.25 Hz.Fig. 7 shows the block EFD of E x : ∂A E x /(A E x ∂ν H, Block ), ∂A E x /(A E x ∂ν V, Block ), ∂θ E x /∂ν H, Block , and ∂θ E x /∂ν V, Block by the RDA with respect to horizontal and vertical conductivity of the block abnormal body in sediments with three different anisotropic coefficients.From the results, we can see that the block EFD of E x is much sensitive to the change in anisotropic coefficients of sediment.The larger anisotropic coefficient of sediment, the slower the decay of block EFD with the increase of distance between transmitter and receiver.That is because in the sediment with a large anisotropic coefficient, EM fields decay more slowly and travel farther.However, we also see that the size of EFD near the abnormal body clearly decreases with the increase of anisotropic coefficients of sediment because induction currents in formation become smaller with the increase of the anisotropic coefficients of sediment.

C. PIXEL SENSITIVITY IN ANISOTROPIC MODEL
Finally, we will give a spatial distribution of the pixel EFD of amplitude and phase of E x to understand the pixels at different positions on the influence of the E x .For brevity, the relevant results of the horizontal component H y will be ignored.
We still apply the model in Fig. 2. The resistivity of sediments is assumed as R H, Sedm = 1 • m with two different anisotropic coefficients of λ Sedm = 1 and 2, and operating frequencies of 0.25 Hz.The computation domain in x, y, and z directions is 100 km×100 km×100 km, and divided with (N x , N y , N z ) = (74, 74, 90) grid nodes.The distribution and spacing of grids are given in Table 3. Uniform grids with a spacing of 250 m in the x and y directions range from −8 km to 8 km, and grids with spacings of 100 m, 50 m, and 25 m in the z-direction range from −1.1 km to 3 km.Near the interface between air and seawater, seawater and the seafloor, the mesh spacing is reduced to 25 m, and around the abnormal body, the spacing is reduced to 50 m.The outermost several gradient grids are selected by Lebedev grids [23].The total of discrete vector and scalar potentials is N = 1990156.The receiver and transmitter are located at (−3,0,0) km and (3,0,0.05)km, respectively.
The sub-domain (−5, 5)km × (−5, 5)km × (0, 2)km of the computation domain is divided into M = 40 × 40 × 40 = 64000 pixels.The lengths in the x-and y-directions per pixel are 250 m, and that in the z-direction is 50 m.We then compute the pixel EFD of E x : ∂A E x /(A E x ∂ν H, Pixel ), ∂A E x /(A E x ∂ν V, Pixel ), ∂θ E x /∂ν H, Pixel , and ∂θ E x /∂ν V, Pixel with respect to horizontal and vertical conductivity per pixel with λ Sedm = 1 and 2. Fig. 9 gives the distribution of pixel FD ∂A E x /(A E x ∂ν H, Pixel ) and ∂A E x /(A E x ∂ν V, Pixel ) on the vertical cross sections of y = 0, ±2 km, ±4 km, and ±6 km.Here, Fig. 9(a) and 9(b) show the results of ∂A E x /(A E x ∂ν H, Pixel ) and ∂A E x /(A E x ∂ν V, Pixel ) with λ Sedm = 1, while Fig. 9(c) and 9(d) do that with λ Sedm = 2. From the results, we can see that FD ∂A E x /(A E x ∂ν H; Pixel ) is much less than ∂A E x /(A E x ∂ν V; Pixel ).The sizes of pixel FD on the vertical cross sections of y = 0 km are obviously larger than those on the other vertical cross sections.
Furthermore, the pixel FD with λ Sedm = 1 is larger than that with λ Sedm = 2. Fig. 10 shows the distribution of pixel FD ∂θ E x /∂ν H, Pixel and ∂θ E x /∂ν V, Pixel on the vertical cross sections of y = 0, ±2 km, ±4 km, and ±6 km.We can see that FD of the phase of E x have similar characteristics to that of the magnitude of E x .The total CPU consumption of explicit sensitivity and forward results of the model is 61 minutes, and the maximum memory occupied is 152 GB.

IV. CONCLUSION
In this paper, we apply the 3-D FVM of coupled potentials in Yee's staggered grids to establish an efficient algorithm of numerical modeling and EFD for marine controlled-source EM measurements in a 3-D TI formation.The interpolation operator and projection operator per receiver obtained by the direct solver PARDISO and 3-D Newtown interpolation can greatly enhance the efficiency of forward modeling of MCSEM with a mass of moving electric current sources.
Furthermore, the goal conductivity can be expressed as a piece-wise constant function according to the block or pixel model.The scattered electric currents caused by small perturbations in the conductivity can be decomposed into the superposition of a series of electric current elements distributed on Yee's grids.Through 3-D FVM on Yee's grid, the new discretized system about the scattered potentials can be acquired in terms of block or pixel models.Similar to the forward modeling of MCSEM, by using a projection operator per receiver we can rapidly determine the linear relationship between the changes in the EM fields and the relative perturbation in conductivities per block or pixel, and ultimately obtain the EFD of EM responses.
The numerical results show that the EFD can better evaluate the characteristics and detection ability of the MCSEM with different offsets and operating frequencies.For a single anomalous body model, the maximum offset of the effective sensitivity of horizontal EM fields can reach about 10 km at the operating frequency of 0.25 Hz.The EFD are more sensitive to the change in vertical conductivity of abnormal body and influenced by the changes in anisotropic coefficients of sediment.The larger the anisotropic coefficient of sediment, the slower the decay of block EFD with the increase of distance between transmitter and receiver.

FIGURE 1 .
FIGURE 1. Formation model and two different model spaces: (a) Inhomogeneous formation model and grids, (b) Pixel model, (c) Block model.

Fig. 1
Fig. 1 is a schematic diagram of the inhomogeneous formations model.Here, the formation is assumed to include several anomalous bodies γ with known conductivity σ Block,γ (r), (γ = 1, 2, • • • , ) in known background media σ bk (r) in Fig. 1(a).The spatial distribution of conductivity can be described by two different models: the pixel model shown in Fig. 1(b) and the block model shown in Fig. 1(c)[20],[21].In the pixel model, we suppose the integer grid V i,j,k on Yee's staggered grid as a basic pixel element.Its conductivity

FIGURE 3 .
FIGURE 3. The comparison of the block EFD of amplitude and phase of E x obtained by RDA and FDA: (a) ∂A E x /(A E x ∂ν H, Block ) (b) ∂θ E x /∂ν H, Block ,

FIGURE 4 .
FIGURE 4. The relative error of the block EFD of amplitude and phase of E x obtained by RDA and FDA: (a) ∂A E x /(A E x ∂ν H, Block ) (b) ∂θ E x /

FIGURE 5 .
FIGURE 5.The comparison of the block EFD of amplitude and phase of H y obtained by RDA and FDA: (a) ∂A H y /(A H y ∂ν H,Block ), (b) ∂θ H y / ∂ν H,Block (c) ∂A H y /(A H y ∂ν V ,Block ), (d) ∂θ H y /∂ν V ,Block .

FIGURE 6 .
FIGURE 6.The relative error of the block EFD of amplitude and phase of H y obtained by RDA and FDA: (a) ∂A H y /(A H y ∂ν H, Block ), (b) ∂θ H y / ∂ν H,Block (c) ∂A H y /(A H y ∂ν V, Block ), (d) ∂θ H y /∂ν V, Block .

TABLE 1 .
The distribution and spacing of grids.

Table 2
gives the comparison of computing time spent by RDA and FDA.From Table2, we can see the computation of F−1 by the directive solver PARDISO takes 11.03 min, which costs the longest time of all operators.When we use RDA to calculate the responses and FD of MCSEM, we need to call PARDISO only one time.However, by using FDA to calculate the responses and FD of MCSEM, we need to call PARDISO three times to determine EM responses at three different

TABLE 2 .
Comparison of computing time by RDA and FDA.
shows the block EFD of amplitude and phase of H y : ∂A H y /(A H y ∂ν H, Block ), ∂A H y /(A H y ∂ν V, Block ), ∂θ H y /∂ν H, Block , and ∂θ H y /∂ν V, Block with respect to horizontal and vertical conductivity of the abnormal body in the same FIGURE 7. The block EFD of amplitude and phase of Ex with respect to horizontal and vertical conductivity of the block abnormal body in sediments with R H,Sedm = 1 • m and different anisotropic coefficients.(a)∂AEx /(A E x ∂ν H, Block ) (b) ∂θ E x /∂ν H, Block , (c) ∂A E x /(A E x ∂ν V, Block ), (d) ∂θ E x /∂ν V, Block .The block EFD of amplitude and phase of Hy with respect to horizontal and vertical conductivity of the block abnormal body in sediments with R H,Sedm = 1 • m and different anisotropic coefficients.(a)∂AHy /(A H y ∂ν H, Block ), (b) ∂θ H y /∂ν H, Block , (c) ∂A H y /(A H y ∂ν V, Block ), (d) ∂θ H y /∂ν V, Block .modelasFig.7.From the results, we can see that the block EFD of H y have similar characteristics to that of E x .

TABLE 3 .
The distribution and spacing of grids.