Method and device for constructing three-dimensional fault plane according to seismic directory

Abstract

The invention discloses a method and device for constructing a three-dimensional fault plane according to an earthquake catalog, and the method comprises the steps: extracting the position information of an earthquake source according to the earthquake catalog of a target region, obtaining the density information of each data point, and generating a density map; based on the density information and the density map, utilizing an OPTICS algorithm to carry out preliminary clustering analysis on the seismic source position information, outputting an accessibility map, and determining a clustering number N; calculating the local density and the center offset distance of each seismic source point based on the Euclidean distance and the truncation distance between every two seismic source points, generating a decision diagram, and selecting the first N seismic source points with the maximum product of the local density and the center offset distance as a clustering center; and PCA analysis and PCA plane determination are carried out on the DPC classification result, and plane features of the fault plane are obtained when set conditions are satisfied. According to the method, an OPTICS algorithm, a DPC algorithm and PCA analysis are combined, and end-to-end parameterization-free three-dimensional fault plane modeling is achieved only by means of seismic source position information.

Description

Technical Field

[0001] This invention relates to the field of earthquake three-dimensional fault plane construction technology, and in particular to a method and apparatus for constructing three-dimensional fault planes based on an earthquake catalog. Background Technology

[0002] Three-dimensional geological modeling is a core technology for earthquake disaster research, land spatial planning, and site selection for major projects. One of its key tasks is to construct high-precision three-dimensional fault models by integrating multi-source data such as surface fault traces, geophysical profiles, and earthquake catalogs. However, traditional methods are highly dependent on data, specifically in the following aspects:

[0003] (1) Limitations of data sources: Surface fault surveys rely on manual mapping or high-resolution remote sensing technology, but in some areas, due to complex terrain, vegetation cover, or lack of systematic surveys, surface fault trace data are missing; geophysical profiles (such as seismic reflection, gravity, and magnetic exploration) can reveal underground structures, but data collection costs are high and it is difficult to popularize them across the entire region.

[0004] (2) Shortcomings of existing modeling methods: Existing three-dimensional fault modeling techniques usually require joint inversion of multiple data sources (such as surface traces and geophysical profiles). When data is missing or economic conditions are limited, the accuracy of the model is significantly reduced. In addition, although some studies have attempted to use source location data (such as epicenter location, source depth, magnitude, etc.) in the earthquake catalog to assist in fault identification, the fault activity information (such as source distribution patterns, spatial correlation of earthquake swarms, etc.) implicit in the earthquake catalog has not been systematically extracted and modeled. Existing methods lack efficient algorithms to directly convert massive source data into three-dimensional fault planes, resulting in the inability to independently construct a complete fault model.

[0005] Earthquake catalogs contain hidden information about active faults. Mining fault information from these catalogs can aid in earthquake disaster research and provide a scientific basis for decision-making in major project site selection, land spatial planning, and the management of major earthquakes and disasters. Therefore, there is an urgent need to develop a new workflow that can directly construct three-dimensional fault planes based on a large amount of seismic source location data. Summary of the Invention

[0006] The purpose of this invention is to provide a method and apparatus for constructing three-dimensional fault planes based on an earthquake catalog. For the first time, it combines the OPTICS (Ordering points to identify the clustering structure) method, density peak clustering (DPC) algorithm, and principal component analysis (PCA) to obtain three-dimensional fault planes solely based on the epicenter location, without relying on other information. This invention can process data fully automatically, reducing inaccuracies caused by manual intervention and offering greater convenience.

[0007] To achieve the above objectives, the present invention provides the following solution:

[0008] A method for constructing three-dimensional fault planes based on a seismic catalog, the method comprising the following steps:

[0009] S1, Data Preprocessing: Based on the earthquake catalog of the target area, extract the source location information, obtain the density information of each data point, and generate a density map; where the source location information includes longitude coordinates, latitude coordinates, and depth coordinates in the UTM coordinate system;

[0010] S2, Preliminary classification using the OPTICS algorithm: Based on density information and density map, determine the input parameters of the OPTICS algorithm, use the OPTICS algorithm to perform preliminary cluster analysis on the source location information, output the accessibility map, and determine the number of clusters N based on the valley-peak structure in the accessibility map, where N is an integer greater than 1;

[0011] S3, Clustering optimization based on DPC algorithm: Calculate the Euclidean distance between all source points, calculate the local density and center offset distance of each source point based on the selected cutoff distance, and use them as the horizontal and vertical axes respectively to generate a decision map. Select the top N source points with the largest product of local density and center offset distance as cluster centers, and assign the remaining points to the nearest neighbor cluster with greater local density, and output the DPC classification results.

[0012] S4, PCA analysis: Perform PCA analysis on the DPC classification results to extract feature values ​​λ1, λ2, λ3 and corresponding feature vectors, where λ1 > λ2 > λ3. λ1, λ2, and λ3 represent the feature values ​​of the first, second, and third principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system.

