Geological modeling data spatial interpolation method based on anisotropic Huber robust reweighting

Abstract

The invention discloses a geologic modeling data spatial interpolation method based on anisotropic Huber robust reweighting, and the method comprises the steps: obtaining stratigraphic scatter data, obtaining a stratigraphic elevation gradient through local neighborhood construction and gradient estimation, expanding the stratigraphic elevation gradient into a continuous gradient field, and constructing a trend-inclination direction field; selecting neighborhood sample points of the query points, projecting the neighborhood sample points to a constructed coordinate system, and defining an anisotropic distance to construct a space kernel weight coefficient; after initial weighted least square fitting, a Huber robust reweighting strategy is introduced, the comprehensive weight is optimized through residual segmentation processing, and finally an interpolation result is obtained through refitting. The method breaks through the traditional isotropic limitation, effectively inhibits outlier and noise interference, does not need an artificial prior geological model, maintains high stability and high precision in a complex scene, and provides reliable data support for geological modeling.

Description

Technical Field

[0001] This invention relates to the field of geological modeling data interpolation technology, and in particular to a spatial interpolation method for geological modeling data based on anisotropic Huber robust reweighting. Background Technology

[0002] As oil and gas exploration and development continue to advance into deeper and more complex structural areas, the structural morphology of target subsurface layers is becoming increasingly complex. High-precision subsurface structural geological models are crucial for the deployment and decision-making of exploration and development plans. However, raw stratigraphic data is often affected by factors such as seismic acquisition conditions and interpretation accuracy, resulting in problems such as data gaps and sparse grid density, making it difficult to meet the needs of high-precision geological modeling. Furthermore, stratigraphic structural changes exhibit significant directional characteristics, typically changing gently along the strike direction and rapidly along the dip direction. Additionally, raw data often contains noise points and outliers. These factors pose a severe challenge to the accuracy of spatial interpolation of stratigraphic data, necessitating a high-precision interpolation method that can adapt to geological structural characteristics and effectively suppress interference from anomalous data.

[0003] Existing traditional interpolation methods have two major flaws: First, they assume spatial isotropy by default, assigning equal weight to stratigraphic variations in different directions, failing to fully consider the differentiated variation patterns of strata along strike and dip directions. This results in interpolation results that cannot accurately reproduce the actual morphology of underground structures, thereby reducing the accuracy of geological modeling. Second, they do not design effective suppression mechanisms for noise and outliers during data processing. When using methods such as ordinary weighted least squares for fitting, extreme residuals can significantly distort the position of the fitting plane, making the interpolation results severely affected by abnormal data. This makes it difficult to maintain stable interpolation accuracy under complex data conditions and cannot provide reliable data support for high-precision geological modeling. Summary of the Invention

[0004] To overcome the shortcomings and deficiencies of existing technologies, this invention provides a spatial interpolation method for geological modeling data based on anisotropic Huber robust reweighting.

[0005] The technical solution adopted in this invention is a spatial interpolation method for geological modeling data based on anisotropic Huber robust reweighting, comprising the following steps: S1, obtaining a scatter dataset of stratigraphic points in the target structural region, which includes information on the correlation between planar coordinates and stratigraphic elevation; S2, constructing a local neighborhood for each original sample point, selecting several nearest neighboring points to form a neighborhood set, assuming that the stratigraphy is approximately planar within a small range, constructing a design matrix and solving for the planar parameters using least squares, and using the correlated components in the planar parameters as the stratigraphic elevation gradient at that sample point; S3, performing scatter interpolation on the discrete gradient data, constructing a continuous gradient function, obtaining the gradient parameters at any query point using this function, and defining the stratigraphic maximum... S4. The descending direction is used as the plane projection of the dip direction, and an orthogonal strike direction vector is constructed; S5. For any query point to be interpolated, several nearest original sample points are selected to form a neighborhood set, and the plane coordinate difference between the neighborhood sample points and the query point is calculated. This difference vector is projected onto the strike and dip directions; S6. The strike association length and dip association length are introduced to define the anisotropic distance. Based on this distance, a spatial kernel weight coefficient is constructed to express the characteristics of slower weight decay along the strike direction and faster weight decay along the dip direction; S7. Assuming that the strata in the neighborhood of the query point are approximately a local plane, a design matrix and observation vector are constructed. A diagonal weight matrix is ​​constructed based on the spatial kernel weight coefficient, and the plane parameter estimation is solved by the weighted least squares method.

[0006] Furthermore, the model for solving the plane parameters using least squares in S2 is as follows: ,in, The fitted plane parameter column vector includes the rate of change in the x-direction, the rate of change in the y-direction, and a plane constant term. To design the transpose of the matrix, The design matrix for local plane fitting. It is the column vector of elevations of the neighborhood sample points.

[0007] Furthermore, the calculation model for the unit vector of the tendency direction at the query point in S3 is as follows: ,in, The unit vector is the direction of inclination. The gradient in the x-direction at the query point. The gradient in the y-direction at the query point is normalized by the square root of the sum of squares of the gradient parameters in the denominator.

[0008] Furthermore, the calculation model for the anisotropic distance in S5 is as follows: ,in, The anisotropic distance from the neighboring sample points to the query point. This represents the component of the planar difference between neighboring sample points and the query point along the direction. The component along the direction of inclination, The length is associated with the direction of travel. The length associated with the direction of inclination.

[0009] Furthermore, the construction model for the spatial weighting coefficients in S5 is as follows: ,in, For spatial anisotropy weights, For anisotropic distances, This is a bandwidth parameter used to adjust the rate at which the weight decays with distance.

