Multi-source geological data processing method and system for three-dimensional geological model

Abstract

The invention relates to the technical field of multi-source data fusion, in particular to a multi-source geological data processing method and system of a three-dimensional geological model.The method comprises the following steps that mountain landform and river valley images are obtained, gray frequency characteristics are extracted, a frequency energy gradient layer is constructed, a frequency continuous response area is screened to generate a structure boundary set, and a structure boundary set is constructed; the method comprises the following steps of: extracting a boundary normal vector by utilizing principal component analysis, identifying boundary sections with consistent directions, estimating a physical property parameter gradient direction, judging an included angle screening blocking region, generating a space attribute limiting layer, carrying out space alignment analysis on an overlapping region vector included angle, updating a boundary label, and generating an available attribute path structure set in three-dimensional geological modeling through a Dijkstra algorithm. According to the method, a conduction model is constructed through frequency domain decomposition and logarithmic transformation enhanced recognition, frequency window analysis noise reduction, principal component extraction vector analysis direction and center difference estimation, dynamic matching is promoted through alignment, a Dijkstra algorithm optimizes a path, and the geological model bedding characterization and conduction simulation precision is improved through cooperation of a multi-dimensional technology.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of multi-source data fusion, in particular to a multi-source geological data processing method and system for three-dimensional geological models. BACKGROUND

[0002] The technical field of multi-source data fusion includes obtaining, analyzing, converting and uniformly processing data from multiple heterogeneous sources to form a comprehensive information expression form that can support high-level decision-making or model construction. The core content mainly includes data collection standardization, coordinate unification between data, semantic matching, format standardization, time and space alignment, etc. The key steps aim to effectively integrate data with different formats, different accuracies and different sources into a unified data set with consistency and comparability. Systematically, multi-source data fusion not only covers data preprocessing, data correlation and integration, but also involves data consistency verification, data fusion strategy selection and execution, and is widely used in remote sensing and mapping, geographic information systems, intelligent transportation, urban planning and other industries. In particular, in the field of three-dimensional modeling and geological exploration, through the fusion processing of multiple types of data sources, a more realistic and accurate spatial information expression system can be constructed.

[0003] Among them, the multi-source geological data processing method and system for three-dimensional geological models refers to the problems of inconsistent spatial resolution, different coordinate reference, heterogeneous attribute fields, etc. between different types of geological data. Through the construction of data uniform conversion rules, the establishment of geological attribute mapping relationship, the uniform transformation of coordinate system, and the combination of specific sub-class data such as stratigraphic drilling data, geological profile image, geophysical survey data, etc., the original geological data is standardized in format, spatially aligned and semantically reconstructed, and then a unified data set that can be used to construct three-dimensional geological models is generated. Through the coordinate conversion function, the spatial reference is unified, the attribute data integration is completed by using the field matching logic and attribute classification coding, and the image analysis algorithm is combined to extract the profile and identify the horizon of the geological image information, realizing the multi-source geological information aggregation processing required for three-dimensional modeling.

[0004] In the processing of multi-source geological data, there is a lack of effective means for feature extraction of micro-fabric structure, and the boundary of layer is determined by artificial experience, which leads to insufficient precision of micro-fabric structure representation in three-dimensional modeling. In the coordinate system transformation process, the influence of local structure trend on spatial alignment is not considered, which may cause the direction deviation of interlayer property conduction path and real geological structure. In the attribute data integration process, the static field matching method is used, which is difficult to adapt to the spatial gradient change characteristics of physical parameters in different survey areas, leading to deviation of interlayer conduction simulation results from the actual physical law. In the data preprocessing stage, the correlation model between frequency characteristics and fabric structure is not established, which cannot effectively distinguish noise signals from real geological structure characteristics, affecting the reliability of layer division under complex geological conditions. In the spatial alignment process, there is a lack of dynamic matching mechanism of direction vector, and when there are local structure direction differences in multi-source data, the traditional method may cause topological error of interlayer conduction path. In the existing method, the global algorithm is used in the analysis of physical parameters, and the constraint effect of local blocking area on the conduction path is ignored, which limits the engineering application value of three-dimensional model under complex geological conditions. For example, the traditional coordinate transformation does not combine with structure trend analysis, which may cause more than 15-degree deviation between shale fabric direction and permeability simulation path, directly affecting the effect of oil and gas exploration and development plan. SUMMARY

[0005] To solve the technical problems existing in the prior art, the embodiments of the present application provide a multi-source geological data processing method and system of a three-dimensional geological model. The technical scheme is as follows:

[0006] In order to achieve the above-mentioned purpose, the present application adopts the following technical scheme, a multi-source geological data processing method of a three-dimensional geological model, comprising the following steps:

[0007] S1: Obtain mountain landform and valley area image data, extract the gray sequence of the area image and decompose the frequency, construct a frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude by using a DFT function and align, and generate a frequency energy gradient layer;

[0008] S2: Screen the frequency continuous response area of mountain slope area and river alluvial zone through the frequency energy gradient layer, sort the frequency energy gradient value and judge the boundary position, and generate a structure boundary set;

[0009] S3: Based on the structure boundary set, construct a normal vector group for the boundary of the mountain slope area and the river alluvial zone by using a principal component analysis method, analyze the included angle relationship between the vectors and identify the direction consistent boundary segment, and generate direction classification information;

[0010] S4: Estimate the gradient direction of the physical property parameters of multiple points based on the direction classification information, judge the direction included angle area and screen the blocking area, and generate a spatial attribute restriction layer;

[0011] S5: the spatial attribute restriction layer is spatially aligned with the structure boundary set, the direction vector included angle of the mountainous topography and the valley region overlapping region is analyzed, the boundary label is updated, the rationality of the execution direction is judged through the Dijkstra algorithm, and a three-dimensional geological modeling available attribute path structure set is generated.

[0012] As a further scheme of the present application, the frequency energy gradient layer includes a main frequency response value distribution, a gray frequency variation section, and a frequency amplitude logarithmic transformation, the structure boundary set specifically refers to a boundary line segment spatial sequence, a slope response range section, and a continuous boundary frequency window label, the direction classification information includes a boundary normal vector set, a main direction included angle interval, and a direction consistency grouping label, the spatial attribute restriction layer specifically refers to an attribute gradient direction value distribution, a direction included angle partition boundary line, and a blocking area boundary point number, and the attribute path structure set includes a shortest path length value, a cumulative turning included angle degree, and a continuous boundary section node sequence.

[0013] As a further scheme of the present application, the mountainous topography and the valley region image data are acquired, the gray sequence of the region image is extracted and the frequency is decomposed, the main frequency amplitude is constructed by constructing a frequency window, the main frequency amplitude is logarithmically transformed and aligned through a DFT function, and the steps of generating the frequency energy gradient layer are specifically as follows:

