Quantitative Assessment of Two Skeletonization Algorithms Adapted to Rectangular Grids

The goal of this paper is to improve the capability of gravity inversion assessment of complex structures using a multiple-source model as indicated by a study of a foothill belt of northwestern Taiwan. In this study, a gravity inversion computer program is applied based on research developed by Tsai, which contains a Marquardt inversion algorithm for mathematical calculations and incorporates constraints on geological parameters. The computed and modified geological parameters are transformed into coordinates. The response of a proper geological model is calculated using the Talwani technique. The method is applied to a field example in a foothills belt of Taiwan that possesses complex subsurface structures. The density profile from 2D gravity data, obtained using our inversion computer program, reveals a good correlation with the geological model obtained from seismic, borehole, and geologic data. Furthermore, potential hydrocarbon traps associated with some of the interpreted geological structures (where no seismic data are available) in the eastern and western seismic sections are being evaluated and assessed for possible future development. This profile highlights features that are important to a good understanding of an area with complex subsurface structures. It also helps to define subsurface geology and assess potential hydrocarbon traps. This


INTRODUCTION
Two-dimensional (2-D) gravity inversion has been used for many years to automatically interpret gravity data (Al-Chalabi 1972;Corbato 1965;Dyrelius and Vogel 1972;Negi and Garde 1969;Tanner 1967;Bott 1967;Pedersen 1975;Enmark 1981;Garci'a-Abdeslem 2003).In a conventional approach, the forward modeling is based on a set of 2-D prisms of arbitrary shape (Talwani et al. 1959;Bhattacharyya and Navolio 1976;Isaacs 1966;Martin-Atienza and Garci'a-Abdeslem 1999) and uses density contrast varying in depth.It works well for gravity data collected across a sedimentary basin or an ore body.In general, seismic data provides a more precise subsurface picture than gravity data.
However, in an area of complex geological structures, such as the Chutung and Guansi foothills in northwestern Taiwan -areas with access difficulties in their higher elevations --seismic data quality is not good or unavailable.Gravity survey may be one of the most important techniques to aid seismic interpretation.Geological models obtained from conventional inversion of gravity anomalies are unable to satisfy current geological demands.A more powerful gravity inversion technique is required.
As pointed out by Roy et al. (2005), an observed gravity anomaly can be exactly reproduced by each of an infinite number of subsurface distributions other than the one that originally produced it.
The key to reducing this gravity ambiguity is to use multiple-source models for forward modeling and incorporate constraints while inverting gravity data.A reasonable gravity model can then be obtained.Tsai (1992) has studied the Colorado plateau west of the Albuquerque basin.In his research, he extensively modified a conventional two-dimensional gravity inversion program using multiple bodies incorporating real geological constraint such as depth, dip angles, fault throws, and densities.During implementation, the inversion and geological constraints work together.To demonstrate the validity of the method, Tsai (1992) selected the synthetic observed data for graben and syncline models to demonstrate the basic facilities of program.All inverted model responses calculated with or without constraints are compared to a given gravity anomaly.The Root Mean Square Error (RMSE) limits were less than 1 mGal.In addition, the method applied to an inverted density model from the field data of the Albuquerque basin in New Mexico, USA was geologically reasonable.
A gravity inversion-modeling program is developed here based on Tsai's 1992 study.This program is used to study complex subsurface geological structures in a foothills belt of northwestern Taiwan.The features of the algorithm used include an allowable variety of geological constraints such as bed densities, depths and dips, and fault locations and dips.Also included are the continuity of layering densities and boundaries, constant fault throws, and fixing of some model parameters.
The inversion program computes and corrects the geological parameters first, then transforms them into vertex coordinates and calculates the response of a model using a polygon method (Talwani et al. 1959).
This paper describes a gravity inverse program based on Tsai's algorithm for two-dimensional (2-D) gravity data and demonstrates its effectiveness with complex field data in Taiwan.Its results indicate that reasonable subsurface geological pictures can be obtained in the study area.The new algorithm has many possible applications, including for future oil exploration in the foothills belt of Taiwan.(Marquardt, 1963) to compute and correct the parameters of a model.In doing so, they apply a polygon method (Talwani et al. 1959) or rectangular prism (e.g., Tanner 1967;Negi and Garde 1969;Dyrelius and Vogel 1972;Braile et al. 1974).The 2-D algorithm takes the vertex coordinates of a model to calculate the gravity response of a model.In most cases, this technique only deals with comparatively simple subsurface structures.It is difficult to incorporate geological constraints.Tsai (1992) proposed a new parameterization algorithm to overcome the limitations imposed by conventional parameterization.

