SUMMARY

When performing a cooperative inversion using a structural constraint, extracting an accurate average direction from the high-resolution model is important because two models with different resolutions should be used in the procedure. However, when the average direction of the high-resolution model, which indicates the major change direction of the model at each inversion block location, is calculated, the conventional gradient methods such as cross-gradient have a limitation on components having opposite directions. Therefore, in this study, an effective average-direction extraction algorithm was developed by introducing the concept of the structure tensor to accurately calculate average-direction information. And finally, based on the extracted average-direction information, structure-tensor-constrained cooperative inversion algorithm was proposed. To verify the effectiveness of the proposed inversion method, the method was tested using data obtained from the synthetic model containing an anomalous body with a complicated shape and the result compared with results of individual EM inversion and the conventional cross-gradient-constrained cooperative inversion. Lastly, to evaluate the performance with a realistic mode, the proposed cooperative inversion was applied to the data acquired using the complex SEG Advanced Modelling Program model. In all experiments, the cooperative inversion with the structure-tensor constraint provided better location estimation results as well as better estimations of the shape of the anomaly. In addition, the resistivity distribution of the anomaly was estimated to be closer to the truth in the inversion result.

INTRODUCTION

Seismic exploration data are generally used in oil exploration because these enable high-resolution structural information to be obtained. However, since the velocity is not directly related to oil saturation, the acquisition of accurate quantitative estimations of a target oil reservoir is limited using seismic exploration data. As an alternative, controlled-source electromagnetic (CSEM) methods have recently been introduced and studied to determine the spatial distribution of electrical resistivity in oil fields (Constable 2010). Electrical resistivity obtained from the acquisition of EM data is sensitive, among others, to oil saturations, and therefore it can be used for the quantitative estimation of the target oil reservoir. However, the resolution provided by the CSEM method is often not sufficient to image deep reservoir. Therefore, the joint interpretation using seismic data and CSEM data has been actively researched to acquire reliable inversion results (Gallardo & Meju 2003; Harris & MacGregor 2006; Hoversten et al. 2006; Gao et al. 2012; Um et al. 2014). We can classify these joint interpretation methods as the joint inversion and the cooperative inversion. The first term typically refers to simultaneous inversion of multiple parameters, while the second term refers to the inversion of one model parameter subject to constraints provided by another parameter distribution which is not inverted for.

Joint interpretation methods using the cross-gradient constraint can be easily applied in many circumstances in a robust manner. For this reason, many researchers have been studying joint interpretation methods using the cross-gradient constraint with various geophysical data (Gallardo & Meju 2003, 2004; Tryggvason & Linde 2006; Colombo & DeStefano 2007; Hu et al. 2009; Jeong et al. 2016; Colombo & Rovetta 2018). Among them, Hu et al. (2009) proposed the joint inversion method with seismic and EM data using the cross-gradient constraint. Jeong et al. (2016) presented various examples of a cooperative inversion with seismic and EM data using the cross-gradient constraint. In these studies, to acquire high-resolution structural information, a velocity model acquired using a seismic full-waveform inversion (seismic FWI) was used. Seismic FWI is well known as a promising method to acquire potentially high-resolution velocity images. However, when the method applied to real data, it is difficult to acquire velocity models that can be used as structural information in the joint or cooperative inversion, and also there exists a strong demand for modelling (Margrave et al. 2012). Furthermore, it is difficult to delineate the boundaries of the target reservoir using only the seismic velocity model, when there is a low degree of contrast between the seismic velocities of the target and non-target areas (Hu et al. 2009). Therefore, there are limitations in using the velocity model as the high-resolution structural information when preforming the joint or cooperative inversion.

On the other hand, the cross-gradient technique is known to be an effective structural constraint for the cooperative inversion of different types of geophysical data sets and is based on structural similarity. When performing a cooperative inversion using the cross-gradient constraint, it is important to match the gradient information of two inverted models with different resolutions. Therefore, the average orientation of the high-resolution model, which indicates the major change in direction of the model at the location of each inversion block, should be carefully calculated. However, when calculating the average orientation within a calculation range using the gradient method, if gradient components points opposite directions, orientations are cancelled and the average orientation is calculated as having no directionality. To overcome the limitation of the gradient method, the concept of a structure tensor has been actively researched in the image processing field (Lucas & Kanade 1981; Forstner & Gulch 1987; Weickert 1999; Zhang et al. 2009; Baghaie & Yu 2015; Sheng et al. 2016). Recently, various studies using a structure tensor also have been conducted in the geophysical field. For example, Zhou et al. (2014) proposed a structural constrained inversion of the 2-D apparent resistivity data. They extracted the structural features using a structure tensor from a migrated section of ground penetrating radar image and performed the structural constrained inversion by weighting the four-direction smoothing matrix to smooth along structural features. Scholl et al. (2017) applied the structure tensor to the 2-D depth migrated seismic section and performed the 2-D structurally guided EM inversion with directional information. Additionally, for detecting faults and channels, Wu (2017) introduced the directional structure-tensor-based coherence method.

In this paper, the limitation of gradient will be briefly explained and the process of the proposed extraction scheme using the structure tensor will be explained step by step. Based on the extracted average orientation from the high-resolution structural information, an improved cooperative inversion using the structure-tensor constraint will be finally introduced. Lastly, the effects of the structure-tensor constraint will be demonstrated on synthetic models and on the SEG Advanced Modelling (SEAM) Program model.

AVERAGE ORIENTATION EXTRACTION SCHEME USING A STRUCTURE TENSOR

There are many ways to calculate the orientation of data. The easiest and fastest way is to calculate the gradient vector of the data. This is identical to calculating the first-order derivative of the data. However, calculating the gradient vector has limitations when calculating the average orientation of complex or thin structures such as salt domes, gas chimneys or interbedded layers. For example, a model containing an anomalous body with a complicated shape was chosen. Figs 1(a) and (b) show the model containing the anomalous body and YZ cross-sections of the model, which were extracted at intervals of 80 m from y = 1660 m, respectively. The average orientation of the model containing the anomalous body was calculated using the gradient method. Many vectors cancelled out and a discontinuous and unstable extraction result was acquired (as shown in Fig. 2). This limitation is crucial when extracting a complex structure from a high-resolution structural information.

(a) Model containing an anomalous body with a complicated shape, (b) YZ cross-sections taken at intervals of 80 m from x = 1660 m.
Figure 1.

(a) Model containing an anomalous body with a complicated shape, (b) YZ cross-sections taken at intervals of 80 m from x = 1660 m.