[0014] S101: mountainous topography and valley region image data in a three-dimensional geological model are acquired and divided into a plurality of sub-regions, for each sub-region image, pixel values are extracted and converted into a gray sequence, the range of the gray values is adjusted through linear mapping, and a sub-region gray sequence is generated;

[0015] S102: according to the sub-region gray sequence, the frequency components in the gray sequence are decomposed and the main frequency amplitude is extracted through a fast Fourier transform function, the logarithmic amplitude of the main frequency amplitude is calculated through a DFT function and aligned with the position parameter, and a main frequency amplitude logarithmic transformation result is generated;

[0016] S103: the main frequency amplitude logarithmic transformation result is called, the energy gradient values corresponding to the multiple position frequencies are calculated and arranged, and a frequency energy gradient layer is acquired.

[0017] As a further scheme of the present application, the mountainous topography and the valley region image data are acquired, the gray sequence of the region image is extracted and the frequency is decomposed, the main frequency amplitude is constructed by constructing a frequency window, the main frequency amplitude is logarithmically transformed and aligned through a DFT function, and the steps of generating the frequency energy gradient layer are specifically as follows:

[0018] S201: the frequency energy gradient values of multiple pixel positions in the frequency energy gradient layer are acquired, the frequencies of the mountainous slope region and the river alluvial zone layer are screened and position aggregated, and a continuous response aggregation value interval is generated;

[0019] S202: Calculating the decreasing amplitude of the edge position frequency of the mountain slope area and the river alluvial zone area according to the continuous response aggregation value interval, constructing a boundary contour map based on the geometric shape of the edge of the response area, determining the boundary direction, and obtaining a continuous response boundary direction coefficient;

[0020] S203: calling the continuous response boundary trend coefficient, screening boundary line segments with a gentle change rate for the trend change rate of all boundary points, numbering and dividing the selected boundary line segments, and establishing a structural boundary set.

[0021] As a further solution of the present invention, based on the structural boundary set, a normal vector group is constructed for the boundaries of the mountain slope area and the river alluvial zone through principal component analysis, the angle relationship between the vectors is analyzed, and the boundary segments with consistent directions are identified. The specific steps of generating direction classification information are as follows:

[0022] S301: extracting a boundary point set within the mountain slope area and the river alluvial zone based on the structural boundary set, constructing boundary direction vectors for the continuous boundary point set and summarizing them according to regional affiliation to generate a regional boundary direction vector set;

[0023] S302: Calling the region boundary direction vector set, extracting the first principal axis direction vector corresponding to the principal component, calculating the angle between the first principal axis direction vector and the boundary direction vector by principal component analysis, determining the vector consistency relationship, and generating a discrete value interval of the boundary direction angle;

[0024] S303: Calling the discrete value interval of the boundary direction angle to screen the angle consistency of multiple boundary segments, classifying the boundary segments with angles smaller than the consistency judgment reference value into the same direction classification number and marking them respectively, thereby generating direction classification information;

[0025] The consistency judgment reference value is set based on the directional consistency feature exhibited when the cosine value of the angle between the main axis direction vectors is greater than 0.97.

[0026] As a further solution of the present invention, the angle value θ between the first principal axis direction vector and the boundary direction vector is calculated by principal component analysis. i , using the formula:

[0027]

[0028] Among them, u k represents the kth principal axis vector extracted by principal component analysis, v irepresents the i-th original direction vector in the set of regional boundary direction vectors, σ represents a projection variance coefficient of the set of boundary vectors, and takes a value of 0.2-0.5, η represents an orthogonal component adjustment factor, and takes a value of 0.1-1.2, p represents the number of principal axes selected by principal component analysis, m represents the total number of boundary direction vectors, and v j represents the j-th original direction vector in the set of regional boundary direction vectors.

[0029] As a further scheme of the present application, the step of estimating the property parameter gradient direction of the multi-point based on the direction classification information, judging the direction angle region and screening the blocking area, and generating the spatial attribute restriction layer is specifically:

[0030] S401: Based on the direction classification information, the spatial angle between the multi-point and the corresponding boundary segment principal axis direction vector is calculated, the point direction projection vector is obtained and recorded based on the angle value and the corresponding position relationship of the point on the boundary segment, and a point direction angle sequence is generated;

[0031] S402: Based on the point direction angle sequence, the continuous point group is screened and the direction gradient change trend between adjacent points is judged, the gradient direction change rate between continuous direction abrupt change points is calculated, and a direction gradient change rate interval is generated;

[0032] S403: According to the direction gradient change rate interval, the point group with severe direction abrupt change and discontinuous direction distribution is spatially clustered, the point concentration distribution region is located, and a spatial attribute restriction layer is generated.

[0033] As a further scheme of the present application, the spatial attribute restriction layer and the structure boundary set are spatially aligned, the direction vector angle of the mountainous landform and the valley region overlap region is analyzed, the boundary label is updated, the execution direction rationality is judged by Dijkstra algorithm, and the step of generating the available attribute path structure set in three-dimensional geological modeling is specifically:

[0034] S501: Based on the spatial attribute restriction layer and the structure boundary set, the spatial alignment coordinate offset value is calculated, and the spatial position of the boundary line segment is adjusted, and the spatial alignment coordinate offset value is generated;

[0035] S502: The spatial alignment coordinate offset value is called to extract the boundary segment principal axis direction vector of the mountainous landform region and the valley region overlap region after alignment, calculate the angle cosine value between adjacent boundary segment principal axis vectors and screen the region, and obtain the direction angle difference value interval;

[0036] S503: According to the direction angle difference value interval, the shortest path length of multiple nodes in the direction difference section is calculated by Dijkstra algorithm, the priority path is constructed, the boundary section direction label is updated through the path node sequence, and the available attribute path structure set in three-dimensional geological modeling is generated.

[0037] As a further scheme of the present application, the shortest path length of multiple nodes in the direction difference section is calculated by Dijkstra algorithm The formula is:

[0038]

[0039] Wherein, represents the three-dimensional space Euclidean distance module length of node a and node b, The subscripts x and y represent the longitude and latitude direction grid numbers, the subscript z represents the vertical direction geological horizon number, ΔH is the elevation difference between nodes a and b, θ ab represents the instantaneous angle radian value, which is calculated after obtaining the node coordinate data by the three-dimensional laser scanner, represents the direction angle sliding average, which is set to 5-10 meters of local geological area.

[0040] On the other hand, a multi-source geological data processing system of three-dimensional geological model is provided, which is applied to the multi-source geological data processing method of three-dimensional geological model, and the system comprises:

[0041] The layer construction module is used to obtain mountainous landform and valley area image data, extract the gray sequence of the area image and decompose the frequency, construct the frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude by DFT function and align, generate the frequency energy gradient layer and pass to the boundary extraction module;

