Drilling construction intelligent decision support and risk early warning method, device and equipment

Abstract

The invention discloses a drilling construction intelligent decision support and risk early warning method, device and equipment, and the method comprises the steps: collecting drilling sound wave response data, carrying out harmonic analysis to recognize a reservoir resonance node, and constructing a reservoir energy conduction network; a fluctuation propagation path is tracked based on an energy conduction network, and a construction risk area is identified through stress wave spectrum features; porosity and permeability data are extracted for phase difference analysis, a seepage-construction mismatch factor is generated, and a construction optimization window is actively generated; constructing a geological potential energy gradient field to determine an overflow influence range, and performing tensor coupling on multi-dimensional risk factors to generate a three-dimensional construction decision matrix; an optimal construction scheme is generated through energy minimization and risk-efficiency dynamic balance processing; construction parameter response changes are monitored in real time, abnormal parameter combination symptoms are recognized, a grading risk early warning instruction is generated, and intelligent decision support and risk early warning of drilling construction are achieved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of drilling engineering technology, in particular to a drilling construction intelligent decision support and risk early warning method, device and equipment. BACKGROUND

[0002] Drilling construction is an important technical means of underground engineering construction, and its construction quality and safety directly affect the success of the project. In actual construction, the internal state of the stratum changes complexly, and the existing monitoring means mainly relies on single-point measurement and local observation, which is difficult to obtain the overall response characteristics and internal connection of the stratum. At the same time, there is a complex interaction between the natural state of the stratum and the construction disturbance, and when they are not coordinated, engineering accidents are easily caused, but the existing technology lacks an effective matching degree evaluation method, and mainly relies on experience to arrange construction.

[0003] The various monitoring data generated during the construction process are independent of each other, lacking comprehensive analysis and fusion processing mechanism, leading to one-sided decision basis and difficulty in achieving global optimization. The existing early warning system mainly uses single parameter threshold judgment, which cannot identify the correlation between parameters and abnormal conduction mode, and often cannot find problems until after the accident occurs. In addition, the construction path planning mainly considers geometric constraints and geological conditions, lacking support for multi-objective optimization and dynamic adjustment, and is difficult to adapt to the complex underground environment. Therefore, a method is needed to solve at least one of the above problems. SUMMARY

[0004] The present application provides a drilling construction intelligent decision support and risk early warning method, device and equipment, which aims to build a reservoir energy conduction network model, deeply analyze the dynamic matching relationship between seepage and construction, establish a multi-dimensional tensor coupled decision matrix, realize intelligent optimization of construction trajectory and accurate early warning of risk, and form a complete drilling construction intelligent management system through multi-well collaborative control technology.

[0005] The present application provides a drilling construction intelligent decision support and risk early warning method, device and equipment, which aims to build a reservoir energy conduction network model, deeply analyze the dynamic matching relationship between seepage and construction, establish a multi-dimensional tensor coupled decision matrix, realize intelligent optimization of construction trajectory and accurate early warning of risk, and form a complete drilling construction intelligent management system through multi-well collaborative control technology. Acquire node coupling strength based on the reservoir resonance node, and construct a reservoir energy conduction network based on the node coupling strength; Track the wave propagation path based on the reservoir energy conduction network, and determine the construction risk area through the wave propagation path; Extract the porosity distribution density and permeability variation curve of each drilling block, identify the high permeability zone based on the porosity distribution density, perform phase difference analysis on the high permeability zone and the permeability variation curve to generate a seepage-construction mismatch factor, and automatically generate a construction optimization window according to the seepage-construction mismatch factor; construct a geological potential energy gradient field by using the construction optimization window, determine an overflow influence range by performing seepage diffusion analysis through the geological potential energy gradient field, perform tensor coupling on the overflow influence range, the construction risk area and the seepage-construction mismatch factor to form a coupling tensor, and generate a three-dimensional construction decision matrix based on the coupling tensor; determine an energy optimal point by performing energy minimization processing on the three-dimensional construction decision matrix, generate a first construction trajectory set according to the energy optimal point, generate a second construction trajectory set by performing risk-efficiency dynamic balance processing on the first construction trajectory set, and determine an optimal construction scheme through the second construction trajectory set; execute the optimal construction scheme and monitor parameter response changes of mud pressure, drilling speed and drilling tool torque, identify abnormal parameter combination signs based on the parameter response changes, and generate a graded risk early warning instruction based on the abnormal parameter combination signs.

[0006] The second aspect of the present application proposes a drilling construction intelligent decision support and risk early warning device, comprising: A data acquisition module is configured to acquire drilling acoustic response data by using an acoustic sensor array, identify reservoir resonance nodes by performing harmonic analysis on the drilling acoustic response data, obtain node coupling strength based on the reservoir resonance nodes, and construct a reservoir energy conduction network based on the node coupling strength. A risk identification module is configured to track wave propagation paths based on the reservoir energy conduction network and determine a construction risk area through the wave propagation paths. A parameter optimization module is configured to extract porosity distribution density and permeability variation curves of each drilling block, identify high permeability zones based on the porosity distribution density, generate a seepage-construction mismatch factor by performing phase difference analysis on the high permeability zones and the permeability variation curves, and automatically generate a construction optimization window according to the seepage-construction mismatch factor. A decision generation module is configured to construct a geological potential energy gradient field by using the construction optimization window, determine an overflow influence range by performing seepage diffusion analysis through the geological potential energy gradient field, perform tensor coupling on the overflow influence range, the construction risk area and the seepage-construction mismatch factor to form a coupling tensor, and generate a three-dimensional construction decision matrix based on the coupling tensor. A trajectory optimization module is configured to determine an energy optimal point by performing energy minimization processing on the three-dimensional construction decision matrix, generate a first construction trajectory set according to the energy optimal point, generate a second construction trajectory set by performing risk-efficiency dynamic balance processing on the first construction trajectory set, and determine an optimal construction scheme through the second construction trajectory set. A parameter early warning module is configured to execute the optimal construction scheme and monitor parameter response changes of mud pressure, drilling speed and drilling tool torque, identify abnormal parameter combination signs based on the parameter response changes, and generate graded risk early warning instructions based on the abnormal parameter combination signs.

[0007] A third aspect of the present application provides a computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor implements the steps of the drilling construction intelligent decision support and risk early warning method disclosed in the first aspect when executing the program.

[0008] The beneficial effects of the present application are embodied in the following aspects: first, the drilling acoustic response data is collected by the acoustic sensor array and harmonic analysis is performed to identify the reservoir resonance nodes and construct the energy transmission network, thereby realizing the visualization analysis of the energy transmission path in the reservoir. Based on the network, the wave propagation path is tracked, the construction risk area is identified through the analysis of the stress diffusion range and stress wave spectrum characteristics of the formation, the traditional single-point monitoring is upgraded to a networked monitoring system, and the correlation mechanism between wave propagation and risk area is established, thereby improving the comprehensiveness of reservoir feature identification and the accuracy of risk prediction. Second, the phase difference analysis is performed by extracting the high permeability zone and the permeability variation curve to establish a quantitative evaluation method of the seepage-construction mismatch factor, which can accurately reflect the matching degree of the seepage state of the reservoir and the construction rhythm. By coupling the mismatch factor with the construction risk area and the overflow influence range to generate a three-dimensional construction decision matrix, intelligent construction decision based on multi-dimensional information fusion is realized, thereby improving the safety and efficiency of drilling construction. Finally, the risk-efficiency dynamic balance processing technology is used to optimize the construction trajectory, and a multi-level risk early warning system is established in combination with the abnormal parameter combination sign identification mechanism. By monitoring the response changes of key parameters such as mud pressure, drilling speed and drilling tool torque in real time, abnormal propagation paths and combination signs are identified, thereby realizing the transition from post-processing to early warning.

[0009] It should be understood that the above general description and the following detailed description are only exemplary and explanatory, and cannot limit the present application. BRIEF DESCRIPTION OF DRAWINGS

[0010] The drawings herein show specific examples of the technical solutions described in the present application, and constitute part of the specification together with the specific embodiments, for explaining the technical solutions, principles and effects of the present application.

[0011] Unless specifically stated, the same reference signs in different drawings represent the same or similar technical features, and different reference signs may also be used to represent the same or similar technical features.

[0012] Figure 1It is a flowchart of a drilling construction intelligent decision support and risk early warning method.

[0013] Figure 2 It is a structural block diagram of a drilling construction intelligent decision support and risk early warning device.

[0014] Figure 3 It is a structural diagram of a computer device. DETAILED DESCRIPTION

[0015] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the scope of protection of the present application.

[0016] It should be noted that all directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of the present application are only used to explain the relative positional relationship, movement condition, etc. between components in a certain posture (as shown in the drawings), and if the certain posture changes, the directional indications also change accordingly.

[0017] It should also be noted that when an element is referred to as being “fixed to” or “set on” another element, it can be directly on the other element or can have a middle element. When an element is referred to as being “connected” to another element, it can be directly connected to the other element or can have a middle element.

[0018] In addition, the description of “first”, “second” and the like in the present application is only for the purpose of description, and cannot be understood as indicating or implying the relative importance of the indicated technical features or implicitly indicating the number of the indicated technical features. Therefore, the features defined as “first”, “second” can explicitly or implicitly include at least one of the features. In addition, the technical solutions of each embodiment can be combined with each other, but it must be based on the realization of a person of ordinary skill in the art, and when the combination of technical solutions contradicts each other or cannot be realized, it should be considered that the combination of technical solutions does not exist, and is not within the scope of protection claimed by the present application.