Each component of the average orientation calculated using the gradient vector. (a) x-directional gradient, (b) y-directional gradient and (c) z-directional gradient.
Figure 2.

Each component of the average orientation calculated using the gradient vector. (a) x-directional gradient, (b) y-directional gradient and (c) z-directional gradient.

To overcome these limitations, we focused on the structure-tensor method. The structure tensor has been studied in various fields, such as data processing or computing, for analysing the orientation of an image in a stable manner. Because the structure tensor is expressed as a gradient-squared tensor, it can be considered the correct way to calculate the average orientation. The average orientation extraction scheme based on the structure tensor concept consists of the following five steps: (1) calculate the three-components gradient vector in each grid location; (2) construct a 3-D gradient-squared matrix for each grid location; (3) build a block-wise 3-D structure-tensor based on the inversion block location; (4) extract 3-D directional information from the block-wise 3-D structure tensor and (5) thin the extracted directional-based boundary. A detailed description of each step follows.

Step 1. A three-component gradient vector of a 3-D model is calculated in each grid location x, and can be expressed as follows:
(1)
where gx, gy and gz are the x-, y- and z-directional gradient components, respectively, and are identical to the first-order spatial derivatives of the 3-D model. Figs 3(b)–(d) show the x-, y- and z-directional gradient components, respectively, of a model containing an anomalous body with a complicated shape.
(a) Original model containing an anomalous body with a complicated shape in the YZ cross-section taken at intervals of 80 m from x = 1660 m and the corresponding (b) x-directional, (c) y-directional and (d) z-directional gradient components.
Figure 3.

(a) Original model containing an anomalous body with a complicated shape in the YZ cross-section taken at intervals of 80 m from x = 1660 m and the corresponding (b) x-directional, (c) y-directional and (d) z-directional gradient components.

Step 2. A 3-D gradient-squared matrix S(x) is constructed using the calculated gradient vector components (Fig. 4) and can be written as follows:
(2)
Each component of the constructed gradient-squared matrices of a 3-D model containing an anomalous body with a complicated shape in YZ cross-sections (extracted at intervals of 80 m from x = 1660 m): (a) Sxx component, (b) Syy components, (c) Szz components, (d) Sxz components, (e) Sxy components and (f) Syz components.
Figure 4.

Each component of the constructed gradient-squared matrices of a 3-D model containing an anomalous body with a complicated shape in YZ cross-sections (extracted at intervals of 80 m from x = 1660 m): (a) Sxx component, (b) Syy components, (c) Szz components, (d) Sxz components, (e) Sxy components and (f) Syz components.

Step 3. A block-wise 3-D structure tensor is constructed based on the inversion block location. To match the directional information of the high-resolution structure information with the relatively low-resolution EM inversion result, the resolution of the structure information must be effectively reduced. To do so, a Gaussian weighting kernel is applied to the structure tensor in each block coordinate, which is equal to those of the EM inversion. The block-wise 3-D structure tensor T at a block coordinate |$( {{x_0},{y_0},{z_0}} )$| can be expressed as follows:
(3)
(4)
(5)
where |${K_\sigma }({x_i},{y_i},{z_i})$| is the normalized Gaussian weighted kernel at a certain distance d from a block location |$( {{x_0},{y_0},{z_0}} )$|⁠. |${{\rm{\sigma }}_{\rm{x}}},{\rm{\ }}{{\rm{\sigma }}_{\rm{y}}}$| and |${\rm{\ }}{{\rm{\sigma }}_{\rm{z}}}$| are standard deviations of the Gaussian kernel in each direction. Based on this equation, the block-wise 3-D structure tensor can be acquired, as shown in Fig. 5.
Each component of the block-wise 3-D structure tensor in YZ cross-sections: (a) txx component, (b) tyy components, (c) tzz components, (d) txz components, (e) txy components and (f) tyz components.
Figure 5.

Each component of the block-wise 3-D structure tensor in YZ cross-sections: (a) txx component, (b) tyy components, (c) tzz components, (d) txz components, (e) txy components and (f) tyz components.

Step 4. Directional information is finally extracted from the established block-wise 3-D structure tensor. To achieve this, an eigenvalue decomposition is applied to each established structure tensor. The structure tensor T(x) is the real symmetric positive-definite matrix, and therefore the eigenvalue decomposition of the structure tensor can be expressed as follows:
(6)
where u, v and w are unit eigenvectors and |${e_u}$|⁠, |${e_v}$| and |${e_w}$| are their corresponding eigenvalues, respectively. The principal eigenvector u indicates the direction of major change in the model. Additionally, eigenvalues can be ordered as follows: |${e_u} \ge {e_v} \ge {e_w} = 0$|⁠; the magnitude of the eigenvalues indicates the degree of change in the model. In this study, the principal eigenvector u and its corresponding eigenvalue |${e_u}$| were only used as the extracted average directional information. Fig. 6(a) shows the major eigenvalues extracted from the established block-wise structure tensor in the previous step and Fig. 6(b) is a vector plot of the principal eigenvector, with a magnitude equal to the corresponding eigenvalue (⁠|${{\rm{e}}_{\rm{u}}}( {{\bf x}} ){{\bf u}}( {{\bf x}} )$|⁠).
(a) Principle eigenvalue and (b) principle eigenvector with a magnitude equal to the corresponding eigenvalue extracted from the established block-wise structure tensor shown in Fig. 5.
Figure 6.

(a) Principle eigenvalue and (b) principle eigenvector with a magnitude equal to the corresponding eigenvalue extracted from the established block-wise structure tensor shown in Fig. 5.