[0013] S5, PCA plane determination: If the set conditions are met: λ1>λ2>e and λ2>>λ3, then a fault plane is formed based on the DPC classification result, where e is the positioning error threshold;

[0014] S6, obtain the planar features of the fault plane.

[0015] Furthermore, in S2, the input parameters of the OPTICS algorithm include the neighborhood radius ε and the minimum number of points M;

[0016] The minimum number of points M is determined based on the density information obtained from S1: half of the maximum density value in the density map is selected as the value of the minimum number of points M.

[0017] The neighborhood radius ε is set to infinity.

[0018] Furthermore, in S3, the method for selecting the cutoff distance is as follows:

[0019] Calculate the local density information entropy at different cutoff distances;

[0020] The cutoff distance that minimizes the local density information entropy is chosen as the DPC clustering parameter.

[0021] Furthermore, S3, the clustering optimization based on the DPC algorithm, also includes:

[0022] Using marginal density, the epicenter is divided into core points and halo points, with halo points considered as noise.

[0023] In the DPC algorithm, each cluster has a marginal region where all points are less than the cutoff distance from other clusters. Within this marginal region, the highest density value is the marginal density of the cluster. When the local density of a source point is greater than the marginal density, it is considered the core point of the cluster; otherwise, it is considered a halo point.

[0024] Furthermore, in S3, the local density ρ i The calculation formula is as follows:

[0025]

[0026] Where, d ij d represents the distance between focal point i and focal point j in the earthquake catalog. c Represents the cutoff distance; ρ i This indicates the cutoff distance d c The set of all earthquake points within the region; χ(x) represents the unit step function;

[0027] Alternatively, the local density ρ can be calculated using the Gaussian kernel density function. i :

[0028]

[0029] Where, d c Used to control the weight degradation rate, ρi The local density of earthquake focal points in the earthquake catalog is a weighted sum of the distances between all focal points.

[0030] Furthermore, in S3, the center offset distance δ i The calculation formula is as follows:

[0031]

[0032] In the formula, the center offset distance δ i This represents the distance from the epicenter i to its nearest higher density point.

[0033] Furthermore, step S5, PCA plane determination, also includes:

[0034] If the set conditions are not met, it means that the fault group cannot form a fault plane and is composed of multiple fault planes. Then repeat steps S2-S4 to re-cluster the fault group using the OPTICS and DPC algorithms, and re-perform PCA analysis until each cluster meets the set conditions.

[0035] Furthermore, in S6, the planar features of the fault plane include: the dip angle and direction of the fault plane, its length and width, the x-axis coordinate of the centroid, the y-axis coordinate of the centroid, and the z-axis coordinate of the centroid;

[0036] The method for obtaining the planar features of the fault plane specifically includes:

[0037] Among the three eigenvalues ​​λ1, λ2, λ3 and their corresponding eigenvectors u1, u2, u3, the fault plane is drawn based on the two larger eigenvalues ​​λ1, λ2 and the eigenvectors u1, u2, and the smallest eigenvector λ3 is perpendicular to the fault plane.

[0038] Among them, the largest eigenvalue λ1 determines the length information of the cluster, the eigenvector u1 determines the direction of the length, the second largest eigenvalue λ2 determines the width information of the cluster, and the eigenvector u2 determines the direction of the width, thus obtaining the length and width of the cross-section.

[0039] The minimum eigenvalue λ3 and the corresponding eigenvector u3 determine the thickness information of the cluster;