[0042] The boundary extraction module is used to screen the mountain slope area and the river alluvial zone frequency continuous response area through the frequency energy gradient layer, sort the frequency energy gradient value and judge the boundary position, generate the structure boundary set and pass to the direction identification module;

[0043] The direction identification module is used to construct the normal vector group of the boundary of the mountain slope area and the river alluvial zone through the principal component analysis method through the structure boundary set, analyze the angle relationship between the vectors and identify the direction consistent boundary section, generate the direction classification information and pass to the gradient estimation module;

[0044] The gradient estimation module is used to estimate the physical property parameter gradient direction of multiple points through the direction classification information, judge the direction angle area and screen the blocking area, generate the space attribute restriction layer and pass to the path generation module;

[0045] A path generation module is configured to spatially align the spatial attribute restriction layer with the structure boundary set, analyze the direction vector angle of the mountainous landform and the valley area overlapping area, update the boundary label, determine the rationality of the execution direction through the Dijkstra algorithm, and generate a three-dimensional geological modeling available attribute path structure set.

[0046] The technical scheme provided by the embodiment of the present application has at least the following beneficial effects:

[0047] Through the frequency decomposition of the gray sequence and the logarithmic transformation processing, the energy gradient feature recognition capability of the micro-layer structure in the geological image is effectively enhanced, and the geological profile boundary detection accuracy is improved. The frequency window main frequency amplitude analysis combined with the continuous response area screening mechanism reduces the influence of noise interference on the horizon division, and ensures the spatial continuity of the micro-layer structure profile division boundary. The principal component analysis method extracts the boundary normal vector group and analyzes the direction consistency, accurately identifies the geological structure trend, and enhances the explainability of the interlayer structure direction feature. The central difference method estimates the physical property parameter gradient direction, combined with the direction angle area blocking judgment, establishes the correlation model of the conduction restriction and the structure trend, and improves the physical rationality of the interlayer attribute conduction path analysis. The spatial alignment and direction vector angle analysis mechanism realize the dynamic matching of the structure boundary and the conduction restriction, and optimize the topology construction efficiency of the attribute path in the three-dimensional modeling process. The application of the Dijkstra algorithm in the direction restriction label graph effectively avoids the direction conflict problem in the traditional path planning, and ensures that the attribute conduction path conforms to the geological structure mechanical characteristics. The whole set of processing logic cooperates through multi-dimensional technologies such as frequency feature enhancement, direction consistency analysis and physical gradient direction constraint, significantly improves the accuracy of the three-dimensional geological model in the micro-layer structure representation and attribute conduction simulation. BRIEF DESCRIPTION OF DRAWINGS

[0048] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following will briefly introduce the drawings needed to be used in the embodiment description. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0049] Figure 1 It is a workflow diagram of the present application.

[0050] Figure 2 It is a system flowchart of the present application. DETAILED DESCRIPTION

[0051] The technical solutions in the present application will be described below with reference to the drawings.

[0052] In the embodiments of the present invention, words such as "exemplarily" and "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design described as an "exemplary" in the present invention should not be interpreted as being preferred or advantageous over other embodiments or designs. Rather, the use of the word "exemplary" is intended to present concepts in a concrete manner. Furthermore, in the embodiments of the present invention, "and / or" can mean both or either of the two.

[0053] In the embodiments of the present invention, the terms "image" and "picture" may be used interchangeably. It should be noted that, when the distinction between them is not emphasized, their intended meanings are the same. The terms "of," "corresponding," and "corresponding" may be used interchangeably. It should be noted that, when the distinction between them is not emphasized, their intended meanings are the same.

[0054] In the embodiments of the present invention, sometimes a subscript such as W1 may be written as a non-subscript such as W1. When the difference is not emphasized, the meanings to be expressed are the same.

[0055] In order to make the technical problems, technical solutions and advantages to be solved by the present invention clearer, a detailed description will be given below with reference to the accompanying drawings and specific embodiments.

[0056] See also Figure 1 The present invention provides a technical solution, a method for processing multi-source geological data of a three-dimensional geological model, comprising the following steps:

[0057] S1: Acquire image data of mountainous landforms and river valley areas, extract the grayscale sequence of the regional image and decompose the frequency, construct a frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude through the DFT function and align them to generate a frequency energy gradient layer;

[0058] S2: Use the frequency energy gradient layer to screen the frequency continuous response areas of mountain slopes and river alluvial zones, sort the frequency energy gradient values, determine the boundary positions, and generate a structural boundary set;

[0059] S3: Based on the structural boundary set, the principal component analysis method is used to construct a normal vector group for the boundaries of the mountain slope area and the river alluvial zone. The angle relationship between the vectors is analyzed and the boundary segments with consistent directions are identified to generate direction classification information.

[0060] S4: Estimate the physical parameter gradient direction of multiple points based on the direction classification information, determine the direction angle area and filter the blocking area to generate a spatial attribute restriction layer;

[0061] S5: Spatially align the spatial attribute restriction layer with the structure boundary set, analyze the direction vector angle of the overlapping area of the mountainous terrain and the valley region, update the boundary label, determine the rationality of the execution direction by Dijkstra algorithm, and generate the attribute path structure set available in three-dimensional geological modeling.

[0062] The frequency energy gradient layer includes the main frequency response value distribution, the gray frequency change section, and the frequency amplitude logarithmic transformation. The structure boundary set specifically refers to the boundary line segment spatial sequence, the slope response range section, and the continuous boundary frequency window label. The direction classification information includes the boundary normal vector set, the main direction angle interval, and the direction consistency grouping label. The spatial attribute restriction layer specifically refers to the attribute gradient direction value distribution, the direction angle partition boundary line, and the blocking area boundary point number. The attribute path structure set includes the shortest path length value, the cumulative turning angle degree, and the continuous boundary section node sequence.

[0063] Referring to Figure 1 , obtain the image data of the mountainous terrain and the valley region, extract the gray sequence of the region image and decompose the frequency, construct the frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude by DFT function and alignment, and the steps of generating the frequency energy gradient layer are specifically as follows:

[0064] S101: Obtain the image data of the mountainous terrain and the valley region in the three-dimensional geological model and divide it into multiple sub-regions. For each sub-region image, extract the pixel value and convert it into a gray sequence. Adjust the range of gray values by linear mapping to generate a sub-region gray sequence.