Step 5. Using a Gaussian kernel with surrounding values is important to increase the continuity when calculating a block-wise structure tensor. However, heavily smoothed results also make it difficult to yield precise directional information and can generate thicker boundaries than is actually the case when performing a cooperative inversion. Therefore, in this last step, the modified non-maximum suppression method is additionally applied to the eigenvalue acquired from the structure tensor. The thinning method used in this study was initially developed in the image processing field for thinning the edge detection results in 2-D images (Canny 1986). The principle of this method is to use the interpolation that locates the pixels where the norms of the gradient are at the local maximum and eliminate those pixels that are not local maxima of the magnitude in the direction of the gradient (Harbi & Syed 2014). The detailed procedure of this method in 2-D imaging is as follows. First, gradient components in the x-direction (gx) and the z-direction (gz) are calculated at each pixel. Figs 7(b) and (c) show a map of x- and z-directional gradient components, respectively. Second, the norm of the gradient (⁠|${\rm{Ng}}( {i,j} )$|⁠) in any pixel location (i, j) is calculated. Fig. 7(d) shows the norm of the gradient section calculated using gradient components shown in Figs 7(b) and (c). Finally, if the value of the norm of the gradient (C0) at the selected pixel location is higher than the minimum magnitude of the norm of the gradient, it is compared with two comparison values, C1 and C2, which are located at a unit distance from the selected pixel location to the gradient directions, as shown in Fig. 8. If the value C0 at a selected pixel location is smaller than the two comparison values (C1 and C2), the value at the pixel location is replaced by the minimum magnitude of the norm of the gradient (min (Ng)). Otherwise, the value remains unchanged. This procedure can be expressed as follows:
(7)
where Ng is the norm of the gradient section calculated in all pixel locations. C0 is the standard value at the selected pixel point and C1 and C2 are the comparison values at a normalized gradient vector location. In this process, to estimate two comparison values (C1 and C2), the bilinear interpolation method is used with four adjacent pixel values. Then, finally, we could acquire the thinned image as shown in Fig. 7(f)
An example of non-maximum suppression. (a) Original image. (b) Map of the x-directional normalized gradient. (c) Map of the z-directional normalized gradient. (d) The norm of the gradient. (e) The norm of the gradient after noise suppression. (f) Final thinned result. Source code is available at https://github.com/ajdecon/gradschool_matlab/blob/master/canny_edge.m.
Figure 7.

An example of non-maximum suppression. (a) Original image. (b) Map of the x-directional normalized gradient. (c) Map of the z-directional normalized gradient. (d) The norm of the gradient. (e) The norm of the gradient after noise suppression. (f) Final thinned result. Source code is available at https://github.com/ajdecon/gradschool_matlab/blob/master/canny_edge.m.

Schematic representation of the maximum value selection principle in the maximum suppression method (modified from Robert 2007).
Figure 8.

Schematic representation of the maximum value selection principle in the maximum suppression method (modified from Robert 2007).

To apply the method to thin the extracted 3-D principal eigenvalue volume, we changed the procedures as follows. First, instead of calculating the gradient and norm of gradients, the eigenvector and eigenvalue at each position were used. Second, to apply the thinning method to the 3-D objects, the bilinear interpolation method was extended to the trilinear interpolation method. Through these modifications, we could effectively thin the 3-D principal eigenvalue volume as shown in Fig. 9(b). If we compare the results with an average orientation acquired using gradient (Fig. 2), we can realize that the continuity of boundary are effectively increased.

Principal eigenvalue (a) before applying the thinning method and (b) after applying the thinning method.
Figure 9.

Principal eigenvalue (a) before applying the thinning method and (b) after applying the thinning method.

COOPERATIVE INVERSION USING A STRUCTURE-TENSOR CONSTRAINT

The basic goal of an inversion is to minimize the objective function (⁠|${\phi _{\rm{d}}}( {{\bf m}} )$|⁠), which consists of the data misfit between observed data (⁠|${{{\bf d}}^{{\rm{obs}}}}$|⁠) and predicted data (⁠|${\rm{F}}( {{\bf m}} )$|⁠), by forward modelling based on the estimated model parameter (⁠|${{\bf m}}$|⁠). To estimate an ideal model parameter in the iterative nonlinear inverse problem based on the least-squares method, the updated model parameter that minimizes the objective function must be calculated. However, ill-posed problems occur frequently in many practical applications when faced with limited and insufficient data. To enhance the resolving power and stabilize the ill-posed problem, the objective function can be normally modified by adding model constraints as follows:
(8)
where |${\phi _{\rm{d}}}( {{\bf m}} ){\rm{\ }}$| and |${\phi _{\rm{m}}}( {{\bf m}} )$| are the data misfit terms and model structure terms, respectively, and |${{\bf \Lambda }}$| is the trade-off parameter that controls the effect of each term in the objective function. Therefore, an appropriate |${{\bf \Lambda }}$| value must be carefully selected by compromising the resolution and stability (Yi et al. 2003). The data misfit terms can be expressed as follows:
(9)
where Wd is a diagonal matrix, in which each component is determined depending on the quality of each data point, and the operator |${\| \cdot \|^2}$| denotes the usual Euclidean norm. In general, Occam's method is used as the model structure term (⁠|${\phi _{\rm{m}}}( {{\bf m}} )$|⁠) to obtain a smooth solution and it can be written as follows:
(10)
where |${{{\bf W}}_{\rm{m}}}$| is a diagonal weighting matrix that results in a similar resolution in the inversion pertaining to different model parameters of the model and C is the roughness matrix consisting of the second-order finite-difference operator. Because we can acquire a smooth solution using Occam's method, the numerical stability of the inversion can be improved (Kim et al. 1999). |${{{\bf m}}^{( {{\boldsymbol{k}} - 1} )}}$| is the model parameter vector at the |${( {k - 1} )^{{\rm{th}}}}$| iteration step and δm is the updated model parameter to be estimated.
The basic 3-D EM inversion algorithm used in this study was originally developed by Chung (2011) and the equation of the updated model parameter vector for a regularized Gauss–Newton inversion at the kth iteration can be expressed as follows:
(11)
where J is the Jacobian matrix, which is defined as the derivative of the output error vector with respect to the model parameters. The lower and upper bounding constraints on model parameters, using the unified transformation function, are used to exclude unrealistic values in the predicted model parameters (Kim & Kim 2011). Therefore, the Jacobian matrix can be expressed as follows (Noh et al. 2014):
(12)
where diag |$( \cdot )$| means the diagonal matrix.

Additionally, in eq. (11), Rd and Rm are general measurements made using an iteratively reweighted least-squared (IRLS) algorithm. Therefore, depending on the data and target model distribution, these measurements can not only be the l2-measure of data misfit and model structure, but also their non-l2 measures, such as l1 or Huber (Farquharson & Oldenburg 1998).