[0040] The plane corresponding to the fault plane passes through the data center μ D =(x d y d , z d According to the eigenvector u3=(u 3，1 u 3，2 u 3，3 Construct the plane equations:

[0041] u 3，1 xd +u 3，2 y d +u 3，3 z d =O (6)

[0042] Among them, u 31 u 32 and u 33 These are the components of the eigenvector u3, μ D Let x represent the centroid of the plane. d y d and z d These represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the centroid, respectively.

[0043] The strike and dip angle of the fault plane are calculated based on the eigenvector u3.

[0044] Furthermore, the method also includes:

[0045] S7, Second-order fault segment analysis: For the seismic source points within the fault plane formed in step S5, the DPC algorithm is used for secondary classification to obtain multiple second-order sub-clustering results. PCA analysis and PCA plane determination of step VS5 are performed on each second-order sub-clustering result. If the set conditions are met, the fault segments that make up the fault plane are obtained based on the second-order sub-clustering results.

[0046] The present invention also provides an apparatus for constructing a three-dimensional fault plane based on a seismic catalog, applied to performing the above-described method for constructing a three-dimensional fault plane based on a seismic catalog, comprising:

[0047] The data preprocessing module is used to extract the source location information, obtain the density information of each data point, and generate a density map based on the earthquake catalog of the target area; the source location information includes longitude coordinates, latitude coordinates, and depth coordinates in the UTM coordinate system.

[0048] The OPTICS algorithm preliminary classification module is used to determine the input parameters of the OPTICS algorithm based on density information and density map, perform preliminary cluster analysis on the source location information using the OPTICS algorithm, output an accessibility map, and determine the number of clusters N based on the valley-peak structure in the accessibility map, where N is an integer greater than 1.

[0049] The DPC algorithm clustering module is used to calculate the Euclidean distance between all source points. Based on the selected cutoff distance, it calculates the local density and center offset distance of each source point and uses them as the horizontal and vertical axes respectively to generate a decision map. The top N source points with the largest product of local density and center offset distance are selected as cluster centers, and the remaining points are assigned to the nearest neighbor cluster with the larger local density. The DPC classification results are then output.

[0050] The PCA analysis module is used to perform PCA analysis on the DPC classification results and extract feature values ​​λ1, λ2, λ3 and their corresponding feature vectors. Here, λ1 > λ2 > λ3, and λ1, λ2, and λ3 represent the feature values ​​of the first, second, and third principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system.

[0051] The PCA plane determination module is used to determine whether the set conditions are met: λ1>λ2>e and λ2>>λ3. If they are met, a fault plane is formed based on the DPC classification result, where e is the positioning error threshold.

[0052] The feature acquisition module is used to acquire planar features of the fault plane.

[0053] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed, implements the method described above for constructing three-dimensional fault planes based on an earthquake catalog.

[0054] According to specific embodiments provided by the present invention, the following technical effects are disclosed: The method and apparatus for constructing three-dimensional fault planes based on an earthquake catalog provided by the present invention, for the first time, combine the OPTICS algorithm (accessibility clustering), the DPC algorithm (density peak clustering), and PCA (principal component analysis). By mining fault information through the earthquake catalog, that is, identifying fault planes solely based on epicenter location information, it breaks through the dependence of traditional modeling on multi-source data. In particular, by introducing the accessibility map of the OPTICS algorithm to assist the selection of cluster centers in the DPC algorithm, the problem of the high sensitivity of the DPC algorithm to human intervention is solved, human error is reduced, and the clustering results are ensured to be closer to the actual fault distribution, significantly improving the clustering accuracy and robustness. Then, PCA analysis is applied to each cluster to construct fault planes. PCA analysis can accurately extract the spatial features of each cluster (such as the strike, dip, dip angle, and geometric dimensions of the fault plane), making the fault plane modeling closer to the real geological structure.

[0055] This invention requires only basic earthquake catalog information such as the location of the seismic source to construct a three-dimensional fault plane model in areas lacking surface surveys or geophysical data, greatly expanding the scope of application of the technology. It is suitable for fault activity research in earthquake-prone areas around the world, and provides a low-cost, high-efficiency fault modeling solution, especially for areas with scarce data. Attached Figure Description

[0056] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0057] Figure 1 This is a flowchart of the method for constructing three-dimensional fault planes based on an earthquake catalog according to the present invention;

[0058] Figure 2 This is a three-dimensional display effect diagram of the earthquake source location information of the present invention;

[0059] Figure 3 This is an accessibility diagram of the present invention;

[0060] Figure 4 This is the decision graph for the present invention, wherein (a) shows the selected point, and (b) shows ρ. i With δ i Data points with values ​​greater than the selected point,

[0061] Figure 5 In a specific embodiment of the present invention, the clustering results obtained by the DPC algorithm are shown, wherein (a) is a decision graph with four cluster centers selected; and (b) is the four clustering results based on the cluster centers.

[0062] Figure 6 This is a schematic diagram of the fracture surface obtained by performing PCA analysis on four clustering results in a specific embodiment of the present invention. Detailed Implementation

[0063] The embodiments of the present invention are described in detail below. These embodiments are intended to explain the present invention and should not be construed as limiting the present invention. Where specific techniques or conditions are not specified in the embodiments, they are performed according to the techniques or conditions described in the literature in the art or according to the product instructions. Materials or instruments whose manufacturers are not specified are all conventional products that can be obtained commercially.

[0064] This invention combines the OPTICS algorithm with the DPC algorithm. The OPTICS algorithm is used to assist in selecting cluster points in the decision graph of the DPC algorithm, making the obtained clusters more accurate. Then, PCA analysis is applied to each cluster to construct the cross-sectional plane. This method consists of three parts: First, the OPTICS algorithm is used to obtain the number of clusters. The OPTICS algorithm requires two input parameters: a minimum number of clusters and a neighborhood radius. It also generates an accessibility map. Second, the DPC algorithm is used to classify the source location information, providing quantitative information about them. The DPC algorithm incorporates local density entropy to assist in selecting the cutoff distance. It generates a decision map, requiring manual selection of cluster centers. During this process, the accessibility map obtained from the OPTICS algorithm is used as an aid to more accurately select cluster centers. Third, each cluster identified by DPC undergoes PCA analysis. PCA analysis obtains the eigenvalues ​​and eigenvectors of each cluster. Based on the relationships between the eigenvalues, each cluster is interpreted to obtain the plane formed by the cluster, including its spatial position, orientation, length, and height. By applying PCA analysis, the planar features obtained from the clusters can be interpreted as fault planes. This invention is the first to combine the OPTICS algorithm, PCA analysis, and DPC algorithm to construct a workflow for three-dimensional fault plane determination. It can obtain the three-dimensional fault plane solely based on the seismic source location, without relying on other information. This invention can process data fully automatically, reducing inaccuracies caused by manual intervention and offering greater convenience.

[0065] like Figure 1 As shown, the method for constructing three-dimensional fault planes based on an earthquake catalog provided by the present invention includes the following steps:

[0066] S1, Data Preprocessing: Based on the earthquake catalog of the target area, extract the focal location information, obtain the density information of each data point, and generate a density map; the focal location information includes longitude, latitude, and depth coordinates in the UTM coordinate system; the earthquake catalog is processed to remove unnecessary information, retaining only information related to earthquake location. The three-dimensional display of the focal location information is shown below. Figure 2 As shown, it is only necessary to collect earthquake location information for a certain area, obtain the density information of each data point, and display it in three dimensions. Figure 2 The different colors on the right-hand side of the color bar represent different density information.

[0067] S2, Preliminary classification using the OPTICS algorithm: Based on density information and density map, determine the input parameters of the OPTICS algorithm, use the OPTICS algorithm to perform preliminary cluster analysis on the source location information, output an accessibility map, and determine the number of clusters N based on the valley-peak structure in the accessibility map, where N is an integer greater than 1; the input parameters of the OPTICS algorithm include the neighborhood radius ε and the minimum number of points M;

[0068] The minimum number of points M is determined based on the density information obtained from S1: half of the maximum density value in the density map is selected as the minimum number of points M; for example, based on... Figure 2 The value of the minimum dot number M is 225, which is half of the maximum value of the medium color stick, which is 450.

[0069] This algorithm is insensitive to the neighborhood radius input (generally ε is fixed to infinity), and can obtain an reachability graph, such as... Figure 3 As shown, for example, Figure 3 The maximum neighborhood radius is 1.6. The output of the OPTICS algorithm is a graph called an reachability graph, which visualizes the reachability distance of each point in the dataset. Points belonging to a cluster have low reachability distances to their nearest neighbors and are identified as valleys in the reachability graph. The deeper the valley, the denser the cluster. Valleys are separated by peaks, which are points in the dataset that are relatively far from their neighbors. In the reachability graph, a point on the vertical axis can be selected as a threshold to set a horizontal line. The value of the horizontal line is the value of the neighborhood radius. This line cuts the data; for example, if it intersects with a peak at four points, it can be divided into four clusters.

[0070] S3, Clustering optimization based on DPC algorithm: Calculate the Euclidean distance between all source points, calculate the local density and center offset distance of each source point based on the selected cutoff distance, and use them as the horizontal and vertical axes respectively to generate a decision map. Select the top N source points with the largest product of local density and center offset distance as cluster centers, and assign the remaining points to the nearest neighbor cluster with greater local density, and output the DPC classification results.

[0071] The method for selecting the cutoff distance is as follows:

[0072] Calculate the local density information entropy at different cutoff distances:

[0073]

[0074] In the formula, ρ represents the density of all data points i The sum of the local function information entropy of the data points. If the local function information entropy of the data points is minimized, then the local density difference of the data points is maximized, which is beneficial for determining the cluster center. Then a cutoff distance d can be found. c This minimizes the local density information entropy, at which point dc For the final cutoff distance:

[0075]

[0076] Local density ρ i The formula for calculating the x-axis of the decision graph is as follows:

[0077]

[0078] Where, d ij d represents the distance between focal point i and focal point j in the earthquake catalog. c Represents the cutoff distance; ρ i This indicates the cutoff distance d c The set of all earthquake points within the region; χ(x) represents the unit step function;

[0079] Alternatively, the local density ρ can be calculated using the Gaussian kernel density function. i :

[0080]

[0081] Where, d c Used to control the weight degradation rate, ρ i The local density of earthquake focal points in the earthquake catalog is a weighted sum of the distances between all focal points.

[0082] In the specific implementation process, first determine the cutoff distance d. c A range of values ​​is given, within which a test d is first provided. c After calculating the local density using formula (3) or (4), the calculated local density is substituted into (1) to obtain an H(d) value. c ), and then d c Within the range of values, repeat this calculation process to obtain multiple H(d) values. c The value of ) is used to find the minimum local density information entropy H(d) according to formula (2). c ), the corresponding d c The value is the final cutoff distance, which is then substituted into formula (3) or (4) to obtain the learned local density ρ. i .