[0019] The technical solutions of the embodiments of the present application will be introduced below.

[0020] As shown in Figure 1 The embodiments of the present application provide a drilling construction intelligent decision support and risk early warning method, which includes the following steps S110-S160: Step S110, collecting borehole acoustic response data by using the acoustic wave sensor array, identifying reservoir resonance nodes by harmonic analysis on the borehole acoustic response data, obtaining node coupling strength based on the reservoir resonance nodes, and constructing a reservoir energy conduction network based on the node coupling strength.

[0021] Specifically, the acoustic response data of the reservoir around the borehole is collected by the multi-channel acoustic wave sensor array arranged downhole. The sensor array adopts a ring distribution structure, and a high-precision piezoelectric sensor is arranged every 60 degrees, forming a 6-channel synchronous acquisition system, with a working frequency range of 0.1-20 kHz and a sampling rate of 50 kHz to ensure the accuracy of frequency domain analysis. The acoustic excitation source adopts a controllable source technology, and a vibrator driven by hydraulic pressure generates a frequency-adjustable sine wave signal, with a starting frequency of 1 Hz and a step sweep frequency of 0.5 Hz to 15 Hz, and each frequency point lasts for 3 seconds to ensure signal stability. The data acquisition system is configured with a 24-bit A / D converter, with a dynamic range of 120 dB, which can capture weak reservoir response signals. The vertical spacing of the sensor array is set to 2 meters, covering a borehole depth range of 100-1500 meters, forming a three-dimensional data grid of 700 acquisition points. Signal preprocessing includes 50 Hz power frequency filtering, high-pass filtering (cutoff frequency 0.05 Hz) and low-pass filtering (cutoff frequency 25 Hz) to eliminate environmental noise interference. Real-time compression algorithm is used for data storage, with a compression ratio of 8:1, and the daily data volume of a single well is controlled within 500 GB.

[0022] The borehole acoustic response data is analyzed by harmonic analysis to identify the reservoir resonance nodes. Fast Fourier Transform (FFT) technology is used to convert time domain signals to frequency domain for analysis, with an FFT window length of 8192 points and an overlap rate of 50% to ensure a frequency resolution of 6.1 mHz. The reservoir resonance characteristics mainly manifest as energy concentration at specific frequencies, and the resonance frequencies are usually distributed in the range of 2-12 Hz, closely related to the elastic modulus and density characteristics of the reservoir. The resonance node identification algorithm uses peak detection technology, with an amplitude threshold of 3 times the standard deviation of the background noise and an adjacent peak spacing of no less than 0.2 Hz to avoid false peak interference. The center frequency, amplitude and quality factor Q value of each detected resonance peak are determined by Gaussian fitting, and the Q value calculation formula is Q=f0 / Δf, where f0 is the resonance frequency and Δf is the half-power bandwidth. The reservoir type discrimination is based on the resonance characteristic parameters: the resonance frequency of sandstone reservoirs is concentrated in the range of 3-6 Hz, and the Q value is usually between 15-25; the resonance frequency of carbonate reservoirs is high, distributed in the range of 6-10 Hz, and the Q value can reach 30-40; the resonance characteristics of shale layers are weak, with frequencies scattered in the range of 2-8 Hz.

[0023] In some embodiments, the obtaining the node coupling strength based on the reservoir resonance node comprises: establishing an energy conduction path based on the reservoir resonance node; monitoring acoustic impedance along the energy conduction path to obtain impedance data; inversely deducing the energy exchange efficiency between nodes according to the impedance data; and converting the energy exchange efficiency into the node coupling strength.

[0024] An energy conduction path is established based on the reservoir resonance node. The resonance node is taken as a network vertex, and the acoustic connection between nodes is taken as a network edge. The path weight is initially set as the Euclidean distance between nodes. The minimum spanning tree algorithm adopts the Kruskal method, and the edges are sorted from small to large according to the edge weight. The edges that do not form a loop are gradually added to ensure that all nodes are connected and the total path length is the shortest. The path validity verification is filtered by the geological constraint condition to delete the paths that pass through the disconnected areas such as faults and mudstone interlayers. The multi-path redundancy design establishes a backup path in the main conduction direction to improve the network robustness. The main conduction direction is determined based on the reservoir dip angle and the permeability anisotropy. The reservoir with horizontal permeability greater than vertical permeability mainly conducts along the bedding direction. The conduction path geometric parameters include path length, tortuosity, and average cross-sectional area, which affect the energy conduction efficiency. The path classification is set according to the importance of the connected nodes. The path connecting the strong resonance nodes is the first-level path, and the path connecting the medium nodes is the second-level path.

[0025] The acoustic impedance is monitored along the energy conduction path to obtain impedance data. The acoustic logging method is used to monitor the impedance. The measurement points are arranged along the conduction path by using a multi-channel acoustic sensor. The sensor spacing is determined according to the complexity of the path, and the measurement frequency range is 0.1 Hz-15 kHz, covering the reservoir response frequency band. The acoustic impedance is calculated by the sound propagation time and amplitude attenuation. The impedance formula is Z=ρ×v, where ρ is the rock density and v is the acoustic velocity. The monitoring point spacing is determined according to the complexity of the path. The spacing between straight sections is 20 meters, and the spacing between curved sections is 10 meters. Multi-frequency impedance measurement reveals the frequency dispersion characteristics of the reservoir. Low-frequency impedance reflects the macrostructure, and high-frequency impedance reflects the micro-pore characteristics. The data filtering processing adopts the median filtering to eliminate abnormal values, and the filtering window length is 5 data points. The impedance gradient calculation identifies the impedance mutation position. The positions with a gradient greater than the threshold value of 0.5 may exist reservoir boundaries or fractures. The statistical analysis calculates the mean, standard deviation, and coefficient of variation of the impedance on the path to quantify the distribution characteristics of the impedance.

[0026] The energy exchange efficiency between nodes is inversely deduced according to the obtained impedance data. Based on the transmission line theory, the reservoir conduction path is equivalent to a distributed parameter transmission line. The characteristic impedance Z0=√(R+jωL) / (G+jωC), where j is the imaginary unit, ω is the angular frequency, R is the resistance, L is the inductance, G is the conductance, and C is the capacitance. The transmission efficiency η=4×Z_in×Z_out / (Z_in+Z_out) ²where Z_in and Z_out are the input and output impedance, respectively. The frequency-dependent analysis reveals the transmission characteristics at different frequencies, with low-frequency transmission efficiency mainly affected by resistive loss, and high-frequency transmission also considering dielectric loss. The energy exchange efficiency between node pairs is calculated by path integration, with the integration path being the shortest conduction path connecting the two nodes. The multi-path effect analysis considers the contribution of parallel paths to the total transmission efficiency, with the total efficiency being the parallel combination of each path efficiency. The influence of reservoir heterogeneity on transmission efficiency is handled by segmentation, with the assumption of uniform medium within each segment and connection between segments through interface conditions. The least squares method is used to invert the optimal transmission efficiency parameters with the measured impedance data as the constraint condition.

[0027] The energy exchange efficiency is converted to the node coupling strength. Through a nonlinear mapping function, the energy exchange efficiency η (range 0-1) is mapped to the coupling strength C, with the mapping formula C = η^α, where α is the adjustment parameter, determined according to the reservoir type: sandstone α = 0.8, carbonate α = 1.0, shale α = 1.2. The coupling strength classification standard is set as: C > 0.8 for strong coupling, 0.5 < C ≤ 0.8 for moderate coupling, 0.2 < C ≤ 0.5 for weak coupling, and C ≤ 0.2 for extremely weak coupling. Normalization processing ensures the comparability of the coupling strength of all node pairs, with the maximum-minimum normalization method C_norm = (C - C_min) / (C_max - C_min) used. Symmetry verification C ij ji ensures the symmetry of the coupling strength matrix. Sparse processing deletes weak connections with coupling strength less than 0.1, reducing storage requirements and computational complexity. Spatial distribution analysis of coupling strength identifies strong coupling areas and weak coupling areas within the reservoir, with strong coupling areas usually corresponding to high permeability reservoirs. Time stability evaluation verifies the reproducibility of coupling strength through multi-period measurement, with a coefficient of variation less than 20% considered stable.

[0028] ​The reservoir energy conduction network is constructed based on the coupling strength of nodes. A weighted undirected graph G=(V, E, W) is used to represent the network, where V is the set of resonance nodes, E is the set of node connection edges, and W is the set of edge weights (i.e., coupling strength). During the network construction process, a coupling strength threshold of 0.2 is set, and only connections with a strength greater than the threshold are retained to form a sparse network structure. Network topology parameter calculation includes node degree, clustering coefficient, shortest path length, and centrality index, which reflect the structural characteristics and conduction performance of the network. Key node identification uses betweenness centrality analysis, and nodes with high betweenness centrality control multiple conduction paths and play a key role in the network. Community structure detection uses a modularity optimization algorithm to identify relatively independent conduction units within the reservoir. Network robustness analysis evaluates the network's resistance to damage by randomly deleting nodes and edges. Conduction path optimization uses the shortest path algorithm to calculate the optimal conduction path between any two nodes. Network hierarchical structure is divided into core layer, intermediate layer, and edge layer based on node strength and connection degree. Based on network modeling and topology analysis, the construction of the reservoir energy conduction network is finally completed.

[0029] In step S120, the wave propagation path is tracked based on the reservoir energy conduction network, and the construction risk area is determined through the wave propagation path.