However, if the basic EM inversion algorithm is used, a smoothed inversion result is usually obtained. Therefore, for an accurate interpretation based on EM inversion results, an additional method is required to improve the resolution of EM inversion results. For this, Jeong et al. (2016) developed a seismic-constrained EM inversion algorithm using the cross-gradient constraint. This method uses the structural similarity between the P-wave velocity and resistivity model. The objective function of the method is as follows:
(13)
where |${{{\bf m}}_\rho }$| and |${{{\bf m}}_v}$| are the resistivity and velocity model parameter vectors, respectively, |${\phi _{{\rm{ref}}}}( {{{\bf m}}_\rho ^{( k )}} )$| and |${\phi _{{\rm{cg}}}}( {{{\bf m}}_\rho ^{( k )},{{{\bf m}}_v}} )$| are the reference model constraint term and the cross-gradient constraint term, respectively, and β and γ are the trade-off parameters and play a role similar to |${{\bf \Lambda }}$|⁠, which balances each term in the objective function. The general model reference term |$( {{\phi _{{\rm{ref}}}}( {{{\bf m}}_\rho ^{( k )}} )} )$| can be expressed as follows:
(14)
where P is a diagonal matrix with values of 0 or 1. If we have prior information that can be used as a reference model (⁠|${{{\bf m}}_{{\rm{ref}}}}$|⁠) in a specific region of the inversion block, the information can be applied by assigning a value of one to the corresponding diagonal component of P. To calculate the cross-gradient constraint term |$( {{\phi _{{\rm{cg}}}}( {{{\bf m}}_\rho ^{( k )},{{{\bf m}}_v}} )} )$|⁠, Jeong et al. (2016) used the estimated resistivity model-parameter vector |$( {{{\bf m}}_\rho ^{( k )}} )$| at the |${k^{{\rm{th}}}}$| iteration step and the final velocity model parameter vector |$( {{{{\bf m}}_v}} )$| obtained from a seismic FWI. The cross-gradient constraint term can be defined as follows:
(15)
where ∇ represents the gradient and × is the cross-product. This term means that the constraint forces two different models to be a similar structure. When the velocity model and resistivity model are structurally related, the seismic velocity model, which has a relatively high resolution, can help to improve the resolution of the resistivity model.
Considering the cross-gradient and model reference terms, the equation of the updated model parameter vector in a seismic-guided EM inversion at the |${k^{{\rm{th}}}}$| iteration can be written as
(16)
where GCG indicates the derivative operator matrix calculated from the seismic velocity obtained by an FWI and the GCG matrix for the 2-D velocity section can be written as follows:
(17)
Jeong et al. (2016) developed a seismic-constrained EM inversion algorithm with a cross-gradient constraint based on the 2.5-D EM inversion algorithm and 2-D seismic FWI algorithm. Therefore, we developed the 3-D version of the inversion algorithm, with the cross-gradient constraint based on the basic EM inversion algorithm introduced above. The G matrix for the 3-D velocity model can be written as follows:
(18)
With the rapid development of seismic interpretation technology, more detailed and complex information can be used as a structural constraint and the importance of an accurate average orientation extraction scheme has increased. Therefore, for the effective integration of the EM inversion method with the structural information, a cooperative inversion algorithm, with a structure tensor, using seismic interpretation results was developed. A flowchart of the proposed algorithm is given in Fig. 10, and the objective function can be written as follows:
(19)
where |${\phi _{{\rm{ST}}}}( {{{\bf m}}_\rho ^{( k )},{{{\bf E}}_{\rm{U}}},{{{\bf U}}_{{\rm{ST}}}}} )$| is the structure-tensor constraint term, |${\rm{\ }}{{{\bf E}}_{\rm{U}}}$| is the principal eigenvalue acquired from the structure tensor and its corresponding eigenvector |$( {{{{\bf U}}_{{\rm{ST}}}}} )\ $| can be expressed as follows:
(20)
Flowchart of the proposed cooperative inversion method using a structure tensor.
Figure 10.

Flowchart of the proposed cooperative inversion method using a structure tensor.

The structure-tensor constraint term can be expressed in a similar way to eq. (15) as follows:
(21)
Considering the structure-tensor constraint term, the equation of the updated model parameter vector in a seismic-guided EM inversion at the |${k^{{\rm{th}}}}$| iteration can be written as
(22)
where GST is a derivative operator matrix calculated using principal eigenvectors and eigenvalues obtained from the structure tensor and it can be written as follows:
(23)

The developed algorithms will be tested on two synthetic models. The first model is the case of a model containing an anomalous body and the second model is a realistic case from the SEAM Project.

SYNTHETIC TEST

To verify the effectiveness of the proposed cooperative inversion algorithm with the structure tensor, we compared three different inversion results with the data acquired from the model containing an anomalous body with a complicated shape (see Fig. 1). In the case of CSEM data used in this study, we suppose that the horizontal electric dipole (HED) source towed about 60 m above the seafloor with frequencies of 0.1, 0.5, 1 and 2 Hz and 15 receiver nodes placed on the seafloor. To confirm only the effectiveness of the proposed inversion method, we first supposed that the structural information is perfectly extracted from seismic data. The model, assuming that it was acquired through interpretation results of seismic data, consisted of 61, 81 and 61 grid numbers in the x-, y- and z-directions, respectively, while the EM inversion results consisted of 15, 20 and 16 block numbers in the x-, y- and z-directions, respectively. To match the directional information with EM inversion results, the average directional information was extracted from the model with structurally high resolution.

During the inversion process, we set the 1 ohm-m half-space initial model and a reference model was not used. Selecting the proper trade-off parameter |${{\bf \Lambda }}$| is important when performing the inversion process. However, it is difficult to find and formulate rules. In this research, the optimum parameter was found through a line search method and for reducing the effect of smoothing, the parameter was set to exponentially decrease at each iteration. Additionally, when selecting parameter γ expressed in eqs (13) and (19), it was set so that the structure-tensor constraint term |$( {{\phi _{{\rm{ST}}}}( {{{\bf m}}_\rho ^{( k )},{{{\bf E}}_{\rm{U}}},{{{\bf U}}_{{\rm{ST}}}}} )} )$| is 10 to 100 times smaller than the model structure constraint term |${\rm{\ }}( {{\phi _{\rm{m}}}( {{\bf m}} )} )$|⁠.

Fig. 11 shows the inversion results in the YZ cross-section. Fig. 11(a) shows cross-sections of the true resistive model. Fig. 11(b) shows inverted cross-sections obtained from the individual EM inversion. Fig. 11(c) displays inverted cross-sections obtained from the cooperative inversion with a cross-gradient constraint. Finally, Fig. 11(d) shows inverted cross-sections obtained from the proposed cooperative inversion with the structure-tensor constraint. Each inversion was continued until normalized data misfit is sufficiently small or until the difference between successive estimates is small enough. Fig. 12 shows the comparison of normalized RMS error curves between inversion method with cross-gradient constraint and structure-tensor constraint. Both inversion methods yielded stable convergence results. The results obtained from the cooperative inversion with the structure-tensor constraint had the highest resolution among the three different EM inversion results. The inverted model shown in Fig. 11(d) provided better location estimation results as well as better estimations of the shape of the anomaly. In addition, the resistivity distribution of the anomaly was estimated to be closer to the truth in the inversion result acquired using the cooperative inversion with the structure-tensor constraint.