[0065] Obtain the terrain data by DEM digital elevation model, determine the terrain boundary by using Canny edge detection algorithm, set the sliding window step to 50 pixels to implement grid division, including dividing a 1024x768 pixel image into 256 sub-regions, for the valley region, use the self-matching segmentation algorithm to automatically divide the sub-region boundary when the elevation difference of adjacent pixels exceeds 5 meters, when extracting the RGB color space data of multiple sub-regions, use the red channel as the key analysis object (because the iron component has high reflectivity in the R band), normalize the sub-region R channel [123, 135], [118, 129], set the lower limit of the geological feature interval to 100 (the critical value of bedrock reflection) and the upper limit to 200 (the saturation value of vegetation coverage), when the pixel value P satisfies 100≤P≤200, perform linear transformation, assign P<100 pixels to 0 (including pixel 95 to 0 after processing), and assign P>200 pixels to 1 (including pixel 205 to 1), including the conversion process of the four feature points shown in Table 1, and the processed gray sequence is stored in a GeoTIFF format file. Each file header includes coordinate system parameters (longitude 112.35°, latitude 28.47°, and elevation reference surface WGS84).

[0066] Table 1 Typical pixel gray scale mapping table

[0067] Original coordinates Pixel value Mapping result (15，28) 123 0.23 (32，45) 88 0.00 (57，63) 192 0.92 (89，12) 210 1.00

[0068] As shown in Table 1, the threshold setting is based on the spectral analysis of 50 borehole core samples in the area, and the statistical average of the bedrock area is 98±15, and the average of the vegetation coverage area is 203±22. When the mapping value of the continuous 5 pixel points is >0.8, it is marked as a potential landslide body.

[0069] S102: According to the gray sequence of the sub-area, the frequency components in the gray sequence are decomposed and the main frequency amplitude is extracted by the fast Fourier transform function, the DFT function is called to calculate the logarithmic amplitude of the main frequency amplitude and is aligned with the position parameter, and the main frequency amplitude logarithmic transformation result is generated;

[0070] The FFTW library is used for fast Fourier transform, the sampling frequency is set to 100Hz (corresponding to the unmanned aerial vehicle surveying speed of 5m / s), the 512-point sequence is windowed (Hanning window coefficient 0.08), and the direct current component (0Hz component) is excluded when calculating the frequency domain modulus. The amplitude of the single sub-area spectrum increases sharply in the 50-75Hz interval, the arithmetic average μ=0.32 and the standard deviation σ=0.12 are calculated, the main frequency screening threshold is set to μ+3σ=0.68, the amplitude corresponding to the frequency 55Hz is extracted 0.73, the natural logarithm is calculated ln(0.73)=-0.3147, the result is bound and stored with the three-dimensional coordinates (X:102345, Y:283746, Z:356), when the elevation change rate ΔZ / ΔX>8% (from 356m to 385m across 50m horizontal distance), ΔZ represents the vertical coordinate, ΔX represents the longitude coordinate, the logarithmic amplitude weighting coefficient increases by 0.15, and the 256×256 amplitude distribution is generated by the spatial interpolation algorithm.

[0071] S103: Call the main frequency amplitude logarithmic transformation result, calculate the energy gradient value corresponding to the multi-position frequency and arrange it, and obtain the frequency energy gradient layer;

[0072] The Sobel operator is used to calculate the energy gradient, the horizontal direction convolution kernel is [-1, 0, 1; -2, 0, 2; -1, 0, 1], and the vertical direction is [-1, -2, -1; 0, 0, 0; 1, 2, 1], the non-maximum suppression is performed on the gradient amplitude, when the gradient difference of adjacent pixels exceeds 0.15 (including the second point difference of [0.12, 0.28, 0.17] sequence reaches 0.11 and 0.11), a Gaussian filter with a standard deviation σ = 1.2 is applied for smoothing, and the processed gradient is classified at an interval of 0.05, wherein the 0.25-0.30 level corresponds to the fault fracture zone (accounting for 12.3%), and the >0.30 level is the active tectonic zone (accounting for 6.8%), and the output layer uses RGBA encoding, the red channel represents the gradient intensity (0-255 corresponds to 0-0.35 gradient value), and the alpha channel marks the confidence (regions with a confidence of more than 80% are verified by three times of sampling).

[0073] 4. The multi-source geological data processing method of the three-dimensional geological model according to claim 1, characterized in that the step of generating the structural boundary set by sorting the frequency energy gradient values and judging the boundary position is specifically:

[0074] S201: Obtain the frequency energy gradient values of multiple pixel positions in the frequency energy gradient layer, filter the frequency of the mountain slope area and the river alluvial zone layer, and generate a continuous response aggregated value interval by position aggregation;

[0075] Based on the frequency energy gradient layer data, set the gradient threshold value to 0.25 (according to the statistics of the first 1000 sampling points, the area exceeding 12.3%), and use the sliding window method (window size 15x15 pixels) to scan the mountain slope area, and mark it as a candidate area when the gradient point number ratio in the window is >60%, including 23 points detected by a single window with gradient values between 0.28 and 0.35 (total 25 points), after meeting the conditions, extract all gradient values in the region, calculate the weighted aggregate value A, when a single point d i = 30m w i = 0.5488, d i is the distance from the pixel to the center of the window, if the three gradient values are 0.32, 0.28, and 0.30, the corresponding weights are 0.5488, 0.3012, and 0.2231, then the aggregate value A = (0.32x0.5488+0.28x0.3012+0.30x0.2231) / (0.5488+0.3012+0.2231) = 0.302, as shown in Table 2, the aggregate calculation results of three typical regions, generating an aggregate value interval of 0.18-0.42.

[0076] Table 2 Aggregate Value Calculation Table

[0077] Region number Effective point number Maximum gradient value Aggregation result SLOPE-05 18 0.38 0.33 RIVER-12 22 0.41 0.29 SLOPE-09 15 0.35 0.25

[0078] As shown in Table 2, the aggregation radius of 50 meters is set according to the UAV aerial survey resolution (0.2 meters / pixel), and when the difference between the aggregation values of three consecutive windows is <0.05, they are merged into the same response area.

[0079] S202: According to the interval of the continuous response aggregation value, the decreasing amplitude of the frequency of the edge position of the mountain slope area and the river alluvial zone is calculated, the boundary contour map is constructed combined with the geometric shape of the edge of the response area, and the boundary trend is judged to obtain the continuous response boundary trend coefficient;

[0080] The curvature change amount is calculated every 5 meters along the boundary Wherein represents the i-th boundary vector x i , represents the i-th boundary vector y i : The numerator part in the coordinate dot product operation is the dot product of the two vectors, and the denominator is the modulus length product of the two vectors, including the boundary point sequence (102, 356)→(105, 358)→(108, 359), the vectors are (3, 2) and (3, 1), the dot product is 3×3+2×1=11, the modulus length and Then the curvature change amount Δθ=arccos(11 / (3.605×3.162))=24.1°, when the curvature change amount of three consecutive points is greater than 15°, it is marked as a trend mutation point, and the trend coefficient of a single river bank boundary is calculated N smooth The number of segments of the smoothed boundary, θ max The maximum estimated angle change sum, N total The total number of segments, ∑|Δθ| The sum of the absolute values of all angle changes, when the proportion of smoothed segments is 70% (N smooth =14 / 20), the angle sum is 85° (the maximum estimated value is 20×180°=3600°), then C=0.7×0.7+0.3×85 / 3600=0.490+0.007=0.497, greater than 0.45 is determined as a tectonic activity boundary.

