New method for estimating strike and dip based on structural expansion orientation for 3D geological modeling

Abstract

Strike and dip are essential to the description of geological features and therefore play important roles in 3D geological modeling. Unevenly and sparsely measured orientations from geological field mapping pose problems for the geological modeling, especially for covered and deep areas. This study developed a new method for estimating strike and dip based on structural expansion orientation, which can be automatically extracted from both geological and geophysical maps or profiles. Specifically, strike and dip can be estimated by minimizing an objective function composed of the included angle between the strike and dip and the leave-one-out cross-validation strike and dip. We used angle parameterization to reduce dimensionality and proposed a quasi-gradient descent (QGD) method to rapidly obtain a near-optimal solution, improving the time-efficiency and accuracy of objective function optimization with the particle swarm method. A synthetic basin fold model was subsequently used to test the proposed method, and the results showed that the strike and dip estimates were close to the true values. Finally, the proposed method was applied to a real fold structure largely covered by Cainozoic sediments in Australia. The strikes and dips estimated by the proposed method conformed to the actual geological structures more than those of the vector interpolation method did. As expected, the results of 3D geological implicit interface modeling and the strike and dip vector field were much improved by the addition of estimated strikes and dips.

Graphical abstract

Keywords

strike and dip / structural expansion orientation / leave-one-out cross-validation / covered area

1 Introduction

Strike and dip are essential to the description, analysis, and modeling of subsurface geological structures. Three-dimensional geological modeling methods are widely used in oil and gas fields, groundwater aquifers, and ore deposits (Perrin and Rainaud, 2013). In most existing geomodeling methods, strike and dip play an important role as, for example, the local anisotropy direction in attribute modeling (Xu, 1996; Mariethoz and Kelly, 2011; Machuca-Mory and Deutsch, 2013) and the gradient constraint in implicit surface modeling (Lajaunie et al., 1997; Hillier et al., 2014; Guo et al., 2018). However, the spatial distribution of measured strikes and dips is often sparse and uneven, especially with limited funds and in areas with inaccessible topography. Therefore, it is often necessary to estimate and add extra strikes and dips in 3D geological modeling.

From a computational perspective, existing methods for adding strike and dip can be broadly categorized into geometric and interpolation approaches. Geometric approaches perform geometric operations on points, lines, and surfaces that are located at the same interface to calculate the additional strike and dip. The three-point method (Martinez-Torres et al., 2012), which is the most popular geometric approach, uses the elevation of three adjacent points in the same rock interface to calculate an approximate strike and dip. The three-point method was used to compute subsurface strike and dip via adjacent boreholes in densely drilled areas (Bucher, 1943). Geological maps, orthophotographs, and Light Detection and Ranging are being used in combination with digital terrain models to compute strikes and dips of the ground using the three-point method (Chorowicz et al., 1991; Banerjee and Mitra, 2004; Yeh et al., 2014; Thiele et al., 2019); however, the three-point method is unusable in covered areas. Subsurface strike and dip can also be calculated from seismic signals by geometric approaches. For example, Brown et al. (1971) calculated the strike and dip of crustal boundaries using seismic cross-profiles. Dalley et al. (2007) acquired the strike and dip field through curvature calculations using 3D seismic data. However, drilling or seismic data that can be used to add subsurface strike and dip are often lacking.

Interpolation approaches, which are based on the spatial correlation between the measured strike and dip, can theoretically add the strike and dip at any position. These approaches are often used to construct the strike and dip field (also known as the locally varying anisotropy field), which is used in geological stochastic simulation modeling (Xu, 1996; Mariethoz and Kelly, 2011; Martin et al., 2019) and structure visualization (Hillier et al., 2013). Since the incoherence of arithmetic means for periodic data cannot be avoided, the interpolation of angle data has always been difficult (Grancher et al., 2012). Therefore, the 3D unit vector instead of dip and strike angles is often used when interpolating the strike and dip (Young, 1987a, 1987b; Hillier et al., 2013). However, it is difficult to match the interpolated strike and dip with the actual geological structure in complex areas of the structure or where the measured strike and dip are sparse (Grose et al., 2019).

The aim of this study was to rapidly and automatically estimate the strike and dip based on structural expansion orientation (SEO) for 3D geological modeling. SEO is commonly used in numerous geological applications, such as local anisotropy direction in 2D interpolation (Boisvert et al., 2009) and tangent constraint in implicit interface modeling (Lajaunie et al., 1997; Hillier et al., 2014). SEO is orthogonal to the strike and dip, widely distributed, and easy to extract in practice (Lillah and Boisvert, 2015). We proposed a new method to obtain the strike and dip at SEO points (abbreviated as “SEOSD”) based on two types of constraints: The first is to guarantee the orthogonality of SEO and SEOSD through angle parameterization, and the second is to use the leave-one-out cross-validation (LOOCV) method, which can apply the spatial correlation between SEOSDs and measured strikes and dips when the objective function is optimized.