[0010] Furthermore, the model for solving the plane parameters using weighted least squares in S6 is as follows: ,in, For the parameter estimation of the 0th plane fitting, Let be the coordinate matrix of the neighborhood sample points. Let be the first spatial weight diagonal matrix. Let the elevation column vector of the neighborhood sample points be denoted as . It is the transpose of the neighborhood sample point coordinate matrix.

[0011] Further, S2 includes the following sub-steps: S21, selecting a predetermined number of neighboring points closest to the target sample point from all sample points based on Euclidean distance, and constructing a local neighborhood set for the sample point; S22, assuming that the strata are distributed in a planar manner within the local neighborhood, constructing a design matrix including the planar coordinates of the neighboring points and constant terms, and simultaneously organizing the elevation data corresponding to the neighboring points to form an elevation column vector; S23, performing matrix multiplication on the design matrix and its transpose, inverting it, and then performing matrix multiplication sequentially with the transpose of the design matrix and the elevation column vector to obtain a planar parameter column vector; S24, extracting the first two components from the planar parameter column vector as the gradients of the strata elevation at the target sample point in the x and y directions.

[0012] Further, step S3 includes the following sub-steps: S31, performing scatter interpolation on the planar coordinates of the original sample points and the corresponding x-direction gradient and y-direction gradient respectively to construct two continuous gradient functions to characterize the gradient parameters at any planar position; S32, substituting the planar coordinates of the query point into the two continuous gradient functions to obtain the x-direction gradient and y-direction gradient corresponding to the query point; S33, obtaining the unit vector of the maximum descent direction of the stratum, i.e., the dip direction, by inverting the gradient parameters and dividing by the square root of the sum of the squares of the gradient parameters; S34, constructing a strike direction unit vector orthogonal to the dip direction by exchanging the components and changing the sign of one of the components based on the components of the dip direction unit vector.

[0013] Further, step S4 includes the following sub-steps: S41, according to the Euclidean distance metric, a preset number of sample points closest to the query point are selected from the original scattered data to form a local neighborhood set of the query point; S42, the difference between the planar coordinates of each sample point in the neighborhood set and the planar coordinates of the query point is calculated one by one to obtain the x-direction coordinate difference and y-direction coordinate difference of each neighboring point relative to the query point; S43, the coordinate difference vector of each neighboring point is multiplied by the unit vector of the direction of movement to obtain the component of the neighboring point relative to the query point along the direction of movement; S44, the coordinate difference vector of each neighboring point is multiplied by the unit vector of the direction of movement to obtain the component of the neighboring point relative to the query point along the direction of movement.

[0014] Further, S5 includes the following sub-steps: S51, setting the orientation association length representing the extension characteristics of the orientation direction and the tendency association length representing the change characteristics of the tendency direction, and clarifying the difference in length parameters between the two directions; S52, based on the components of the neighboring point along the orientation direction and the tendency direction, performing ratio calculations with the corresponding association lengths respectively, and then squaring the two ratio results; S53, summing the two squared results to obtain the squared value of the anisotropic distance from the neighboring point to the query point; S54, substituting the squared value of the anisotropic distance into the exponential function, and controlling the attenuation rate by setting the bandwidth parameter to construct the spatial weight coefficient of each neighboring point.

[0015] Beneficial Effects: This invention proposes a spatial interpolation method for geological modeling data based on anisotropic Huber robust reweighting. By constructing a strike-dip bidirectional structural feature coordinate system, it quantifies the differential influence of stratigraphic strike and dip on the interpolation weight coefficients, overcoming the limitations of spatial isotropy in traditional methods. This ensures that the interpolation process closely matches the actual structural characteristics of strata that change gently along the strike and rapidly along the dip, accurately restoring the subsurface structural morphology and effectively overcoming the problem of insufficient modeling accuracy caused by neglecting structural directionality in traditional methods. Simultaneously, the invention introduces a Huber robust reweighting strategy to segment the fitting residuals, maintaining high weights for low residual points and automatically reducing the weights of outliers and noise points. This significantly suppresses the interference of local anomalies on the fitting trend, solving the defects of traditional weighted least squares interpolation, which is affected by extreme residuals and suffers from unstable interpolation accuracy. This method does not rely on artificial prior geological models and still maintains high stability and accuracy in complex scenarios such as sparse seismic interpretation data, abrupt changes in stratigraphic structure, and data containing anomalous noise. Through joint control of anisotropy and robustness, it achieves accurate reconstruction of stratigraphic boundary sharpness, accurate reconstruction of fault gradient, and maintenance of strike continuity, providing reliable data support for high-precision geological modeling, structural interpretation, and reservoir structure description, which is significantly better than traditional interpolation methods. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the overall steps of the method of the present invention; Figure 2 This is a flowchart of method step S2 of the present invention; Figure 3 This is a flowchart of method step S3 of the present invention; Figure 4 This is a flowchart of method step S4 of the present invention; Figure 5 This is a flowchart of step S5 of the method of the present invention. Detailed Implementation

[0017] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. The application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0018] like Figure 1 As shown, the spatial interpolation method for geological modeling data based on anisotropic Huber robust reweighting is characterized by the following steps: S1, Obtain the scatter dataset of stratigraphic points of the target structural area, which includes information on the relationship between planar coordinates and stratigraphic elevation; Specifically, step S1 involves acquiring a complete and valid stratigraphic scatter dataset of the target structural area to provide fundamental data support for subsequent interpolation calculations. During implementation, scatter data reflecting stratigraphic distribution within the target area should be comprehensively collected, covering key structural components and boundary regions to ensure representativeness and completeness. This dataset includes the planar coordinates and corresponding stratigraphic elevation information for each sample point. The total number of sample points should be reasonably determined based on the target area's scope, structural complexity, and interpolation accuracy requirements. Typically, the number of sample points should not be less than 1000; for extremely complex structural areas, the number should be increased to over 3000. During the data collection process, data errors must be strictly controlled. Planar coordinate errors should not exceed 0.5 units, and stratigraphic elevation errors should be controlled within 0.3 units to avoid affecting the accuracy of subsequent interpolation results due to errors in the original data. After collection, the data undergoes preliminary screening to remove abnormal data that significantly exceeds reasonable limits, forming a standardized original stratigraphic scatter dataset. This dataset directly determines the fundamental quality of the interpolation calculations and is a prerequisite for ensuring the smooth progress of subsequent steps and the reliability of the final interpolation results.

