CSEM inversion was performed using the open-source inversion package MARE2DEM (Modeling with Adaptively Refined Elements 2D EM) (Key, 2016). This package was developed for 2D anisotropic modeling and inversion of both offshore and onshore CSEM and magnetotelluric (MT) data. MARE2DEM is based on the Occam approach (Constable et al., 1987), which is a variant of Gauss-Newton minimization. The starting point of the inversion scheme is a nonlinear problem formulation, which is solved iteratively by minimizing a cost function (Key, 2016; Ren and Kalscheuer, 2020). ϕ m , α = d − F m † W d † W d d − F m + α m † W m † W m m , (1) where d of size i = 1,2, … N is the measured complex field data (i.e., frequency domain), F m is the corresponding model response, W d is the weighting matrix for the data misfit, α is the Langrangian weight factor for the regularization term, and W m is the regularization matrix. While dealing with complex fields, the Hermitian † (i.e., matrix transpose and complex conjugation) notation should be adopted for the matrices involved. In MARE2DEM, W d is the diagonal matrix composed of the inverse of the standard error δ for each sample and W m is the weighting matrix that forces smoothness on the model. The latter is obtained by use of a gradient roughness operator. In the case of anisotropic Earth models, the roughness is implemented by partitioning the model vector into anisotropic subsets (Key, 2016). In MARE2DEM, the model parameter m represents the logarithm of resistivity log ρ (bounded to a user-defined interval).

In practice and due to the nonlinearity of the inverse problem, the forward (modeling) operator F in Eq.1 is quasi-linearized using a Taylor series expansion. This leads to an iterative formulation where the (k+1) ^th update is given as ϕ l i n m k + 1 , α = [ d − F m k − J m k + 1 − m k † W d † W d d − F m k − J m k + 1 − m k ] + α m k + 1 † W m † W m m k + 1 , (2) where J is the model Jacobian matrix with entries ∂ F i m k / ∂ m j where m j = log ρ j and ρ j is the resistivity in cell j. Finally, after differentiating the cost function (2) with respect to the current model and setting ∂ ϕ l i n m k + 1 , α / ∂ m k + 1 = 0 , a least-squares solution is obtained after rearrangement: m k + 1 = J w − g W d d k , (3) with d k = d − F m k + J m k being the modified data vector and J w − g being the generalized inverse matrix defined as J † W d † W d J + α W m † W m − 1 J † W d † .

In MARE2DEM, Eq. 3 is solved iteratively by applying the Occam approach. This implies that the Langrangian multiplier α is optimized as part of the inversion. For more details, the reader is referred to Key (2016) and Constable et al. (1987).

In general, for an EM problem, there is a total of six different data components, which correspond to the three different directions of the magnetic and electric fields. However, in this study, we only used the embedded inline horizontal electric field (Ey), which is the most important component for marine CSEM.

If we assume a noise-free case and that the true model has been obtained from the inversion, the modified data vector can be written as d k = d t r u e = J m t r u e . (4)

The combination of Eqs 3, 4 gives then m k + 1 = R M m t r u e , (5) with R M being the resolution matrix, which can be written explicitly as (Ren and Kalscheuer, 2020) R M = R J † W d † W d J + α W m † W m − 1 J † W d † W d J , (6) and where R implies taking the real part.

In a practical inversion case m t r u e is unobtainable. The model resolution matrix reveals how close the preferred inversion model is to the true model, which relies on the Lagrangian multiplier α. By letting α → 0, the model resolution matrix approximates the unity matrix. In this case, the inverse problem is perfectly solved if no noise is present. As a pragmatic approach, we assume that m k + 1 represents the preferred inversion model if the inversion is terminated after iteration number k.

The resolution matrix is not calculated as part of the output from MARE2DEM. We have therefore developed an extension to the inversion package where this quantity is computed.

To gain further insight, we consider a 1D case first and decompose the corresponding resolution matrix into its column vectors: R M = r 1 ⋯ r j ⋯ r M , (7) where r j is the jth column vector ( j = 1,2 , ⋯ , M ) and M represents the total number of 1D image points. Each column vector in Eq. 7 represents now a point-spread function (PSF) associated with a corresponding fixed image point (cf. Figure 1A).

The concept of a PSF is well known from imaging theory (Rossmann, 1969) and describes how much a point or pixel in an image (i.e., a model parameter) is blurred due to the imaging system (a combination of acquisition system and choice of inversion algorithm in our case). A 2D image, as considered in this paper, is represented by a lexicographical ordering as illustrated in Figure 2. The resolution (blur) matrix then takes a more complex form as discussed in Section 3.2 (cf. Eq. 9). A perfectly resolved case exhibits a PSF with the value of 1 at the location of the image point and 0 elsewhere. Figures 1B, C show examples of a well-resolved and a poorly resolved case, respectively, for a 2D image. The PSF in Figure 1B is characterized by a small spread centered on the corresponding model parameter. However, in Figure 1C, the PSF is characterized by a large spread.