Since many SEOs can be extracted from a general scenario, high dimensional optimization problems may occur when the objective function is optimized. Instead of using nonlinear optimization methods (such as particle swarm optimization, PSO) to directly optimize the objective function, we proposed a comprehensive strategy that can quickly and accurately optimize the objective function. In this optimization strategy, the angle parameterization method was used to reduce dimensionality and the QGD method to quickly obtain a near-optimal solution that can be used as initial solution for further nonlinear optimization. Finally, the effectiveness and practicability of the SEO-based addition of strike and dip was validated using synthetic data generated by the simulation of a basin fold model and real complex fold-structure in Australia, of which approximately 28% is covered by Cainozoic sediments.

2 Methods

From the perspective of optimization problem, the SEOSD is the parameter to be optimized. As a popular parameter optimization method, the cross-validation method is often applied to the selection of parameters, e.g., selecting a good shape parameter for radial basis interpolation (Fasshauer and Zhang, 2007; Scheuerer, 2011), super-parameter in machine learning (Rodriguez et al., 2010), and variogram models (Marcotte, 1995). We used the LOOCV method (Rippa, 1999) to select an “optimal” SEOSD because of its high accuracy. Our method was divided into the following steps: 1) extract the SEO from multi-source geological data, 2) build the objective function based on the LOOCV method, and 3) propose an optimization strategy to rapidly and accurately calculate the SEOSD.

2.1 SEO extraction

SEO can be extracted from multi-source geological data. The common sources and methods for obtaining SEO are shown in Table 1. In the non-covered surface area of a geological map, the moment-of-inertia method (Lajevardi et al., 2015) can be used to obtain the SEO within a layer, and tangential sampling of a conformable contact boundary line can be used to obtain the SEO at an interface. Although it is impossible to obtain the SEO of bedrock from a geological map, the geophysical signals with high penetration, such as gravity and magnetic signals with high horizontal resolution, in the covered surface area can be used to automatically compute the SEO using the boundary detection method (Archibald et al., 1999), moment-of-inertia method (Lajevardi et al., 2015), or automatic feature interpolation method (Boisvert, 2010). Subsurface SEO can be automatically extracted from geological or geophysical profiles using the above methods (Scholl et al., 2017).

2.2 Objective function

In this study, the 3D unit vector was used to express the strike and dip (Hillier et al., 2013) and was called strike and dip vector \overrightarrow{\mathit{V}} (Fig. 1(a)). For any strike and dip vector,

(1)\overrightarrow{V}=\left[\begin{array}{l}{v}_x\hfill \\ v_y\hfill \\ v_z\hfill \end{array}\right]={\left[\begin{array}{l}\sin (\varphi -\mathrm{\pi }/\mathrm{2})\cos (\mathrm{\pi }/\mathrm{2}-\theta )\hfill \\ \cos (\varphi -\mathrm{\pi }/\mathrm{2})\cos (\mathrm{\pi }/\mathrm{2}-\theta )\hfill \\ -\sin (\mathrm{\pi }/\mathrm{2}-\theta )\hfill \end{array}\right]}^{\mathrm{T}},

where \varphi and \theta are the strike and dip angles, respectively, and {[•]}^{\mathrm{T}} is the transposed matrix of [•]. The 3D unit vector was used to express the SEO, which was named SEO vector \overrightarrow{\mathit{E}} (Fig. 1(a)). The SEO vector is orthogonal to the strike and dip vector. Supposing that there are m measured strike and dip vectors \overrightarrow{{\mathit{V}}^M}=\left[\overrightarrow{v_{\mathrm{1}}^M},\overrightarrow{v_{\mathrm{2}}^M},\mathrm{...},\overrightarrow{v_m^M}\right] and n SEOSD vectors to be solved \overrightarrow{{\mathit{V}}^E}=\left[\overrightarrow{v_{\mathrm{1}}^E},\overrightarrow{v_{\mathrm{2}}^E},\mathrm{...},\overrightarrow{v_n^E}\right] at the SEO points, then all strike and dip vectors are \overrightarrow{\mathit{V}}=\left[\overrightarrow{v_{\mathrm{1}}^M},\overrightarrow{v_{\mathrm{2}}^M},\mathrm{...},\overrightarrow{v_m^M},\overrightarrow{v_{\mathrm{1}}^E},\overrightarrow{v_{\mathrm{2}}^E},\mathrm{...},\overrightarrow{v_n^E}\right].

According to the LOOCV method, the LOOCV vector \overrightarrow{{\mathit{V}}_{\left[i\right]}} can be defined as

\overrightarrow{{\mathit{V}}_{\left[i\right]}}=\left[\overrightarrow{v_{\mathrm{1}}},\mathrm{...},\overrightarrow{v_{i-\mathrm{1}}},\overrightarrow{v_{i+\mathrm{1}}},\mathrm{...},\overrightarrow{v_{m+n}}\right],

(2)\mathrm{if}i\in [\mathrm{1},\mathrm{2},\mathrm{...},m+n].

The LOOCV strike and dip vector \overrightarrow{v_{\left[i\right]}^{CV}} interpolated using the LOOCV vector \overrightarrow{{\mathit{V}}_{\left[i\right]}} can then be written as

(3)\overrightarrow{v_{\left[i\right]}^{CV}}=\mathrm{interpolater}\left({\omega }_i^{\mathrm{u}\mathrm{n}\mathrm{c}}\overrightarrow{{\mathit{V}}_{\left[i\right]}}\right),