[0030] Specifically, the wave propagation path is tracked using the topological structure and node connection relationship of the reservoir energy conduction network. Based on the node degree distribution and edge weight distribution in the network, high-value nodes serve as the main hub for wave propagation, and high-weight edges serve as the preferred propagation channel. The propagation path search uses an extended version of the Dijkstra shortest path algorithm, which considers both propagation distance and energy attenuation. The path cost function is Cost=Distance×exp(-α×CouplingStrength), where Distance is the distance between nodes, CouplingStrength is the average coupling strength of edges on the path (with a value range of 0-1), and α=0.8 is the attenuation coefficient. The wave source point is set to the drill bit position, and the mechanical vibration during drilling is considered as a new excitation source, generating initial disturbance energy according to the vibration excitation during drilling. The multi-path propagation analysis considers the branching and converging phenomena of waves in the network, and uses a parallel search strategy to track the top K optimal paths simultaneously. The value of K is set to 5-15 according to the complexity of the network. The propagation time calculation is based on the propagation speed of each network edge, which is determined by the elastic parameters of the reservoir medium, v=√(E / ρ), where E is the elastic modulus and ρ is the density. The energy attenuation on the path follows an exponential decay law, with E=E0×exp(-β×L) after passing through each edge, where E0 is the initial energy, β is the attenuation constant, and L is the edge length. The effectiveness of the propagation path is determined by setting an energy threshold, and the tracking of the path is terminated when the propagation energy is lower than 5% of the initial energy.

[0031] In some embodiments, the determining the construction risk area through the wave propagation path comprises: A formation stress diffusion range is determined through the wave propagation path; a stress wave frequency spectrum feature is extracted through the formation stress diffusion range; and a construction risk area is identified according to the stress wave frequency spectrum feature.

[0032] A formation stress diffusion range is determined through the wave propagation path. The spatial distribution characteristics of the propagation path are utilized: the area with high path density is set as the main stress diffusion zone, the path convergence point is marked as the stress concentration zone, and the path termination point is determined as the diffusion boundary. The formation stress diffusion modeling adopts a weighted diffusion algorithm based on the propagation path. Each propagation path corresponds to a diffusion sub-area, and the stress intensity of the sub-area is determined according to the residual energy on the path. The multi-path superposition effect is handled through the energy superposition principle. When multiple paths pass through the same area, the total stress of the area is the vector sum of the stress contributed by each path. The stress diffusion directionality is determined according to the main direction of the propagation path. An elliptical diffusion model is adopted, with the long axis direction along the main direction of the path and the short axis direction perpendicular to the main direction. The diffusion intensity decay function is established based on the path energy distribution. The farther away from the path, the faster the stress attenuation, and the attenuation coefficient is determined according to the reservoir medium characteristics. The anisotropic effect is considered by introducing the directional diffusion coefficient. The horizontal propagation velocity v h =√(E_h / ρ), and the vertical propagation velocity v_v=√(E_v / ρ), where E_h is the horizontal elastic modulus, E_v is the vertical elastic modulus, and ρ is the density. The diffusion range boundary is determined by the path termination position and the energy threshold double criteria, forming an irregular diffusion area profile.

[0033] The stress wave frequency spectrum characteristics are extracted by the stress diffusion range of the stratum. The short-time Fourier transform (STFT) technique is adopted to analyze the frequency component variation of the stress signal in the time-frequency domain synchronously. The time window length is set to 0.1 second, and the overlap rate is 75%, so as to ensure the balance of the time-frequency resolution. The spectrum feature extraction focuses on the following parameters: the main frequency fp (the frequency component with the maximum energy), the frequency band width Δf (the frequency range containing 80% of the energy), and the spectral centroid fc (the center of gravity frequency of the spectrum energy). The main frequency identification is realized by searching the peak value of the power spectral density, and the search range is limited to 2-50 Hz to avoid high-frequency noise interference. The spectral centroid calculation formula is fc=∫f×S(f)df / ∫S(f)df, where f is the frequency, and S(f) is the power spectral density function. The spatial variability analysis calculates the difference of the spectrum characteristics at different positions, and the spatial variability degree is quantified by the coefficient of variation CV=σ / μ, where σ is the standard deviation, and μ is the mean. The multi-band energy distribution analysis divides the spectrum into low frequency band (2-8 Hz), middle frequency band (8-20 Hz), and high frequency band (20-50 Hz), calculates the energy proportion of each frequency band, and the energy distribution mode of different frequency bands reflects different reservoir characteristics. The spectrum characteristic vector construction organizes each parameter into a characteristic vector F=[fp,Δf,fc,E_low,E_mid,E_high].

[0034] The construction risk area is identified according to the stress wave frequency spectrum characteristics. A multi-level discriminant model is adopted to compare and analyze the spectrum characteristics and the feature mode of the known risk area. The risk type classification includes: collapse risk (low frequency energy proportion >60%, main frequency <5 Hz), water gushing risk (middle frequency energy proportion >50%, spectral centroid >15 Hz), gas outburst risk (high frequency energy proportion >40%, frequency band width >15 Hz), and stratum crushing risk (spectral centroid deviation >20%). The risk degree quantification adopts the fuzzy comprehensive evaluation method, establishes the risk membership function, and maps the spectrum characteristic value to the risk probability. The collapse risk membership function is μ1=1 / (1+exp(5×(fp-4))), and when the main frequency fp is lower than 4 Hz, the risk probability rises rapidly. The water gushing risk is judged by the joint discrimination of the spectral centroid and the middle frequency energy, and the discrimination function is R2=w1×fc+w2×E_mid, where fc is the spectral centroid, E_mid is the middle frequency band energy proportion, the weight coefficients w1=0.6 and w2=0.4. The risk level is divided into five levels: extremely high risk (risk value >0.8), high risk (0.6-0.8), medium risk (0.4-0.6), low risk (0.2-0.4), and extremely low risk (<0.2). The comprehensive risk assessment considers the superposition effect of multiple risk types, and the total risk value R_total=√(R1 ² +R2 ² +R3 ² +R4 ²), where R1-R4 are risk values of four risk types respectively. The risk area boundary is determined by the gradient change rate of the risk, and the position where the gradient change rate is greater than the threshold 0.1 / m is marked as the risk boundary.

[0035] In step S130, the porosity distribution density and the permeability variation curve of each drilling block are extracted, the high permeability zone is identified based on the porosity distribution density, the phase difference analysis is performed by using the high permeability zone and the permeability variation curve to generate the seepage-construction mismatch factor, and the construction optimization window is automatically generated according to the seepage-construction mismatch factor.

[0036] Specifically, the porosity distribution density and the permeability variation curve of each drilling block are extracted. The porosity data is extracted by using the method of combining the nuclear magnetic resonance logging and the acoustic logging. The nuclear magnetic resonance logging provides the total porosity and the effective porosity. The acoustic logging calculates the acoustic porosity by using the Willey formula φ=(Δt-Δt_ma) / (Δt_f-Δt_ma), where Δt is the formation acoustic time difference, Δt_ma is the rock skeleton acoustic time difference, and Δt_f is the pore fluid acoustic time difference. The porosity distribution density is calculated by using the kernel density estimation method. The porosity values of each measuring point are taken as sample points, and the density distribution function is determined by statistical analysis. The permeability variation curve is extracted based on the mercury injection test and the fluid pressure test data. The relationship between the porosity and the permeability is established by using the Kozeny-Carman equation K=(φ ³ / (1-φ) ² )×(d ² / 180), where K is the permeability, φ is the porosity, and d is the average diameter of the particles. The vertical permeability profile is obtained by stratified statistics, and the horizontal direction is extended to two-dimensional distribution by using the interpolation algorithm. The permeability variation curve is fitted by using the piecewise linear method, and the time sequence variation characteristics are recorded.

[0037] Based on the extracted porosity distribution density data, the distribution of high permeability zones in the reservoir is identified. First, the threshold value is determined by statistical analysis of the porosity distribution density, and the P95 quantile (the value below 95% of the data points) is used as the high porosity threshold. The area with porosity higher than the threshold is marked as a potential high permeability zone. The porosity-permeability relationship correction is achieved by establishing a local regression model, and the model form is logK=a×logφ+b×φ+c, where a, b, and c are regression coefficients determined by least squares fitting. The permeability threshold adopts a double standard: the absolute threshold is set to 50 mD (millidarcy), and the relative threshold is twice the average permeability of the block. The area that meets both conditions is determined as the high permeability candidate zone. Spatial connectivity analysis uses the connected component labeling algorithm to merge spatially adjacent and similar permeability pixels into the same connected component, and the connectivity criterion uses the 8-neighborhood criterion. The geometric characteristics of high permeability zones include length, width, thickness, and inclination. The length is calculated by the skeleton extraction algorithm to obtain the principal axis length, the width is the average thickness perpendicular to the principal axis direction, and the inclination is determined by the principal component analysis to obtain the principal direction vector. The strength classification of high permeability zones is based on the average permeability: strong high permeability zone (>200 mD), medium high permeability zone (100-200 mD), and weak high permeability zone (50-100 mD).

[0038] In some embodiments, the phase difference analysis using the high permeability zone and the permeability variation curve generates a seepage-construction mismatch factor, including: extracting a seepage dominant frequency from the high permeability zone; extracting a construction rhythm frequency from the permeability variation curve; based on the phase deviation of the seepage dominant frequency and the construction rhythm frequency, calculating a deviation angle; and performing spatial energy mapping on the deviation angle to form a seepage-construction mismatch factor.