(a) YZ cross-section of the true resistive model. (b) Individual electromagnetic (EM) inversion result. (c) Cooperative inversion result with the cross-gradient constraint. (d) Cooperative inversion result with the structure-tensor constraint. The white dotted line indicates the structure information used in the cooperative inversion.
Figure 11.

(a) YZ cross-section of the true resistive model. (b) Individual electromagnetic (EM) inversion result. (c) Cooperative inversion result with the cross-gradient constraint. (d) Cooperative inversion result with the structure-tensor constraint. The white dotted line indicates the structure information used in the cooperative inversion.

Comparison of normalized RMS error curves.
Figure 12.

Comparison of normalized RMS error curves.

For a more detailed comparison, the model misfit errors were calculated using the mean absolute percentage error equation (MAPE):
(24)
where |${m_i}$| and |$m_i^{\mathrm{ true}}$| are the inverted and true model parameters at the ith inversion block, respectively, and n is the number of inversion blocks. The average model misfit calculated using the individual EM inversion result is shown in Table 1. The inversion results acquired from the proposed cooperative inversion had the smallest error with the true model. The estimation accuracy when using the proposed cooperative inversion was improved by 2.77 per cent compared to the accuracy when using the cooperative inversion with cross-gradient constraints.
Table 1.

Comparison of model misfits between true and inverted results.

Individual EM inversion resultCooperative inversion result + Cross-gradient constraintCooperative inversion result + Structure-tensor-based constraint
MAPE*22.67 per cent21.62 per cent18.85 per cent
Individual EM inversion resultCooperative inversion result + Cross-gradient constraintCooperative inversion result + Structure-tensor-based constraint
MAPE*22.67 per cent21.62 per cent18.85 per cent

*MAPE: Mean absolute percentage error.

Table 1.

Comparison of model misfits between true and inverted results.

Individual EM inversion resultCooperative inversion result + Cross-gradient constraintCooperative inversion result + Structure-tensor-based constraint
MAPE*22.67 per cent21.62 per cent18.85 per cent
Individual EM inversion resultCooperative inversion result + Cross-gradient constraintCooperative inversion result + Structure-tensor-based constraint
MAPE*22.67 per cent21.62 per cent18.85 per cent

*MAPE: Mean absolute percentage error.

To investigate how sensitive the proposed method is to noise, we added 5 per cent independent random noise to the synthetic data relative to the amplitudes of the CSEM data. Fig. 13(a) shows the true resistive model. Figs 13(b) and (c) display the cooperative inversion result of noise-free data and noise-added data, respectively. Compared with the case without noise, the anomaly region is imaged similarly, but the boundary of anomaly is smoothed due to the noise. Additionally, it is confirmed that the background noise in the inversion result increases due to the noise in the data. Therefore, when using the method, if a large amount of noise is presented in the data, it is difficult to obtain an accurate quantitative interpretation result.

(a) YZ cross-section of the true resistive model. (b) Cooperative inversion result of the noise-free data with the structure-tensor constraint. (c) Cooperative inversion result of noisy data with the structure-tensor constraint (5 per cent random noise added to the original data).
Figure 13.

(a) YZ cross-section of the true resistive model. (b) Cooperative inversion result of the noise-free data with the structure-tensor constraint. (c) Cooperative inversion result of noisy data with the structure-tensor constraint (5 per cent random noise added to the original data).

To consider the result in presence of errors in the model used as a reference, the model that is moved in the y-direction was used as a reference model. Fig. 14 shows the cooperative inversion result using the moved structural information. The black dotted region shows the moved structural information and the white dotted region shows the original structural information of the true model. Even if we use the erroneous structural information, we can see that the location of high-resistivity distribution is located in original dotted region. However, when compared with the inversion result using the correct structural information, high-resistivity distribution is only located in the regions where two structures overlap. Therefore, for more precise quantitative interpretation, it is important to obtain and use accurate structural information using this inversion method.

Cooperative inversion result using shifted structural information. The white dotted region indicates the structure of true model and the black dotted region indicates the shifted structure model used as structural constraint in this test.
Figure 14.

Cooperative inversion result using shifted structural information. The white dotted region indicates the structure of true model and the black dotted region indicates the shifted structure model used as structural constraint in this test.

APPLICATION TO A REALISTIC MODEL

To evaluate the performance with a realistic model, the proposed cooperative inversion method was applied to the data acquired using the SEAM Subsalt Earth Model. The SEAM model used in this study was developed by the SEAM corporation to provide the geoscience exploration community with a useful model for important geophysical challenges that provide a high business value to the petroleum resource industry (Fehler & Keliher 2011). The SEAM models have shown the 3-D representation of a deepwater Gulf of Mexico salt domain, including an oil and gas reservoir. There are not only various models of geophysical properties, such as density, P- and S-wave velocities, and electrical resistivity, but also various models of reservoir properties, such as shale volume and porosity. Fig. 15 shows a 3-D oblique view of the top of the salt body in the SEAM model. The dimensions of the original model are 40 × 35 km laterally and 15 km in depth and it was defined on a 20 × 20 × 10 m grid. In this study, considering the memory and processing capability limitations, the small region marked by the yellow dashed window in Fig. 15 was used as the model.

3-D oblique view of the SEG Advanced Modelling Program (SEAM) model. The region marked by the yellow-dashed window was used in this study.
Figure 15.

3-D oblique view of the SEG Advanced Modelling Program (SEAM) model. The region marked by the yellow-dashed window was used in this study.

In this study, a new sequential PI (Poisson impedance) method using TCCA (Target correlation coefficient analysis) which is proposed by Kim et al. (2016) was used to effectively extract the target anomaly distribution. PI, originally proposed by Quakenbush et al. (2006), can be defined as follows
(25)
where AI and SI are the acoustic and shear impedance, respectively, and c is the scaling or rotation factor.