where \mathrm{interpolater}(•) is the 3D vector interpolation function, such as the inverse distance weighting interpolation method (Hillier et al., 2013), the kriging-based interpolation method (Young, 1987b), and the radial basis functions interpolation method (Smolik and Skala, 2016), and {\omega }^{\mathrm{u}\mathrm{n}\mathrm{c}} is the weight coefficient corresponding to the data uncertainty (\mathrm{0}\le {\omega }^{\mathrm{u}\mathrm{n}\mathrm{c}}\le \mathrm{1}). The greater the uncertainty about the strike and dip is, the closer {\omega }^{\mathrm{u}\mathrm{n}\mathrm{c}} is to 0. Where the uncertainty of the SEOSD is unknown, the SEO uncertainty can be used instead. It should be emphasized that when the LOOCV strike and dip vector is calculated by Eq. (3), both the measured strike and dip and the SEOSD is used. Moreover, the LOOCV strike and dip vector needs to be normalized to the 3D unit vector.

In this study, the estimation error was represented by the included angle between the strike and dip vector \overrightarrow{\mathit{V}} and the LOOCV vector \overrightarrow{{\mathit{V}}^{CV}}:

(4){\mathrm{IA}}_i=\arccos (\overrightarrow{v_i}\cdot \overrightarrow{v_{\left[i\right]}^{CV}}),\mathrm{if}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}i\in [\mathrm{1},\mathrm{2},\mathrm{...},m+n].

However, since the true SEOSD vector \overrightarrow{v_i}(i\in [m+\mathrm{1},m+\mathrm{2},\mathrm{...},m+n\left]\right) was unknown, we chose the approximation of the true SEOSD vector as an alternative. For example, the projection vector of the LOOCV strike and dip vector on the plane perpendicular to the SEO vector \overrightarrow{e_{i-m}}(i\in [m+\mathrm{1},m+\mathrm{2},\mathrm{...},m+n\left]\right) can be regarded as an approximation of the true SEOSD vector:

\overrightarrow{v_i}=\overrightarrow{v_{\left[i\right]}^{CV}}-(\overrightarrow{v_{\left[i\right]}^{CV}}\cdot \overrightarrow{e_{i-m}})\cdot \overrightarrow{e_{i-m}},

(5)\mathrm{if}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}i\in [m+\mathrm{1},m+\mathrm{2},\mathrm{...},m+n].

The objective function constructed in this study is given by

(6)\mathrm{ObjFunc}={\displaystyle \sum _{i=\mathrm{1}}^{m+n}\left({\omega }_i^{\mathrm{u}\mathrm{n}\mathrm{c}}{\mathrm{IA}}_i\right)}/(m+n)\to \min ,

where, the weight coefficient {\omega }^{\mathrm{u}\mathrm{n}\mathrm{c}} corresponding to the data uncertainty in Eq. (3) is also used.

2.3 Optimization strategy

Numerous SEOs can be extracted from many scenarios because SEO has a wide range of sources. However, the optimization algorithm faces problems (Wang et al., 2011) when there are many parameters to be solved. For example, many nonlinear algorithms show premature and slow convergence in high-dimensional complex optimization questions (Patnaik et al., 2017). Furthermore, because the computational complexity of the LOOCV operation is O\left(N^{\mathrm{4}}\right) (Rippa, 1999), the computational cost is high when a large number of parameters needs to be optimized.

To improve the efficiency and accuracy of the traditional nonlinear optimization method, the proposed optimization strategy was as follows (Fig. 2(a)): First, to reduce the dimension of the optimization problem, the geometric properties of the SEO and SEOSD were used to parameterize the SEOSD vector into the angle parameter. Secondly, the QGD method was proposed to rapidly obtain a near-optimal solution (Fig. 2(b)). Finally, if a high-accuracy solution is required, traditional nonlinear algorithms can be used for further optimization, in which the near-optimal solution of the QGD method can be used as initial solution.

2.3.1 Parameterizing strike and dip vector to the angle

For any SEO point, the SEOSD vector \overrightarrow{v_j^E}(j\in [\mathrm{1},\mathrm{2},\mathrm{...},n\left]\right) is perpendicular to the SEO vector \overrightarrow{e_j}(j\in [\mathrm{1},\mathrm{2},\mathrm{...},n\left]\right) (Fig. 1(b)). The reference SEOSD vector \overrightarrow{v_j^{E,\mathrm{0}}} is defined as