[0019] S2. Construct a local neighborhood for each original sample point, select several nearest neighboring points to form a neighborhood set, assume that the strata are approximately planar in a small range, construct a design matrix and solve the plane parameters by least squares, and use the correlation components in the plane parameters as the stratum elevation gradient at the sample point. Specifically, step S2 obtains the stratigraphic elevation gradient of each original sample point through local neighborhood construction and plane fitting, providing key parameters for subsequent structural orientation field construction. During implementation, for each original sample point, the nearest neighbor is selected from all sample points to form a local neighborhood according to the Euclidean distance metric. The number of neighboring points, i.e., the number of neighborhood sample points used for gradient estimation, needs to be adjusted based on data density and structural characteristics, generally ranging from 15 to 30. The upper limit is used in complex structural regions, and the lower limit is used in densely datad and gently sloping structural regions. Within the constructed local neighborhood, assuming the stratigraphy is approximately planar, a design matrix including the planar coordinates of neighboring points and constant terms is constructed based on the planar coordinates and stratigraphic elevation information of the neighboring points. Simultaneously, the elevation data corresponding to the neighboring points are organized into an elevation column vector. The plane parameters are solved using the ordinary least squares method to obtain the rate of change of the local plane in the x-direction, the rate of change in the y-direction, and the plane constant term. The first two components are directly used as the stratigraphic elevation gradient at that sample point. This step quantifies the local variation trend of the strata through local plane fitting. The gradient data of each sample point will become the core basis for constructing a continuous gradient field and determining the strike and dip of the strata. Its calculation accuracy directly affects the ability of the entire interpolation method to characterize the directionality of geological structures.

[0020] S3 performs scatter interpolation on discrete gradient data to construct a continuous gradient function. At any query point, the gradient parameters are obtained through this function. The maximum descent direction of the strata is defined as the plane projection of the dip direction, and an orthogonal strike direction vector is constructed. Specifically, step S3 expands the gradient data of discrete sample points into a continuous gradient field and constructs the strike-dip direction field at the query point, providing directional basis for anisotropic interpolation. In practice, the planar coordinates of the original sample points and their corresponding x-direction and y-direction gradients are first processed using scattered interpolation. Radial basis function interpolation is used to ensure that the two constructed continuous gradient functions have good smoothness and continuity, accurately representing the gradient parameters at any planar location. For any query point to be interpolated, its planar coordinates are substituted into the constructed continuous gradient function to accurately obtain the x-direction and y-direction gradients corresponding to that point. Based on the obtained gradient parameters, the maximum descent direction of the formation is calculated by inverting the gradient components and dividing by the square root of the sum of squares of the gradient parameters. This direction serves as the planar projection of the dip direction, while a small constant is introduced to avoid a zero denominator. Subsequently, by exchanging the components of the dip direction unit vector and changing the sign of one of the components, a strike direction unit vector orthogonal to the dip direction is constructed, ultimately forming a set of normalized dip and strike direction vectors at each query point. This step enables the transformation from discrete gradients to continuous directional fields, allowing the interpolation process to adapt to the directional characteristics of stratigraphic structures and laying the foundation for anisotropic distance metrics and weight construction.

[0021] S4. For any query point to be interpolated, select several nearest original sample points to form a neighborhood set, calculate the difference in planar coordinates between the neighborhood sample points and the query point, and project the difference vector onto the direction and dip directions. Specifically, step S4 completes the selection of neighborhood sample points for the query point and constructs coordinate projection, converting the planar coordinate difference into components in the directional-dip coordinate system, preparing for anisotropic distance calculation. During implementation, for each query point to be interpolated, based on the Euclidean distance metric, the nearest neighborhood sample points to the query point are selected from the original scattered data. The number of neighborhood sample points ranges from 20 to 40, adjusted according to the original data density. In sparse data regions, the number of neighborhood sample points is appropriately increased to ensure fitting reliability. For each selected neighborhood sample point, the difference between its planar coordinates and the query point's planar coordinates is calculated, yielding the x-direction and y-direction coordinate differences. Sufficient decimal places must be retained in the coordinate difference calculation to ensure accuracy. The coordinate difference vector of each neighborhood sample point is multiplied by the directional unit vector at the query point to obtain the component of the neighborhood sample point relative to the query point along the directional direction; simultaneously, the coordinate difference vector is multiplied by the dip unit vector to obtain the component along the dip direction. This step achieves the transformation from a conventional planar coordinate system to a structural coordinate system through coordinate projection, enabling subsequent distance measurement and weight allocation to fully consider the strike and dip characteristics of the stratigraphic structure. It is a key step in realizing anisotropic interpolation.