[0083] Center offset distance δ i The formula for calculating the vertical axis of the decision diagram is as follows:

[0084]

[0085] In the formula, the center offset distance δ i This represents the distance from the epicenter i to its nearest higher density point.

[0086] A decision map is drawn with local density on the horizontal axis and center offset distance on the vertical axis, as shown below. Figure 4 As shown.

[0087] Select a point on the decision map, typically choosing a large local density and center offset distance, and then automatically obtain all ρ values. i and δ i Points whose values ​​are all greater than the selected point. In this invention, based on the reachability graph obtained in step S2, it is understood how many clusters the data can be divided into, for example, it can be divided into 4 clusters. Then, in the decision graph, the value ρ of each point is calculated. i ×δ i The values ​​are used to determine the cluster centers, with the four largest values ​​being selected as the cluster values. The number of points obtained represents the number of clusters. Figure 5 As shown in (a), the data is divided into four categories. This step is the first to utilize OPTICS to assist DPC in selecting cluster points, thus improving clustering accuracy.

[0088] Assign the remaining points to the nearest neighbor clusters whose central local density is greater than theirs, and obtain the final DPC classification result, such as... Figure 5 As shown in (b).

[0089] In addition, to identify anomalies in the earthquake catalog, marginal density ρ can be used. b The seismic catalog is divided into core points and clusters called HALO. Halo points can be considered as noise. In the original DPC algorithm, each cluster has a marginal region where all points are less than the cutoff distance dc from other clusters. Within this marginal region, the highest density value is the marginal density ρ of that cluster. b When ρ b <ρ i If the data point is true, it is considered the core point of the cluster; otherwise, it is considered a HALO cluster.