(7)\overrightarrow{v_j^{E,\mathrm{0}}}=\{\begin{array}{l}\left[\begin{array}{lll}{e}_{j,y}/\sqrt{e_{j,x}{}^{\mathrm{2}}+e_{j,y}{}^{\mathrm{2}}}\hfill & -e_{j,x}/\sqrt{e_{j,x}{}^{\mathrm{2}}+e_{j,y}{}^{\mathrm{2}}}\hfill & \mathrm{0}\hfill \end{array}\right]\mathrm{,\; if}{e}_{j,z}\ne \mathrm{1}\hfill \\ \left[\begin{array}{lll}\mathrm{0}\hfill & \mathrm{1}\hfill & \mathrm{0}\hfill \end{array}\right]\mathrm{,\; if}{e}_{j,z}=\mathrm{1}\hfill \end{array}.

Consequently, \overrightarrow{v_j^E} can be regarded as rotating the reference SEOSD vector \overrightarrow{v_j^{E,\mathrm{0}}} about the SEO vector \overrightarrow{e_j}=\left[\begin{array}{lll}{e}_{j,x}\hfill & e_{j,y}\hfill & e_{j,z}\hfill \end{array}\right] by the angle {\theta }_j (Taylor and Kriegman, 1994) using

(8)\overrightarrow{v_j^E}=\overrightarrow{v_j^{E,\mathrm{0}}}\mathrm{R}(\overrightarrow{e_j},{\theta }_j)\mathrm{,if}j\in [\mathrm{1},\mathrm{2},\mathrm{...},n],

where R(•) is the rotation matrix, written as

(9)R(\overrightarrow{e_j},{\theta }_j)=\left[\begin{array}{lll}{e}_{j,x}{}^{\mathrm{2}}(\mathrm{1}-\cos {\theta }_j)+\cos {\theta }_j\hfill & e_{j,x}{e}_{j,y}(\mathrm{1}-\cos {\theta }_j)-e_{j,z}\sin {\theta }_j\hfill & e_{j,x}{e}_{j,z}(\mathrm{1}-\cos {\theta }_j)+e_{j,y}\sin {\theta }_j\hfill \\ e_{j,x}{e}_{j,y}(\mathrm{1}-\cos {\theta }_j)+e_{j,z}\sin {\theta }_j\hfill & e_{j,y}{}^{\mathrm{2}}(\mathrm{1}-\cos {\theta }_j)+\cos {\theta }_j\hfill & e_{j,y}{e}_{j,z}(\mathrm{1}-\cos {\theta }_j)-e_{j,x}\sin {\theta }_j\hfill \\ e_{j,x}{e}_{j,z}(\mathrm{1}-\cos {\theta }_j)-e_{j,y}\sin {\theta }_j\hfill & e_{j,y}{e}_{j,z}(\mathrm{1}-\cos {\theta }_j)+e_{j,x}\sin {\theta }_j\hfill & e_{j,z}{}^{\mathrm{2}}(\mathrm{1}-\cos {\theta }_j)+\cos {\theta }_j\hfill \end{array}\right],

and {\theta }_j is the angle parameter defined as

(10)

2.3.2 The QGD method

The LOOCV SEOSD vector, which is interpolated by the surrounding strike and dip vectors, can be regarded as an approximation of the true SEOSD vector. Similar to the gradient descent method (Snyman, 2005), the difference between the SEOSD and LOOCV SEOSD vectors was used as the “approximate” gradient. Accordingly, we proposed a fast algorithm, namely the QGD method. The steps of the method are as follows:

1) Initialize the SEOSD vector and its angle parameters {\theta }_{\mathrm{0}} using Eq. (10).

2) Calculate the LOOCV strike and dip vector \overrightarrow{V^{CV}} at all strike and dip vector points using Eq. (3).

3) Compute the projection vector of the LOOCV strike and dip vector based on Eq. (5) at the SEO points and compute the angle parameter {\theta }^{CV} between the projection and reference SEOSD vectors using Eq. (10).

4) Update the angle parameter \theta using Eq. (11) and compute the updated SEOSD vector \overrightarrow{{\mathit{V}}^E} according to Eq. (8). When the kth iteration is performed, the angle parameter is updated as

(11){\theta }_k={\theta }_{k-\mathrm{1}}+\alpha ({\theta }^{CV}-{\theta }_{k-\mathrm{1}}),

where, \alpha is the step length (taken as 1 in this study).

5) Evaluate the objective function using Eq. (6). If the objective function decreases and does not reach the stop threshold value, repeat steps (2) to (5); if not, proceed to the next step.

6) Choose the SEOSD vector \overrightarrow{V^E} computed in step (4) of the last loop as the result.

2.3.3 Nonlinear optimization algorithm

Many traditional nonlinear optimization algorithms can be used to solve our objective functions, such as PSO (Sammut and Webb, 2011), genetic algorithms (Haupt and Ellen Haupt, 2004), and simulated annealing (Wang et al., 2011). In this study, PSO was selected because of its simple parameter setting and because it does not easily fall into a local optimum. The rough solution obtained by the QGD method was used as the initial PSO solution. After optimization by PSO, the optimal angle parameter was obtained, and the SEOSD vector was calculated using Eq. (8).

3 Application

3.1 Synthetic example

A synthetic basin fold model was used to test the proposed method (Fig. 3(a)). The simulated deposition time function is \mathrm{T}(x,y,z)=\mathrm{200}-\sqrt{x^{\mathrm{2}}{\mathrm{+}y}^{\mathrm{2}}\mathrm{+}{(z-\mathrm{500})}^{\mathrm{2}}}\mathrm{Ma/m}. In this synthetic example, the weight coefficient {\omega }^{\mathrm{u}\mathrm{n}\mathrm{c}} in Eqs. (3) and (6) was set to 1, which means that the data uncertainty was not considered. In addition, the inverse distance weighted interpolation method, which has a high calculation speed and is easily implemented (Hillier et al., 2013; Mariethoz and Kelly, 2011), was chosen to interpolator the vector in this synthetic example.

3.1.1 Estimate the SEOSD