From now on, we will use the notation A for the resolution matrix corresponding to a lexicographic ordering of the 2D image or model. The following general relationship between the true image m and its blurred counterpart b (i.e., the output from the CSEM inversion) holds (forward model) A m = b + n , (8) where A is the blurring (resolution) matrix and n represents an additive noise term. Since m (and b ) is organized in lexicographical order (cf. Figure 2), the structure of the blur matrix takes a special form as discussed in the next section.

We start by considering the case of a space-invariant PSF. Assume that the number of image points is M = 2N + 1 along each direction (and with indices running from -N to N). Assume also that the PSF has the same size as the image (cf. upper part of Figure 3). The first step in forming the blur matrix A is to organize each column of the PSF in a Toeplitz matrix as shown in the lower part of Figure 3 for row n.

The blur matrix A can now be constructed as a block Toeplitz matrix with Toeplitz blocks (BTTB) with zero boundary conditions (Hansen et al., 2006) and where each block element a n n = − N , . . . . . − 1,0,1 . . . . . N is defined in Figure 3: A = a 0 a − 1 ⋯ a − N 0 ⋯ 0 a 1 ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ ⋮ ⋱ a 0 a − 1 ⋱ a − N 0 a N ⋱ a 1 a 0 a − 1 ⋱ a − N 0 ⋱ ⋱ a 1 a 0 ⋱ ⋮ ⋮ ⋱ a N ⋱ ⋱ ⋱ a − 1 0 ⋯ 0 a N ⋯ a 1 a 0 . (9)

The blur matrix A has dimensions M 2   x   M 2   w h e r e   M = 2 N + 1 .

To discuss a more general case characterized by space-variant PSFs, a pragmatic approach would be to subdivide the model space into space-invariant regions and perform deblurring separately for each region. This implies that each region is assigned a deblur matrix of the form given by Eq. 9, but with its own representative PSF. The final image is then constructed by combining the space-invariant regions after deblurring (and with possibly some smoothing applied to avoid edge effects). A more attractive approach, however, is to construct a space-variant A matrix (Nagy and O’Leary, 1997). Let Figure 4 represent an idealized case where the model space is horizontally split into two regions, each of them characterized by distinct and different PSFs. To minimize edge effects, a transition zone, which defines a gradual transition between the two PSFs, has also been introduced.

Before constructing the space-variant blur matrix A, we need to define a corresponding space-invariant blur matrix for each subregion (same form as in Eq. 9). In the idealized case shown in Figure 4, two blur matrices (A ₁ and A ₂ ) need therefore to be constructed. In this demonstration example, we have defined the two PSFs as simple 2D Gaussian functions with different degrees of blurring. More specifically, we chose the PSF of region 1 to introduce less blurring than the corresponding PSF of region 2. This implies that the blur matrix A ₁ has a more narrow band of values concentrated along its diagonal compared to the blur matrix A ₂ (cf. upper row in Figure 5).

The next step is to calculate a weighting matrix for each of the two regions in Figure 4 (D ₁ and D ₂ ). To avoid edge effects, the PSF should vary smoothly between different subregions. This can be achieved by applying linear tapering between neighboring subregions. In such a transition zone, an effective PSF is constructed as the linear combination between the two neighboring PSFs. The two weighting matrices for the idealized case in Figure 4 are shown in the middle row of Figure 5. A zoomed version of a section of the weight matrix D ₂ is also included to better visualize the smooth transition between the two subregions (i.e., no sharp edges). The final blur matrix A can then be constructed as the sum of the Hadamard product of the two space-invariant matrices and the associated weighting matrices A = A 1 ⊙ D 1 + A 2 ⊙ D 2 , (10a) where the effective blur matrix A is shown in the bottom row of Figure 5. A zoomed version of a section of this matrix is also shown to better illustrate the effect of the smooth transition zone introduced between the two subregions in Figure 4.

For a general case with N regions, Eq. 10a generalizes to A = ∑ i = 1 N A i ⊙ D i (10b)

As the inversion can be sensitive to input parameters (e.g., the choice of PSFs and sizes of the transition zones), developing an interactive user interface to assist in the selection of the optimized parameters is crucial. Figure 6 shows the typical plot output from this user interface (i.e., the blurred model output from the CSEM inversion along with three selected PSFs). In an interactive system, a PSF plot (Figures 6B–D) automatically updates when clicking on the corresponding cell location in the blurred model (green dots in Figure 6A). This allows users to interactively evaluate whether a given PSF is suitable, and the user can then add this PSF to a list. After selecting the optimal PSFs, the user can define boundaries and the size of the transition zones (cf. Figure 5). This set of parametric choices is then used to construct the space-variant A -matrix. This blur matrix A is stored in the memory of the user interface, allowing the user to efficiently test different sets of input parameters.