The workflow of the QSI method used in this study is shown in Fig. 16. First, using proper well-log data, we determined the optimal c values and the reference values for determining target properties in every three steps. The three steps consisted of lithology impedance (LI), porosity impedance (ϕI) and fluid impedance (FI). For calculating the LI equation, correlation coefficient analysis was performed using the gamma ray (GR) log and PI values having different c values. The maximum absolute correlation coefficient was 0.9514 when the c value was 1.42 (Fig. 17a). Figs 17(b) and (c) show a crossplot and a joint histogram between the selected LI value and shale volume at well locations. The largest LI value when the shale volume was 0.5 v/v was selected as the reference LI value. The sand and shale areas were determined using the reference LI value, and the shale area was then excluded from the next ϕI analysis. Then, for eliminating low-porosity areas, the optimal ΦI volume was estimated based on the correlation coefficient analysis between effective porosity logs and PI curves with different c values (Fig. 17d). Like in the previous step, considering the crossplot and joint histogram, the largest ΦI value when the porosity was 0.25 was selected as the reference value for discriminating the porous area. Lastly, FI was determined with the resistivity log only for the porous sand area. Referring to the results of the interpretation of porosity and fluid constituents in the Gulf of Mexico reservoir (Ijasan et al. 2013), the log10(resistivity) value that is larger than 0.1416 ohm-m (corresponding to 60 per cent oil saturation) was designated as an oil-bearing reservoir. Finally, after determining c values and reference values using well-log data, this process was repeated using the AI and SI volume, which is calculated using a simultaneous inversion with seismic data. We could then identify the distribution of promising oil-bearing reservoirs in the seismic survey area (Fig. 18).

The workflow for a sequential Poisson impedance (PI) analysis (modified from Kim et al. 2016).
Figure 16.

The workflow for a sequential Poisson impedance (PI) analysis (modified from Kim et al. 2016).

(a) Correlation coefficient between gamma ray (GR) log data and Poisson impedance (PI) values with different c values. (b) Crossplot between the estimated lithology impedance (LI) and shale volume. (c) Joint histogram of the estimated LI and shale volume. (d) Correlation coefficient between effective porosity-log data and PI values with different c values. (e) Crossplot between the estimated porosity impedance ($\phi {\rm{I}}$) and porosity. (f) Joint histogram of the estimated $\phi {\rm{I}}$ and porosity. (g) Correlation coefficient between resistivity-log data and PI values with different c values. (h) Crossplot between the estimated fluid impedance (FI) and resistivity. (i) Joint histogram of the estimated FI and resistivity.
Figure 17.

(a) Correlation coefficient between gamma ray (GR) log data and Poisson impedance (PI) values with different c values. (b) Crossplot between the estimated lithology impedance (LI) and shale volume. (c) Joint histogram of the estimated LI and shale volume. (d) Correlation coefficient between effective porosity-log data and PI values with different c values. (e) Crossplot between the estimated porosity impedance (⁠|$\phi {\rm{I}}$|⁠) and porosity. (f) Joint histogram of the estimated |$\phi {\rm{I}}$| and porosity. (g) Correlation coefficient between resistivity-log data and PI values with different c values. (h) Crossplot between the estimated fluid impedance (FI) and resistivity. (i) Joint histogram of the estimated FI and resistivity.

Sequential Poisson impedance (PI) result. (a) XZ cross-section of the original acoustic impedance model. (b) The acoustic impedance section remaining after performing a lithology impedance (LI). The shale-dominated area was eliminated based on the LI result. (c) The acoustic impedance remaining section after performing a porosity impedance ($\phi {\rm{I}}$). The tight-porous area was eliminated based on the PI result. (d) The acoustic impedance section remaining after performing a fluid impedance (FI). Only the promising target reservoir area remained in the acoustic impedance model. (e) The final target distribution identified using a sequential Poisson impedance (PI) using a target correlation coefficient analysis.
Figure 18.

Sequential Poisson impedance (PI) result. (a) XZ cross-section of the original acoustic impedance model. (b) The acoustic impedance section remaining after performing a lithology impedance (LI). The shale-dominated area was eliminated based on the LI result. (c) The acoustic impedance remaining section after performing a porosity impedance (⁠|$\phi {\rm{I}}$|⁠). The tight-porous area was eliminated based on the PI result. (d) The acoustic impedance section remaining after performing a fluid impedance (FI). Only the promising target reservoir area remained in the acoustic impedance model. (e) The final target distribution identified using a sequential Poisson impedance (PI) using a target correlation coefficient analysis.

Next, to extract the average directional information of the target oil reservoir distribution identified in the previous section, the 3-D extraction scheme developed in this study was used. The extracted oil reservoir distribution model consisted of 96, 95 and 85 grid numbers in the x-, y- and z-directions, respectively. However, the inversion results consisted of 23, 23 and 24 block numbers in the x-, y- and z-directions, respectively. To match the directional information from the final property distribution model with one from the EM inversion results, the average directional inversion should be calculated in each inversion block location. As a result, the average directional information was appropriately calculated using the proposed extraction scheme at the inversion block location as shown in Fig. 19.

Vector plot of the principle eigenvector calculated from the structure tensor in each inversion block location.
Figure 19.

Vector plot of the principle eigenvector calculated from the structure tensor in each inversion block location.

In the case of CSEM data used in this study, we supposed that the HED source towed about 60 m above the seafloor and 20 receiver nodes placed on the seafloor. In general, noise level of mCSEM data is about 10−15 (V Am−2) for an electric field and 10−12 (A m−1Am−1) for a magnetic field (Constable 2010). Therefore, the data above the noise level were selected and used for inversion. In addition, referring the feasibility study of mCSEM data (Kang 2012), the frequencies lower than 0.5 Hz was selected (0.05, 0.1, 0.2 and 0.4 Hz).

To show the effect of the cooperative inversion method developed in this study, the results of an individual EM inversion and a cooperative inversion were compared. The initial resistivity model was set based on the average value of resistive logs in oil-free regions. Figs 20(a) and (b) show the true and initial resistivity models, respectively. Fig. 20(c) displays the individual EM inversion result. The individual inversion result was very different from the true model. Using the individual inversion result, it was not possible to successfully estimate the quantity of the oil reservoir.

(a) True resistivity model. (b) Linearly increasing initial model. (c) Individual electromagnetic (EM) inversion result. (d) Upgraded initial model using extracted target distribution information. (e) Proposed inversion result with the upgraded initial model shown in (d).
Figure 20.

(a) True resistivity model. (b) Linearly increasing initial model. (c) Individual electromagnetic (EM) inversion result. (d) Upgraded initial model using extracted target distribution information. (e) Proposed inversion result with the upgraded initial model shown in (d).