Ten measured strikes and dips on the ground (surface of z=\mathrm{0}\mathrm{m}) and 50 SEOs in the entire space were randomly sampled (Fig. 3(b)). Simulated measured strike and dip vector \overrightarrow{V^M}was computed by \overrightarrow{V^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}=-\nabla T/\Vert \nabla T\Vert. Simulated SEO vector \overrightarrow{E} was randomly selected, which met \overrightarrow{V^{E,\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}·\overrightarrow{E}=\mathrm{0}. The SEOSD was subsequently calculated using traditional PSO and our optimization strategy (Section 2.3). All approaches were based on angle parameterization. Before the optimization, the initial SEOSD was interpolated by the inverse distance weighting interpolation method. The Pearson correlation coefficient between the initial and true SEOSD was small (Fig. 5(a) and 5(b)); that of the dip angle was especially weak (0.49). The QGD method was used to reduce the objective function from 7.37° to 6.19° after only four iterations; the objective function converged to 6.09° after 12 iterations, which approached the optimal value (Fig. 4(b)). However, traditional PSO required 503 iterations to converge to 6.09° (Fig. 4(a)), and the computational cost was about 503 × 100/12≈5192 times that of the QGD method. After optimization with our optimization strategy, the Pearson correlation coefficient of the strike and dip between the estimated and true SEOSD was 0.99 and 1, respectively (Fig. 5(c) and 5(d)).

The strike and dip vector field indicating the anisotropy direction is widely used in stochastic simulation modeling. Another important role of the strike and dip vector field is visualizing the structure to aid geological interpretation, e.g., Hillier et al. (Hillier et al., 2013) used the strike and dip vector field to visualize the structure trend by the form line. Based on these authors’ work, we incorporated the line integral convolution method (Cabral and Leedom, 1993) to calculate the vector field texture of the structural geometric contour corresponding to the strike and dip vector field in the profile. (Appendix A describes the method in detail.) Since this method can provide detailed visualization of the strike and dip vector field, it can be used to visually check whether the strike and dip vector field matches the geological structure. The vector field textures corresponding to the strike and dip vector sets \left\{\overrightarrow{V^M}\right\}, \{\overrightarrow{V^M},\overrightarrow{V^E}\}, and \{\overrightarrow{V^M},\overrightarrow{V^{E,\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}\} of profile AA' are shown in Figs. 6(b)–6(d). As a comparison, the true strike and dip vector field of profile AA' was computed by \overrightarrow{V^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}=-\nabla T/\Vert \nabla T\Vert, and the corresponding vector texture is shown in Fig. 7(a). Since there were no measured subsurface strike and dip and strike and dip vector near the A' side of profile AA' (Fig. 3(b)), the structural vector texture (Fig. 6(b)) obtained using only the surface measured strike and dip vector \left\{\overrightarrow{V^M}\right\} had a significant morphological deviation compared to the theoretical model (Fig. 6(a)). However, the vector field texture of the strike and dip vector set \{\overrightarrow{V^M},\overrightarrow{V^E}\} (Fig. 6(d)) was consistent with the theoretical model (Fig. 6(a)), and Fig. 6(d) is similar to Fig. 6(c), which shows the vector field texture of the strike and dip vector set \{\overrightarrow{V^M},\overrightarrow{V^{E,\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}\}.

3.1.2 Evaluation of the QGD method

Using the QGD method to calculate a rough solution was a critical step in the proposed optimization strategy. To evaluate the accuracy and speed of the QGD method, the SEOSD vector \overrightarrow{V^E} was computed by different SEO vectors \overrightarrow{E} and measured strike and dip vectors \overrightarrow{V^M}. The accuracy of our results was evaluated by the mean included angle between the calculated SEOSD vector \overrightarrow{V^E} and the true SEOSD vector \overrightarrow{V^{E,\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}}:

(12){\mathrm{MIA}}^{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}={\displaystyle \sum _{i=\mathrm{1}}^n\arccos (\overrightarrow{v_i^E}\cdot \overrightarrow{v_i^{E,\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l}}})}/n.

After the SEOSD was computed 100 times with different numbers of SEOs and measured strikes and dips, the mean and variance of the mean included angle were calculated (Fig. 7(a) and 7(b)). When the number of SEOs was constant, increasing the number of measured strikes and dips had a limited effect on reducing the mean included angle because the measured strike and dip points were distributed on the ground. However, when the number of measured strikes and dips was constant, the average of the mean included angle decreased with increasing number of SEOs; the average of the mean included angle was reduced by about half when the number of SEOs increased to 100. The standard deviation of the mean included angle mostly ranged between 0.5° and 3°, which indicates that the QGD method is stable. The stop criterion (where k is the number of iterations) was chosen, and the average number of iterations was computed (Fig. 7(c)). With increasing number of SEOs, the number of iterations increased, and only 5–24 iterations were needed to reach stable convergence.

3.2 Real application

The study area is located in the Bonney Well region, Northern Territory, Australia (Fig. 8(a)) and extends from 134°21′E–134°42′E and 20°23′S–20°50′S. The geological map (Fig. 8(a)), aeromagnetic map (Fig. 8(b)), and geological profile AA' (Fig. 9(a)) used in the case study are collected from geoscience australia (Carr et al., 2016). Several sets of Proterozoic rock strata are exposed, and approximately 28% of the area is covered by Cainozoic sediments. The bedrock is folded with SE–NW-trending axes (Macdonald et al., 2005; Grose et al., 2019). Ninety ground strikes and dips were unevenly distributed along the exposed bedrock. Because the faults are mostly small, this study did not consider the influence of faults, for the sake of simplicity. However, the effect of the fault can be considered later by adding the drift function during the vector interpolation process (Verly et al., 1984). In addition to the geological map, an expert profile AA' in the SW–NE direction diagonally traversed the center of the study area (Fig. 9(a)).

Tangent sampling of the depositional contact boundary of the bedrock strata in the geological map and profile AA' was the most direct way to determine SEO. However, since the quality standards of geological mapping are unknown, the corresponding uncertainty weight of SEO obtained from the geological map was uniformly selected as 0.8. Considering that the uncertainty of expert cognition increases with depth and there is no borehole to constrain geological profile AA', the corresponding uncertainty weight of the SEO extracted by this profile was selected as {\omega }_{\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{f}\mathrm{i}\mathrm{l}\mathrm{e}}^{\mathrm{u}\mathrm{n}\mathrm{c}}=\mathrm{1}-h/H (where h is depth with a range of \mathrm{0}\le h\le \mathrm{10}\mathrm{km} and H=\mathrm{10}\mathrm{km}). Although SEO could not be extracted from the geological map in the area covered by Cainozoic sediments, the aeromagnetic map (Fig. 8(b)) corresponded well with the bedrock. The moment-of-inertia method was used to extract the SEO from the aeromagnetic map, where the local primary directions represent the SEO. The corresponding uncertainty weight was selected as {\omega }_{\mathrm{m}\mathrm{a}\mathrm{g}\mathrm{n}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{c}}^{\mathrm{u}\mathrm{n}\mathrm{c}}=\mathrm{1}-I_{xx}/I_{xy} (I_{xx} and I_{xy} are the moments of inertia in the primary and secondary directions, respectively). Through the above method, 699 SEOs were extracted from the geological map, aeromagnetic map, and profile AA' (Figs. 8(a), 8(b), and 9(a)).

Since the variation of strike and dip vectors has anisotropic characteristics, the anisotropic inverse distance weighting method (Tomczak, 1998) was used to perform 3D vector interpolation. The initial SEOSD was interpolated by measured strikes and dips (Figs. 9(b) and 10(a)). However, many interpolated SEOSDs did not match the actual structure, especially in the intense bending of the formation (zones I and III in Fig. 10(a)) and Cainozoic coverage (zones II and III in Fig. 10(a)). Considering that many SEOs were extracted, the QGD method was selected to estimate the SEOSD. In the optimization process, the objective function decreased rapidly (Fig. 11(a)), and the mean included angle calculated by Eq. (4) decreased from 9.40° to 4.86° with only 10 iterations. Finally, the number of included angles between 0° and 3° increased from 27.5 to 61.2% (Figs. 11(b) and 11(c)). The SEOSD estimated by the QGD method not only followed the actual geological structure, but the spatial transition of estimated SEOSDs was also more stable (Figs. 9(c) and 10(b)). To determine whether the strike and dip vector field matched the real structure, the vector field texture of three profiles, namely AA', BB', and CC' (Figs. 9(b), 9(c), and 12), was calculated by using measured strikes and dips and all strikes and dips (including measured strikes and dips and SEOSD estimated by the QGD method). In the near surface, the vector field textures of these three profiles were similar and showed small folds. However, in the deeper surface, these small folds disappeared or flipped in the vector field texture of measured strikes and dips (Figs. 9(b), 12(a) and 12(c)), which is inconsistent with the actual geological structure.

The interfaces between kudinga basalt or frew river formation and coulters sandstone were modeled using the radial basis function-based implicit modeling method (Hillier et al., 2014), where the cube function was chosen as the radial basis function. In the modeling process, on-contact points (Fig. 8(a)) were taken as the same-interface constraint, and measured strikes and dips and all strikes and dips (including estimated SEOSDs) were used as gradient constraints. The interface built using the measured strike and dip and on-contact points showed that the folds in the blue dotted circle did not extend as expected (Fig. 13(a)). However, the interface model was smooth and followed the actual geological structure (Fig. 13(b)) when estimated SEOSDs were added.

4 Discussion

An important consideration for using SEO to supplement the strike and dip is that it can be automatically extracted from multiple types of geological information. In particular, SEOs can be extracted from geophysical signals with high penetration, and the strike and dip that match the actual geological structure in covered and deep areas can thus be added, whereas this is difficult in other methods. In the case study, the SEOSDs estimated by the QGD method matched the actual geological structure better than did the initial SEOSDs interpolated by the vector interpolation method (Figs. 9(b) and 10(a)).

In this study, the objective function was constructed based on the estimated deviation, including the angle between the strike and dip and the LOOCV strike and dip. When the LOOCV strikes and dips were calculated, both the measured strikes and dips and SEOSDs to be estimated were used. Therefore, the objective function ensured that all strikes and dips were spatially related to each other, thus improving the accuracy of the estimated SEOSDs and extrapolation ability of the method. In the simulation test, although there were only 10 measured (and sparse) strikes and dips, the mean included angle continually decreased with increased number of SEOs (Fig. 7(a)). Similarly, although there were not many measured strikes and dips in regions I and III of the case study (Fig. 10(b)), the SEOSD spatial variation obtained by our method was spatially stable and matched the geological structure. However, since the SEO essentially comprises additional information, the measured strikes and dips should be as dense as possible, covering the entire area, to improve the accuracy of the SEOSD.

In theory, the objective function (Eq. (6)) can be directly optimized by traditional nonlinear optimization algorithms. However, when many SEOs are being extracted, these traditional algorithms face problems in terms of computing speed and accuracy caused by high-dimensional optimization. In this study, angle parameterization and the QGD method were used to mitigate high-dimensional optimization. After angle parameterization, solving an SEOSD vector is equivalent to solving scalar parameter \theta, but the number of parameters to be solved is reduced by two-thirds. There are two more advantages to optimization through angle parameterization: 1) The solved SEOSD vector is mathematically guaranteed to be perpendicular to the SEO vector. 2) The parameter θ to be optimized is a cycle angle variable, which is advantageous for parameter searching. In addition, spatial sampling (Carmichael and Ailleres, 2016) of SEO can be used to mitigate the high-dimensional optimization problem, which could be studied in the future.