[0090] S4, PCA Analysis: Perform PCA analysis on the DPC classification results to extract eigenvalues ​​λ1, λ2, λ3 and their corresponding eigenvectors, where λ1 > λ2 > λ3. λ1, λ2, and λ3 represent the eigenvalues ​​of the first (PC1), second (PC2), and third (PC3) principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system. The three eigenvalues ​​represent the lengths of the coordinate axes, and the corresponding eigenvectors represent the directions of the coordinate axes. λ1 (the largest eigenvalue) represents the dispersion of seismic events along the main extension direction on the fault plane (such as fault strike or dip), λ2 (the second eigenvalue) represents the dispersion of the secondary extension direction within the fault plane (the direction perpendicular to PC1), and λ3 (the smallest eigenvalue) represents the dispersion of the fault plane normal direction. In an ideal plane, λ3 should approach 0; in reality, due to positioning errors, λ3 reflects the fault thickness or data noise.

[0091] λ1 and λ2 quantify the distribution intensity within the cross-section (forming a planar coordinate system); λ3 quantifies the degree of deviation from the ideal plane (the smaller the value, the flatter the surface).

[0092] S5, PCA plane determination: If the set conditions are met: λ1>λ2>e and λ2>>λ3, then a fault plane is formed based on the DPC classification result, where e is the positioning error threshold; if the set conditions are not met, it means that the fault group cannot form a fault plane, but is composed of multiple fault planes; then repeat steps S2-S4, re-apply the OPTICS algorithm and DPC algorithm to cluster the fault group, and re-perform PCA analysis until each cluster meets the set conditions.

[0093] S6, obtain the planar features of the fault plane, including: the dip angle and strike, length and width of the fault plane, the x-axis coordinate of the centroid, the y-axis coordinate of the centroid, and the z-axis coordinate of the centroid;

[0094] The method for obtaining the planar features of the fault plane specifically includes:

[0095] Among the three eigenvalues ​​λ1, λ2, λ3 and their corresponding eigenvectors u1, u2, u3, the fault plane is drawn based on the two larger eigenvalues ​​λ1, λ2 and the eigenvectors u1, u2, and the smallest eigenvector λ3 is perpendicular to the fault plane.

[0096] Among them, the largest eigenvalue λ1 determines the length information of the cluster, the eigenvector u1 determines the direction of the length, the second largest eigenvalue λ2 determines the width information of the cluster, and the eigenvector u2 determines the direction of the width, thus obtaining the length and width of the cross-section.

[0097] The minimum eigenvalue λ3 and the corresponding eigenvector u3 determine the thickness information of the cluster;