[0081] S203: Call the continuous response boundary trend coefficient to screen the trend change rate of all boundary points, the boundary line segment with a smooth change rate is selected, the selected boundary line segment is numbered and divided, and a structural boundary set is established;

[0082] The change rate of the adjacent boundary point direction is calculated, where As=5 meters is the sampling interval, when the change rate of the adjacent boundary point direction <5° / m is considered as a flat section, including three point direction angle sequences [45°, 47°, 49°], the change rates are 2° / 5m=0.4° and 2° / 5m=0.4° respectively, after meeting the conditions, it is merged into a line segment numbered BOUNDARY-07, the storage structure is a linked list node {start point coordinates (1024, 768), end point (1056, 792), length 32 meters, average change rate 0.35° / m}, a micro-layer structure set including 28 boundaries is generated, where the longest boundary section reaches 125 meters (BOUNDARY-19), and the average change rate is 2.8° / m.

[0083] Please refer to Figure 1 Based on the structure boundary set, the normal vector set of the boundaries in the mountain slope area and the river alluvial zone is constructed by principal component analysis method, the angle relationship between the vectors is analyzed and the direction consistent boundary section is identified, and the steps of generating the direction classification information are as follows:

[0084] S301: Based on the structure boundary set, the boundary point set in the mountain slope area and the river alluvial zone is extracted, the boundary direction vector is constructed for the continuous boundary point set, and the regional boundary direction vector set is generated by summarizing according to the regional attribution;

[0085] The boundary line sections numbered BOUNDARY-01 to BOUNDARY-28 are extracted from the structure boundary set, the boundary points are sampled at intervals of 5 meters, the coordinate sequence including the sampling points 1024, 768, 1027, 770, 1030, 772 of BOUNDARY-05 is obtained, the direction vector is calculated by calculating the coordinate difference of adjacent points, in the example, the vectors between two points are (3, 2) and (3, 2), the direction angle is calculated by the ratio of x coordinate and y coordinate, when the x coordinate=3.2 meters and the y coordinate=1.8 meters, the direction angle θ=arctan(1.8 / 3.2)×(180 / π)=29.4°, the special horizontal movement is uniformly marked as 0 degree, and is stored according to the geomorphic type, 58 vectors are accumulated in the mountain slope area, and 42 vectors are accumulated in the river alluvial zone, as shown in the typical data in Table 3, the storage format is a JSON object including the coordinate difference, the direction angle and the boundary number field.

[0086] Table 3 Boundary direction vector data table

[0087] Boundary number Horizontal difference Vertical difference Direction angle SL-05-01 3.2 1.8 29.4 RV-12-03 2.1 -0.5 346.6

[0088] As shown in Table 3, the direction angle accuracy is retained to one decimal place, when the angle fluctuation of three consecutive vectors is less than 2 degrees, it is merged into the same direction section.

[0089] S302: Call the regional boundary direction vector set, extract the first principal axis direction vector corresponding to the principal component, calculate the angle measure value between the first principal axis direction vector and the boundary direction vector by principal component analysis, and determine the vector consistency relationship, and generate the boundary direction angle discrete value interval;

[0090] Perform principal component analysis on the 58 vectors of the mountain slope area, m represents the total number of vectors, and 58 effective vectors are generated for the boundary segments of the mountain slope area, u k Obtained by covariance characteristic decomposition, calculate the vector mean When constructing the covariance, take the i-th vector The deviation from the mean is (-0.014, 0.282), and the covariance elements C 11 = 0.084, C 12 = 0.019, C 22 = 0.023, and the first principal axis vector u k = (0.894, 0.447) is obtained by characteristic decomposition. When calculating the angle measure, take the projection variance coefficient σ = 0.35, which is obtained by calculating the projection value standard deviation of 20 groups of samples, and the orthogonal adjustment factor η = 0.8. Through 50 times of iteration test, it is determined that it can minimize the within-class variance. Take vector as an example, calculate the dot product with the principal axis 0.6x0.894+0.8x0.447=0.895, and the cross product modulus |0.894x0.8-0.447x0.6|=0.429, v j represent the j-th original direction vector in the regional boundary direction vector set, and when the modulus of the 58 vectors is 1, the denominator item · represents the vector dot product operation, and × represents the vector cross product operation. Substituting into the formula gives After traversing all vectors, the discrete value interval [1.85, 3.12] is obtained.

[0091] S303: Call the boundary direction angle discrete value interval, and screen the angle consistency of the multi-boundary segment. The boundary segments with an angle less than the consistency judgment reference value are classified into the same direction classification number and are respectively marked, and the direction classification information is generated;

[0092] The consistency judgment reference value is set based on the direction consistency feature that the angle cosine value between the principal axis direction vectors is above 0.97;

[0093] The consistency judgment reference value is set to 14 degrees (corresponding to a cosine value of 0.97), and the threshold is set according to the angle distribution of 100 known stable boundary segments (95% data points are less than 14 degrees). The average angle of the angle sequence [12°, 15°, 13°] of the BOUNDARY-07 segment is 13.3 degrees, and the standard deviation is 1.08 degrees, which meets the threshold And the vector standard deviation σ θ When the angle is less than 2°, it is classified as CLASS-03. The boundary segment IDs [07, 09, 15] included in this category are stored. The average angle is 13.2 degrees and the length interval is [35 meters, 82 meters]. Six direction categories are generated, of which the longest category includes 9 boundary segments with an average angle deviation of 1.8 degrees.

[0094] See also Figure 1 Based on the directional classification information, the physical parameter gradient direction of multiple points is estimated, the directional angle area is determined, and the blocking area is screened. The specific steps for generating the spatial attribute restriction layer are as follows:

[0095] S401: Based on the direction classification information, calculate the spatial angles between the multiple points and the main axis direction vectors of the corresponding boundary segments. Based on the angle values ​​and the relative position relationship between the points on the boundary segments, obtain and record the direction projection vectors of the points, and generate a point direction angle sequence.

[0096] Extract multi-classification principal axis vectors based on direction classification information CLASS-01 to CLASS-06 (such as CLASS-03 spindle ), perform spatial angle calculation on the 15 points of boundary segment BOUNDARY-07, take the coordinates of point 5 (x5, y5) = (1032, 775), calculate the relative position ratio to the starting point of the boundary 0.333, and obtain the principal axis vector of the boundary segment to which the point belongs Calculate the point direction projection vector in is the local direction vector of the point, and the calculated Projection vector The recorded angle θ5 = 12.3°, as shown in Table 4, is calculated for three typical points, generating a sequence of 128 angles in the range of [8.5°, 23.7°].

[0097] Table 4 Example of point angle calculation