Since the true SEOSD is unknown in the iterative process of the QGD method, we believe that the LOOCV SEOSD can be considered an approximation of the true SEOSD. In the simulation model test, the results obtained by the QGD method and PSO were close (Fig. 4). Moreover, the SEOSD obtained by the QGD method was consistent with the actual geological structure in the real application (Fig. 10(b)). When the number of SEOs is too large or an approximate result that can be quickly obtained is acceptable, we suggest that the results of the QGD method can be selected.

5 Conclusions

We proposed a new method for automatically estimating the strike and dip by using SEOs, which are widely distributed and easily extracted from geological and/or geophysical maps. This method addresses an important problem in 3D geological modeling in terms of estimating strikes and dips in covered and deep areas when using geometric approaches (e.g., the three-point method). The synthetic example and real application suggested that the SEOSD estimated by the proposed method was more consistent with the actual geological structure than that estimated by vector interpolation approaches because of SEO constraints in the former. The proposed method cannot only be widely applied to 3D modeling scenes, such as oil, gas, and mineral resources exploration, which is closely related to fold structure, but it also has practical significance for supplementing the strike and dip of fault surfaces when modeling faults.

References

[1] Archibald N, Gow P, Boschetti F (1999). Multiscale edge analysis of potential field data. Explor Geophys, 30(1–2): 38–44