[0039] The dominant frequency is extracted from the identified high permeability zone. The geometric parameters of the high permeability zone are used: the length L of the zone determines the characteristic time scale of seepage, the width W of the zone affects the lateral diffusion speed, the thickness H of the zone controls the vertical seepage capacity, and the inclination α of the zone affects the three-dimensional seepage path. The seepage dominant frequency calculation formula is fp=D / (2πL ² ), where the diffusion coefficient D is determined according to the comprehensive geometric characteristics of the high permeability zone: D=(k_avg×H×cos ²where k avg is the average permeability of the zone, φ avg is the average porosity, and μ is the fluid viscosity. The frequency correction for different shapes of high permeability zones: the frequency correction factor for linear high permeability zones (W / L < 0.2) is 1.2, the frequency correction factor for sheet-like high permeability zones (H / W < 0.1) is 0.8, and the frequency correction factor for block-like high permeability zones (H / W > 0.5) is 1.0. The composite frequency of multi-layer high permeability zones is calculated by thickness-weighted average: fp_total =∑(Hi x fp_i) / ∑Hi, where Hi is the thickness of the i-th layer, and fp_i is the dominant frequency of the i-th layer. The effect of connectivity of high permeability zones on frequency is corrected by a connectivity coefficient: the frequency response of zones with high connectivity is fast, and the frequency response of zones with low connectivity is slow.

[0040] The construction rhythm frequency is identified from the permeability variation curve. The autocorrelation analysis method is used, and the autocorrelation function R(τ) =∫k(t)k(t+τ)dt, where k(t) is the permeability time series, and τ is the time delay. The peak position of the autocorrelation function corresponds to the periodic characteristics of the construction activity, and the main periods include the drilling period, the drill tool changing period, the cycle period, etc. The frequency spectrum analysis is realized by Fourier transform of the autocorrelation function, and the main peak frequency in the frequency spectrum corresponds to the main construction rhythm frequency fs. The drilling rhythm identification is realized by amplitude analysis of the permeability variation. The mechanical vibration in the drilling process will cause periodic fluctuation of the permeability measurement. The stop drilling rhythm is identified by the stable period of the permeability. During the stop drilling period, the permeability tends to be stable, forming a platform-like feature. The drill tool changing rhythm is identified by the step change of the permeability. The permeability usually changes abruptly before and after the drill tool changing. The construction efficiency evaluation is realized by stability analysis of the rhythm frequency. Large frequency fluctuation indicates irregular construction, and small frequency fluctuation indicates good construction rhythm. The multi-well coordination analysis compares the construction rhythm frequencies of different drill holes, and identifies the coordination of the overall construction plan.

[0041] The synchronization deviation angle between the seepage dominant frequency and the construction rhythm frequency is calculated. The cross-spectrum analysis method is adopted, and the cross-power spectral density function is Sxy(ω)=X(ω)Y*(ω), wherein X(ω) is the frequency spectrum of the seepage signal, Y*(ω) is the conjugate of the frequency spectrum of the construction signal, the phase difference is φ(ω)=arg[Sxy(ω)], and X(ω) and Y*(ω) are the frequency spectrum of the seepage signal and the conjugate of the frequency spectrum of the construction signal, respectively. The deviation angle θ is defined as the phase difference value at the dominant frequency, θ=φ(fp)-φ(fs), and the angle range is-π to π. The frequency matching degree is evaluated by the frequency ratio fp / fs, and the closer the ratio to 1, the better the frequency matching, and the farther the ratio deviates from 1, the poorer the frequency matching. The phase-locked analysis calculates the time derivative of the phase difference, and the closer the derivative to 0, the better the phase locking, and the greater the fluctuation of the derivative, the greater the phase drift. The synchronization window identification analyzes the local stability of the phase difference through a sliding window, and the interval in which the phase difference variation coefficient is less than 0.1 in the window is marked as a synchronization window. The phase lead-lag discrimination: a positive deviation angle indicates that the construction rhythm lags behind the seepage change, a negative deviation angle indicates that the construction rhythm leads the seepage change, and a zero deviation indicates complete synchronization. The spatial variability of the deviation angle reflects the differences in the matching degree of the seepage-construction at different positions.

[0042] For example, the spatial energy mapping of the deviation angle forms a seepage-construction mismatch factor, including: determining a mismatch intensity level based on the deviation angle; constructing a mismatch energy field according to the spatial energy mapping of the mismatch intensity level; obtaining an energy distribution feature based on the mismatch energy field, the energy distribution feature including a mismatch sensitive area and a propagation attenuation coefficient; and establishing a seepage-construction mismatch factor using the energy distribution feature.

[0043] The deviation angle is used to determine the mismatch intensity level. A five-level classification method is adopted: extremely strong mismatch (|θ|>3π / 4, intensity value 0.9-1.0), strong mismatch (π / 2<|θ|≤3π / 4, intensity value 0.7-0.9), moderate mismatch (π / 4<|θ|≤π / 2, intensity value 0.4-0.7), weak mismatch (π / 8<|θ|≤π / 4, intensity value 0.1-0.4), and extremely weak mismatch (|θ|≤π / 8, intensity value 0-0.1). The intensity value calculation formula is S=(2|θ| / π)^β, wherein β is a nonlinear adjustment parameter, and is set to 1.5 to highlight the influence of strong mismatch. The mismatch type identification is based on the positivity of the deviation angle: positive deviation indicates that the construction lags behind the seepage change, and negative deviation indicates that the construction leads the seepage change. The intensity spatial distribution is extended to the entire reservoir space through the Kriging interpolation method. The time evolution analysis of the mismatch intensity tracks the dynamic changes of the intensity level, and identifies the time periods of mismatch aggravation and mismatch relief. The multi-parameter coupling analysis considers the correlation between the deviation angle and other parameters, such as the influence of permeability difference and construction parameter change on the intensity level.

[0044] The mismatch energy field is constructed by spatial energy mapping according to the mismatch intensity level. The potential field theory is adopted to regard the mismatch point as an energy source to generate potential energy distribution around it. The energy potential function of single-point mismatch is V(r)=S×exp(-r / λ), where S is the mismatch intensity, r is the distance from the mismatch point, and λ is the energy attenuation length constant. The total energy field of multi-point mismatch is obtained by linear superposition of the potential energy of each point: V_total=∑V i (r). The energy gradient calculation determines the energy flow direction, and the gradient vector points to the direction of the fastest energy increase. The energy streamline tracking solves the differential equation using the Euler method to obtain the trajectory line of energy propagation. The boundary condition processing considers the reflection and absorption effects of the reservoir boundary on energy propagation. The mirror method is used to process the rigid boundary, and the decay boundary condition is used for the absorbing boundary. The field strength grading divides the energy field strength into high intensity area (V>0.7V_max), medium intensity area (0.3V_max<V≤0.7V_max), and low intensity area (V≤0.3V_max).

[0045] The energy distribution characteristics are obtained based on the mismatch energy field, including the mismatch sensitive area and the propagation attenuation coefficient. The mismatch sensitive area is determined by energy gradient analysis, and the area with gradient greater than the threshold is marked as the sensitive area, which is most sensitive to mismatch changes. The geometric characteristics of the sensitive area include area, shape factor, and directionality. The sensitive area A is obtained by pixel counting, the shape factor is represented by the ratio of perimeter square to area P ² / A, and the directionality is determined by the principal component analysis to obtain the principal axis direction. The propagation attenuation coefficient λ is obtained by fitting the energy decay curve with distance, using the nonlinear least squares method to fit the exponential decay model E(r)=E0×exp(-r / λ), where E(r) is the energy value at distance r, E0 is the initial energy value, and λ is the propagation attenuation coefficient. The attenuation coefficients in different directions may be different, reflecting the anisotropic characteristics of the reservoir, and the horizontal attenuation coefficient λ h and the vertical attenuation coefficient λ v are calculated respectively. The energy aggregation degree G is obtained by calculating the concentration degree of energy distribution, and the energy dispersion degree is quantified by the ratio of standard deviation to mean value, and the aggregation degree G is inversely proportional to the dispersion degree. The spatial correlation analysis calculates the correlation coefficient of energy values at different positions, and identifies the spatial correlation length L of energy distribution through the attenuation characteristics of the correlation function, which reflects the spatial continuity of mismatch influence.

[0046] The extracted energy distribution characteristics are used to establish a seepage-construction mismatch factor. The seepage-construction mismatch factor quantitatively represents the degree of incoordination between the reservoir seepage state and the construction activity, and is constructed by fusing multi-dimensional characteristic parameters of the mismatch energy field. A multi-dimensional feature fusion method is used to comprehensively consider parameters such as mismatch sensitive area A, propagation attenuation coefficient λ, energy aggregation degree G, and spatial correlation length L. The mismatch sensitive area A reflects the spatial range of the mismatch influence, and the larger the area, the more extensive the mismatch influence. The propagation attenuation coefficient λ reflects the propagation ability of the mismatch energy, and the smaller λ, the farther the mismatch influence propagates. The energy aggregation degree G reflects the concentration degree of the mismatch, and the larger G, the more concentrated the mismatch in the local area. The spatial correlation length L reflects the spatial continuity of the mismatch, and the larger L, the stronger the spatial correlation of the mismatch. The mismatch factor calculation formula is F = w1 × (A / A_ref) + w2 × (λ_ref / λ) + w3 × G + w4 × (L_ref / L), where w1, w2, w3, and w4 are weight coefficients, and A_ref, λ_ref, and L_ref are reference values. The numerical range of the mismatch factor is normalized to the interval of 0-1, 0 represents complete matching, and 1 represents complete mismatch. The mismatch type subdivision is based on the combination mode of the energy distribution characteristics: diffusion type mismatch (λ is large and G is small), aggregation type mismatch (λ is small and G is large), and anisotropic mismatch (λ h ≠λ v ). The mismatch severity evaluation is determined by threshold: slight mismatch (F < 0.2), general mismatch (0.2 ≤ F < 0.5), severe mismatch (0.5 ≤ F < 0.8), and extremely severe mismatch (F ≥ 0.8). Time series analysis tracks the evolution trend of the mismatch factor, and identifies the time nodes of mismatch deterioration and improvement. Multi-well correlation analysis compares the mismatch factors of different boreholes, and identifies regional mismatch patterns.

[0047] The construction optimization window is automatically generated according to the seepage-construction mismatch factor. The seepage-construction mismatch factor reflects the degree of mismatch between the natural seepage law of the reservoir and the construction rhythm. When the two are seriously out of sync, engineering risks such as water inrush and collapse are easily caused. The generation process of the construction optimization window first analyzes the spatiotemporal distribution characteristics of the mismatch factor and identifies time periods and spatial regions with low mismatch degree. In the time dimension, by tracking the dynamic change curve of the mismatch factor, time periods with low values and stable trends are found, which indicate that the seepage state and the construction rhythm are relatively coordinated, and are suitable for key construction operations. In the spatial dimension, by analyzing the distribution of mismatch sensitive regions in reverse, high sensitive zones that respond dramatically to construction disturbance are avoided, and stable regions with relatively slow mismatch response are selected as the preferred construction positions. The determination of the optimization window also needs to consider various constraints, including construction efficiency requirements, safety risk control and economic cost limitations. By weighing these factors comprehensively, the best construction time and space range is determined under the guidance of the mismatch factor, forming a construction optimization window with clear time start and end points and spatial boundaries. The generation strategy of the optimization window is classified according to the severity of the mismatch: in the slight mismatch region (F<0.2), normal construction is required, only routine monitoring is needed; in the general mismatch region (0.2≤F<0.5), monitoring frequency needs to be strengthened, and construction parameters need to be adjusted appropriately; in the serious mismatch region (0.5≤F<0.8), high-risk periods should be avoided, and the time window with the lowest mismatch factor should be selected for construction; in the extremely serious mismatch region (F≥0.8), the construction plan needs to be reevaluated, and if necessary, the drilling trajectory needs to be adjusted or special construction technology needs to be used.

[0048] In step S140, a geological potential energy gradient field is constructed using the construction optimization window, and seepage diffusion analysis is performed through the geological potential energy gradient field to determine the overflow influence range. The overflow influence range, the construction risk region and the seepage-construction mismatch factor are coupled to form a coupling tensor, and a three-dimensional construction decision matrix is generated based on the coupling tensor.

[0049] Specifically, a geological potential energy gradient field is constructed using the construction optimization window. Based on the geometric parameters and energy characteristics of the construction optimization window: the center position of the window is the lowest point of potential energy, the boundary of the window is the demarcation line of the change of potential energy gradient, and the size of the window determines the spatial range of the gradient field. The potential energy function adopts a multipole expansion form, and the potential energy distribution of a single optimization window is V(r)=-A / r+B / r ² +C / r ³where A, B, C are multipole coefficients, r is the distance to the center of the window. The multipole coefficients are determined according to the characteristic parameters of the optimization window: the A coefficient is proportional to the mismatch improvement potential of the window, the B coefficient reflects the shape factor of the window, and the C coefficient considers the directional feature of the window. The multi-window coupling effect is processed by potential superposition, and the total potential V_total = ∑V_i(r), where V_i is the contribution of the i-th window. The potential gradient calculation adopts the finite difference method, and the gradient vector points to the direction of the fastest potential increase. The gradient field intensity classification is divided according to the gradient module length: strong gradient area , medium gradient area , weak gradient area . The boundary condition processing considers the constraint of the gradient field by the geological structure, and the fault boundary adopts the discontinuous boundary condition and the rock layer interface adopts the continuous boundary condition.

[0050] The seepage diffusion analysis is performed by the geopotential gradient field to determine the overflow influence range. The gradient field is used to drive the seepage to diffuse preferentially along the direction of the potential gradient, and the gradient size controls the diffusion speed and the gradient direction determines the diffusion path. The diffusion equation is established as where C is the concentration field, D is the diffusion tensor, and S is the source term. The diffusion tensor D is related to the potential gradient field: where D0 is the baseline diffusion coefficient, and a is the gradient enhancement coefficient. The source term S is determined according to the position and intensity of the construction optimization window, and the center of the window is the main source point, and the source intensity is proportional to the improvement potential of the window. The diffusion boundary condition is set according to the boundary characteristics of the gradient field: the strong gradient boundary adopts the reflection boundary, and the weak gradient boundary adopts the absorption boundary. The diffusion front tracking is realized by numerical method, and the main path and arrival time of the diffusion are identified. The overflow influence range is determined by using the concentration threshold criterion: when the concentration C decreases to 10% of the initial concentration, the position is defined as the influence boundary. The range geometric characteristics include the maximum diffusion distance, the main diffusion direction and the influence area, the maximum distance is obtained by tracking the longest diffusion path, and the main direction corresponds to the main characteristic vector of the gradient field.

[0051] In some embodiments, the tensor coupling of the overflow influence range, the construction risk area and the seepage-construction mismatch factor to form a coupled tensor includes: converting the overflow influence range into a spatial dimension tensor; mapping the construction risk area into a risk dimension tensor; reconstructing the seepage-construction mismatch factor into a time dimension tensor; and cross-fusing the spatial dimension tensor, the risk dimension tensor and the time dimension tensor to form a coupled tensor.

[0052] Convert the spill impact range into a spatial dimension tensor. Based on the geometric characteristics and intensity distribution of the spill impact range: the boundary coordinates of the impact range are converted into the spatial indices of the tensor, and the spill intensity values are taken as the numerical elements of the tensor. Spatial discretization uses a regular grid, and the grid size is determined according to the minimum characteristic scale of the impact range, ensuring that the spatial resolution is fine enough. The tensor dimension is designed as a three-dimensional tensor of n x m x p, where n, m, and p correspond to the number of grids in the x, y, and z directions respectively. The spatial distribution of the spill intensity is extended to the entire tensor space through an interpolation algorithm, and the interpolation method uses cubic spline interpolation to maintain the smoothness of the data. Boundary effects are handled by boundary conditions, and tensor elements outside the impact range are set to zero or a minimum value. Intensity normalization maps the spill intensity to the 0-1 interval, facilitating coupled operations with other dimensional tensors. Spatial gradient information is extracted through the difference operator of the tensor, preserving the spatial variation characteristics. Symmetry testing ensures that the spatial tensor meets physical constraints such as mass and momentum conservation.

[0053] Map the construction risk area to a risk dimension tensor. According to the type classification and intensity level of the risk area: different risk types (collapse, water inrush, gas outburst, and strata fracture) correspond to different components of the tensor, and the risk intensity level is taken as the numerical value of the tensor element. The tensor structure is designed as a multi-component form R^{k x l}, where k is the number of risk types and l is the number of risk levels. Risk intensity mapping uses a nonlinear transformation, which transforms the original risk value R through the function f(R)=R^α, where α is a nonlinear parameter that highlights the importance of high-risk areas. Spatial positioning information is maintained through the spatial index of the risk tensor, with each tensor element corresponding to a specific spatial location and risk type combination. Risk superposition processing considers the coexistence of multiple risk types at the same location, and uses the joint distribution method of probability theory to calculate the total risk. Time evolution information is represented by a sequence of dynamic risk tensors, recording the trend of risk over time. Weight allocation is determined according to the severity and probability of occurrence of different risk types, with the highest weight for life safety-related risks.

[0054] Reconstruct the seepage-construction mismatch factor into a time dimension tensor. The time dimension tensor is designed as T^{s x t}, where s is the number of spatial positions and t is the number of time nodes. The time resolution is determined according to the frequency of change of the mismatch factor, ensuring that the main time dynamic characteristics can be captured. Periodic analysis extracts the main periodic components of the mismatch factor through Fourier transform, and the main period corresponds to the main time mode of the tensor. Trend analysis extracts the change trend of the mismatch factor through numerical difference method, identifying growth, decay and oscillation patterns. Time correlation analysis calculates the correlation degree of mismatch states at different times, reflecting the time dependence characteristics. Multi-time scale representation is achieved through wavelet transform, preserving both short-term fluctuations and long-term trend information. Time window division divides the continuous time axis into decision-making time windows, and the mismatch characteristics within each window are relatively stable.

[0055] The cross-fusion process employs tensor outer product operation to combine three tensors of different dimensions into a high-order coupled tensor: wherein denotes tensor outer product. The fusion weight assignment considers the relative importance and contribution of the three dimensions: spatial dimension weight ws reflects the importance of spatial distribution, risk dimension weight wr reflects the priority of safety factors, and time dimension weight wt reflects the urgency of time constraints. Weighted fusion is achieved by adjusting the weight coefficients of the dimension tensors, and the spatial weight ws, risk weight wr, and time weight wt act on the corresponding tensor components. Tensor contraction reduces the order of the tensor by summing or averaging along a specific dimension: partial contraction preserves the main dimension information, and full contraction obtains a scalar indicator. Symmetrization ensures that the coupled tensor meets the physical symmetry requirement, and is achieved by standard symmetrization methods. Singular value decomposition extracts the main modes of the coupled tensor, which reflect the dominant characteristics of the system, and the secondary modes reflect the detailed characteristics.

[0056] A three-dimensional construction decision matrix is generated based on the coupled tensor. First, the main coupled modes are extracted by modal decomposition of the coupled tensor; then, the three-dimensional tensor information is projected into a two-dimensional matrix space without losing key information. The matrix element M_ij is defined as the comprehensive decision indicator of position i at time j, and the calculation formula is M_ij = ∑_kw_kT_ijk, where w_k is the weight of the kth risk level, and T_ijk is the coupled tensor element. The matrix structure design adopts a block form: the spatial dimension corresponds to the row of the matrix, the time dimension corresponds to the column of the matrix, and the risk dimension is reflected by color coding or numerical size. Decision indicator standardization processing unifies indicators of different dimensions to the 0-1 interval, facilitating comparison and decision-making. Matrix eigenvalue analysis identifies the main decision-making mode of the system, and the largest eigenvalue corresponds to the most important decision direction. Decision priority ranking is determined by the order of matrix element values, and regions with larger element values are given priority in construction. Spatial correlation analysis calculates the correlation coefficient between matrix rows to identify spatially correlated decision areas. Time dependency analysis calculates the correlation between matrix columns to identify temporal decision constraint relationships. Decision constraint conditions are represented by zero elements of the matrix, indicating that the space-time combination is infeasible or not recommended.

[0057] In step S150, energy minimization processing is performed on the three-dimensional construction decision matrix to determine the energy optimal point, a first construction trajectory set is generated according to the energy optimal point, risk-efficiency dynamic balance processing is performed on the first construction trajectory set to generate a second construction trajectory set, and the optimal construction scheme is determined through the second construction trajectory set.

[0058] Specifically, the energy minimization processing is performed on the three-dimensional construction decision matrix to determine the energy optimal point. According to the numerical distribution of the matrix elements M_ij, the matrix elements are regarded as discrete sampling points of the energy potential function, and the extremum point of the energy function is found by a numerical optimization method. The energy function is constructed by using the negative value of the matrix element as the potential energy, E(i,j)=-M_ij, so that the area with high decision index corresponds to the energy low, which conforms to the energy minimization principle. The gradient calculation adopts the finite difference method, and the energy gradient points to the direction of the fastest energy increase, and the point with zero gradient is the energy extremum point. The global optimal point search adopts the simulated annealing algorithm to avoid falling into local optimal solution, and the algorithm parameters include the initial temperature T0, the cooling rate a and the termination temperature T_min. The local optimal point identification is performed by calculating the Hessian matrix of the energy function, and the positive definite Hessian matrix corresponds to the local minimum point and the negative definite matrix corresponds to the local maximum point.

[0059] The first construction trajectory set is generated according to the energy optimal point. The energy optimal point is taken as the key node of the trajectory, and the adjacent nodes are connected to form the trajectory by a path search algorithm. The trajectory planning adopts the A* search algorithm to ensure the optimality and operability of the path. The trajectory risk value evaluation is the core performance index, which is obtained by integrating the construction risk area intensity along the trajectory, and the calculation formula is R=(∫ρ(s)ds) / L, the integration is performed along the trajectory L, wherein R is the trajectory risk value, ρ(s) is the risk density of the position s on the trajectory, and L is the total length of the trajectory. The risk density ρ(s) considers three dimensions of geological risk, construction risk and environmental risk, and the comprehensive risk level of each position is determined by weighted superposition. The trajectory efficiency value evaluation reflects the economy and time of construction, which is calculated based on the path length and geometric complexity, and the formula is E=L0 / (L×(1+κ_avg)), wherein E is the efficiency value, L0 is the straight line distance, L is the actual trajectory length, and κ_avg is the average curvature. The efficiency value considers the influence of path detour loss and construction difficulty increase, and the greater the curvature, the more frequent the steering, and the lower the construction efficiency. The multi-trajectory generation strategy considers different construction schemes and emergency plans, generates multiple candidate trajectories for each pair of start-end points, and ensures the robustness of the scheme.

[0060] In some embodiments, the risk-efficiency dynamic balance processing on the first construction trajectory set generates a second construction trajectory set, including: performing risk-efficiency potential energy gradient field construction on the first construction trajectory set to obtain potential energy distribution characteristics and gradient direction vectors; performing energy flow direction tracking based on the potential energy distribution characteristics and the gradient direction vectors to identify energy convergence points and energy divergence points; performing potential energy redistribution on the energy convergence points and the energy divergence points to establish a gradient field balance regulation strategy; and performing potential energy guidance optimization on the first construction trajectory set by using the gradient field balance regulation strategy to generate a second construction trajectory set.

[0061] Risk-efficiency potential gradient field construction based on the first construction trajectory set: convert the risk value and efficiency value of each trajectory into potential density, high-risk low-efficiency area corresponds to high potential, low-risk high-efficiency area corresponds to low potential. The potential function is designed as V(r, e) = w1 x r ² +w2 / e ² , where r is the risk value, e is the efficiency value, w1, w2 are weight coefficients. Spatial interpolation extends the potential of discrete trajectory points to continuous space, and radial basis function interpolation method is adopted to maintain the smoothness of the potential field. Potential distribution feature extraction is realized through statistical analysis, including mean, variance, skewness and kurtosis of potential, reflecting the overall distribution characteristics of potential field. Gradient calculation uses central difference format, and the gradient direction vector points to the direction of the fastest increase of potential. The spatial variability of the gradient field is analyzed by the variogram function, which identifies the correlation length and anisotropy characteristics of the gradient. Potential contour line drawing reveals the spatial distribution pattern of potential, and the dense contour line area corresponds to the area with large gradient. The boundary effect processing considers the influence of trajectory set boundary on potential field, and adopts reflective boundary condition or periodic boundary condition. Multi-scale analysis extracts potential features at different scales through wavelet transform, and large scale reflects the overall trend and small scale reflects local changes.

[0062] Based on the potential distribution characteristics and gradient direction vector, energy flow direction tracking is carried out to identify energy convergence points and energy divergence points. Streamline tracking method is adopted, and the gradient direction vector is regarded as the velocity field, and the energy flow path is tracked by integrating the streamline equation . The starting point of the streamline is determined according to the potential distribution characteristics, and is set in the boundary area and high potential area of the potential field. The streamline trajectory is solved by numerical integration method. Energy convergence point identification is realized by analyzing the convergence of streamlines, and the area where multiple streamlines converge is marked as energy convergence point. According to the potential distribution characteristics, these points correspond to the local minimum value area of potential. Energy divergence point identification is realized by analyzing the divergence of streamlines, and the area where streamlines diverge from a point to all directions is marked as divergence point. According to the potential distribution characteristics, these points correspond to the local maximum value area of potential. Convergence intensity quantification is realized by calculating the streamline density around the convergence point, and the convergence intensity level is determined by combining the potential gradient size. The divergence intensity is determined by analyzing the streamline distribution density and potential gradient characteristics around the divergence point. Streamline topology analysis constructs the connection relationship between convergence points and divergence points, and identifies the main energy transmission channel.

[0063] The identified energy sink and energy source are subjected to potential energy redistribution, and a gradient field balance control strategy is established. The specific information of the sink and source is used: according to the position coordinates (x_c, y_c) of the sink and the sink strength S_c (the sink strength is defined as the ability of the point to absorb energy), the spatial position and absorption capacity of energy absorption are determined; according to the position coordinates (x_d, y_d) of the source and the divergence strength S_d (the divergence strength is defined as the ability of the point to release energy), the spatial position and release strength of energy release are determined. The energy transfer path is planned to establish the shortest distance connection between the source and the sink, and the transmission efficiency is inversely proportional to the distance. The energy of the strong source (S_d>0.8S_max) is preferentially transferred to the strong sink (S_d>0.8×maximum divergence strength), realizing efficient energy redistribution. Energy balance constraints ensure that the output energy of each source is equal to the input energy of the corresponding sink, maintaining local energy conservation. The gradient field balance control strategy includes three levels: global control is based on the overall distribution of sinks and sources, and macroscopic balance is achieved by adjusting the strength of the main divergence source; local control is aimed at specific sink-source pairs, and local balance is achieved by accurate energy transmission; microscopic control is in the neighborhood of the sink and the source, and local imbalance is eliminated by gradient fine-tuning. The control effect is evaluated by the change of the strength of the sink and the source: the decrease of the strength of the source and the increase of the strength of the sink indicate that the control effect is good.

[0064] The gradient field balance control strategy is used to guide the optimization of the first construction trajectory set, and a second construction trajectory set with better performance is generated. Based on the distribution of the regulated potential energy field: the direction and amplitude of the trajectory adjustment are determined by the potential energy gradient, and the trajectory moves to the low potential energy area to avoid the high potential energy area. The trajectory adjustment algorithm is based on the potential energy gradient guiding principle, and the trajectory is gradually adjusted from the high-risk high-cost area to the low-risk low-cost area. The adjustment process is carried out along the direction of the fastest potential energy decrease, and the step size of each adjustment is controlled by the learning rate to ensure smooth transition of the trajectory without mutation. The core significance of this optimization is to improve the local defects of the initial trajectory. Although the first construction trajectory set connects the energy optimal point, there may be local high-risk or low-efficiency sections. Through potential energy guided optimization, the trajectory will automatically avoid high-risk areas and tend to low-risk high-efficiency areas, realizing the dynamic balance of risk and efficiency. The adaptive step size control adjusts the learning rate according to the gradient size to avoid oscillation in steep areas and stagnation in flat areas. Trajectory constraint processing ensures that the adjusted trajectory meets the geometric and physical constraints, such as minimum turning radius and maximum slope. The trajectory quality evaluation is based on comprehensive performance indicators, including risk reduction rate, efficiency improvement rate and path smoothness.

[0065] The optimal construction scheme is determined through a second set of construction trajectories. First, each optimized trajectory in the second set is transformed into an executable construction scheme, including determining the drilling path, setting construction parameters for each segment, and establishing a timeline and resource allocation. Each trajectory corresponds to a complete construction scheme, encompassing technical elements such as the trajectory's three-dimensional coordinate sequence, recommended drilling speed values ​​for each segment, mud pressure control parameters, and drill string selection recommendations. The scheme evaluation system employs the analytic hierarchy process (AHP), establishing an evaluation framework that includes four primary indicators: safety, economy, technology, and environment. Safety indicators include accident probability, personnel risk, and equipment risk, with a weight of 0.4 reflecting the safety priority principle. Economic indicators include cost control, schedule control, and resource utilization, with a weight of 0.3 reflecting economic efficiency requirements. Technical indicators cover construction difficulty, quality assurance, and technical reliability, with a weight of 0.2 reflecting technical level requirements. Environmental indicators consider environmental impact, resource consumption, and sustainability, with a weight of 0.1 reflecting environmental protection requirements. The TOPSIS method is used for scheme optimization, calculating the relative closeness of each construction scheme to the ideal scheme; the scheme with the highest closeness is determined as the optimal construction scheme.

[0066] Step S160: Execute the optimal construction plan and monitor the parameter response changes of mud pressure, drilling speed and drill bit torque. Identify abnormal parameter combination symptoms based on parameter response changes and generate graded risk warning instructions based on abnormal parameter combination symptoms.

[0067] Specifically, the optimal construction plan is implemented and the response changes of key parameters are monitored in real time. Mud pressure monitoring uses a high-precision pressure sensor with a measurement range of 0-50 MPa, an accuracy of ±0.1 MPa, and a sampling frequency of 100 Hz. Installation locations include three key points: mud pump outlet, drill pipe internal pressure, and wellhead back pressure. Drilling speed monitoring is achieved through an encoder and displacement sensor, measuring the instantaneous speed and cumulative footage of the drill bit. The speed measurement range is 0-10 m / min, with an accuracy of ±0.01 m / min and a time resolution of 1 second. Drill string torque monitoring uses a torque sensor with a measurement range of 0-5000 N·m and an accuracy of ±1 N·m. It is installed at key nodes of the drill string transmission system to record torque fluctuations in real time. The parameter data acquisition system uses a 24-bit A / D converter with multi-channel synchronous acquisition capability. Data storage uses a circular buffer to store complete data from the most recent 24 hours. Data preprocessing includes filtering and noise reduction, outlier removal, and data synchronization to ensure a unified time reference for the three types of parameters. Parameter response feature extraction is performed using sliding window analysis with a window length of 60 seconds and an overlap rate of 50%. The mean, variance, peak value, and rate of change of the parameters within each window are calculated.

[0068] In some embodiments, the identifying the abnormal parameter combination symptom based on the parameter response changes comprises: identifying an abnormal source feature and an abnormal diffusion feature from the parameter response changes; establishing an abnormal correlation node based on the abnormal source feature; establishing a conduction link between the abnormal correlation nodes using the abnormal diffusion feature; and determining an abnormal parameter combination symptom based on the conduction link.

[0069] The abnormal source feature and the abnormal diffusion feature are identified from the parameter response changes. The abnormal source feature identifies the starting time and the initial parameter of the abnormality through time series analysis: the time series of three types of parameters, mud pressure, drilling speed and drilling tool torque, are analyzed to identify the parameter type and the time point at which the abnormal deviation first occurs. The source feature parameters include the abnormal starting time t0, the first abnormal parameter type P0, the initial deviation amplitude ΔP0 and the deviation direction (positive deviation or negative deviation). The mud pressure abnormal source is usually manifested as a sudden rise or drop in pressure, and the starting feature is the pressure change rate dp / dt>threshold value; the drilling speed abnormal source is manifested as a sudden change in speed, and the feature is that the speed change rate dv / dt exceeds the normal range; the drilling tool torque abnormal source is manifested as torque fluctuation or sticking. The abnormal diffusion feature describes how the abnormality propagates from the source parameter to other parameters: the propagation time delay Δt_ij of the abnormality between different parameters is calculated, i is the source parameter, and j is the affected parameter; the attenuation or amplification characteristics of the abnormal intensity are analyzed, and the propagation effect is quantified by the propagation coefficient α_ij=ΔP_j / ΔP_i. The diffusion mode classification includes cascade diffusion (abnormality propagates to related parameters in turn), concurrent diffusion (abnormality affects multiple parameters at the same time), and oscillation diffusion (abnormality propagates back and forth between parameters).

[0070] Based on the extracted abnormal source feature information, an abnormal correlation node reflecting the origin and nature of the abnormality is established. Each abnormal source feature is mapped to a node in the network: the node identification contains attributes such as abnormal type, starting time, impact strength and spatial location. The node weight is determined according to the severity and duration of the abnormality, and the weight calculation formula is w=α×|ΔP|×√t_duration, where α is the parameter importance coefficient, |ΔP| is the deviation amplitude, and t_duration is the duration. The spatial location of the node corresponds to the three-dimensional coordinates of the drilling tool when the abnormality occurs, and the accurate position is determined through construction trajectory recording. The time attribute of the node records the starting time, peak time and end time of the abnormality, forming a life cycle description of the node. Multi-source abnormality processing considers the case of multiple abnormal sources in the same time period, and combines spatially adjacent abnormal sources into a composite node through node clustering.

[0071] The abnormal diffusion characteristic information obtained is used to build a conduction link between the established abnormal correlation nodes. According to the space-time relationship of abnormal diffusion: the causal relationship between nodes is determined according to the time delay Δt_ij of abnormal propagation, and a short delay time indicates a strong correlation, and a long delay time indicates a weak correlation. The link weight is determined according to the propagation coefficient α_ij, and the weight formula is W_ij=α_ij×p_ij, where α_ij is the propagation coefficient and p_ij is the propagation probability. The link direction reflects the direction of abnormal propagation: from the source node to the affected node, forming a directed edge structure. The conduction path type classification includes direct conduction (direct connection between nodes), indirect conduction (propagation through intermediate nodes), and bidirectional conduction (mutual influence between nodes). The link strength classification is divided according to the propagation effect: strong link (α_ij>0.8) indicates strong abnormal propagation, medium link (0.3<α_ij≤0.8) indicates medium propagation, and weak link (α_ij≤0.3) indicates weak propagation. The multi-path conduction handles the case where the anomaly propagates through multiple paths at the same time, and is modeled by parallel links and path merging.

[0072] Based on the abnormal correlation node network and the conduction link relationship constructed, the systemic abnormal parameter combination signs are determined through network analysis. The abnormal parameter combination sign determination adopts a graph theory analysis method: the network composed of abnormal correlation nodes and conduction links is analyzed as a directed weighted graph G=(V,E,W), where V is the node set, E is the edge set, and W is the weight set. The key path identification finds the main channel of abnormal propagation through the shortest path algorithm, and the path length considers the conduction delay and attenuation effect. Strongly connected component analysis identifies the tightly connected area in the network, and the nodes in the strongly connected component affect each other to form an abnormal feedback loop. The abnormal sign type is classified according to the network topology characteristics: chain sign (abnormal propagation along linear path), star sign (abnormal propagation from center node to outside), ring sign (abnormal forms closed propagation loop), and network sign (abnormal propagation in complex network). The sign intensity is quantified by normalizing the network connectivity and weight sum, and the sign intensity S=(C / C_max)×(W_sum / W_max), where C is the network connectivity, W_sum is the weight sum, C_max and W_max are the corresponding maximum values, and the numerical range of the sign intensity S is 0-1.

[0073] The abnormal parameter combination symptom-based hierarchical risk early warning instruction is generated. Based on the topological characteristics and intensity quantization results of the abnormal symptom: chain-shaped symptom corresponds to linear risk propagation, star-shaped symptom corresponds to central diffusion risk, ring-shaped symptom corresponds to cyclic amplification risk, and network-shaped symptom corresponds to complex coupling risk. The early warning level is divided into four levels: blue early warning (symptom intensity <0.3, single chain-shaped or weak star-shaped symptom), yellow early warning (symptom intensity 0.3-0.6, multiple chain-shaped or strong star-shaped symptom), orange early warning (symptom intensity 0.6-0.8, ring-shaped symptom or local network-shaped symptom), and red early warning (symptom intensity >0.8, strong ring-shaped symptom or complex network-shaped symptom). The blue early warning instruction includes strengthening monitoring and parameter fine-tuning; the yellow early warning instruction includes equipment inspection and process optimization; the orange early warning instruction includes equipment shutdown inspection and expert analysis; and the red early warning instruction includes immediate shutdown and starting of emergency response. The early warning instruction content includes early warning level, symptom type, intensity value, influence range and disposal measures. The early warning information is released through multiple channels: on-site sound and light alarm, mobile terminal push, monitoring center display and management layer notification, finally completing the drilling construction intelligent decision support and risk early warning, and realizing the intelligentization, safety and collaboration management of drilling construction.

[0074] In order to perform the drilling construction intelligent decision support and risk early warning method corresponding to the above-mentioned method embodiment, to realize the corresponding functions and technical effects. Referring to Figure 2 , Figure 2 The structure block diagram of the drilling construction intelligent decision support and risk early warning device 200 provided by the embodiment of the application is shown. For ease of illustration, only the parts related to the embodiment are shown. The drilling construction intelligent decision support and risk early warning device 200 provided by the embodiment of the application comprises: The data acquisition module 201 is configured to collect drilling acoustic wave response data by using an acoustic wave sensor array, identify reservoir resonance nodes through harmonic analysis on the drilling acoustic wave response data, obtain node coupling strength based on the reservoir resonance nodes, and construct a reservoir energy conduction network based on the node coupling strength. The risk identification module 202 is configured to track wave propagation paths based on the reservoir energy conduction network, and determine a construction risk area through the wave propagation paths. The parameter optimization module 203 is configured to extract the porosity distribution density and the permeability variation curve of each drilling block, identify a high permeability zone based on the porosity distribution density, perform phase difference analysis by using the high permeability zone and the permeability variation curve to generate a seepage-construction mismatch factor, and automatically generate a construction optimization window according to the seepage-construction mismatch factor. The decision generation module 204 is configured to construct a geological potential energy gradient field by using the construction optimization window, perform seepage diffusion analysis by using the geological potential energy gradient field to determine an overflow influence range, perform tensor coupling on the overflow influence range, the construction risk area and the seepage-construction mismatch factor to form a coupling tensor, and generate a three-dimensional construction decision matrix based on the coupling tensor. The trajectory optimization module 205 is configured to perform energy minimization processing on the three-dimensional construction decision matrix to determine an energy optimal point, generate a first construction trajectory set according to the energy optimal point, perform risk-efficiency dynamic balance processing on the first construction trajectory set to generate a second construction trajectory set, and determine an optimal construction scheme by using the second construction trajectory set. The parameter early warning module 206 is configured to execute the optimal construction scheme and monitor parameter response changes of mud pressure, drilling speed and drilling tool torque, identify abnormal parameter combination signs based on the parameter response changes, and generate a graded risk early warning instruction based on the abnormal parameter combination signs.

[0075] The drilling construction intelligent decision support and risk early warning device 200 described above can implement a drilling construction intelligent decision support and risk early warning method described in the method embodiment. The optional items in the method embodiment are also applicable to the present embodiment, and will not be described in detail herein. The remaining content of the present embodiment can be referred to the content of the method embodiment, and will not be described in detail herein.

[0076] As shown in Figure 3 The third embodiment of the present application further provides a computer device, which comprises a memory 301, a processor 302 and a computer program stored in the memory 301 and capable of running on the processor 302, and characterized in that the processor 302 implements the steps of the drilling construction intelligent decision support and risk early warning method according to the first embodiment of the present application when executing the program.

[0077] The above only describes some or preferred embodiments of the present application, and neither the text nor the drawings can limit the scope of protection of the present application. Any equivalent structural transformation based on the content of the present application and the drawings, or direct / indirect application in other related technical fields is included in the scope of protection of the present application.

Claims

1. A method for intelligent decision support and risk early warning for drilling construction, characterized in that, The method comprises the following steps: acquiring borehole acoustic response data by using an acoustic sensor array, identifying reservoir resonance nodes by harmonic analysis on the borehole acoustic response data, obtaining node coupling strength based on the reservoir resonance nodes, and constructing a reservoir energy conduction network based on the node coupling strength; tracking wave propagation paths based on the reservoir energy conduction network, and determining a construction risk area through the wave propagation paths; extracting porosity distribution density and permeability variation curves of each borehole block, identifying a high-permeability zone based on the porosity distribution density, and generating a seepage-construction mismatch factor by phase difference analysis using the high-permeability zone and the permeability variation curves, and automatically generating a construction optimization window according to the seepage-construction mismatch factor; constructing a geological potential energy gradient field by using the construction optimization window, determining an overflow influence range by seepage diffusion analysis through the geological potential energy gradient field, and forming a coupling tensor by tensor coupling of the overflow influence range, the construction risk area and the seepage-construction mismatch factor, and generating a three-dimensional construction decision matrix based on the coupling tensor; determining an energy optimal point by energy minimization processing on the three-dimensional construction decision matrix, generating a first construction trajectory set according to the energy optimal point, generating a second construction trajectory set by risk-efficiency dynamic balance processing on the first construction trajectory set, and determining an optimal construction scheme through the second construction trajectory set; executing the optimal construction scheme and monitoring parameter response changes of mud pressure, drilling speed and drilling tool torque, identifying abnormal parameter combination signs based on the parameter response changes, and generating a hierarchical risk early warning instruction based on the abnormal parameter combination signs.

2. The method of claim 1, wherein, The method comprises the following steps: establishing an energy conduction path based on the reservoir resonance nodes; monitoring acoustic impedance along the energy conduction path to obtain impedance data; backstepping energy exchange efficiency between nodes according to the impedance data; converting the energy exchange efficiency into node coupling strength.

3. The method of claim 1, wherein, The method comprises the following steps: determining a formation stress diffusion range through the wave propagation paths; extracting stress wave spectral features through the formation stress diffusion range; identifying a construction risk area according to the stress wave spectral features.

4. The method of claim 1, wherein, The method comprises the following steps: extracting a seepage dominant frequency from the high-permeability zone; extracting a construction rhythm frequency from the permeability variation curve; calculating a deviation angle based on the phase deviation between the seepage dominant frequency and the construction rhythm frequency; forming a seepage-construction mismatch factor by spatial energy mapping on the deviation angle.

5. The method of claim 4, wherein, The method comprises the following steps: determining a mismatch strength level based on the deviation angle; constructing a mismatch energy field by spatial energy mapping according to the mismatch strength level; obtaining energy distribution features based on the mismatch energy field, the energy distribution features including a mismatch sensitive area and a propagation attenuation coefficient; establishing a seepage-construction mismatch factor by using the energy distribution features.

6. The method of claim 1, wherein, The tensor coupling of the overflow influence range, the construction risk area and the seepage-construction mismatch factor includes: Converting the overflow influence range into a spatial dimension tensor; Mapping the construction risk area into a risk dimension tensor; Reconstructing the seepage-construction mismatch factor into a time dimension tensor; Cross-fusing the spatial dimension tensor, the risk dimension tensor and the time dimension tensor to form a coupling tensor.

7. The method of claim 1, wherein, The risk-efficiency dynamic balance processing of the first construction trajectory set includes: Risk-efficiency potential energy gradient field construction of the first construction trajectory set to obtain potential energy distribution characteristics and gradient direction vectors; Energy flow direction tracking based on the potential energy distribution characteristics and the gradient direction vectors to identify energy convergence points and energy divergence points; Potential energy redistribution of the energy convergence points and the energy divergence points to establish a gradient field balance regulation strategy; Potential energy orientation optimization of the first construction trajectory set by using the gradient field balance regulation strategy to generate a second construction trajectory set.

8. The method of claim 1, wherein, The abnormal parameter combination symptom based on the parameter response changes includes: Identifying abnormal source characteristics and abnormal diffusion characteristics from the parameter response changes; Establishing abnormal association nodes based on the abnormal source characteristics; Establishing conduction links between the abnormal association nodes by using the abnormal diffusion characteristics; Determining the abnormal parameter combination symptom based on the conduction links.

9. A drilling construction intelligent decision support and risk early warning device, characterized in that, It includes: A data acquisition module for collecting drilling acoustic response data using an acoustic sensor array, identifying reservoir resonance nodes by harmonic analysis of the drilling acoustic response data, obtaining node coupling strength based on the reservoir resonance nodes, and constructing a reservoir energy conduction network based on the node coupling strength; A risk identification module for tracking wave propagation paths based on the reservoir energy conduction network and determining construction risk areas through the wave propagation paths; A parameter optimization module for extracting porosity distribution density and permeability variation curves of each drilling block, identifying high permeability zones based on the porosity distribution density, generating a seepage-construction mismatch factor by phase difference analysis using the high permeability zones and the permeability variation curves, and automatically generating a construction optimization window according to the seepage-construction mismatch factor; A decision generation module for constructing a geological potential energy gradient field using the construction optimization window, determining an overflow influence range by seepage diffusion analysis through the geological potential energy gradient field, tensor coupling of the overflow influence range, the construction risk area and the seepage-construction mismatch factor to form a coupling tensor, and generating a three-dimensional construction decision matrix based on the coupling tensor; A trajectory optimization module for energy minimization processing of the three-dimensional construction decision matrix to determine an energy optimal point, generating a first construction trajectory set according to the energy optimal point, generating a second construction trajectory set by risk-efficiency dynamic balance processing of the first construction trajectory set, and determining an optimal construction scheme through the second construction trajectory set. A parameter early warning module is configured to execute the optimal construction scheme and monitor parameter response changes of mud pressure, drilling speed and drilling tool torque, identify abnormal parameter combination signs based on the parameter response changes, and generate graded risk early warning instructions based on the abnormal parameter combination signs.

10. A computer device, comprising: A computer program product comprising a memory and a processor; the memory is configured to store a computer program; the processor is configured to execute the computer program and realize the method according to any one of claims 1 to 8 when executing the computer program.

Images