Nevertheless, the inversion is sensitive to the choice of PSFs. In the example shown in Figure 6, one of the PSFs corresponds to a poorly resolved model parameter (Figure 6B), whereas the PSFs in Figures 6C, D both describe a fairly well-resolved model parameter. These observations stress the important role our interactive user interface plays in controlling the quality of the selected parameters. The choice of an anomalous PSF (Figure 6B) can thus be easily avoided. In general, PSFs located (well) outside the target area should not be employed. Relevant PSFs are those near and inside the target area or structure. The PSFs in Figures 6C, D are examples of proper selections. It is well known that the output from a CSEM inversion is poorly resolved along the vertical direction. If the two PSFs in Figures 6C, D are used to construct the blur matrix, the corresponding deblurred image is expected to show improved vertical resolution. For a practical application, each PSF is delimited to a smaller area with tapering and normalization that ensures that the sum of its values inside the tapered area adds to 1.

In this Section, we discuss how to deblur the output from the CSEM inversion. This step represents a new inversion problem to be solved, namely, the one with a forward model given by Eq. 8. Several solution alternatives exist, and in this study, we used the nonnegative flexible conjugate gradient least-squares (NN-FCGLS) algorithm (Gazzola et al., 2017), which is implemented as an inner-outer scheme. The updated model of the inner iteration can be written as m k + 1 = m k + α k p k , (11) where α k is a bounded step size and p k is a vector of direction. To fulfill the condition of nonnegativity in the solution space, the step size α k is reduced so that m k + 1 ≥ 0   i f   m k > 0 . In NN-FCGLS, this is accomplished using a bounded step size calculated with the following expression: α ¯ k = α k   i f   p k ≥ 0 o t h e r w i s e   α ¯ k = min α k , min ⁡ − m k j / p k j j ∈ K , (12) where K is a set of indices j such that p k j < 0 and the direction p k is obtained by a linear combination of at most κ k previously computed p j with j varying in { max{0, k − κ k },. . ., k−1}. If the maximum number of iterations k _max is assigned to the inner cycle, the choice κ k = k _max corresponds to a full recursion, while a lower κ k corresponds to a truncated recursion and with κ k = 1, only the last computed vector p k − 1 is used. The outer cycle relies on suitable restarts to avoid stagnation. For further details about the algorithm, the reader is referred to Gazzola et al. (2017). In this study, we employed a code taken from the MATLAB library IR Tools (Gazzola et al., 2019).

Because the NN-FCGLS method enforces a nonnegativity constraint at each iteration, we consider that this algorithm will produce a more accurate approximate solution when the output from the CSEM inversion is a true nonnegative (i.e., log ( ρ ) and the resistivity ρ is bounded by ρ ≥ 1 Ω-m) like in this case. The proposed deblurring approach is based on PSFs extracted from the resolution matrix associated with a linearized approximation of the original nonlinear problem. Thus, this procedure does not represent an exact solution to the blurring problem and the results obtained should always be treated with caution.

No precalculated PSFs are needed for blind deconvolution. This implies that the blur matrix A is now unknown, and the blind deconvolution technique estimates both the PSFs and the unknown “true” image m (Ayers and Dainty, 1988; Krishnamurthy et al., 1995; Holmes et al., 2006). To solve this joint problem in an efficient manner, the maximum likelihood estimation (MLE) principle is employed (van Trees, 1968). The main idea is to search for a statistical solution that maximizes the likelihood function, given some observations. Thus, the parameter values obtained (i.e., PSF and m ) are most likely to lead to the observed data. For more details, the reader is referred to Krishnamurthy et al. (1995) and Biggs and Andrews (1997).

In this paper, we applied blind deconvolution to benchmark our proposed approach, using the MATLAB routine deconvblind. It restores the image and the point-spread function (PSF) simultaneously employing the accelerated, damped Richardson-Lucy algorithm in each iteration (Richardson, 1972; Lucy, 1974).

The Wisting oil field is located in the southwestern Barents Sea. The reservoir is highly segmented and very shallow (approximately 250 m below the seafloor), and contains oil in sandstone from the late Triassic (Fruholmen Formation) and early Jurassic (Nordmela and Stø Formations) periods.

In this Chapter, we start by considering an idealized model of the Wisting field that can serve as a complex controlled-data case. Such a synthetic data set is vital to perform quality control in the proposed deblurring approach. The controlled-data example is then followed by a case where CSEM field data acquired across the Wisting field is employed.