[0022] S5 introduces the orientation association length and the dip association length to define the anisotropic distance. Based on this distance, a spatial kernel weight coefficient is constructed to express the feature that the weight decays more slowly along the orientation direction and the weight decays more quickly along the dip direction. Specifically, step S5 defines anisotropic distances and constructs spatial weighting coefficients to achieve differentiated weight allocation along different structural directions. In implementation, strike association length and dip association length are first set, with the strike association length being 3 to 5 times the dip association length. The specific ratio is determined based on the stratigraphic structural characteristics; the gentler the change in strata along the strike direction, the larger the ratio, ensuring a feature expression where the weight decays more slowly along the strike direction and more rapidly along the dip direction. Based on the components of neighboring samples along the strike and dip directions, ratios are calculated with the corresponding association lengths. The two ratios are then squared, and the sum of the squared results yields the squared anisotropic distance from the neighboring sample to the query point. Based on the calculated anisotropic distances, spatial weighting coefficients are constructed using a Gaussian kernel function with a bandwidth parameter ranging from 5 to 15. This parameter controls the overall rate of weight decay with distance; a smaller bandwidth parameter results in faster weight decay, and vice versa. Through the above process, each neighborhood sample point obtains a corresponding spatial weight coefficient. This weight can reflect the correlation between the sample point and the query point in the structural direction, so that distant neighbor points in the strike direction can still maintain a certain weight, while the weight of nearby neighbor points in the dip direction decays rapidly, fully reflecting the anisotropic characteristics of stratigraphic structure.

[0023] S6. Assuming the strata are approximately local planes in the neighborhood of the query point, construct the design matrix and observation vector, construct the diagonal weight matrix based on the spatial kernel weight coefficient, and solve the plane parameter estimation by weighted least squares method.

[0024] Specifically, step S6 performs the 0th weighted least squares fitting, obtaining initial plane parameter estimates based on the spatial weighting coefficients, providing an initial solution for subsequent robust reweighting. During implementation, within the neighborhood of the query point, assuming the strata are approximately a local plane, a design matrix is ​​constructed based on the plane coordinates of the neighborhood sample points and the stratum elevation information, including the sample point's x-coordinate, y-coordinate, and constant terms. Simultaneously, the elevation data of the neighborhood sample points are organized into observation vectors. The construction of the design matrix and observation vectors must strictly correspond to the sample point order to avoid data misalignment. Based on the anisotropic spatial weighting coefficients obtained in step S5, a diagonal weight matrix is ​​constructed, where the diagonal elements are the spatial weighting coefficients of each neighborhood sample point, and the off-diagonal elements are all 0. The weighted least squares method is employed. This involves calculating the product of the transpose of the design matrix and the diagonal weight matrix, multiplying this product by the design matrix, inverting the result, and then multiplying the product by the transpose of the design matrix, the diagonal weight matrix, and the observation vector. This yields the zeroth plane parameter estimate, which includes the rate of change of the plane in the x-direction, the rate of change in the y-direction, and a plane constant term. This step utilizes only spatial anisotropic weights for smooth fitting without introducing a robust weighting strategy. The resulting initial plane parameters serve as the basis for subsequent residual calculations and the construction of Huber robust weights, providing an initial reference for obtaining the final high-precision interpolation results.

[0025] Preferably, the model for solving the plane parameters using least squares in S2 is as follows: ,in, The fitted plane parameter column vector includes the rate of change in the x-direction, the rate of change in the y-direction, and a plane constant term. To design the transpose of the matrix, The design matrix for local plane fitting. It is the column vector of elevations of the neighborhood sample points.

[0026] Specifically, in step S2, the plane parameter model is solved using least squares. This model serves as the basis for calculating the elevation gradient of the sample points. During implementation, the plane parameter column vector in the model includes three key components, corresponding to the rate of change of the local plane in the x-direction, the rate of change in the y-direction, and a plane constant term. These three parameters collectively characterize the approximate planar morphology of the strata in the local neighborhood. The design matrix must be constructed strictly according to the planar coordinates of the neighborhood sample points. Each row of the matrix sequentially arranges the x-coordinate, y-coordinate, and constant 1 of a single neighborhood sample point. The number of columns in the matrix is ​​fixed at 3, and the number of rows is consistent with the number of neighborhood sample points used for gradient estimation. The number of neighborhood sample points is typically between 15 and 30 to ensure the reliability of the plane fitting. The transpose of the design matrix must be accurately calculated according to matrix operation rules. After multiplying it by the original design matrix, a square matrix is ​​obtained. This square matrix is ​​then inverted, and subsequently multiplied sequentially with the transpose of the design matrix and the elevation column vector of the neighborhood sample points to finally obtain the plane parameter column vector. The elevation column vector is composed of the stratigraphic elevations of each sample point in the neighborhood, arranged sequentially, with the number of elements matching the number of rows in the design matrix. This model minimizes the deviation between the actual elevation of the sample points and the plane-fitted elevation using the least squares principle, ensuring that the solved plane parameters accurately reflect the local stratigraphic trends. This provides accurate foundational data for subsequent gradient extraction and directly impacts the accuracy of the strike-dip direction field construction.

[0027] Preferably, the calculation model for the unit vector of the tendency direction at the query point in S3 is as follows: ,in, The unit vector is the direction of inclination. The gradient in the x-direction at the query point. The gradient in the y-direction at the query point is normalized by the square root of the sum of squares of the gradient parameters in the denominator.

[0028] Specifically, the calculation model for the dip direction unit vector at the query point in step S3 is a crucial step in constructing the strike-dip direction field. During implementation, the x-direction gradient and y-direction gradient of the query point in the model are obtained through the continuous gradient function constructed in step S3. This function is obtained by interpolating the discrete gradient data of the original sample points, ensuring accurate characterization of the gradient features of any query point. In the calculation process, the x-direction gradient and y-direction gradient of the query point are first inverted to obtain the original vector of the maximum descent direction of the formation, which approximately corresponds to the planar projection of the dip direction. To convert it into a unit vector, the sum of the squares of the x-direction gradient and the y-direction gradient is calculated. The square root of this sum is then used as the denominator. The two components of the original vector are divided by this denominator to obtain the normalized dip direction unit vector. To avoid the special case of a zero denominator, a small constant is introduced into the denominator, typically not exceeding 0.0001, to ensure the stability of the calculation process. The dip direction unit vector obtained by this model has a magnitude of 1, which can accurately indicate the dip direction of the strata at the query point. It provides an orthogonal reference for the construction of subsequent strike direction vectors and is an important directional basis for realizing anisotropic interpolation. It directly determines the degree of adaptation of the interpolation process to the directionality of the stratigraphic structure.