[2] Banerjee S, Mitra S (2004). Remote surface mapping using orthophotos and geologic maps draped over digital elevation models: application to the Sheep Mountain anticline, Wyoming. AAPG Bull, 88(9): 1227–1237

[3] Boisvert J (2010). Geostatistics with Locally Varying Anisotropy. Dissertation for the Doctor Degree. Alberta: University of Alberta

[4] Boisvert J B, Manchuk J G, Deutsch C V (2009). Kriging in the presence of locally varying anisotropy using non-Euclidean distances. Math Geosci, 41(5): 585–601

[5] Brown R J, Borg H, Båth M (1971). Strike and dip of crustal boundaries — a method and its application. Pure Appl Geophys, 88(1): 60–74

[6] Bucher W H (1943). Dip and strike from three not parallel drill cores lacking key beds (stereographic method). Econ Geol, 38(8): 648–657

[7] Cabral B, Leedom L C (1993). Imaging vector fields using line integral convolution. In: Kajiya J T, ed. SIGGRAPH ’93: Proceedings of the 20th annual conference on Computer graphics and interactive techniques. Anaheim, CA, 263–270

[8] Carr L, Ayling B, Morris R, Higgins K, Edwards D, Khider K, Smith T, Lawrie J, Morse M, Boreham C, Tenthorey E, Buckler T, Cox W, Hatch L, Woods M, Cortese A, Southby C, Chirinos A (2016). Georgina Basin geoscience data package, Various repository formats

[9] Carmichael T, Ailleres L (2016). Method and analysis for the upscaling of structural data. J Struct Geol, 83: 121–133

[10] Chorowicz J, Breard J, Guillande R, Morasse C, Prudon D, Rudant J (1991). Dip and strike measured systematically on digitized three-dimensional geological maps. Photogramm Eng Remote Sensing, 57(4): 431–436

[11] Dalley R M, Gevers E, Stampfli G M, Davies D J, Gastaldi C N, Ruijtenberg P A, Vermeer G (2007). Dip and azimuth displays for 3D seismic interpretation. First Break, 25(12): 101–108

[12] Fasshauer G E, Zhang J G (2007). On choosing “optimal” shape parameters for RBF approximation. Numer Algorithms, 45: 345–368

[13] Grancher D, Bar-Hen A, Paris R, Lavigne F, Brunstein D (2012). Spatial interpolation of circular data: application to tsunami of December 2004. Adv Appl Stat, 30(1): 19–29

[14] Grose L, Ailleres L, Laurent G, Armit R, Jessell M (2019). Inversion of geological knowledge for fold geometry. J Struct Geol, 119: 1–14