[0098] Point number Local vector Projection vector Angle P-07-05 (2，1) (2.02，0.97) 12.3 P-09-12 (1，3) (1.15，0.55) 19.8 P-15-08 (-1，2) (-0.45，0.22) 28.5

[0099] As shown in Table 4, the standard deviation of the included angles of points with relative position ratio r>0.5 is 37% higher than that of points with r≤0.5.

[0100] S402: Based on the point direction angle sequence, continuous point groups are screened and the directional gradient change trend between adjacent points is determined. The gradient direction change rate between continuous directional mutation points is calculated to generate a square gradient change rate interval.

[0101] In the angle sequence, screen more than 5 consecutive points, calculate the gradient change of adjacent points, when the gradient change of adjacent points of three consecutive points > 2° / m, mark the mutation section, including the gradient change rate of point sequence [12.3°, 15.6°, 19.1°, 24.5°, 28.2°] is [0.66, 0.70, 1.08, 0.74]° / m, wherein the gradient change of adjacent points of 3-4 points = 1.08 > 1.0 threshold, record the interval as the mutation section, the change rate interval of all data is [0.15, 3.42]° / m, wherein the abnormal section > 2.0° / m accounts for 9.7%.

[0102] S403: According to the direction gradient change rate interval, the point group with severe direction mutation and discontinuous direction distribution is spatially clustered, the point concentration distribution area is located, and a spatial attribute restriction layer is generated;

[0103] Using DBSCAN clustering algorithm, 24 point positions with mutation rate > 2.0° / m are spatially analyzed, the Euclidean distance between point positions is calculated, including the distance between point P-07-05(1032, 775) and P-07-06(1035, 778) m, when 4 points are included in the neighborhood, a cluster CLUSTER-02 is formed, the center coordinates (1038, 780) are recorded, the coverage radius is 8.5 meters, 6 spatial clustering areas are generated, and the maximum cluster includes 9 point positions, and the average mutation rate is 2.8° / m.

[0104] Please refer to Figure 1 The spatial attribute restriction layer and the structure boundary set are spatially aligned, the direction vector included angle of the mountainous topography and the valley area overlapping area is analyzed, the boundary label is updated, the execution direction rationality is judged by Dijkstra algorithm, and the steps of generating the attribute path structure set available in three-dimensional geological modeling are as follows:

[0105] S501: Based on the spatial attribute restriction layer and the structure boundary set, the spatial alignment is performed, the difference value of the start and end points of the boundary line segment in the corresponding space of the two is calculated, and the spatial position of the boundary line segment is adjusted, the spatial alignment coordinate offset value is generated;

[0106] The spatial attribute restriction layer CLUSTER-02 (center point coordinates 1038, 780, 152) and the boundary set BOUNDARY-07 (start point 1024, 768, 148) are aligned in three-dimensional coordinates, the offset amounts of X / Y / Z axes are calculated as 14 m, 12 m and 4 m respectively, and the boundary points are gradually displaced. Taking point P-07-03 (original coordinates 1027, 770, 149) as an example, the point distance start path length s = 9.8 m (total length L = 32 m), and the adjustment amount is calculated: Δx, Δy represent the longitude and latitude direction adjustment amount, Δz represents the vertical direction adjustment amount, the coordinates become (1027+4.29, 770+3.68, 149+1.23) = (1031.29, 773.68, 150.23), and the offset table (table 5) is generated by traversing all boundary points, wherein the maximum height adjustment occurs at point P-07-15 (Δz = 2.81m), and statistics show that 82% of the point positions have a plane displacement amount in the range of 5-15m.

[0107] Table 5: Three-dimensional coordinate offset statistical table

[0108] Coordinate axis Minimum offset (m) Maximum offset (m) Average offset (m) X 3.2 14.0 8.7 Y 2.1 12.0 6.3 Z 0.5 2.8 1.4

[0109] As shown in table 5, the Z-axis offset is positively correlated with the terrain slope (correlation coefficient r = 0.78), and the slope increases by 1°, the height adjustment amount increases by 0.23m through least square fitting.

[0110] S502: Call the spatial alignment coordinate offset value, extract the boundary segment main axis direction vector of the coincident section of the mountainous terrain area and the valley area after alignment, calculate the angle cosine value between adjacent boundary segment main axis vectors and screen the area, and obtain the direction angle difference value interval;

[0111] Extract the main axis vector of BOUNDARY-07 segment after alignment And its adjacent segment BOUNDARY-19 main axis Calculate the three-dimensional direction similarity: Set the threshold cosθ th = 0.950, when the fitting degree is lower than, mark it as a direction abnormal segment, and statistics show that among the 128 boundary segments, 23 segments (accounting for 18.0%) are marked, the average length of the abnormal segments is 15.3m, and the key is distributed in the area with a slope change > 8°.

[0112] S503: According to the direction angle difference value interval, calculate the shortest path length of multiple nodes in the direction difference section by Dijkstra algorithm, construct the priority path and update the boundary segment direction label through the path node sequence, and generate the attribute path structure set available in three-dimensional geological modeling;

[0113] Construct a graph structure including 28 nodes, and take node 3 (1031, 771, 150) to node 7 (1035, 775, 151) as an example: ΔH is extracted through DEM data, and the sampling interval is 5m, ΔH = 1m, ΔH = 151-150 = 1m, The three-dimensional space Euclidean distance module length of representative node a and node b, Representing the sliding average direction angle, sliding window: 7-meter radius sphere, including points: ≥5 valid nodes, selecting the direction angle of 6 nodes within 7 meters around node a [0.175, 0.182, 0.178, 0.190, 0.185, 0.180], θ ab Representing the instantaneous direction angle, Calculate the dot product: At the same time, calculate the radian value At the same time, substitute the formula, Through Dijkstra algorithm iteration calculation, determine that the node 3→7→12→15 is the optimal path, the total cost is 34.7m, update the direction label to NE35°, the path passes through 3 abnormal sections, the average direction deviation is 4.2°, update the boundary section direction label, and generate the available attribute path structure set in three-dimensional geological modeling.

[0114] Please refer to Figure 2 A multi-source geological data processing system of a three-dimensional geological model, the multi-source geological data processing system of the three-dimensional geological model is used to execute the multi-source geological data processing method of the three-dimensional geological model, and the system comprises:

[0115] A layer construction module, which is used to acquire mountainous topography and valley region image data, extract a gray sequence of the region image and decompose a frequency, construct a frequency window to extract a main frequency amplitude, perform logarithmic transformation on the main frequency amplitude through a DFT function and align, generate a frequency energy gradient layer and deliver to a boundary extraction module;

[0116] A boundary extraction module, which is used to screen a mountain slope area and a river alluvial zone frequency continuous response region through the frequency energy gradient layer, sort the frequency energy gradient values and judge a boundary position, generate a structure boundary set and deliver to a direction identification module;