Figure 7 shows a plot of the synthetic model color-coded with vertical resistivity. Although the model includes anisotropy in most layers, we focus here on vertical resistivity because the CSEM method is known to be most sensitive to this polarization. For more details about the construction of the synthetic model, the reader is referred to Thorkildsen and Gelius (2023). Synthetic marine CSEM data were generated using the model shown in Figure 7, employing the MARE2DEM package in a forward-modeling mode. Five receivers evenly deployed on the seafloor at a 3000 m interval were used for the calculations. On the transmitter side, approximately 180 source positions were simulated using a spatial sample interval of around 200 m. Moreover, a total of 11 frequencies ranging from 0.2 to 12 Hz were included. Random noise with a standard deviation equal to 1% of the data amplitude was added to each recording.

The blurred model obtained from the CSEM inversion is displayed in Figure 8A. The inversion was characterized by a smooth transition from the low resistivity background to the high resistivity inside the reservoir. During the inversion, we constrained the model update to a rectangular area enclosing the main reservoir, while keeping the rest of the model unchanged. We also employed a vertical sample interval of 5 m, which is denser than that normally employed in CSEM inversion. This was performed to enhance the deblurring effect from a visualization point of view. When compared to the “true” model in Figure 7, the inversion was able to reproduce the main features of the reservoir. However, the main resistive structure was too shallow for both the left and (parts of) the middle compartment. The white vertical lines introduced in Figure 8A indicate the borders between the different subregions employed when constructing the space-variant blur matrix A (transition zones not shown). In addition, ideal PSFs were introduced along the red frame surrounding the target region to constrain the outer parts of the image. An ideal PSF is a point-spread function that takes the value of 1 at the corresponding parameter location and 0 elsewhere.

If we use the proposed deblurring approach, we generally expect to observe a better-resolved reservoir along the vertical direction (i.e., thinner). A direct comparison between the output from the CSEM inversion (Figure 8A) and the corresponding deblurred image obtained from NN-FCGLS (Figure 8B) supports this assumption. A total of six PSFs were employed during deconvolution (their actual locations are represented by the superimposed green dots in Figure 8). By using our interactive interface software the PSFs were carefully chosen based on visual inspection and extensive testing. After deblurring, the target structure appeared thinner overall and the right compartment showed the greatest change in resolution. This was also expected because that part of the reservoir is known to be the most poorly resolved. On the other hand, the result of blind deconvolution (Figure 8C) exhibited no resolution enhancement but rather the opposite. Note, however, that the blind deconvolution technique applied can only handle the case of a space-invariant (and unknown) PSF. The blind deconvolution result in Figure 8C was initialized employing a flat PSF. We also tested a Gaussian PSF as initialization but the results were poor. Use of multiple windows for both initializations did not improve the quality of the result.

To demonstrate the effectiveness of the proposed method on field data, the proposed deblurring approach was applied to a 2D CSEM line extracted from the BSMC08W 3D survey acquired across the Wisting oil field during the summer of 2008 (with similar acquisition parameters as in the synthetic model case). We used data corresponding to a frequency range of 0.2–12 Hz (total of 23 frequencies) as the input for the CSEM inversion.

Figure 9A shows the blurred model obtained from an unconstrained CSEM inversion. This result was achieved after cubic interpolation of the originally coarser output grid resulting from MARE2DEM. This coarser inversion grid was sampled at 200 m laterally and 30 m vertically and was chosen in collaboration with the industry to ensure a reasonable computational time. However, the interpolated grid was sampled at 50 m laterally and 5 m vertically. Resampling of the image (and also selected PSFs) was performed to ensure that the deblurring step would be stable (otherwise, very few grid points would define both the blurred image and the corresponding PSFs). Figure 9A through C show superimposed arrows indicating an estimate of the average reservoir thickness (about 40 m).

Figure 9B displays the deblurred image obtained after six iterations of the NN-FCGLS method. A total of three carefully selected PSFs were employed during deconvolution (their actual locations are represented by the superimposed green dots in Figure 9). The overall effect of deconvolution is that the reservoir has become thinner and with sharper boundaries. In addition, the right portion of the reservoir, which was more poorly resolved, appears slightly uplifted. Due to the coarse sampling of the original inversion grid and the fact that Wisting is a thin reservoir, this field data example represents a limiting case. However, the vertical resolution of the deblurred image has still improved and is now closer to the known average reservoir thickness. Thus demonstrating the potential of the proposed method. By directly comparing the deblurred or deconvolved result (Figure 9B) with the result obtained through blind deconvolution (Figure 9C), it is clear that the latter technique does a poor job of deblurring. Employing blind deconvolution actually produces an image exhibiting lower resolution as compared to the original input.