[0029] Preferably, the calculation model for the anisotropic distance in S5 is as follows: ,in, The anisotropic distance from the neighboring sample points to the query point. This represents the component of the planar difference between neighboring sample points and the query point along the direction. The component along the direction of inclination, The length is associated with the direction of travel. The length associated with the direction of inclination.

[0030] Specifically, the anisotropic distance calculation model in step S5 provides a foundation for constructing spatial weighting coefficients by quantifying the relative positional relationship between neighboring sample points and the query point in the structural direction. During implementation, the strike-direction and dip-direction components of the neighboring sample points in the model are obtained from the coordinate projection calculation in step S4, accurately reflecting the position of the sample points relative to the query point in the structural coordinate system. The strike association length and dip association length are key parameters of the model and need to be reasonably set according to the stratigraphic structural characteristics. The strike association length is taken as 3 to 5 times the dip association length; the gentler the change in strata along the strike direction, the larger this ratio should be, ensuring that the weight decay along the strike direction is relatively slow and the weight decay along the dip direction is relatively fast. In the calculation process, the strike-direction component is first compared to the strike association length, and the dip-direction component is also compared to the dip association length. Then, the two ratios are squared separately, and finally, the two squared values ​​are summed to obtain the squared value of the anisotropic distance. This model breaks through the limitation of the isotropic nature of traditional Euclidean distance. By introducing directional correlation length differences, the distance measurement can fully reflect the actual variation law of stratigraphic structure. It allows samples that are far apart in the strike direction to still maintain a certain correlation, while samples that are close together in the dip direction have a rapidly enhanced correlation, providing a scientific distance basis for subsequent differentiated weight allocation.

[0031] Preferably, the construction model for the spatial weighting coefficients in S5 is as follows: ,in, For spatial anisotropy weights, For anisotropic distances, This is a bandwidth parameter used to adjust the rate at which the weight decays with distance.

[0032] Specifically, in step S5, the spatial weight coefficient construction model is used to allocate differentiated weights to neighboring sample points based on anisotropic distance. During implementation, the anisotropic distance in the model is calculated using the model corresponding to weight 4, comprehensively reflecting the relative positional relationship between the sample point and the query point in the construction direction. The bandwidth parameter is key to controlling the weight decay rate, ranging from 5 to 15. The smaller the bandwidth parameter, the faster the weight decays with increasing anisotropic distance, and vice versa. It needs to be reasonably selected based on the original data density and construction complexity. For densely dataed and gently constructed regions, the bandwidth parameter can be appropriately increased, while for sparsely dataed and complexly constructed regions, the bandwidth parameter should be decreased. In the calculation process, the square of the anisotropic distance is first divided by the product of 2 and the square of the bandwidth parameter, then the result is negative, and finally substituted into the exponential function to obtain the spatial weight coefficient of the neighboring sample point. The model uses a Gaussian kernel function to construct weights, ensuring that the weight values ​​are between 0 and 1. Spatially, the closer the sample point is to the query point and along the strike direction, the closer the weight value is to 1, and the greater its contribution to the interpolation result. The farther the sample point is and along the dip direction, the closer the weight value is to 0, and the smaller its contribution. This achieves differentiated weight allocation along the structural direction, allowing the interpolation process to fully fit the changing characteristics of the strata, which are gentle in strike and steep in dip.

[0033] Preferably, the model for solving the plane parameters using weighted least squares in S6 is as follows: ,in, For the parameter estimation of the 0th plane fitting, Let be the coordinate matrix of the neighborhood sample points. Let be the first spatial weight diagonal matrix. Let the elevation column vector of the neighborhood sample points be denoted as . It is the transpose of the neighborhood sample point coordinate matrix.

[0034] Specifically, the weighted least squares model for solving the plane parameters in step S6 is used to obtain the 0th plane fitting parameter estimate, providing an initial solution for subsequent robust reweighting. In implementation, the neighborhood sample point coordinate matrix in the model is constructed by sequentially arranging the x-coordinates, y-coordinates, and constant 1 of the samples within the neighborhood of the query point. The number of rows in the matrix is ​​consistent with the number of neighborhood samples, which ranges from 20 to 40 to ensure the reliability of the fitting results. The diagonal elements of the initial spatial weight diagonal matrix are the spatial weight coefficients of each neighborhood sample point, while the off-diagonal elements are all 0. The spatial weight coefficients are calculated using the model corresponding to weight 5 and reflect the spatial correlation of the samples. The neighborhood sample point elevation column vector is composed of the stratigraphic elevations of the samples within the neighborhood, arranged sequentially, and is consistent with the number of rows in the neighborhood sample point coordinate matrix. In the calculation process, the transpose of the neighborhood sample point coordinate matrix is ​​first multiplied by the initial spatial weight diagonal matrix, and then multiplied by the neighborhood sample point coordinate matrix to obtain a square matrix. After inverting this square matrix, matrix multiplication is performed sequentially with the transpose of the neighborhood sample point coordinate matrix, the initial spatial weight diagonal matrix, and the neighborhood sample point elevation column vector. Finally, the parameter estimates for the 0th plane fitting are obtained, including the rate of change of the plane in the x-direction, the rate of change in the y-direction, and the plane constant term. This model only fits based on spatial anisotropic weights and does not introduce a robust weight adjustment strategy. Its results serve as a benchmark for subsequent residual calculations, providing an initial reference for the construction of Huber robust weights and representing an important starting point for realizing "anisotropic + robust" joint control.