[0117] A direction identification module, which is used to construct a normal vector group of the boundaries of the mountain slope area and the river alluvial zone through the structure boundary set through a principal component analysis method, analyze the included angle relationship between the vectors and identify a direction consistent boundary section, generate direction classification information and deliver to a gradient estimation module;

[0118] A gradient estimation module, which is used to estimate a property parameter gradient direction of multiple points through the direction classification information, judge a direction included angle region and screen a blocking area, generate a spatial attribute restriction layer and deliver to a path generation module;

[0119] A path generation module, which is used to perform spatial alignment on the spatial attribute restriction layer and the structure boundary set, analyze a direction vector included angle of a mountainous topography and valley region overlapping region, update a boundary label, judge the direction rationality through a Dijkstra algorithm, and generate an available attribute path structure set in three-dimensional geological modeling.

[0120] It should be understood that the term "and / or" in this document is only used to describe the association relationship of the associated objects, which means that there can be three relationships, for example, A and / or B can represent the following three cases: A exists alone, A and B exist together, and B exists alone, where A and B can be singular or plural. In addition, the character " / " in this document generally represents an "or" relationship between the associated objects before and after it, but it can also represent an "and / or" relationship. The specific meaning can be understood according to the context before and after it.

[0121] In this application, "at least one" means one or more, and "multiple" means two or more. "At least one of the following" or similar expressions means any combination of these items, including any combination of single or multiple items. For example, at least one of a, b, or c can mean a, b, c, a-b, a-c, b-c, or a-b-c, where a, b, and c can be single or multiple.

[0122] It should be understood that in various embodiments of the present application, the size of the sequence number of the above-mentioned multiple processes does not mean the order of execution, and the execution order of the multiple processes should be determined according to its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present application.

[0123] Those of ordinary skill in the art can realize that the units and algorithm steps of the examples described in combination with the embodiments disclosed herein can be realized in electronic hardware or a combination of computer software and electronic hardware. Whether the functions are executed in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present application.

[0124] Those skilled in the art can clearly understand that, for the convenience and brevity of description, the specific working process of the above-mentioned devices, apparatuses and units can refer to the corresponding process in the foregoing method embodiments, which will not be repeated here.

[0125] In several embodiments provided by the present application, it should be understood that the disclosed devices, apparatuses and methods can be implemented in other manners. For example, the embodiments of the apparatus described above are merely schematic. For example, the division of the units is only a logical function division. There can be another division manner for the actual implementation, for example, multiple units or components can be combined or integrated into another device, or some features can be ignored or not executed. In addition, the displayed or discussed mutual couplings or direct couplings or communication connections between different units, can be indirect couplings or communication connections through some interfaces, devices or units, and can be in electrical, mechanical or other forms.

[0126] The units described as separated components can or can not be physically separated, and the components displayed as units can or can not be physical units, i.e., can be located in one place, or can be distributed on a plurality of network units. Some or all of the units can be selected according to actual needs to achieve the purposes of the embodiments of the present application.

[0127] In addition, the multifunctional units in the embodiments of the present application can be integrated in one processing unit, or can be physically separated units, or two or more units can be integrated in one unit.

[0128] If the functions are realized in the form of software function units and sold or used as independent products, they can be stored in a computer readable storage medium. Based on this understanding, the technical solutions of the present application essentially or the parts of the prior art that make contributions or parts of the technical solutions can be embodied in the form of a software product. The computer software product is stored in a storage medium, and includes a plurality of instructions for causing a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods described in the embodiments of the present application. The aforementioned storage medium includes: a U disk, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and various media that can store program codes.

[0129] The above description is only a specific implementation of the present application, but the protection scope of the present application is not limited thereto. Any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the present application, which should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A method for processing multi-source geological data of a three-dimensional geological model, characterized in that: The method comprises: S1: Acquire image data of mountainous landforms and river valley areas, extract the grayscale sequence of the regional image and decompose the frequency, construct a frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude through the DFT function and align them to generate a frequency energy gradient layer; S2: screening the frequency continuous response areas of mountain slope areas and river alluvial belts through the frequency energy gradient layer, sorting the frequency energy gradient values ​​and determining the boundary positions to generate a structural boundary set; S3: Based on the structural boundary set, a normal vector group is constructed for the boundaries of the mountain slope area and the river alluvial zone through principal component analysis, the angle relationship between the vectors is analyzed, and the boundary segments with consistent directions are identified to generate direction classification information; S4: estimating the physical parameter gradient direction of multiple points based on the direction classification information, determining the direction angle area and screening the blocking area, and generating a spatial attribute restriction layer; S5: Spatially align the spatial attribute restriction layer with the structural boundary set, analyze the direction vector angle of the overlapping area of ​​the mountain landform and the river valley area, update the boundary label, judge the rationality of the execution direction through the Dijkstra algorithm, and generate the attribute path structure set available in 3D geological modeling.

2. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: The frequency energy gradient layer includes the main frequency response value distribution, grayscale frequency change segment, and frequency amplitude logarithmic transformation. The structural boundary set specifically includes the boundary line segment space sequence, slope response range segment, and continuous boundary frequency window label. The direction classification information includes the boundary normal vector set, the main direction angle interval, and the direction consistency grouping label. The spatial attribute restriction layer specifically refers to the attribute gradient direction value distribution, the direction angle partition boundary, and the blocking area boundary point number. The attribute path structure set includes the shortest path length value, the cumulative turning angle, and the continuous boundary segment node sequence.

3. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: Obtain image data of mountainous landforms and river valley areas, extract the grayscale sequence of the regional image and decompose the frequency, construct a frequency window to extract the main frequency amplitude, perform logarithmic transformation on the main frequency amplitude through the DFT function and align it, and generate the frequency energy gradient layer in the following steps: S101: Obtain image data of mountain landforms and river valley regions in a three-dimensional geological model and divide the data into multiple sub-regions. For each sub-region image, extract pixel values ​​and convert them into a grayscale sequence. Adjust the grayscale value range through linear mapping to generate a sub-region grayscale sequence. S102: Decomposing the frequency components in the grayscale sequence of the sub-region by using a fast Fourier transform function and extracting the main frequency amplitude, calling a DFT function to calculate the logarithmic amplitude of the main frequency amplitude and aligning it with the position parameter to generate a logarithmic transformation result of the main frequency amplitude; S103: Calling the logarithmic transformation result of the main frequency amplitude, calculating the energy gradient values ​​corresponding to the frequencies at multiple positions and arranging them to obtain a frequency energy gradient layer.

4. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: The frequency energy gradient layer is used to screen the frequency continuous response areas of mountain slope areas and river alluvial belts, sort the frequency energy gradient values ​​and determine the boundary positions, and the steps of generating the structural boundary set are as follows: S201: Obtain frequency energy gradient values ​​at multiple pixel positions in the frequency energy gradient layer, filter the frequencies of the mountain slope area and river alluvial zone layers, perform position aggregation, and generate a continuous response aggregation value interval; S202: Calculating the decreasing amplitude of the edge position frequency of the mountain slope area and the river alluvial zone area according to the continuous response aggregation value interval, constructing a boundary contour map based on the geometric shape of the edge of the response area, and determining the boundary direction to obtain a continuous response boundary direction coefficient; S203: calling the continuous response boundary trend coefficient, screening boundary line segments with a gentle change rate for the trend change rate of all boundary points, numbering and dividing the selected boundary line segments, and establishing a structural boundary set.

5. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: Based on the structural boundary set, the principal component analysis method is used to construct a normal vector group for the boundaries of the mountain slope area and the river alluvial zone. The angle relationship between the vectors is analyzed and the boundary segments with consistent directions are identified. The specific steps for generating direction classification information are as follows: S301: extracting a boundary point set within the mountain slope area and the river alluvial zone based on the structural boundary set, constructing boundary direction vectors for the continuous boundary point set and summarizing them according to regional affiliation to generate a regional boundary direction vector set; S302: Calling the region boundary direction vector set, extracting the first principal axis direction vector corresponding to the principal component, calculating the angle between the first principal axis direction vector and the boundary direction vector by principal component analysis, determining the vector consistency relationship, and generating a discrete value interval of the boundary direction angle; S303: Calling the boundary direction angle discrete value interval to screen the angle consistency of multiple boundary segments, classifying the boundary segments with angles smaller than the consistency judgment reference value into the same direction classification number and marking them respectively, thereby generating direction classification information; The consistency judgment reference value is set based on the directional consistency feature exhibited when the cosine value of the angle between the main axis direction vectors is greater than 0.

97.

6. The multi-source geological data processing method for a three-dimensional geological model according to claim 5, characterized in that: The angle value θ between the first principal axis direction vector and the boundary direction vector is calculated by the principal component analysis method. i , using the formula: Among them, u k represents the kth principal axis vector extracted by principal component analysis, v i represents the i-th original direction vector in the regional boundary direction vector set, σ represents the projection variance coefficient of the boundary vector set, which is 0.2-0.5, η represents the orthogonal component adjustment factor, which is 0.1-1.2, p represents the number of principal axes selected by principal component analysis, m represents the total number of boundary direction vectors, and v j Represents the jth original direction vector in the region boundary direction vector set.

7. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: The steps of estimating the physical parameter gradient directions of multiple points based on the direction classification information, determining the direction angle area and screening the blocking area, and generating the spatial attribute restriction layer are as follows: S401: Based on the direction classification information, calculate the spatial angles between multiple points and the main axis direction vectors of the corresponding boundary segments, obtain and record the direction projection vectors of the points based on the angle values ​​and the corresponding positions of the points on the boundary segments, and generate a point direction angle sequence; S402: Based on the point direction angle sequence, continuous point groups are screened and the directional gradient change trend between adjacent points is determined, the gradient direction change rate between continuous directional mutation points is calculated, and a directional gradient change rate interval is generated; S403: Based on the directional gradient change rate interval, spatially cluster the point groups with drastic directional changes and discontinuous directional distribution, locate the concentrated distribution area of ​​the points, and generate a spatial attribute restriction layer.

8. The multi-source geological data processing method for a three-dimensional geological model according to claim 1, characterized in that: The steps of spatially aligning the spatial attribute restriction layer with the structural boundary set, analyzing the direction vector angle of the overlapping area of ​​the mountain landform and the river valley area, updating the boundary label, and judging the rationality of the execution direction by the Dijkstra algorithm to generate the attribute path structure set available in the three-dimensional geological modeling are as follows: S501: spatially aligning the layer and the structure boundary set based on the spatial attribute restriction, calculating the difference in the three-dimensional coordinates of the starting and ending points of the boundary segments in the corresponding spaces of the two, adjusting the spatial positions of the boundary segments, and generating a spatial alignment coordinate offset value; S502: Calling the spatial alignment coordinate offset value, extracting the principal axis direction vector of the boundary segment of the overlapped section between the mountainous landform area and the river valley area after alignment, calculating the cosine value of the angle between the principal axis vectors of adjacent boundary segments and filtering the area to obtain the direction angle difference value interval; S503: According to the direction angle difference value interval, the shortest path length of multiple nodes in the direction difference segment is calculated by Dijkstra algorithm, a priority path is constructed, and the boundary segment direction label is updated through the path node sequence to generate an attribute path structure set available in three-dimensional geological modeling.

9. The multi-source geological data processing method for a three-dimensional geological model according to claim 8, characterized in that: The shortest path length of multiple nodes in the direction difference section is calculated by the Dijkstra algorithm Using the formula: in, Represents the three-dimensional Euclidean distance modulus between node a and node b, Subscripts x and y represent the grid numbers in the longitude and latitude directions, subscript z represents the geological layer number in the vertical direction, ΔH is the elevation difference between nodes a and b, θ ab Represents the instantaneous angle in radians, which is obtained by obtaining the node coordinate data through a 3D laser scanner. It represents the sliding mean of the azimuth angle and is set to 5-10 meters for the geological local area.

10. A multi-source geological data processing system for a three-dimensional geological model, characterized in that: The system is used to implement the multi-source geological data processing method for a three-dimensional geological model according to any one of claims 1 to 9, and the system comprises: The layer construction module is used to obtain image data of mountainous landforms and river valley areas, extract the grayscale sequence of the regional image and decompose the frequency, construct a frequency window to extract the main frequency amplitude, perform logarithmic transformation and alignment on the main frequency amplitude through the DFT function, generate a frequency energy gradient layer and pass it to the boundary extraction module; A boundary extraction module is used to screen the frequency continuous response areas of mountain slope areas and river alluvial belts through the frequency energy gradient layer, sort the frequency energy gradient values ​​and determine the boundary positions, generate a structural boundary set and pass it to the direction recognition module; A direction recognition module is used to construct a normal vector group for the boundaries of the mountain slope area and the river alluvial zone through the principal component analysis method using the structural boundary set, analyze the angle relationship between the vectors and identify the boundary segments with consistent directions, generate direction classification information and transmit it to the gradient estimation module; A gradient estimation module is used to estimate the gradient direction of the physical property parameters of multiple points based on the direction classification information, determine the direction angle area and filter the blocking area, generate a spatial attribute restriction layer and pass it to the path generation module; The path generation module is used to spatially align the spatial attribute restriction layer with the structural boundary set, analyze the direction vector angle of the overlapping area of ​​mountain landforms and river valley areas, update the boundary label, determine the rationality of the execution direction through the Dijkstra algorithm, and generate an attribute path structure set available in three-dimensional geological modeling.

Images