[15] Guo J, Wu L, Zhou W, Li C, Li F (2018). Section-constrained local geological interface dynamic updating method based on the HRBF surface. J Struct Geol, 107: 64–72

[16] Haupt R L, Ellen Haupt S (2004). Practical Genetic Algorithms. 2nd ed. New Jersey: John Wiley & Sons

[17] Hillier M, Kemp E, Schetselaar E (2013). 3D form line construction by structural field interpolation (SFI) of geologic strike and dip observations. J Struct Geol, 51: 167–179

[18] Hillier M J, Schetselaar E M, de Kemp E A, Perron G (2014). Three-dimensional modelling of geological surfaces using generalized interpolation with radial basis functions. Math Geosci, 46(8): 931–953

[19] Lajaunie C, Courrioux G, Manuel L (1997). Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation. Math Geol, 29(4): 571–584

[20] Lajevardi S, Babak O, Deutsch C V (2015). Estimating barrier shale extent and optimizing well placement in heavy oil reservoirs. Petrol Geosci, 21(4): 322–332

[21] Lillah M, Boisvert J B (2015). Inference of locally varying anisotropy fields from diverse data sources. Comput Geosci, 82: 170–182

[22] Macdonald F A, Mitchell K, Stewart A J (2005). Amelia Creek: a proterozoic impact structure in the Davenport Ranges, Northern Territory. Aust J Earth Sci, 52(4-5): 631–640

[23] Machuca-Mory D F, Deutsch C V (2013). Non-stationary geostatistical modeling based on distance weighted statistics and distributions. Math Geosci, 45(1): 31–48

[24] Marcotte D (1995). Generalized cross-validation for covariance model selection. Math Geol, 27(5): 659–672

[25] Mariethoz G, Kelly B F (2011). Modeling complex geological structures with elementary training images and transform-invariant distances. Water Resour Res, 47(7): W07527

[26] Martin R, Machuca-Mory D, Leuangthong O, Boisvert J B (2019). Non-stationary geostatistical modeling: a case study comparing LVA estimation frameworks. Nat Resour Res, 28(2): 291–307

[27] Martinez-Torres L M, Lopetegui A, Eguiluz L (2012). Automatic resolution of the three-points geological problem. Comput Geosci, 42: 200–202

[28] Patnaik S, Yang X, Nakamatsu K (2017). Nature-Inspired Computing and Optimization: Theory and Applications. Switzerland: Springer Nature

[29] Perrin M, Rainaud J F (2013). Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models. Paris: Editions Technip

[30] Rippa S (1999). An algorithm for selecting a good value for the parameter c in radial basis function interpolation. Adv Comput Math, 11: 193–210

[31] Rodriguez J D, Perez A, Lozano J A (2010). Sensitivity analysis of kappa-fold cross validation in prediction error estimation. IEEE Trans Pattern Anal Mach Intell, 32(3): 569–575

[32] Sammut C, Webb G I (2011). Encyclopedia of Machine Learning. New York: Springer Science & Business Media

[33] Scheuerer M (2011). An alternative procedure for selecting a good value for the parameter C in RBF-interpolation. Adv Comput Math, 34(1): 105–126

[34] Scholl C, Hallinan S, Miorelli F, Watts M D (2017). Geological consistency from inversions of geophysical data. In: 79th EAGE Conference and Exhibition, 2017, Paris, France

[35] Smolik M, Skala V (2016). Vector field interpolation with radial basis functions. In: Proceedings of SIGRAD 2016, Visby, Sweden, 15–21

[36] Snyman J A (2005). Practical Mathematical Optimization. New York: Springer Science+Business Media

[37] Taylor C J, Kriegman D J (1994). Minimization on the Lie group SO (3) and related manifolds. Connecticut‎: Yale University

[38] Thiele S T, Grose L, Cui T, Cruden A R, Micklethwaite S (2019). Extraction of high-resolution structural orientations from digital data: a Bayesian approach. J Struct Geol, 122: 106–115

[39] Tomczak M (1998). Spatial interpolation and its uncertainty using automated anisotropic inverse distance weighting (IDW)-cross-validation/jackknife approach. J Geogr Inf Decis Anal, 2: 18–30

[40] Verly G, David M, Journel A G, Marechal A (1984). Geostatistics for Natural Resources Characterization. New York: Springer Science+Business Media

[41] Wang Y, Yagola A G, Yang C (2011). Optimization and Regularization for Computational Inverse Problems and Applications. Heidelberg, Dordrecht, London, New York: Springer

[42] Xu W (1996). Conditional curvilinear stochastic simulation using pixel-based algorithms. Math Geol, 28(7): 937–949

[43] Yeh C, Chan Y, Chang K, Lin M, Hsieh Y (2014). Derivation of strike and dip in sedimentary terrain using 3D image interpretation based on airborne LiDAR data. Terr Atmos Ocean Sci, 25(6): 775–790

[44] Young D S (1987a). Indicator kriging for unit vectors: rock joint orientations. Math Geol, 19(6): 481–501

[45] Young‚D S (1987b). Random vectors and spatial analysis by geostatistics for geotechnical applications. Math Geol, 19(6): 467–479

RIGHTS & PERMISSIONS

Higher Education Press