In case of the SEAM model used in this study, the depth of the reservoir is so deep that it is less sensitive to inversion and it is strongly influenced by the initial model when performing inversion. Therefore, to enhance the inversion results, the initial resistivity model was modified using the oil reservoir distribution identified during the PI analysis. The resistivity value at the location of the oil reservoir was replaced with the resistivity value higher than background resistivity and lower than the oil reservoir and was then smoothed. Fig. 20(d) shows the improved initial model based on the physical property distribution. Fig. 20(e) shows the cooperative inversion results. Compared with the previous inversion results shown in Fig. 20(c), the inversion results were notably improved. The cooperative inversion result provided a better location estimation as well as better estimation of the resistivity. The average model misfit calculated using the inversion results are shown in Table 2. Based on these average model misfits, a 28.13 per cent improvement in the inversion result was achieved using the proposed cooperative inversion method.

Table 2.

Comparison of model misfits between true and inverted results.

Individual EM inversion result (shown in Fig. 20c)Cooperative inversion result + modified initial model using property distribution (shown in Fig. 20e)
MAPE*68.12 per cent39.99 per cent
Individual EM inversion result (shown in Fig. 20c)Cooperative inversion result + modified initial model using property distribution (shown in Fig. 20e)
MAPE*68.12 per cent39.99 per cent

*MAPE: Mean absolute percentage error.

Table 2.

Comparison of model misfits between true and inverted results.

Individual EM inversion result (shown in Fig. 20c)Cooperative inversion result + modified initial model using property distribution (shown in Fig. 20e)
MAPE*68.12 per cent39.99 per cent
Individual EM inversion result (shown in Fig. 20c)Cooperative inversion result + modified initial model using property distribution (shown in Fig. 20e)
MAPE*68.12 per cent39.99 per cent

*MAPE: Mean absolute percentage error.

Additionally, we compared the effect of cross-gradient constraints and structure-tensor constraints on the thin-layer structure. Figs 21(a) and (b) show the diagonal components of GCG and GST explained in eqs (18) and (23), respectively. From these results, it was clear that when the structure-tensor constraint was used, the upper and lower boundaries of the thin layer were calculated effectively. On the other hand, when the cross-gradient constraint was used, the boundary information was calculated incorrectly or disappeared. Therefore, if we use the cross-gradient constraint on thin-layer models, we cannot obtain reliable inversion results (Fig. 21c).

Comparison between the cross-gradient constraint and structure-tensor constraint. (a) Diagonal components of the derivative operator matrix GCG explained in eq. (18). (b) Diagonal components of the derivative operator matrix GST explained in eq. (23). (c) Cooperative inversion result using the cross-gradient constraint. (d) Proposed cooperative inversion result using the structure-tensor constraint.
Figure 21.

Comparison between the cross-gradient constraint and structure-tensor constraint. (a) Diagonal components of the derivative operator matrix GCG explained in eq. (18). (b) Diagonal components of the derivative operator matrix GST explained in eq. (23). (c) Cooperative inversion result using the cross-gradient constraint. (d) Proposed cooperative inversion result using the structure-tensor constraint.

CONCLUSIONS

When performing a cooperative inversion with two different sets of geophysical data using the structural constraint, it is crucial to accurately extract a high-resolution model of the target reservoir and calculate the appropriate average orientation from the extracted high-resolution model. Therefore, in this study, to efficiently link the structural information between two models with different resolutions, the structure-tensor approach was applied instead of the conventional gradient approach for effectively extracting the average orientation from high-resolution information. At this scheme, we added the thinning process to compensate for the effect of the Gaussian weight, which increases continuity but blurs the results. Finally, using the extracted average orientation, the previous cooperative inversion algorithm was improved. To analyse the effectiveness of the proposed cooperative inversion algorithm, three different inversion algorithms (individual EM inversion, cooperative inversion with a cross-gradient and cooperative inversion with a structure tensor) were compared with the data acquired from the model, including an anomalous body with a complicated shape. The result of the proposed cooperative inversion showed the highest resolution close to actual resistivity values. Additionally, to evaluate the performance of the cooperative inversion algorithm developed in this study, the proposed method was applied to the challenging SEAM model. It should be noted here that, because there is already information regarding the target reservoir distribution, it can be helpful to set up an initial model that has a substantial impact on the inversion results. Through this process, the relative model misfit of the proposed cooperative inversion was improved by 28.13 per cent compared with the individual EM inversion result.

ACKNOWLEDGEMENTS

This work was supported by the Human Resources Development Program of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grants funded by the Korea government Ministry of Trade, Industry and Energy (No. 20164010201120). This research was also supported by the basic research project of the Korea Institute of Geoscience and Mineral Resources, Daejeon, South Korea, funded by the Ministry of Science and ICT of Korea (GP2017-022).

REFERENCES

Baghaie
A.
,
Yu
Z.
,
2015
.
Structure tensor based image interpolation method
,
AEU-Int. J. Electron. Commun.
,
69
,
515
522
.

Canny
J.
,
1986
.
A computational approach to edge detection
,
IEEE Trans. Pattern Anal. Mach. Intell.
,
8
,
679
698
.

Chung
Y.
,
2011
.
Three-Dimensional Imaging of Subsurface Structure Using Controlled-source Electromagnetic Surveys
,
PhD thesis
,
Seoul National University
,
Seoul, South Korea
.

Colombo
D.
,
Rovetta
D.
,
2018
.
Coupling strategies in multiparameter geophysical joint inversion
,
Geophys. J. Int.
,
215
,
1171
1184
.

Colombo
D.
,
Stefano
M.D.
,
2007
.
Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: application to prestack depth imaging
,
Leading Edge
,
26
,
326
331
.

Constable
S.
,
2010
.
Ten years of marine CSEM for hydrocarbon exploration
,
Geophysics
,
75
,
75A67
75A81
.

Farquharson
C.G.
,
Oldenburg
D.W.
,
1998
.
Non-linear inversion using general measures of data misfit and model structure
,
Geophys. J. Int.
,
134
,
213
227
.

Fehler
M.
,
Keliher
P.J.
,
2011
.

SEAM Phase 1: Challenges of Subsalt Imaging in Tertiary Basins, with Emphasis on Deepwater Gulf of Mexico
,
Soc. Explor. Geophys
.

Forstner
W.
,
Gulch
E.
,
1987
.
A fast operator for detection and precise location of distinct points, corners and centres of circular features
, in
ISPRS Inter-Commission Conference on Fast Processing of Photogrammetric Data
, pp.
281
305
.,
Sci. Res
,
Interlaken
.