Most conventional inversion programs use the Marquardt inversion technique
First, fault locations and dips, bed depths and dips, and densities replace the vertex parameters used in conventional gravity modeling.Then, parameterization incorporating geological constraints and working with general multiple-source models is applied to improve the quality of the inversion results.

The Basic Model
As shown in Fig. 1, there are two types of parameters: fault parameters and bed parameters.The fault parameters include fault location ( ) The forward modeling program does not calculate the vertices of the source bodies directly.
Instead, it computes and corrects the geological parameters first, then transforms them into vertex coordinates and calculates the response of a model using the Talwani technique (Talwani et al. 1959).
The transformation of a model is straightforward as any source body is confined by several lines, which are determined by geological parameters.The vertices are obtained from the intersections of lines.
Calculating the response of one source body at one time, each block's contribution is accumulated from the bottom bed to the top.The total response of the model is obtained by summing up the contribution of each source body (Tsai 1992).

The Least Squares Problem
The inversion scheme is based on a least squares approach.A starting model, This model should not be too far from a model that fits the measured data, is not linear, we expand it in a Taylor series around p r , and ignoring second and higher order terms we obtained: To minimize the error i E , a loss function i Q is defined: Where D r is the column vector: ) , ( ) ( A generalized inverse of matrix A, the first derivatives of the model responses at each observed point on the surface, with respect to the parameters of the model. The properties of s H can be explored using the generalized inverse of A as suggested by e.g.Lanczos (1961), and A is decomposed using the SVD (Singular Value Decomposition) routine given by Golub and Renisch (1970).
Furthermore, weighting factors controlled the correction of the parameters within a reasonable range for each element of A, as well as the correction term of each parameter.Unreasonable corrections are thereby avoided.An iterative second-order Marquardt Least Squares scheme (Marquardt 1963) was used to stabilize the inversion process.

FIELD EXAMPLE
Here the program has been tested by using a set of gravity data provided by the Chinese Petroleum Corporation that covers the Hsinchu and Wufeng areas in northwestern Taiwan.The Chutung, with the nearby Paoshan and Chingtasohu areas, are the main gas production regions of this northwestern foothills belt of Taiwan.Seismic results and well information indicate that the Hsinchu and Wufeng areas have complex structures.As shown in Figs. 2 and 3, the folded thrust fault can be grouped into two zones: the folded zone and en echelon fault zone.In general, an anticline underneath the footwall of a thrust fault or a high angle fault is not easily distinguished from a seismic section.In addition, if the foothills have an irregular topography and complex subsurface structures, the seismic data obtained from this site is poor; sometimes no data is available.
Gravity interpretation of subsurfaces is generally less accurate than seismic interpretations.
However, a gravity survey has the advantages of being rapidly implemented in foothill belts.It may aid interpretation of a seismic section.Specifically, it is possible to refer to the seismic section to set the skeleton of the initial gravity model then invert the measured gravity data to obtained a proper interpreted subsurface geological model.

Geological Setting
The study area is located at the Hsinchu and Taoshan villages.Late Cenozoic to Middle Miocene rocks are exposed to the east of the foothill belt.The clastic deposits in the Late Cenozoic Basin of the Western Foothills belt range in age from late Oligocene to early Pleistocene (Huang 1978 and1980).The Oligocene to Miocene sediment in the Western Foothills belt can be divided into two types of facies on the basis of inferred depositional environments: shelf type and basin type (Ho 1975 and1979).The strata of the shelf type are characterized by near-shore deposits of mixed marine and continental origin.They are mainly white to light gray orthoquartzitic or arkosic sandstones, thin coal beds, dark shale and thin interlaminations of silt, clay, and sand.The strata of the basin type were deposited under conditions of greater subsidence and rapid accumulation and are exclusively marine.
They are represented by poorly sorted clastic sediment, either as thick sandstone or shale members or else in a monotonous sequence of alternating fine-grained sandstones, subgraywackes, and dark shales or claystones.The small density contrast between each layer makes identification of faulted layers difficult.

Integrated 2-D Seismic Profiling and Geology
A 2-D seismic line runs roughly parallel to the geological model section, showing basement fault blocks and tying key wells.The seismic image covers only the area from the Hsincheng Fault to the Chutung Fault (Fig. 3a).The extension of the west and east sides of the interpreted geological model were based on geological and well information.With the aid of this geological and well data, a geological model was established.On the surface, these structures are, from west to east, the Hsinchu Fault, Chingtasohu-Chiting Anticline, Hsincheng Fault, Paoshan Anticline, Chutung-Peipu Fault, Chutung Anticline, Juanchiao Fault, Taping Fault, Taping Anticline, and Okungchi Fault (Fig. 3b).
Some of these structures can be observed on the surface and others are very likely due to blind faulting.
The Hsinchu Fault divides into two high angle thrust faults.One of these is at a lower thrust angle and extends downward to the bottom of Peiliao Formation or the top of Shidi Formation, where the fault plane turns flat.The Chingtasohu-Chiting Anticline on the hanging wall of the Hsinchu Fault has a steep front limb and a gentle dip along its back limb.
The footwall of Hsincheng Fault, located between Hsinchu Fault and Paoshan Anticline, is a normal fault in an early stage of evolution.The Hsincheng thrust fault can be discerned on the seismic section (Fig. 3a).This low angle thrust fault is located between a set of tilted reflectors over the fault plane and another set of horizontal reflections underlying the fault plane.The convex Hsincheng Fault divides into two low branches, which become decollement at a depth; the convex branch turns into a bedding surface within the shallower and the other within Kueichulin Formation.The eastwardly extension of these two branches faults within the Hopai Formation and Kueichulin Formation merge to the Chutung fault after dislocating two normal faults (Yang et al. 1996).

2-D Gravity Modeling and Interpretation
The gravity data was collected using a Worden gravity meter from the Chinese Petroleum Corporation (CPC).The gravimeter looped back to the base station every two hours for drift correction.
The data was tied to the reference station at Chu-Ming Building, Miaoli.Afterwards, drift and tide corrections (reoccupy some of the stations periodically during a gravity survey), free-air, Bouguer corrections, and terrain effect were applied.The gravity data was digitized into 6637 grid points with grid spacing of about 300 m.The resulting Bouguer anomaly map is shown on Fig. 2. We selected one profile (red solid line on Fig. 2) for 2-D modeling, utilizing the algorithm adapted for this research.
Inverse structural modeling uses the iterative feedback approach.Seismic Line 3 and its formal interpretation by CPC (Fig. 3) were used to set the initial model for our gravity inversions.Before quantitative modeling could be done, representative values for the density of the different rock types in the study area were determined using density-log data at various locations throughout the study area.
Density-log data was used not only to support assumptions about the constant density contrast between the layers, but also to estimate its value.The density of the sedimentary section was obtained from wells PS-14, EL-1 and CTU-26 and also the reports of Pan (1967), Hsieh and Hu (1972) and Liu and Tsai (1999).The average density obtained from the density-log measurements of boreholes in the study area and its surroundings are summarized in Table 1.During the inverted the Bouguer gravity anomaly data, the constraints of the thickness of layers and the position and dip of faults from seismic and geological surveys were applied.However, the seismic section covering the range is less than that of gravity data, therefore the gravity data may supply the subsurface information in the eastern and western sides of the study area where no seismic data available.The densities of layers obtained from wells insitu were also used to constrain the range of density of each layer.We used different weighting to control the model parameters during inversion.
The NW-SE direction of the Bouguer gravity data profile (Fig. 4a) is about 30 km long

CONCLUSIONS
This study illustrates the application of gravity data to in-depth detection and delineation of subtle faulting and anticlines in a foothill belt.The results were achieved by applying a multiple-source model inversion method to the gravity data of the foothill belt portion of Hsinchu to Taoshan villages.
In general, the seismic reflection method gives a more accurate subsurface structure than other geophysical methods.
The results of this study are in close agreement with previous structure imaging using seismic interpretations.Furthermore, potential hydrocarbon traps associated with some of the interpreted geologic structures (where no seismic data are available) in the eastern and western seismic sections are being evaluated and assessed for possible future development.
This method confirms subtle discontinuities in the located subsurface structure via an iterative approach.The results clearly show several structural features (such as faults and anticlines).It also provides more reliable geological information about the structure from Hsincheng Fault to Chutung Fault than that provided by seismic imaging.This study demonstrates that gravity data combined with previous seismic interpretation can be a useful, cost-effective tool for imaging complex structures in the foothill belts of Taiwan.Table 1.The average Formation densities in the study area (Liu and Tsai, 1999).

FORMATION DENSITY (g cm
e., the intersection of the fault plane (line) with the earth's surface, and fault dip angle j θ for the jth fault.The bed parameters include the bed location angle ij θ of the bottom surface of the bed, and the density ij ρ for the jth bed in the ith block.Confining the outline of the model to a rectangular shape, the top and bottom of the model is flat, and ij X can be taken as the location of a vertical borehole, to be constant for each block (not shown in the figure).
an m by n system of equations, which are solved simultaneously to estimate p r δ :

(
dashed line).The NW edge is a distance of 0.693 km from Hsinchu Fault, and extends southeastward to the edge of Taoshan Village for about 1.75 km.The observed gravity data generated from a geological model containing folding and faulting with nonuniform thickness is shown in the dotted line.Applying the above gravity inversion procedure in this gravity anomaly profile, the final inversion density section is plotted in Fig.4b.The RMSE between the observed data (dashed line) and the model response (solid line) is 0.569 mGal.As shown on the figure, the west of the footwall of Hsincheng thrust fault continues to have the characters of a normal fault system but its east has a lie low anticline.Between the normal fault and anticline are Paoshan Anticline and reversal structures.Comparison made between Fig.4band Fig.4creveals that this method can clearly detect and locate faults.In addition, the most significant results are shown in Fig.4b:(1) the density of layers underneath the region between Chutung and Taping anticline is lower than its average density.It is possible these layers were fractured when the overlying units slid over them or when they were in the presence of a gas reservoir.(2) The depth of the Toukoshan Formation and Cholan Formation between the Paoshan Anticline and Chutung Fault are shallower than we expected.(3) The density of Tungkeng Formation is lower than its average density in the region between the Taping Fault and Okungchi Fault, which may be associated with faulting.The reversal of layers in the southeast side of the density model indicates high angle thrust fault structures at a shallow depth.

Fig. 2 .
Fig.2.Geological sketch and Bouguer gravity anomalies maps of the study area, showing the known surface evidence of the main faults.Locality of selected gravity profile (red straight line) is also shown on the figure.

Fig. 3 .
Fig.3.(a) Seismic section of Line 3 from the Hsincheng Fault to Chutung Fault, (b) interpreted geological model of Line 3 based on the seismic section shown on (a), boreholes data and geological data.

Fig. 4 .
Fig.4.(a) Gravity profile obtained from the Hsinchu to Taoshan villages.The seismic L3 closes to the center of the gravity profile, (b) 2D density model section obtained from the inversion of the gravity profile data shown on Fig. 4(a), 4(c) location of seismic section of Line 3 corresponds to the 2D density model .