[0098] The plane corresponding to the fault plane passes through the data center μ D =(x d y d , z d According to the eigenvector u3=(u 3，1 u 3，2 u 3，3 Construct the plane equations:

[0099] u 3，1 x d +u 3,2 y d +u 3，3 z d =O (6)

[0100] Among them, u 31 u 32 and u 33 These are the components of the eigenvector u3, μ D Let x represent the centroid of the plane. d y d and z d These represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the centroid, respectively.

[0101] The strike and dip angle of the fault plane are calculated based on the eigenvector u3.

[0102] According to the method of the present invention, if the two largest feature values ​​λ1 and λ2 are greater than the positioning error threshold e extracted directly from the dataset, and are much greater than the third feature value λ3, then the clustering geometrically defines a planar feature. The present invention can automatically analyze the four clusters obtained in step S3, such as... Figure 6 As shown, three of them satisfy the planar conditions, automatically forming a plane. This plane can be considered a fault plane, and the parameters of the three fault planes are automatically obtained, as shown in Table 1; while Figure 6 In the middle, the area where green dots are concentrated does not meet the planar conditions, therefore, it cannot form a plane.

[0103] Table 1

[0104]

[0105] In addition, the method also includes:

[0106] S7, Second-order fault segment analysis: For the seismic source points within the fault plane formed in step S5, the DPC algorithm is used for secondary classification to obtain multiple second-order sub-clustering results. PCA analysis and PCA plane determination of step S5 are performed on each second-order sub-clustering result. If the set conditions are met, the fault segments that make up the fault plane are obtained based on the second-order sub-clustering results.

[0107] Specifically, regarding PCA analysis, it's important to understand that: Principal Component Analysis transforms the original data through a linear transformation into a new set of features, called "principal components." Principal components primarily rely on two mathematical tools: the covariance matrix and eigenvalue decomposition. The covariance matrix reveals the correlations between the features. In other words, it creates a new coordinate basis for the original data. Since the dataset X used in this study contains three variables, the covariance matrix is ​​a 3x3 matrix.

[0108]

[0109] Eigenvalue decomposition of the covariance matrix yields a set of three eigenvalues ​​and corresponding eigenvectors.

[0110] Cu i =λ i u i , i = 1, 2, 3 (8)

[0111] Where, λ i Let u represent the eigenvalues, satisfying λ1>λ2>λ3. i This represents the eigenvector. The magnitude of the eigenvalues ​​corresponds to the variance of the principal components, while the eigenvectors indicate the direction of the principal components. Note that in formula (7), cov(x) i x i ), i = 1, 2, 3 represent the variances of different variables. And it is equivalent to the eigenvalues ​​of PCA:

[0112]

[0113] Assuming the fault is planar and the location uncertainty of earthquakes is uniform throughout the earthquake catalog, the fault plane is drawn based on the two largest eigenvalues ​​λ1, λ2, λ3 and their corresponding eigenvectors u1, u2, u3. The smallest eigenvector λ3 is perpendicular to the fault plane. The largest eigenvalue λ1 determines the length information of the cluster, and the eigenvector u1 determines the direction of the length. The second largest eigenvalue λ2 determines the width information of the cluster, and the eigenvector u2 determines the direction of the width, thus obtaining the length and width of the fault plane. Therefore, the length L and width W of the plane corresponding to this fault plane can be expressed as:

[0114]

[0115] The minimum eigenvalue λ3 and the corresponding eigenvector u3 determine the thickness information of the cluster.

[0116] The plane corresponding to the fault plane passes through the data center μ D =(x d y d , zd According to the eigenvector u3=(u 3，1 u 3,2 u 3，3 Construct the plane equations:

[0117] u 3，1 x d +u 3，2 y d +u 3，3 z d =O (11)

[0118] Among them, u 31 u 32 and u 33 These are the components of the eigenvector u3, μ D Let x represent the centroid of the plane. d y d and z d These represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the centroid, respectively.

[0119] Calculate the strike and dip angle of the fault plane based on the eigenvector u3:

[0120]

[0121] Here, strike and dip represent the orientation and angle of inclination of the plane, respectively.

[0122] The specific process of PCA analysis is as follows:

[0123] Based on the clustering results from the earthquake catalog, it is considered a fault group;

[0124] Principal component analysis was performed on each fault group to obtain the corresponding eigenvalues ​​and eigenvectors;

[0125] Determine whether each fault group satisfies λ1>λ2>e. If not, the analysis ends and the final result is obtained. If yes, λ1>λ2>e, then the calculation continues.

[0126] Determine whether λ2 >> λ3 is satisfied. If satisfied, calculate the length and width of the constructed fault plane according to formula (10).

[0127] Calculate the centroid, strike, and dip angle of the fault plane according to formulas (11) and (12);

[0128] If λ2 >> λ3 is not satisfied, it means that the fault group cannot form a fault plane and may be composed of multiple fault planes. In this case, the above steps need to be repeated to re-cluster the fault group and perform PCA analysis again.

[0129] The analysis ends when each cluster satisfies λ1>λ2>e and λ2>>λ3, yielding the final fault plane and fault segment results.

[0130] This invention proposes a process for processing earthquake catalogs based on the OPTICS algorithm, DCP algorithm, and PCA method. It can obtain fault planes and fault segments that make up fault planes by relying solely on the earthquake catalog, and can also obtain the corresponding plane parameters, thereby constructing a three-dimensional fault model.

[0131] Example 2

[0132] The present invention also provides an apparatus for constructing a three-dimensional fault plane based on a seismic catalog, applied to performing the above-described method for constructing a three-dimensional fault plane based on a seismic catalog, comprising:

[0133] The data preprocessing module is used to extract the source location information, obtain the density information of each data point, and generate a density map based on the earthquake catalog of the target area; the source location information includes longitude coordinates, latitude coordinates, and depth coordinates in the UTM coordinate system.

[0134] The OPTICS algorithm preliminary classification module is used to determine the input parameters of the OPTICS algorithm based on density information and density map, perform preliminary cluster analysis on the source location information using the OPTICS algorithm, output an accessibility map, and determine the number of clusters N based on the valley-peak structure in the accessibility map, where N is an integer greater than 1.

[0135] The DPC algorithm clustering module is used to calculate the Euclidean distance between all source points. Based on the selected cutoff distance, it calculates the local density and center offset distance of each source point and uses them as the horizontal and vertical axes respectively to generate a decision map. The top N source points with the largest product of local density and center offset distance are selected as cluster centers, and the remaining points are assigned to the nearest neighbor cluster with the larger local density. The DPC classification results are then output.

[0136] The PCA analysis module is used to perform PCA analysis on the DPC classification results and extract feature values ​​λ1, λ2, λ3 and their corresponding feature vectors. Here, λ1 > λ2 > λ3, and λ1, λ2, and λ3 represent the feature values ​​of the first, second, and third principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system.

[0137] The PCA plane determination module is used to determine whether the set conditions are met: λ1>λ2>e and λ2>>λ3. If they are met, a fault plane is formed based on the DPC classification result, where e is the positioning error threshold.

[0138] The feature acquisition module is used to acquire planar features of the fault plane.

[0139] Example 3

[0140] The present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed, implements the method described above for constructing three-dimensional fault planes based on an earthquake catalog.

[0141] Matters not covered in this invention are common knowledge.

[0142] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0143] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or improvement of the technology in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

Claims

1. A method for constructing three-dimensional fault planes based on an earthquake catalog, characterized in that, Includes the following steps: S1, Data Preprocessing: Based on the earthquake catalog of the target area, extract the source location information, obtain the density information of each data point, and generate a density map; where the source location information includes longitude coordinates, latitude coordinates, and depth coordinates in the UTM coordinate system; S2. Preliminary classification using the OPTICS algorithm: Based on density information and density maps, the input parameters of the OPTICS algorithm are determined. The OPTICS algorithm is then used to perform preliminary cluster analysis on the source location information, outputting an accessibility map. The number of clusters N is determined based on the valley-peak structure in the accessibility map, where N is an integer greater than 1. S3, Clustering Optimization Based on DPC Algorithm: Calculate the pairwise Euclidean distances between all source points. Based on the selected cutoff distance, calculate the local density and center offset distance of each source point, and use them as the horizontal and vertical axes respectively to generate a decision map. Select the N source points with the largest product of local density and center offset distance as cluster centers, and assign the remaining points to the nearest neighbor clusters with higher local density. Output the DPC classification results. S4, PCA analysis: Perform PCA analysis on the DPC classification results to extract feature values ​​λ1, λ2, λ3 and corresponding feature vectors, where λ1 > λ2 > λ3. λ1, λ2, and λ3 represent the feature values ​​of the first, second, and third principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system. S5, PCA plane determination: If the set conditions are met: λ1>λ2>e and λ2>>λ3, then a fault plane is formed based on the DPC classification result, where e is the positioning error threshold; S6, obtain the planar features of the fault plane.

2. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, In S2, the input parameters of the OPTICS algorithm include the neighborhood radius ε and the minimum number of points M; The minimum number of points M is determined based on the density information obtained from S1: half of the maximum density value in the density map is selected as the value of the minimum number of points M. The neighborhood radius ε is set to infinity.

3. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, In S3, the cutoff distance is selected as follows: Calculate the local density information entropy at different cutoff distances; The cutoff distance that minimizes the local density information entropy is chosen as the DPC clustering parameter.

4. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, The S3 method, which optimizes clustering based on the DPC algorithm, further includes: Using marginal density, the source points are divided into core points and halo points, with halo points being considered as noise; In the DPC algorithm, each cluster has a marginal region where all points are less than the cutoff distance from other clusters. Within this marginal region, the highest density value is the marginal density of the cluster. When the local density of a source point is greater than the marginal density, it is considered the core point of the cluster; otherwise, it is considered a halo point.

5. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, In S3, the local density ρ i The calculation formula is as follows: Where, d ij d represents the distance between focal point i and focal point j in the earthquake catalog. c Represents the cutoff distance; ρ i This indicates the cutoff distance d c The set of all earthquake points within the region; χ(x) represents the unit step function; Alternatively, the local density ρ can be calculated using the Gaussian kernel density function. i : Where, d c Used to control the weight degradation rate, ρ i The local density of earthquake focal points in the earthquake catalog is a weighted sum of the distances between all focal points.

6. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, In S3, the center offset distance δ i The calculation formula is as follows: In the formula, the center offset distance δ i This represents the distance from the epicenter i to its nearest higher density point.

7. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, Step S5, PCA plane determination, further includes: If the set conditions are not met, it means that the fault group cannot form a fault plane and is composed of multiple fault planes. Then repeat steps S2-S4 to re-cluster the fault group using the OPTICS and DPC algorithms, and re-perform PCA analysis until each cluster meets the set conditions.

8. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, In S6, the planar features of the fault plane include: the dip angle and strike, length and width of the fault plane, the x-axis coordinate of the centroid, the y-axis coordinate of the centroid, and the z-axis coordinate of the centroid. The method for obtaining the planar features of the fault plane specifically includes: Among the three eigenvalues ​​λ1, λ2, and λ3 and their corresponding eigenvectors u1, u2, and u3, the fault plane is drawn based on the two larger eigenvalues ​​λ1 and λ2 and the eigenvectors u1 and u2. The smallest eigenvector λ3 is perpendicular to the fault plane. Among them, the largest eigenvalue λ1 determines the length information of the cluster, the eigenvector u1 determines the direction of the length, the second largest eigenvalue λ2 determines the width information of the cluster, and the eigenvector u2 determines the direction of the width, thus obtaining the length and width of the cross-section. The minimum eigenvalue λ3 and the corresponding eigenvector u3 determine the thickness information of the cluster; The plane corresponding to the fault plane passes through the data center μ D =(x d y d , z d According to the eigenvector u3=(u 3，1 u 3，2 u 3，3 Construct the plane equations: u 3，1 x d +u 3,2 y d +u 3，3 z d =0 (6) Among them, u 31 u 32 and u 33 These are the components of the eigenvector u3, μ D Let x represent the centroid of the plane. d y d and z d These represent the x-axis coordinate, y-axis coordinate, and z-axis coordinate of the centroid, respectively. The strike and dip angle of the fault plane are calculated based on the eigenvector u3.

9. The method for constructing three-dimensional fault planes based on an earthquake catalog according to claim 1, characterized in that, The method further includes: S7, Second-order fault segment analysis: For the seismic source points within the fault plane formed in step S5, the DPC algorithm is used for secondary classification to obtain multiple second-order sub-clustering results. PCA analysis and PCA plane determination of step S5 are performed on each second-order sub-clustering result. If the set conditions are met, the fault segments that make up the fault plane are obtained based on the second-order sub-clustering results.

10. An apparatus for constructing three-dimensional fault planes based on a seismic catalog, applied to the method for constructing three-dimensional fault planes based on a seismic catalog as described in any one of claims 1-9, characterized in that, include: The data preprocessing module is used to extract the source location information, obtain the density information of each data point, and generate a density map based on the earthquake catalog of the target area; the source location information includes longitude coordinates, latitude coordinates, and depth coordinates in the UTM coordinate system. The OPTICS algorithm preliminary classification module is used to determine the input parameters of the OPTICS algorithm based on density information and density map, perform preliminary cluster analysis on the source location information using the OPTICS algorithm, output an accessibility map, and determine the number of clusters N based on the valley-peak structure in the accessibility map, where N is an integer greater than 1. The DPC algorithm clustering module is used to calculate the Euclidean distance between all source points. Based on the selected cutoff distance, it calculates the local density and center offset distance of each source point and uses them as the horizontal and vertical axes respectively to generate a decision map. The top N source points with the largest product of local density and center offset distance are selected as cluster centers, and the remaining points are assigned to the nearest neighbor cluster with the larger local density. The DPC classification results are then output. The PCA analysis module is used to perform PCA analysis on the DPC classification results and extract feature values ​​λ1, λ2, λ3 and their corresponding feature vectors. Here, λ1 > λ2 > λ3, and λ1, λ2, and λ3 represent the feature values ​​of the first, second, and third principal component directions obtained through PCA analysis, respectively. The three principal component directions constitute a new three-dimensional coordinate system. The PCA plane determination module is used to determine whether the set conditions are met: λ1>λ2>e and λ2>>λ3. If they are met, a fault plane is formed based on the DPC classification result, where e is the positioning error threshold. The feature acquisition module is used to acquire planar features of the fault plane.

Images