Gallardo
L.A.
,
Meju
M.A.
,
2003
.
Characterization of heterogeneous near-surface materials by joint 2D inversion of dc resistivity and seismic data
,
Geophys. Res. Lett.
,
30
,

Gallardo
L.A.
,
Meju
M.A.
,
2004
.
Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradient constraints
,
J. geophys. Res.
,
109
,
B03311
,
doi:10.1029/2003JB002716
.

Gao
G.
,
Abubakar
A.
,
Habashy
T.
,
2012
.
Joint petrophysical inversion of electromagnetic and full-waveform seismic data
,
Geophysics
,
77
,
WA3
WA18
.

Harbi
N.S.A
,
Syed
M.P.
,
2014
.
A novel approach to gradient edge detection
,
J. Glob. Res. Comput. Sci.
,
5
,
1
4
.

Harris
P.
,
MacGregor
L.
,
2006
.
Determination of reservoir properties from the integration of CSEM, seismic, and well-log data
,
First Break
,
25
,
53
59
.

Hoversten
G.M.
,
Cassassuce
F.
,
Gasperikova
E.
,
Newman
G.A.
,
Chen
J.
,
Rubin
Y.
,
Hou
Z.
,
Vasco
D.
,
2006
.
Direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data
,
Geophysics
,
71
,
C1
C13
.

Hu
W.
,
Abubakar
A.
,
Habashy
T.M.
,
2009
.
Joint electromagnetic and seismic inversion using structural constraints
,
Geophysics
,
74
,
R99
R109
.

Ijasan
O.
,
Torres-Verdin
C.
,
Preeg
W.E.
,
2013
.
Interpretation of porosity and fluid constituents from well logs using an interactive neutron-density matrix scale
,
Interpretation
,
1
,
T143
T155
.

Jeong
S.
,
Byun
J.
,
Seol
S.J.
,
Yi
M.J.
,
Park
G.
,
2016
.
Integrated inversion of seismic and MCSEM data for the estimation of petrophysical parameters
, in
SEG Technical Program Expanded Abstracts 2016
, pp.
912
916
.,
Soc. Explor. Geophys
,
Tulsa
.

Kang
S.
,
Seol
S.J.
,
Byun
J.
2012
.
A feasibility study of CO2 sequestration monitoring using the mCSEM method at a deep brine aquifer in a shallow sea
,
Geophysics
,
77
(2):
E117
E126
.

Kim
H.J.
,
Kim
Y.H.
,
2011
.
A unified transformation function for lower and upper bounding constraints on model parameters in electrical and electromagnetic inversion
,
J. Geophys. Eng.
,
8
,
21
26
.

Kim
H.J.
,
Song
Y.
,
Lee
K.H.
,
1999
.
Inequality constraint in least-squares inversion of geophysical data
,
Earth Planets Space
,
51
,
255
259
.

Kim
S.
,
Lee
J.
,
Kim
B.
,
Byun
J.
,
2016
.
Effective workflow of Poisson impedance analysis for identifying oil reservoir with similar resistivity log response to neighboring medias
, in
SEG Technical Program Expanded Abstracts 2016
, pp.
2851
2855
.,
Soc. Explor. Geophys
,
Tulsa
.

Lucas
B.
,
Kanade
T.
,
1981
.
An iterative image registration technique with an application to stereo vision
, in
IJCAI'81 Proceedings of the 7th international joint conference on Artificial intelligence - Volume 2
, pp.
674
679
.,
Morgan Kaufmann Publishers Inc
,
San Francisco
.

Margrave
G.F.
,
Innanen
K.
,
Yedlin
M.
,
2012
.
A perspective on Full-Waveform Inversion
,
CREWES Res. Rep.
,
24
,
1
18
.

Noh
K.
,
Chung
Y.
,
Seol
S.J.
,
Byun
J.
,
Uchida
T.
,
2014
.
Three-dimensional inversion of CSEM data: water leak detection using a small-loop EM method
,
J. appl. Geophys.
,
102
,
134
144
.

Quakenbush
M.
,
Shang
B.
,
Tuttle
C.
,
2006
.
Poisson impedance
,
Leading Edge
,
25
,
128
138
.

Robert
C.
,
2007
.
CSE486, Lecture 6: Corner dectection
,
The Pennsylvania State University
,
http://www.cse.psu.edu/∼rtc12/CSE486/lecture05_6pp.pdf
.

Scholl
C.
,
Hallinan
S.
,
Miorelli
F.
,
Watts
M.D.
,
2017
.
Geological consistency from inversions of geophysical data
, in
79th EAGE Conf. and Exhibition 2017
,
EAGE
,
the Netherlands
,
doi:10.3997/2214-4609.201700849
.

Sheng
H.
,
Deng
S.
,
Zhang
S.
,
Li
C.
,
Xiong
Z.
,
2016
.
Segmentation of light field image with the structure tensor
, in
Image processing (ICIP), 2016 IEEE Int. conf.
,
IEEE
,
Phoenix, AZ
,
doi:10.1109/ICIP.2016.7532602
.

Tryggvason
A.
,
Linde
N.
,
2006
,
Local earthquake (LE) tomography with joint inversion for P- and S-wave velocities using structural constraints
,
Geophys. Res. Lett.
,
33
,
doi:10.1029/2005GL025485
.

Um
E.
,
Commer
M.
,
Newmann
G.
,
2014
.
A strategy for coupled 3D imaging of large-scale seismic and electromagnetic data sets: application to subsalt imaging
,
Geophysics
,
79
,
ID1
ID13
.

Weickert
J.
,
1999
.
Coherence-enhancing diffusion filtering
,
Int. J. Comput. Vis.
,
31
,
111
127
.

Wu
X.
,
2017
.
Directional structure-tensor-based coherence to detect seismic faults and channels
,
Geophysics
,
82
,
A13
A17
.

Yi
M.J.
,
Kim
J.H.
,
Chung
S.H.
,
2003
.
Enhancing the resolving power of least-squares inversion with active constraint balancing
,
Geophysics
,
68
,
931
941
.

Zhang
L.
,
Zhang
L.
,
Zhang
D.
,
2009
.
A multi-scale bilateral structure tensor based corner detector
,
Asian Conf. Comput. Vis.
,
618
627
.

Zhou
J.
,
Revil
A.
,
Karaoulis
M.
,
Hale
D.
,
Doetsch
J.
,
Cuttler
S.
,
2014
.
Image-guided inversion of electrical resistivity data
,
Geophys. J. Int.
,
197
,
292
309
.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)