[0035] Preferred, such as Figure 2 The S2 step includes the following steps: S21, selecting a predetermined number of neighboring points closest to the target sample point from all sample points based on Euclidean distance, and forming a local neighborhood set for the sample point; S22, assuming that the strata are distributed in a planar manner within the local neighborhood, constructing a design matrix including the planar coordinates of the neighboring points and constant terms, and simultaneously organizing the elevation data corresponding to the neighboring points to form an elevation column vector; S23, performing matrix multiplication on the design matrix and its transpose, inverting it, and then performing matrix multiplication with the transpose of the design matrix and the elevation column vector in sequence to obtain a planar parameter column vector; S24, extracting the first two components of the planar parameter column vector as the gradient of the strata elevation at the target sample point in the x and y directions.

[0036] Specifically, step S2 involves a step-by-step implementation process for obtaining the elevation gradient of the sample points. In step S21, based on the Euclidean distance metric, 15 to 30 nearest neighboring points to the target sample point are selected from all original sample points to form a local neighborhood set. For complex regions, 30 neighboring points are selected; for densely populated and gently constructed regions, 15 neighboring points are selected, ensuring the neighborhood fully reflects local stratigraphic characteristics. Step S22 assumes the stratigraphy is planar within the local neighborhood. Based on the planar coordinates of the neighboring points and the stratigraphic elevation, a design matrix is ​​constructed, with each row containing the planar coordinates of a single neighboring point and a constant term. Simultaneously, elevation data is organized in neighboring point order to form an elevation column vector, with the number of matrix rows matching the number of selected neighboring sample points. Step S23 strictly follows matrix operation rules, first multiplying the design matrix by its transpose to obtain a square matrix. After inverting this square matrix, matrix multiplication is performed sequentially with the transpose of the design matrix and the elevation column vector to obtain a planar parameter column vector containing the x-direction rate of change, the y-direction rate of change, and a planar constant term. Step S24 directly extracts the first two components from the plane parameter column vector, which are used as the gradients of the stratum elevation in the x and y directions at the target sample point, respectively. This gradient data accurately quantifies the local variation trend of the strata, providing a core foundation for the subsequent construction of continuous gradient fields and determination of strike-dip directions. The accuracy of its calculation directly affects the adaptation effect of the entire interpolation method to the directionality of geological structures. The reasonable selection of the number of neighborhood sample points and the accurate execution of matrix operations are the key to ensuring the accuracy of gradient estimation.

[0037] Preferred, such as Figure 3 The S3 step includes the following sub-steps: S31, performing scatter interpolation on the planar coordinates of the original sample points and the corresponding x-direction gradient and y-direction gradient to construct two continuous gradient functions to characterize the gradient parameters at any planar position; S32, substituting the planar coordinates of the query point into the two continuous gradient functions to obtain the x-direction gradient and y-direction gradient corresponding to the query point; S33, obtaining the unit vector of the maximum descent direction of the formation, i.e., the dip direction, by inverting the gradient parameters and dividing by the square root of the sum of the squares of the gradient parameters; S34, constructing a strike direction unit vector orthogonal to the dip direction by exchanging the components and changing the sign of one of the components based on the components of the dip direction unit vector.

[0038] Specifically, step S3 is implemented in steps to transform discrete gradients into continuous directional fields. In step S31, the radial basis function interpolation method is used to perform scattered interpolation on the planar coordinates of the original sample points and their corresponding x-direction and y-direction gradients, constructing two continuous gradient functions. The interpolation process ensures that the functions have good smoothness and continuity, accurately representing the gradient parameters at any planar position, providing a reliable model for obtaining the gradient of the query point. Step S32 substitutes the planar coordinates of any query point to be interpolated into the two constructed continuous gradient functions, accurately obtaining the x-direction and y-direction gradients corresponding to the query point through function calculation. The gradient values ​​are calculated with four decimal places to avoid error accumulation affecting subsequent direction calculations. Step S33 inverts the x-direction and y-direction gradients of the query point to obtain the original vector of the maximum descent direction of the stratum. The sum of the squares of the two components of this original vector is calculated, and the square root is used as the denominator. The original vector components are divided by this denominator to obtain a normalized dip direction unit vector. A small constant of 0.0001 is introduced to avoid the special case of a zero denominator, ensuring computational stability. Step S34 involves swapping the two components of the dip direction unit vector and changing the sign of one of the components to construct a strike direction unit vector orthogonal to the dip direction. This results in a set of normalized dip and strike direction vectors at each query point. This set of vectors accurately reflects the directional characteristics of the stratigraphic structure and provides crucial directional information for subsequent anisotropic interpolation. The choice of interpolation method and the reasonable setting of small constants directly affect the accuracy and stability of the directional field construction.

[0039] Preferred, such as Figure 4 The S4 step includes the following sub-steps: S41, according to the Euclidean distance metric, a preset number of sample points closest to the query point are selected from the original scattered data to form a local neighborhood set of the query point; S42, the difference between the planar coordinates of each sample point in the neighborhood set and the planar coordinates of the query point are calculated one by one to obtain the x-direction coordinate difference and y-direction coordinate difference of each neighboring point relative to the query point; S43, the coordinate difference vector of each neighboring point is multiplied by the unit vector of the direction of movement to obtain the component of the neighboring point relative to the query point along the direction of movement; S44, the coordinate difference vector of each neighboring point is multiplied by the unit vector of the direction of movement to obtain the component of the neighboring point relative to the query point along the direction of movement.

[0040] Specifically, step S4 involves the selection of neighborhood sample points and the construction of coordinate projection for the query point. In step S41, following the Euclidean distance metric, the 20 to 40 nearest sample points to the query point are selected from the original scattered data to form a local neighborhood set. Forty neighborhood sample points are selected for sparse data areas, and 20 for dense data areas, ensuring the fitting process has sufficient data support while reflecting local stratigraphic characteristics. Step S42 calculates the difference between the planar coordinates of each sample point in the neighborhood set and the planar coordinates of the query point, obtaining the x-coordinate difference and y-coordinate difference of each neighboring point relative to the query point. The coordinate difference calculation is rounded to four decimal places to ensure accuracy and avoid affecting subsequent projection results due to coordinate difference errors. Step S43 performs a dot product operation on the coordinate difference vector of each neighboring point and the strike direction unit vector. The dot product operation strictly follows vector operation rules to obtain the strike direction component of the neighboring point relative to the query point. This component accurately reflects the relative positional relationship between the sample point and the query point in the strike direction. Step S44 involves performing a dot product operation between the coordinate difference vector of each neighboring point and the unit vector of the tilt direction, following the same vector operation rules, to obtain the component of the neighboring point relative to the query point along the tilt direction. This component accurately represents the relative position of the sample point to the query point in the tilt direction. Through these four steps, the transformation from the conventional planar coordinate system to the constructed coordinate system is achieved, providing accurate positional parameters for subsequent anisotropic distance measurement and differentiated weight allocation. Reasonable control of the number of neighboring sample points and strict control of computational accuracy are crucial to ensuring the reliability of the projection results.

[0041] Preferred, such as Figure 5 The S5 step includes the following sub-steps: S51, setting the orientation association length representing the extension characteristics of the orientation direction and the orientation association length representing the change characteristics of the orientation direction, and clarifying the difference in length parameters between the two directions; S52, based on the components of the neighboring point along the orientation direction and the orientation direction, performing ratio calculations with the corresponding association lengths respectively, and then squaring the two ratio results; S53, summing the two squared results to obtain the squared value of the anisotropic distance from the neighboring point to the query point; S54, substituting the squared value of the anisotropic distance into the exponential function, and controlling the attenuation rate by setting the bandwidth parameter to construct the spatial weight coefficient of each neighboring point.

[0042] Specifically, step S5 involves constructing anisotropic distance and spatial weighting coefficients. In step S51, a strike correlation length representing the extension characteristics of the strike direction and a dip correlation length representing the change characteristics of the dip direction are set. The strike correlation length is 3 to 5 times the dip correlation length; the gentler the change in strata along the strike direction, the larger the ratio. Typically, a ratio of 5 times is used in areas with gentle strata and a ratio of 3 times is used in areas with significant strata changes, ensuring that the characteristics of gentle strike and steep dip changes in strata are reflected. Step S52, based on the components of neighboring points along the strike and dip directions, performs ratio calculations with the corresponding correlation lengths, retaining four decimal places. The two ratio results are then squared, which amplifies the distance differences in different directions. Step S53 sums the two squared results to obtain the squared value of the anisotropic distance from the neighboring point to the query point. This squared value fully integrates the relative positional information in different directions, overcoming the limitation of the isotropic nature of traditional Euclidean distance. Step S54 substitutes the squared value of the anisotropic distance into the exponential function, setting the bandwidth parameter to 5 to 15. For dense, gently constructed regions, the bandwidth parameter is set to 15, while for sparse, complex regions, it is set to 5. By adjusting the weight decay rate through the bandwidth parameter, a spatial weight coefficient for each neighboring point is constructed. This coefficient, ranging from 0 to 1, accurately reflects the correlation between the sample point and the query point along the construction direction. These four steps are closely linked. Through the reasonable setting of the association length ratio, strict control of computational precision, and flexible adjustment of the bandwidth parameter, differentiated weight allocation along the construction direction is achieved, providing a scientific weight basis for subsequent high-precision interpolation. The optimized configuration of the association parameter directly affects the accuracy of the anisotropic feature representation and the reliability of the interpolation results.

[0043] This geological modeling data spatial interpolation method, based on anisotropic Huber robust reweighting, achieves dual optimization of geological structural directionality and data robustness, overcoming the technical limitations of traditional interpolation methods. It constructs a strike-dip bidirectional structural feature coordinate system through local gradient estimation, quantifying the differentiated impact of different directions on the interpolation weight coefficients. This ensures the interpolation process closely aligns with the actual structural patterns of gently changing strata along the strike and rapidly changing dip, enabling the interpolation results to accurately reproduce subsurface structural morphology. Simultaneously, a Huber robust reweighting strategy is introduced to segment the fitting residuals, maintaining high weights for low-residual points that conform to geological patterns and automatically reducing the weights of outliers and noise points that deviate from the patterns, significantly improving the data processing's resistance to interference. Furthermore, this method does not rely on manual prior geological models, can automatically adapt to complex structural features, and maintains high stability and accuracy even in complex scenarios such as sparse data, abrupt stratigraphic changes, and the presence of anomalous noise, making it widely applicable.

[0044] This method constructs anisotropic distance metrics and spatial weighting coefficients to achieve differentiated expression where the weight decays more slowly along the strike direction and more quickly along the dip direction. This shifts the interpolation process from a "geometric" to a "structural" one, accurately restoring the true morphology of subsurface structures. Addressing the issues of instability in interpolation accuracy caused by outliers and noise in traditional weighted least squares methods, the method employs a robust Huber reweighting strategy and joint control of spatial weights to effectively suppress the distortion of the fitting plane by extreme residuals. While maintaining the overall structural trend, it achieves accurate reconstruction of stratigraphic boundary sharpness, fault gradient, and strike continuity, providing reliable data support for high-precision geological modeling and significantly outperforming traditional interpolation methods.

[0045] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms "set," "install," "connect," "link," and "fix" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal communication between two components. Those skilled in the art will understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0046] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various equivalent changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A geostatistical modeling data space interpolation method based on anisotropic Huber robust reweighting, characterized in that, The method comprises the following steps: S1, obtaining a stratum scattered point data set of a target construction area, the data set comprising plane coordinate and stratum elevation associated information; S2, constructing a local neighborhood for each original sample point, selecting a plurality of neighboring points closest to the sample point to form a neighborhood set, assuming that the stratum is approximately planar in a small range, constructing a design matrix and solving the plane parameters by least square, taking the associated components in the plane parameters as the stratum elevation gradient at the sample point; S3, performing scattered point interpolation on the discrete gradient data, constructing a continuous gradient function, obtaining the gradient parameters at any query point through the function, defining the maximum stratum descending direction as the planar projection of the tendency direction and constructing an orthogonal strike direction vector; S4, for any query point to be interpolated, selecting a plurality of original sample points closest to the query point to form a neighborhood set, calculating the plane coordinate difference value of the neighborhood sample points relative to the query point, and projecting the difference vector onto the strike direction and the tendency direction; S5, introducing the strike associated length and the tendency associated length to define anisotropic distance, constructing a spatial kernel weight coefficient based on the distance, and performing feature expression with slow weight value attenuation in the strike direction and fast weight value attenuation in the tendency direction; S6, assuming that the stratum is approximately planar in the neighborhood of the query point, constructing a design matrix and an observation vector, constructing a diagonal weight matrix based on the spatial kernel weight coefficient, and solving the plane parameter estimation by the weighted least square method.

2. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The model for solving the plane parameters by least square in S2 is: wherein, is the column vector of plane parameters obtained by fitting, including the rate of change in x direction, the rate of change in y direction and the constant term of the plane, is the transpose of the design matrix, is the design matrix of local plane fitting, is the elevation column vector of the neighborhood sample points.

3. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The calculation model of the tendency direction unit vector at the query point in S3 is: wherein, is the tendency direction unit vector, is the x-direction gradient at the query point, is the y-direction gradient at the query point, and the denominator is the square root of the square of the sum of the gradient parameters for normalization.

4. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The calculation model of the anisotropic distance in S5 is: wherein, is the anisotropic distance of the neighborhood sample to the query point, is the component of the planar difference of the neighborhood sample relative to the query point along the strike direction, is the component along the dip direction, is the correlation length of the strike direction, is the correlation length of the dip direction.

5. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The construction model of the spatial kernel weight coefficient in S5 is: wherein, is a spatial anisotropic weight, is an anisotropic distance, is a bandwidth parameter for regulating the rate of weight decay with distance.

6. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The model for solving the plane parameters by weighted least squares in S6 is: wherein, is the parameter estimation of the 0th plane fitting, is the coordinate matrix of the neighborhood sample points, is the initial space weight diagonal matrix, is the elevation column vector of the neighborhood sample points, is the transpose of the coordinate matrix of the neighborhood sample points.

7. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The S2 comprises the following steps: S21, selecting a preset number of neighboring points closest to the target sample point from all sample points according to the Euclidean distance, and assembling a local neighborhood set of the sample point; S22, assuming that the stratum is planarly distributed in the local neighborhood, constructing a design matrix comprising the plane coordinates of the neighboring points and a constant term, and arranging the elevation data corresponding to the neighboring points to form an elevation column vector; S23, performing matrix multiplication operation on the design matrix and its transpose, inverting, and then performing matrix multiplication on the design matrix transpose and the elevation column vector in sequence to obtain a plane parameter column vector; S24, extracting the first two components in the plane parameter column vector as the gradients of the stratum elevation in the x direction and the y direction at the target sample point.

8. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The S3 comprises the following steps: S31, performing scattered point interpolation processing on the plane coordinates of the original sample points and the corresponding x direction gradient and y direction gradient respectively, constructing two continuous gradient functions for representing the gradient parameters at any plane position; S32, substituting the plane coordinates of the query point into the two continuous gradient functions to obtain the x direction gradient and the y direction gradient corresponding to the query point; S33, obtaining the unit vector of the stratum maximum descending direction, i.e. the tendency direction, by taking the inverse of the gradient parameters and dividing by the square root of the square sum of the gradient parameters; S34, based on the components of the tendency direction unit vector, constructing a strike direction unit vector orthogonal to the tendency direction by exchanging the components and changing the sign of one component.

9. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The S4 comprises the following steps: S41, according to the Euclidean distance measurement rule, screening out a preset number of sample points closest to the query point in the original scattered point data to form a local neighborhood set of the query point; S42, calculating the difference between the planar coordinates of each sample point in the neighborhood set and the planar coordinates of the query point to obtain the x-direction coordinate difference and the y-direction coordinate difference of each neighborhood point relative to the query point; S43, performing dot product operation on the coordinate difference vector of each neighborhood point and the unit vector of the strike direction to obtain the component of the neighborhood point along the strike direction relative to the query point; S44, performing dot product operation on the coordinate difference vector of each neighborhood point and the unit vector of the dip direction to obtain the component of the neighborhood point along the dip direction relative to the query point.

10. The anisotropic Huber robust reweighted based geostatistical modeling data space interpolation method of claim 1, wherein, The S5 comprises the following steps: S51, setting a strike correlation length representing the extension characteristics of the strike direction and a dip correlation length representing the change characteristics of the dip direction to clearly distinguish the length parameter difference between the two directions; S52, performing ratio operation on the components of the neighborhood point along the strike direction and the dip direction and the corresponding correlation lengths respectively, and then performing square processing on the two ratio results; S53, performing summation operation on the two results after square processing to obtain the square value of the anisotropic distance of the neighborhood point to the query point; S54, substituting the square value of the anisotropic distance into the exponential function to control the decay rate by setting a bandwidth parameter, and constructing the spatial kernel weight coefficient of each neighborhood point.

Images