Real-time intelligent analysis method and system for oil and gas well production data

Abstract

This invention provides a real-time intelligent analysis method and system for oil and gas well production data, belonging to the field of intelligent analysis technology. The method includes: Step 1, collecting and preprocessing multi-source time-series data of oil and gas wells to obtain a standard dataset; extracting from the standard dataset the real-time difference rate of change between casing pressure and oil pressure, the short-term fluctuation amplitude of instantaneous production, the centroid offset sequence of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference sequence of wellhead and bottomhole temperatures to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions; obtaining the completion structure parameters and tubing string geometric parameters of the oil and gas well to construct a three-dimensional digital model of the wellbore. This invention can realize dynamic characterization of downhole operating conditions, accurate assessment of tubing string damage, and intelligent control of production systems, ensuring the safe and stable operation of oil and gas well tubing strings and improving the level of intelligent production management and production efficiency.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of intelligent analysis, in particular to a real-time intelligent analysis method and system for oil and gas well production data. BACKGROUND

[0002] In the production process of oil and gas wells, as the development enters the middle and late stages, the stress state and flow environment of the downhole string become increasingly complex, and the traditional data analysis method has certain adaptability when dealing with the dynamic correlation of multiple parameters; the existing technology focuses more on the threshold monitoring of a single production parameter or the offline regression analysis of historical data, for example, independently tracking the casing pressure, oil pressure or liquid production, which may need to be improved in terms of fusing multi-source data and identifying the spatio-temporal evolution characteristics of downhole working conditions, and sometimes it is difficult to make early warning in a timely manner for the global risk caused by local damage accumulation.

[0003] Taking the actual operation of a certain offshore gas field as an example, the monitoring system once recorded that the difference between the casing pressure and the oil pressure showed a slow rise, while the instantaneous liquid production showed a small amplitude high frequency fluctuation; according to the conventional experience, the field technical personnel tended to believe that it was a normal production fluctuation, and no targeted measures were taken; after the fact, it was found that this phenomenon may be related to the local change of the downhole temperature and pressure gradient, which may cause the stress state of a section of the string to deviate, and be accompanied by a change in the frequency spectrum characteristics of the micro-vibration; due to the lack of effective means to analyze the spatial correlation of multiple parameters such as pressure, production, vibration and temperature at that time, the potential connection between these early signals may not have been identified in a timely manner, which to some extent affected the optimization of the control strategy; subsequently, the string of the well appeared fatigue damage at a certain connection, causing unplanned maintenance work, which had a certain impact on the production rate and operation cost. SUMMARY

[0004] The present application provides a real-time intelligent analysis method and system for oil and gas well production data, which realizes the equivalent damage quantification calculation of downhole working conditions and the self-adaptive regulation of production system.

[0005] To solve the above technical problems, the technical scheme of the present application is as follows: In a first aspect, a real-time intelligent analysis method for oil and gas well production data, the method comprising: Step 1, collecting and preprocessing multi-source time series data of oil and gas wells to obtain a standard data set, extracting a real-time difference rate sequence of casing pressure and oil pressure, a short-time fluctuation amplitude sequence of instantaneous liquid production, a frequency spectrum barycenter offset sequence of string axial vibration signals and a real-time gradient difference sequence of wellhead and bottomhole temperature from the standard data set to construct a multi-dimensional feature matrix representing the dynamic change of downhole working conditions; obtaining the completion structure parameters and string geometric parameters of the oil and gas well to construct a three-dimensional digital model of the wellbore; Step 2: Perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; using the first, second, third, and fourth nucleation base points as vertices, construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore; Step 3: Discretize the spatial tetrahedral domain into several grid cells, and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. Step 4: Calculate the local coupling damage index of each mesh element, and perform a volume-weighted integral based on the local coupling damage indices of all mesh elements to obtain the global equivalent damage factor. Step 5: Input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision-making information containing control instructions.

[0006] Secondly, the real-time intelligent analysis system for oil and gas well production data includes: The module is used to collect and preprocess multi-source time-series data from oil and gas wells to obtain a standard dataset. From the standard dataset, the real-time difference rate of change of casing pressure and oil pressure, the short-time fluctuation amplitude of instantaneous production, the centroid offset of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference of wellhead and bottom hole temperature are extracted to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions. The completion structure parameters and tubing string geometric parameters of the oil and gas well are obtained to construct a three-dimensional digital model of the wellbore. The analysis module is used to perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; and to construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore using the first, second, third, and fourth nucleation base points as vertices. The partitioning module is used to discretize the spatial tetrahedral domain into several grid cells and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. The calculation module is used to calculate the local coupling damage index of each grid cell, and to perform a volume-weighted integral based on the local coupling damage indices of all grid cells to obtain the global equivalent damage factor. The adjustment module is used to input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision information containing control instructions.

[0007] Thirdly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0008] The above-described solution of the present invention has at least the following beneficial effects: This method constructs a multi-dimensional feature matrix by integrating various time-series features such as casing pressure differential, production fluctuations, tubing vibration spectrum, and well temperature gradient. This matrix comprehensively depicts the dynamic changes in downhole operating conditions, overcoming the shortcomings of single-parameter monitoring information being one-sided and unable to reflect real operating condition changes, and improving the sensitivity of identifying downhole anomalies. Based on the three-dimensional digital model of the wellbore and four types of nucleation points, a spatial tetrahedral association domain is constructed, realizing the spatial binding of production monitoring features with the location of key wellbore structures. This transforms the originally scattered time-series data into spatially locatable and associative analysis objects, effectively revealing the spatial distribution and evolution of operating condition anomalies within the wellbore.

[0009] By constructing a physical field distribution matrix through grid discretization and physical field parameter mapping, the spatial distribution of physical quantities such as temperature, pressure, and strain inside the wellbore can be characterized. A global equivalent damage factor is obtained by combining a local coupling damage index with volume-weighted integration, which quantifies the comprehensive damage degree of the tubing string under the coupling effect of multiple physical fields. This enables a unified evaluation from local minor damage to overall equivalent damage, improving the accuracy of the assessment of tubing string structural safety risks. Based on the global equivalent damage factor, the tubing string wall thickness loss rate and residual strength coefficient are solved, and the production system is dynamically adjusted accordingly. This allows for targeted control instructions based on real-time downhole operating conditions and structural status, achieving proactive early warning and prevention of tubing string damage risks and reducing unplanned operations. Attached Figure Description

[0010] Figure 1 This is a flowchart illustrating the real-time intelligent analysis method for oil and gas well production data provided in an embodiment of the present invention.

[0011] Figure 2 This is a schematic diagram of a real-time intelligent analysis system for oil and gas well production data provided in an embodiment of the present invention. Detailed Implementation

[0012] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0013] like Figure 1 As shown, embodiments of the present invention propose a real-time intelligent analysis method for oil and gas well production data, the method comprising the following steps: Step 1: Collect and preprocess multi-source time-series data of oil and gas wells to obtain a standard dataset. Extract the real-time difference rate of change of casing pressure and oil pressure, the short-term fluctuation amplitude of instantaneous production, the centroid offset of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference of wellhead and bottom hole temperature from the standard dataset to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions. Obtain the completion structure parameters and tubing string geometric parameters of the oil and gas well to construct a three-dimensional digital model of the wellbore. Step 2: Perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; using the first, second, third, and fourth nucleation base points as vertices, construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore; Step 3: Discretize the spatial tetrahedral domain into several grid cells, and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. Step 4: Calculate the local coupling damage index of each mesh element, and perform a volume-weighted integral based on the local coupling damage indices of all mesh elements to obtain the global equivalent damage factor. Step 5: Input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision-making information containing control instructions.

[0014] In this embodiment of the invention, the method constructs a multi-dimensional feature matrix by integrating multiple time-series features such as casing pressure differential, production fluctuation, tubing vibration spectrum, and well temperature gradient. This matrix comprehensively depicts the dynamic changes in downhole operating conditions, overcoming the shortcomings of single-parameter monitoring information being one-sided and unable to reflect real operating condition changes, and improving the sensitivity of identifying downhole anomalies. Based on the three-dimensional digital model of the wellbore and four types of nucleation points, a spatial tetrahedral association domain is constructed, realizing the spatial binding of production monitoring features with the key structural locations of the wellbore. This transforms the originally scattered time-series data into spatially locatable and associative analysis objects, effectively revealing the spatial distribution and evolution of operating condition anomalies within the wellbore.

[0015] By constructing a physical field distribution matrix through grid discretization and physical field parameter mapping, the spatial distribution of physical quantities such as temperature, pressure, and strain inside the wellbore can be characterized. A global equivalent damage factor is obtained by combining a local coupling damage index with volume-weighted integration, which quantifies the comprehensive damage degree of the tubing string under the coupling effect of multiple physical fields. This enables a unified evaluation from local minor damage to overall equivalent damage, improving the accuracy of the assessment of tubing string structural safety risks. Based on the global equivalent damage factor, the tubing string wall thickness loss rate and residual strength coefficient are solved, and the production system is dynamically adjusted accordingly. This allows for targeted control instructions based on real-time downhole operating conditions and structural status, achieving proactive early warning and prevention of tubing string damage risks and reducing unplanned operations.

[0016] In a preferred embodiment of the present invention, step 1 involves collecting and preprocessing multi-source time-series data of oil and gas wells to obtain a standard dataset. The multi-source time-series data includes: casing pressure time-series values, oil pressure time-series values, instantaneous production time-series values, wellhead temperature time-series values, bottom hole temperature time-series values, and tubing string axial vibration time-series signals. From the standard dataset, the sequence of real-time difference in casing pressure and oil pressure, the sequence of short-term fluctuations in instantaneous production, the sequence of centroid offset of the tubing string axial vibration signal spectrum, and the sequence of real-time gradient difference between wellhead and bottom hole temperatures are extracted to construct a multi-dimensional feature matrix characterizing the dynamic changes in downhole operating conditions. The completion structure parameters and tubing string geometric parameters of the oil and gas well are obtained to construct a three-dimensional digital model of the wellbore. These parameters include: wellhead flange coordinates, casing coupling position coordinates, perforation section top depth coordinates, packer installation depth coordinates, casing outer diameter, casing inner diameter, tubing outer diameter, tubing inner diameter, and the elastic modulus and Poisson's ratio of the tubing string material, specifically including: At a preset acquisition frequency of once per minute, multi-source time-series data acquisition operations are carried out simultaneously at the oil and gas well production site. Through various monitoring sensors deployed at the wellhead and downhole, time-series values ​​of casing pressure, oil pressure, instantaneous fluid production, wellhead temperature, bottom hole temperature, and tubing axial vibration are acquired. During the acquisition process, the timestamps of all parameters are kept consistent, so that various monitoring data remain synchronous in the time dimension. The acquired raw multi-source time-series data were preprocessed sequentially. First, an outlier and missing data segments were identified and removed by using a threshold of 30% exceeding the normal production range, combined with time-series continuity checks, to ensure the validity of the data. Next, filtering algorithms were used to smooth signals such as pressure, temperature, and vibration to reduce noise interference and address environmental noise and electromagnetic interference introduced during sensor acquisition. Subsequently, time-series alignment and normalization were performed on data from different acquisition nodes to eliminate time misalignment caused by sampling delays from different devices. Finally, all data were normalized to unify data with different dimensions and amplitude ranges into a standardized numerical range, ultimately forming a standard dataset that is time-continuous, numerically standardized, and free from significant interference.

[0017] Based on the preprocessed standard dataset, four types of key feature sequences were extracted and calculated. For the time series data of casing pressure and oil pressure, the real-time difference between casing pressure and oil pressure was calculated first, and then the change rate was calculated according to the formula: Change rate of real-time difference between casing pressure and oil pressure = (Difference at the current time - Difference at the previous time) ÷ Time interval. This yielded the sequence of change rate of real-time difference between casing pressure and oil pressure. For the time series value of instantaneous fluid production, a fixed-length sliding window was used to traverse the entire time series data. Within each short window, the maximum and minimum values ​​of fluid production were first found. Then, the fluctuation amplitude of fluid production was characterized by the formula: Fluid production fluctuation amplitude = (Maximum value of fluid production within the window - Minimum value of fluid production within the window) ÷ Average value of fluid production within the window. The larger the value, the more violent the fluctuation of fluid production within this period; the smaller the value, the more stable the fluid production. After calculating window by window, a sequence of short-term fluctuation amplitude of instantaneous fluid production was formed. For the time-series signal of axial vibration of the tubing string, the time-domain vibration signal is first converted into a frequency-domain distribution by fast Fourier transform. Then, the centroid of the spectrum is calculated moment by moment according to the formula: vibration spectrum centroid = (sum of the product of the amplitude of each frequency component and the corresponding frequency) ÷ sum of the amplitudes of each frequency component. The difference between the centroids of the spectrum at adjacent moments is used to obtain the sequence of centroid offset of the axial vibration signal spectrum. For the time-series data of wellhead and bottom hole temperature, the real-time difference between the bottom hole and wellhead temperature is first calculated moment by moment. Then, the temperature gradient with depth is calculated according to the formula: wellbore temperature gradient = (bottom hole temperature - wellhead temperature) ÷ total wellbore depth. The real-time gradient difference sequence between the wellhead and bottom hole temperature is formed by the continuous change of this gradient over time.

[0018] The obtained sequences of real-time difference rate of casing pressure and oil pressure, short-term fluctuation amplitude of instantaneous production, centroid offset of tubing axial vibration signal spectrum, and real-time gradient difference of wellhead and bottom hole temperature are time-aligned using a unified timestamp to ensure that each moment corresponds to the synchronous values ​​of the four types of features. Then, they are arranged in an orderly manner with time as the row dimension and the four types of features as the column dimension, forming a well-organized, data-corresponding multidimensional feature matrix. In this multidimensional feature matrix, each row corresponds to a specific moment, and the four types of feature values ​​in that row correspond to the four core operating condition indicators at that moment: dynamic change of downhole pressure difference, stability of production, tubing vibration state, and wellbore temperature gradient distribution. The continuous change trend of the feature values ​​in each column of the matrix over time (each row) reflects the dynamic evolution of various operating conditions. The fluctuation of the real-time difference rate of casing pressure and oil pressure reflects the evolution of the downhole pressure difference from stable equilibrium to gradual imbalance; a continuous increase or decrease in the value corresponds to the change in the pressure state. The following parameters are considered: Continuous abnormal changes; The numerical changes in the short-term fluctuation amplitude sequence of instantaneous production volume can reflect the evolution of production volume from stable output to increased or decreased fluctuation. A larger value indicates poorer production stability, and a continuous increase in the value corresponds to a gradual intensification of abnormal production fluctuations. The numerical changes in the centroid offset sequence of the tubing string axial vibration signal spectrum can reflect the evolution of tubing string vibration state from normal stability to abnormal offset. The baseline range for the centroid offset of the tubing string vibration spectrum is set to 0 to 0.5 Hz. A value deviating from this range and continuously increasing corresponds to the gradual development of abnormal tubing string vibration. The numerical changes in the real-time temperature gradient difference sequence between the wellhead and bottom hole can reflect the evolution of the wellbore temperature gradient from uniform distribution to local distortion. The normal range for the wellbore temperature gradient is set to 0.02 to 0.05℃ / m. Sudden changes or continuous deviations from this normal range correspond to the formation and development of abnormal temperature gradient distribution. Through the synchronous changes and interrelation of these four types of characteristic values, the overall dynamic change law of downhole operating conditions from normal stability to abnormal fluctuations can be captured.

[0019] While constructing data features, the well completion construction data, tubing factory inspection reports, and downhole operation as-built drawings of the oil and gas well were retrieved. The completion structural parameters and tubing geometric parameters of the oil and gas well were collected and organized simultaneously. The key dimensions, location coordinates, and material performance indicators in the drawings and reports were checked and standardized item by item, and errors and duplicate data were eliminated. Specifically, this included the coordinates of the wellhead flange, the location coordinates of the casing coupling, the depth coordinates of the top of the perforated section, the coordinates of the packer installation depth, the outer diameter of the casing, the inner diameter of the casing, the outer diameter of the tubing, the inner diameter of the tubing, and the elastic modulus and Poisson's ratio of the materials used in the tubing. Based on the actual wellbore space dimensions and structural layout of the oil and gas well, and combined with the above complete structural and geometric parameters, a three-dimensional digital model of the wellbore that perfectly matches the actual wellbore was constructed using a three-dimensional modeling method. The model completely reproduces the accurate spatial position, geometric dimensions, and relative assembly relationship of key downhole components such as casing, tubing, casing coupling, perforated section, and packer.

[0020] This embodiment, by simultaneously collecting multiple types of production monitoring time-series data and completing standardized preprocessing, can ensure data quality and time-series consistency, avoiding interference from original data errors on analysis results. It extracts four types of feature sequences that directly reflect changes in downhole operating conditions and constructs a multi-dimensional feature matrix, which can comprehensively characterize the linkage changes of parameters such as pressure, production, vibration, and temperature. This overcomes the limitation that a single parameter cannot fully represent complex downhole operating conditions and improves the ability to capture abnormal changes in operating conditions. Furthermore, by combining complete well completion and tubing string parameters to construct a three-dimensional digital model of the wellbore, it realizes the correspondence between production data characteristics and the actual spatial structure of the wellbore, effectively solving the problem of data and spatial structure disconnect in analysis.

[0021] In a preferred embodiment of the present invention, step 2 includes: Step 200: Local extreme point detection is performed on the real-time difference rate of change sequence of casing pressure and oil pressure in the multidimensional feature matrix. Extreme point features corresponding to the zero-crossing points of the first derivative are extracted. The first nucleation base point is determined based on the spatial location of the extreme point features in the three-dimensional digital model of the wellbore. The first nucleation base point is located at the center point of the wellhead flange. Specifically, this includes: performing local extreme point detection on the real-time difference rate of change sequence of casing pressure and oil pressure in the multidimensional feature matrix; firstly, performing point-by-point differentiation on the sequence with time as the independent variable to obtain the corresponding first derivative sequence; then, traversing the first derivative sequence to find the positions where the value changes from positive to negative or from negative to positive and crosses zero, thereby determining the time of occurrence of the extreme point. During the determination, the difference between the real-time differential rates of casing pressure and oil pressure at two adjacent moments is calculated to determine the trend of its sign change. When the difference changes from positive to negative or from negative to positive, it can be determined that the moment is the zero point of the first derivative, and the corresponding position is the local extreme point. After extracting the feature information corresponding to the extreme point, it is matched with the spatial position of the wellbore three-dimensional digital model. According to the preset correspondence rules between production monitoring features and wellbore structure, the first nucleation base point is determined. The preset correspondence rule is that the extreme point feature of the real-time differential rate of casing pressure and oil pressure directly corresponds to the wellhead position. Therefore, the spatial coordinates of this base point are fixed as the coordinates of the center point of the wellhead flange in the wellbore three-dimensional digital model.

[0022] Step 201 involves performing joint mutation point analysis on the multidimensional feature matrix to extract synchronous anomaly segments from the real-time difference rate of change sequence between casing pressure and oil pressure and the short-term fluctuation amplitude sequence of instantaneous production volume. Based on the spatial location of the starting time of these synchronous anomaly segments in the wellbore 3D digital model, the second nucleation base point is determined. This second nucleation base point is located at the center point of the top of the downhole perforated section. Specifically, this includes: conducting joint mutation point analysis on the multidimensional feature matrix, first aligning the real-time difference rate of change sequence between casing pressure and oil pressure and the short-term fluctuation amplitude sequence of instantaneous production volume along a unified time axis. The numerical changes of the two sets of sequences are then compared hourly to identify segments where both deviate from the normal fluctuation range and exhibit synchronous anomalies within the same time period. The normal fluctuation range for the real-time difference rate of change between casing pressure and oil pressure is set to -0.02 to 0.02 MPa per minute, and the normal fluctuation range for the short-term fluctuation amplitude of instantaneous production is set to 0 to 0.15. When determining the synchronous anomaly segment, the earliest moment when both sequences simultaneously exceed their respective normal fluctuation ranges is taken as the segment's starting moment. This starting moment reflects the initial node of the combined anomaly in the downhole operating conditions. This starting moment is spatially mapped and associated with the wellbore 3D digital model. Based on the preset correspondence between production parameters and downhole structure—that is, the starting moment of synchronous pressure and production anomalies corresponds to the location of the perforation section where downhole fluid enters—a second nucleation base point is determined. The spatial location of this base point is fixed at the center point of the top of the downhole perforation section in the wellbore 3D digital model.

[0023] Step 202 involves performing cumulative offset analysis on the centroid offset sequence of the axial vibration signal spectrum in the multidimensional feature matrix, extracting the cumulative offset extreme points, and determining the third nucleation base point based on the spatial position of the time corresponding to the cumulative offset extreme point in the three-dimensional digital model of the wellbore. The third nucleation base point is located at the midpoint of the casing coupling. Specifically, this includes performing cumulative offset analysis on the centroid offset sequence of the axial vibration signal spectrum in the multidimensional feature matrix, accumulating the offset time by time starting from the beginning of the sequence. The cumulative offset at the current time = the cumulative offset at the previous time + the centroid offset at the current time. The complete cumulative offset change curve is obtained by iteratively calculating the offset. Traversing the cumulative offset change curve, the point with the largest value is extracted as the cumulative offset extreme point. This extreme point can reflect the moment when the abnormal accumulation of tubing vibration is most significant. The normal fluctuation range of the centroid offset of the tubing axial vibration signal spectrum is set to 0 to 0.5 Hz. The moment corresponding to this extreme point is mapped to the wellbore three-dimensional digital model. Combining the vibration characteristics with the preset correspondence between the tubing structure and the vibration spectrum centroid offset extreme point, the casing coupling position where the tubing stiffness changes abruptly is corresponding to the casing coupling position. The third nucleation base point is determined. The spatial position of this base point is fixed as the midpoint of the casing coupling in the wellbore three-dimensional digital model.

[0024] Step 203 involves detecting inflection points in the real-time gradient difference sequence of wellhead and bottomhole temperatures in the multidimensional feature matrix, extracting gradient change inflection points, and determining the fourth nucleus base point based on the spatial location of the gradient change inflection point in the three-dimensional digital model of the wellbore. The fourth nucleus base point is located at the center point of the packer position. Specifically, this includes: detecting inflection points in the real-time gradient difference sequence of wellhead and bottomhole temperatures in the multidimensional feature matrix; firstly, calculating the second derivative of the sequence; identifying gradient change inflection points through the numerical change and sign reversal features of the second derivative; when the second derivative of the sequence changes from positive to negative or from negative to positive, and the number of inflection points increases... When the value crosses zero, the position is determined to be the inflection point of the gradient change. This inflection point corresponds to the moment when the temperature gradient of the wellbore changes significantly. The normal fluctuation range of the real-time temperature gradient difference between the wellhead and the bottom of the well is set to 0.02 to 0.05℃ / m. The time information corresponding to this inflection point is extracted and spatially matched with the three-dimensional digital model of the wellbore. According to the preset correspondence rule between the temperature gradient characteristics and the position of the downhole tool, the temperature gradient inflection point corresponds to the sealing position of the packer in the wellbore. The fourth nucleation base point is determined. The spatial position of this base point is fixed as the center point of the packer position in the three-dimensional digital model of the wellbore.

[0025] Step 204: Using the first, second, third, and fourth nucleation base points as vertices, construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore. Specifically, after determining the precise spatial coordinates of the first, second, third, and fourth nucleation base points in sequence, use these four non-coplanar spatial points as the four vertices of the spatial tetrahedron. Connect each vertex sequentially inside the three-dimensional digital model of the wellbore to construct a closed and continuous spatial surface, ultimately forming a spatial tetrahedral association domain that encloses the key areas of the wellbore. This association domain completely covers the core stress and fluid flow sections of the wellbore, such as the wellhead flange, perforation section, casing coupling, and packer, binding multidimensional temporal characteristics with the location of key wellbore structures.

[0026] This embodiment extracts key feature points reflecting downhole operating anomalies by performing extreme value detection, joint mutation analysis, cumulative offset analysis, and inflection point detection on multidimensional features, and maps them one-to-one with the locations of key components in the wellbore, thus realizing the transformation of time-series data features into spatial locations. By constructing a spatial tetrahedral association domain with four key locations as vertices, the dispersed monitoring features can be concentrated in the core risk area of ​​the wellbore, effectively solving the problem of data disconnection from spatial structure.

[0027] In a preferred embodiment of the present invention, step 3 includes: Step 300a: Extract twelve edges of the spatial tetrahedral associated domain as boundary curves. Perform equidistant sampling on each boundary curve to obtain a set of sampling points. Calculate the curvature value at each sampling point. Detect local maxima of curvature from the curvature values ​​as inflection points of the boundary curves. Specifically, this includes extracting twelve edges from the constructed spatial tetrahedral associated domain contour, namely: from the center point of the wellhead flange to the top center point of the perforated section; from the center point of the wellhead flange to the midpoint of the casing coupling; from the center point of the wellhead flange to the center point of the packer position; from the top center point of the perforated section to the midpoint of the casing coupling; from the midpoint of the casing coupling to the center point of the packer position; and each corresponding reverse connecting edge (from the top center point of the perforated section to the wellhead flange)... The system comprises twelve boundary edges, including the points from the center of the wellhead flange, the point from the midpoint of the casing coupling to the center of the wellhead flange, the point from the center of the packer position to the center of the wellhead flange, the point from the midpoint of the casing coupling to the center of the top of the perforated section, the point from the center of the packer position to the center of the top of the perforated section, and the point from the center of the packer position to the midpoint of the casing coupling. Each edge is formed by directly connecting the spatial coordinates of two adjacent nucleation base points. Each edge is treated as an independent boundary curve for refined geometric processing. At a fixed sampling interval of 0.1 meters, uniform and equidistant sampling is performed on each boundary curve from the start point to the end point. During the sampling process, the sampling points are ensured to be evenly distributed along the curve without uneven density. Finally, a set of sampling points with consistent number, continuous position, and regular distribution is formed on each boundary curve.

[0028] For each sampling point, taking the current sampling point as the center, take its immediate forward and backward sampling points. The tangent change rate = the difference in coordinates between adjacent sampling points ÷ the sampling interval, and the arc length change = the straight-line distance between the forward and backward sampling points. Then, calculate the curvature of the point according to the curvature value = 2 × tangent change rate ÷ arc length change. After calculating point by point, a complete and continuous curvature distribution sequence of the entire boundary curve is obtained. Traverse the curvature values ​​of all sampling points, and determine the position where the curvature value is greater than the curvature values ​​of both the forward and backward adjacent points as a local maximum point of curvature. Mark all local maximum points of curvature that meet this condition as the inflection points of the current boundary curve.

[0029] Step 300b: Divide each boundary curve into multiple curve segments based on the inflection points, and calculate the length of each curve segment. Based on the length of the curve segment and the spatial position of the inflection points, construct a parametric spline surface covering the four triangular faces of the spatial tetrahedral associated domain, using the boundary curve as a constraint, so that the boundary of the corresponding spline surface fits the boundary curve. Specifically, this includes: dividing the continuous and complete boundary curve into multiple independent, smooth curve segments without abrupt changes in angle, using the inflection points as rigid segmentation nodes, based on the spatial coordinates of the identified inflection points on each boundary curve; for each segmented curve segment, calculate the three-dimensional straight-line distance between each group of adjacent sampling points sequentially from the beginning to the end of the segment, and then sum the distances between all adjacent sampling points to obtain the actual geometric length of the curve segment. Using the length values ​​of each curve segment and the spatial coordinates of the inflection points as rigid constraints, and the overall boundary curve as the shape contour constraint, a cubic parametric spline fitting method is used to construct the spatial surface. The specific process is as follows: First, based on each boundary curve (the edge of the spatial tetrahedral associated domain), the boundary curve is sampled at equal intervals according to a fixed sampling interval of 0.1 meters to obtain uniformly distributed sampling points. Combined with the identified inflection points, i.e. the locations of local bends and curvature changes in the curve, as control points, these sampling points and inflection points are used together as control points for surface fitting, and the three-dimensional spatial coordinates of each control point are determined.

[0030] A cubic parametric spline fitting method is used for surface construction. First, based on the sampling points and inflection points of each boundary curve, cubic spline basis functions are constructed. The expression for the basis functions is as follows: (in , , , (where the coefficients are undetermined), The function value (i.e., the three-dimensional spatial coordinates of that point) represents any point on the boundary curve. The parameter variables representing the boundary curve (corresponding to the length direction of the curve) require the establishment of a system of linear equations based on the endpoint conditions and piecewise connection conditions of the boundary curve to solve for the undetermined coefficients. The endpoint conditions refer to the constraints satisfied by the start and end points of the boundary curve segments, specifically including the function values ​​(i.e., the spatial coordinates of the corresponding control points) and the first derivative values ​​(reflecting the trend of the curve at the endpoints). The piecewise connection conditions refer to the constraints at the inflection points of adjacent curve segments, requiring that the function values, first derivatives, and second derivatives of adjacent curve segments at the inflection points be equal to ensure the continuous and smooth segmentation of the curve. Specifically, assuming that the adjacent curve segments are... (Domain) )and (Domain) The inflection point is ,in Corresponding to the curve segmentation on the left, , , , These are the undetermined coefficients of the piecewise basis function. This is the parameter value for the starting point of this segment. These are the inflection point parameter values ​​for two adjacent curve segments; Corresponding to the right curve segmentation, , , , These are the undetermined coefficients of the piecewise basis function. Here is the endpoint parameter value for this segment; the second derivative of the curve is... (Reflecting the curvature change of the curve, i.e., the degree of bending), the set of conditional equations for piecewise connection is: ; In this system of equations, the first equation corresponds to the function values ​​of adjacent curve segments at the inflection points being equal, ensuring continuity at the inflection points. The second equation corresponds to the first derivatives of adjacent curve segments at the inflection points being equal, ensuring consistent tangent directions at the inflection points. The third equation corresponds to the second derivatives of adjacent curve segments at the inflection points being equal, ensuring a smooth curvature transition without significant abrupt changes at the inflection points. Integrating this system of piecewise connection condition equations with the linear equations established based on the boundary curve endpoint conditions (corresponding to the function values ​​and first derivative constraints at the start and end points of the curve segments) forms a complete linear equation system. By solving this complete linear equation system, the undetermined coefficients in the basis functions of each curve segment can be obtained. The specific values ​​are used to ensure that the basis functions meet the requirements of piecewise smoothness, continuity and differentiability, thereby ensuring that adjacent curve segments are continuous in position at inflection points, have consistent tangent directions, and smooth curvature transitions without obvious abrupt changes.

[0031] Using the boundary curve as the contour constraint, the four triangular surfaces of the spatial tetrahedral associated domain (which are respectively the center points of the wellhead flange) are... Center point of the top of the perforation section The triangular surface formed by the midpoint of the casing coupling and the center point of the wellhead flange Center point of the top of the perforation section The triangular surface formed by the center point of the packer location, and the center point of the wellhead flange. Midpoint of sleeve coupling The triangular surface formed by the center point of the packer position, and the center point of the top of the perforation section. Midpoint of sleeve coupling The triangular surface formed by the center point of the packer position was fitted using a bivariate cubic spline interpolation method. Specifically, the parameter domain of the triangular surface was first divided into a uniform parameter grid with a grid spacing of 0.05 meters (meaning the spatial distance between adjacent grid points was 0.05 meters). Two mutually perpendicular interpolation parameters were then defined. and ( ≥0、 ≥0 and + ≤1), using the coordinates of the control points of the three vertices of the triangle as a reference, and through the bivariate interpolation formula, i.e. ,in , , The spatial coordinates of the three vertices of the triangle. (Using the spatial coordinates of the interpolation points), the spatial coordinates of all interpolation points inside the triangular surface are calculated point by point, and a continuous and smooth surface is generated by fitting. During the fitting process, the fit between the surface boundary and the boundary curve is controlled to ensure that the boundary of each spline surface completely coincides with and closely fits the corresponding boundary curve, without misalignment, gaps, protrusions or depressions; at the same time, the undulation of the surface is adjusted by combining constraints such as the length of the curve segment and the position of the inflection point, so that the surface shape completely matches the contour of the spatial tetrahedral associated domain, and finally forms a smooth surface that covers the entire spatial tetrahedral associated domain, fits the boundary curve, and has no obvious flaws.

[0032] Step 300c involves recursively subdividing the parametric spline surface to generate a boundary surface mesh composed of multiple surface patches. The normal vector of each surface patch on the boundary surface mesh is calculated, and the surface curvature at each grid vertex is estimated. Based on the estimated surface curvature distribution characteristics, regions of drastic surface change are identified. According to the positional distribution of these drastic surface change regions on the boundary surface mesh, tetrahedral mesh elements are generated within the spatial tetrahedral associated domain, ensuring that the generated positions of the tetrahedral mesh elements correspond to the positions of the drastic surface change regions. Specifically, this includes performing multi-level recursive subdivision processing on the constructed parametric spline surface. The subdivision level is set to 3 levels: Level 1 coarse subdivision, which evenly divides the entire triangular surface into 4 basic triangular patches based on the center point; Level 2 medium subdivision, which further subdivides each of the previous level's obtained subdivisions into 4 basic triangular patches. Each basic triangular facet is further subdivided into four smaller triangular facets based on its center point. The third level of subdivision repeats this subdivision rule, performing a final fine subdivision on all subdivided facets. The subdivision ratio is 1:4, meaning each parent facet is evenly divided into four child facets. The final subdivision results in a minimum spatial size of 0.05 meters. Based on the surface's spatial undulations, bending angles, and geometric gradients, the surface is gradually refined from a coarse mesh, subdividing the smooth surface layer by layer into numerous small facets of uniform size, regular shape, and consistent topology. All facets are then sequentially assembled according to their spatial positions to form a closed, continuous boundary surface mesh. The normal vector is calculated for each small facet within the boundary surface mesh. The formula for calculating the normal vector is: , Let the two edge vectors of the surface patch be... The normalized normal vector is used to determine the spatial orientation of each surface patch. Simultaneously, the curvature of each vertex in the boundary surface mesh is estimated by combining the undulation, bending angle, and angle between adjacent surface patches. That is, vertex curvature = sum of the angles between the normal vectors of all adjacent surface patches surrounding that vertex ÷ number of adjacent surface patches. Based on the curvature distribution of all mesh vertices, regions with curvature values ​​greater than 1.5 times the average curvature of the surrounding area are considered to have significantly higher curvature values ​​than the surrounding area. These regions represent areas with drastic changes in spatial morphology, obvious bending and undulation (where the average curvature of the surrounding area is the average curvature of the five adjacent vertices surrounding that vertex).

[0033] Based on the specific location distribution of regions with drastic surface changes on the boundary surface mesh, an adaptive tetrahedral mesh generation operation is performed within the spatial tetrahedral associated domain. Specifically, the baseline mesh unit size is initially set to 0.05 meters. For regions with drastic surface changes, the mesh unit size is reduced to 0.5 times the baseline size, using a dense arrangement to generate one tetrahedral mesh unit every 0.025 meters, ensuring the capture of details related to the drastic surface changes. For regions with gentle surfaces, regular shapes, and no obvious bends, the baseline mesh unit size (0.05 meters) is maintained, and a sparse arrangement is used. The grid is generated in a way that creates a tetrahedral grid cell every 0.05 meters to reduce the number of redundant grid cells. During the generation process, the matching degree between the distribution density of the grid cells and the curvature of the surface is detected in real time. If the curvature value of a certain area drops to less than 1.5 times the average curvature of the surrounding area, the grid arrangement is automatically switched to the reference size. If the curvature value continues to be higher than 1.5 times the average curvature of the surrounding area, a dense grid arrangement is maintained to ensure that the density of the grid matches the degree of change in the surface space. This approach ensures the calculation accuracy in areas with drastic surface changes through dense grids, while reducing the overall computational load through sparse grids, thus balancing computational accuracy and efficiency.

[0034] Step 300d involves determining the spatial coordinates of the four vertices of each tetrahedral mesh element in the 3D digital model of the wellbore, using these coordinates as the node coordinates of the corresponding mesh element, and calculating the element volume of each tetrahedral mesh element based on these spatial coordinates. Specifically, this includes: reading the 3D spatial coordinates of the four vertices of each generated tetrahedral mesh element in the 3D digital model of the wellbore. Each vertex is a key location point within the tetrahedral domain, specifically any four non-coplanar points among the center point of the wellhead flange, the top center point of the perforation section, the midpoint of the casing coupling, and the center point of the packer position (the four vertices of each tetrahedral mesh element are selected from these four nucleation base points and the vertices of the surface patches generated after recursive subdivision, ensuring that the vertices fall within the contour or interior of the tetrahedral domain and correspond to the key structural positions of the wellbore). The 3D spatial coordinates of these four vertices are directly recorded, stored, and indexed as the node coordinates of the corresponding tetrahedral mesh element. Based on the spatial coordinates of the four vertices of the tetrahedral mesh element, the volume of the parallelepiped enclosed by the three sets of edge vectors formed by the four vertices is calculated. For the calculation basis, the element volume of each tetrahedral mesh element is calculated one by one, and the calculation is expressed as V= Where V is the volume of a tetrahedral mesh element. The volume of the parallelepiped formed by the three sets of edge vectors of the four vertices is calculated using the size scale of the three-dimensional digital model of the well shaft to ensure that the volume calculation result corresponds completely to the actual space size, without scaling, offset, or calculation deviation.

[0035] Step 301: Based on the standard dataset and the wellbore 3D digital model, a spatial interpolation method is used to calculate the real-time physical field parameters at the center point of each grid cell. These real-time physical field parameters include temperature, pressure, and strain. Specifically, based on the standard dataset after outlier removal, noise filtering, time-series alignment, and numerical normalization, and the wellbore 3D digital model constructed by combining well completion structure parameters, tubing geometry parameters, and material parameters, a spatial interpolation algorithm is used to perform unified physical field calculations on all tetrahedral grid cells within the spatial tetrahedral associated domain. For each tetrahedral grid cell, firstly... The geometric center point is determined by calculating the coordinates of its four vertices (center point of the wellhead flange, center point of the perforated section top, midpoint of the casing coupling, center point of the packer position, and any four non-coplanar points among the vertices of the subdivided surface patch). Spatial interpolation is then performed using real-time data from four fixed, known monitoring points around the wellbore: the wellhead flange center point, the perforated section top center point, the casing coupling midpoint, and the packer position center point. The real-time physical field data collected at each monitoring point includes real-time temperature, pressure, and strain values. The interpolation calculation uses a distance-weighted average method, i.e., P = In the formula, P is the physical field parameter at the center point (one of temperature, pressure, or strain). For the first Real-time acquisition of physical field parameters at each corresponding monitoring point From the center point of the grid cell to the 1st The three-dimensional spatial distance between monitoring points =4 (4 key wellbore monitoring points are fixed for interpolation); through this interpolation calculation, the single-point measured data of discrete monitoring points can be transformed into continuously distributed spatial physical field values, reflecting the physical field distribution of the local environment downhole. The closer to the monitoring point, the higher the weight of the measured data of that point to the center point parameter; the farther away, the lower the weight. This fully restores the continuous distribution law of temperature, pressure, and strain inside the wellbore as they change with spatial position, ensuring that the physical field parameter values ​​are strictly consistent with the actual spatial position of the wellbore and the current production and operation status, and truly reflecting the physical field distribution of the local environment downhole.

[0036] Step 302: Map the real-time physical field parameters at the center points of all grid cells to the corresponding grid cells. Based on the mapping results, construct the physical field distribution matrix within the spatial tetrahedral domain. Specifically, this includes: calculating the real-time temperature, pressure, and strain parameters at the center points of each tetrahedral grid cell, and mapping them sequentially to the corresponding tetrahedral grid cells according to their unique grid cell numbers. Each grid cell corresponds to four distinct vertices: the center point of the wellhead flange, the top center point of the perforated section, the midpoint of the casing coupling, the center point of the packer location, and the vertex of the subdivision surface. The physical field parameters are mapped one-to-one with and bound to spatial grid cells. The physical field distribution matrix is ​​generated as follows: All tetrahedral grid cells are assigned consecutive unique numbers according to their spatial location. These numbers are used as row numbers. Temperature, pressure, and strain are fixed as three columns. The physical field values ​​of the corresponding grid cells are filled in row by row, while simultaneously recording the spatial coordinates of the four vertices of each grid cell in each row. This forms a well-organized two-dimensional numerical matrix with one-to-one numerical correspondence, thus completing the construction of the physical field distribution matrix. The row number of all tetrahedral grid cells is used as the row dimension, and the values ​​of temperature, pressure, strain, and strain are used as the column dimensions. The three types of physical field parameters are transformed into column dimensions and arranged in a fixed order. Each row corresponds to an independent grid cell, and the coordinates of the four vertices of the corresponding grid cell are marked in the row. Each column corresponds to a type of physical field parameter. The value at each position in the matrix represents the real-time physical field magnitude at the center of the corresponding grid cell. The strength of the physical field is directly determined by the magnitude of the values ​​in the matrix. The larger the temperature value, the higher the local temperature in the wellbore; the larger the pressure value, the higher the local fluid pressure in the wellbore; and the larger the strain value, the more severe the deformation of the wellbore structure. Comparison of values ​​within the same column can directly distinguish the high and low distribution of the physical field. Abrupt changes in values ​​indicate drastic changes in the physical field, while stable values ​​indicate a uniform distribution of the physical field. Through the correspondence between spatial coordinates and physical field values ​​in the matrix, the overall change law of the physical field inside the wellbore can be reflected. That is, the pressure and strain values ​​are generally higher in the region near the center point of the packer, the temperature and pressure values ​​gradually decrease in the region near the center point of the wellhead flange, the region at the top center point of the perforated section shows abrupt changes in the physical field, and the physical field in the region at the midpoint of the casing coupling maintains a stable transition. This fully reflects the continuous change law of the physical field with the position of the key structure of the wellbore, spatial depth, and radial distance. Multiple sets of physical field distribution matrices are generated at different production times. By comparing the numerical changes of the same grid cell at different times, the dynamic evolution of the physical field with production time and downhole conditions can be reflected. That is, the temperature and pressure rise rapidly in the early stage of production, the physical field remains stable in the middle stage of production, and the strain value increases slowly with structural fatigue in the later stage of production. Finally, a physical field distribution matrix is ​​constructed that can completely, accurately and continuously characterize the spatial distribution law and spatiotemporal variation characteristics of the physical field in the entire domain of the spatial tetrahedral associated domain. This enables the complete expression and visualization support of the physical field inside the wellbore from single-point data to the entire spatial domain and from static values ​​to dynamic changes.

[0037] This embodiment achieves adaptive mesh generation through boundary curve sampling, curvature analysis, and recursive subdivision. It can densify the mesh in areas with complex wellbore structures and drastic surface changes, improving local analysis accuracy, while simplifying the mesh in flat areas to balance computational efficiency. Based on spatial interpolation to obtain the physical field parameters of the mesh center point, it can realize the mapping of discrete monitoring data to a continuous spatial field. The constructed physical field distribution matrix can comprehensively and accurately reflect the spatial distribution law of temperature, pressure, and strain in key areas of the wellbore.

[0038] In a preferred embodiment of the present invention, step 4 includes: Step 400: Extract the temperature, pressure, and strain values ​​of all grid cells within the spatial tetrahedral correlation domain from the physical field distribution matrix. Calculate the mean and standard deviation of the temperature, pressure, and strain values ​​respectively. Determine the temperature variation coefficient based on the ratio of the standard deviation to the mean of the temperature value, the pressure variation coefficient based on the ratio of the standard deviation to the mean of the pressure value, and the strain variation coefficient based on the ratio of the standard deviation to the mean of the strain value. Normalize the temperature, pressure, and strain variation coefficients to obtain the temperature weighting coefficient, pressure weighting coefficient, and strain weighting coefficient. The force weighting coefficient and strain weighting coefficient specifically involve: extracting data row by row from the physical field distribution matrix according to the unique numbering of the tetrahedral grid cells, and separately aggregating the real-time temperature, pressure, and strain values ​​of all grid cells within the spatial tetrahedral associated domain, forming three independent and complete physical field datasets. The reason for extracting and independently aggregating the three types of data row by row according to the grid cell number is that each row in the physical field distribution matrix corresponds to an independent grid cell, and temperature, pressure, and strain are listed as independent parameters; independent aggregation avoids confusion between different physical field data. Statistical calculations are performed on the temperature dataset to obtain the arithmetic mean and standard deviation of all temperature values; statistical calculations are performed on the pressure dataset to obtain the arithmetic mean and standard deviation of all pressure values; and statistical calculations are performed on the strain dataset to obtain the arithmetic mean and standard deviation of all strain values. Based on the above statistical results, the coefficient of variation was calculated. The temperature coefficient of variation was calculated as the ratio of the temperature standard deviation to the temperature average, i.e., temperature coefficient of variation = temperature standard deviation ÷ temperature average. The pressure coefficient of variation was calculated as the ratio of the pressure standard deviation to the pressure average, i.e., pressure coefficient of variation = pressure standard deviation ÷ pressure average. The strain coefficient of variation was calculated as the ratio of the strain standard deviation to the strain average, i.e., strain coefficient of variation = strain standard deviation ÷ strain average. The coefficient of variation was calculated by dividing the standard deviation by the average because the dimensions of the average and standard deviation are consistent with the corresponding physical fields. The ratio of the two can eliminate the influence of dimensions and take into account both the average level and the degree of dispersion of the physical field, objectively reflecting the relative fluctuation of various physical fields. That is, the larger the coefficient of variation, the more violent the fluctuation of the physical field, and the more significant the impact on the local damage of the tubular structure.

[0039] The temperature, pressure, and strain coefficients of variation are all normalized to eliminate dimensional and numerical differences among the three parameters. This results in the calculation of temperature, pressure, and strain weighting coefficients. These three weighting coefficients objectively reflect the proportion of each physical field's influence on the tubular structure damage. Normalization maps the three coefficients of variation to the same numerical range, ensuring the sum of the calculated weighting coefficients is 1. This accurately and objectively reflects the degree of influence of each physical field on the tubular structure damage. After normalization, a larger weighting coefficient for a particular physical field indicates a relatively larger coefficient of variation, more drastic spatial distribution fluctuations, a more significant local damage effect on the tubular structure, and a higher proportion of damage contribution. Conversely, a smaller weighting coefficient indicates a more uniform distribution and gentler fluctuations of the physical field, a milder impact on the tubular structure damage, and a lower proportion of damage contribution. Therefore, the specific values ​​of the weighting coefficients clearly quantify the degree of influence of each physical field on the tubular structure damage.

[0040] Step 401: The average temperature value of all grid cells within the spatial tetrahedral domain is used as the temperature reference value; the average pressure value of all grid cells within the spatial tetrahedral domain is used as the pressure reference value; and the average strain value of all grid cells within the spatial tetrahedral domain is used as the strain reference value. Based on the temperature weighting coefficient, pressure weighting coefficient, and strain weighting coefficient, as well as the temperature reference value, pressure reference value, and strain reference value, the local coupling damage index of each grid cell is calculated. The local coupling damage index characterizes the coupled influence of temperature, pressure, and strain at the corresponding grid cell on the damage to the tubing, specifically including: The arithmetic mean of the temperature values ​​of all tetrahedral mesh elements within the spatial tetrahedral associated domain is defined as the temperature reference value of the wellbore string; the arithmetic mean of the pressure values ​​of all tetrahedral mesh elements is defined as the pressure reference value of the wellbore string; and the arithmetic mean of the strain values ​​of all tetrahedral mesh elements is defined as the strain reference value of the wellbore string. Using the temperature weighting coefficient, pressure weighting coefficient, and strain weighting coefficient calculated in step 400 as weighted coupling parameters, and combining them with the aforementioned temperature reference value, pressure reference value, and strain reference value, multiphysics coupling calculations are performed on each independent tetrahedral mesh element to obtain the local coupling damage index of that mesh element. That is, the local coupling damage index = temperature weighting coefficient × temperature reference value + pressure weighting coefficient × pressure reference value + strain weighting coefficient × strain reference value. This index is a quantitative evaluation parameter; its specific value directly characterizes the intensity and degree of coupling damage caused by the combined action of temperature, pressure, and strain on the wellbore string structure at the location of the corresponding tetrahedral mesh element. A value greater than 0.8 is considered high. The damage level indicates that the pipe structure at the location of the grid cell has suffered extremely severe coupling damage, posing a high risk of failure and requiring immediate inspection and maintenance. A value between 0.4 and 0.8 indicates a moderate damage level, meaning the pipe at that location has suffered some degree of coupling damage, requiring close monitoring and regular surveillance. A value between 0.1 and 0.4 indicates a low damage level, meaning the pipe at that location has suffered relatively minor coupling damage and the structure is relatively safe. A value less than 0.1 and close to 0 indicates that the pipe at that location has suffered almost no coupling damage and the structure is in good safety condition. Through a clear numerical grading standard, the degree of damage to the pipes at each grid cell can be determined, and the key areas of damage to the pipes can be located.

[0041] Step 402: Multiply the local coupling damage index of all grid cells within the spatial tetrahedral association domain by the volume of the corresponding grid cell to obtain the damage contribution volume value of each grid cell. Then, sum the damage contribution volume values ​​of all grid cells and calculate the ratio with the total volume of the spatial tetrahedral association domain to obtain the global equivalent damage factor. Specifically, this includes: Damage contribution volume value = Local coupling damage index × Tetrahedral grid cell volume. The core logic of this calculation is to quantify the contribution of a single grid cell to the overall damage of the pipe column by combining the spatial volume weight of the grid cells. The higher the local coupling damage index of the grid cell, the greater the contribution of the individual grid cell to the overall damage. The larger the volume, the greater its contribution to the overall damage of the pipe column. The specific value of the damage contribution volume directly corresponds to the total damage borne by a single mesh element. By successively summing the damage contribution volumes of all tetrahedral mesh elements, the total damage contribution volume of the spatial tetrahedral associated domain is obtained. The magnitude of the total damage contribution volume directly reflects the degree of total damage accumulation of the pipe column structure within the entire associated domain. The larger the value, the greater the total amount of damage accumulated in the pipe column within the associated domain, and the more severe the basis of overall damage. The smaller the value, the less the total amount of damage accumulated in the pipe column within the associated domain, and the less severe the basis of overall damage.

[0042] The ratio of the total damage contribution volume to the total volume of the spatial tetrahedral associated domain is calculated to obtain the global equivalent damage factor, i.e., global equivalent damage factor = total damage contribution volume ÷ total volume of the spatial tetrahedral associated domain. This parameter is a comprehensive quantitative evaluation index, and its specific value can comprehensively and objectively characterize the overall damage degree and overall safety status of the entire wellbore string structure under the coupling effect of multiple physical fields. Combined with the downhole string engineering evaluation standard, it is clear that when the global equivalent damage factor ≥ 0.7, it indicates that the overall comprehensive damage degree of the string is extremely high, and the structural safety redundancy is extremely low. When the system is nearing its safety limit, an immediate shutdown and comprehensive overhaul are required. When the global equivalent damage factor is 0.4 ≤ global equivalent damage factor < 0.7, the overall damage level of the tubing is high, the structural safety redundancy is low, and there are obvious safety hazards. The monitoring cycle needs to be shortened and a targeted maintenance plan needs to be developed. When the global equivalent damage factor is 0.1 ≤ global equivalent damage factor < 0.4, the overall damage level of the tubing is low, the structural safety redundancy is high, and only routine periodic monitoring is required. When the global equivalent damage factor is < 0.1, the overall damage level of the tubing is extremely low, the structural safety redundancy is extremely high, and the system is in a good safety condition.

[0043] This embodiment achieves systematic extraction and quantitative analysis of multiple physical field parameters based on the physical field distribution matrix. By calculating the mean, standard deviation, coefficient of variation, and weighting coefficient step by step, it distinguishes the contribution of temperature, pressure, and strain to tubing damage, avoiding the assessment bias caused by single physical field analysis and improving the accuracy and comprehensiveness of tubing condition evaluation. It constructs a two-layer evaluation system of local coupled damage index and global equivalent damage factor. The local index can accurately locate the damage concentration area and weak point of the wellbore tubing, while the global factor can intuitively reflect the comprehensive damage state of the overall tubing structure, realizing a full-dimensional damage assessment from local details to overall overview. By combining the tetrahedral grid unit volume for weighted calculation, the damage evaluation results are highly matched with spatial geometric features, which can truly reflect the actual stress and deformation state of the tubing structure under complex downhole conditions, providing a reliable quantitative basis for wellbore tubing safety monitoring, damage early warning, and maintenance decisions.

[0044] In a preferred embodiment of the present invention, step 5 includes: Step 500: Determine the cumulative effect coefficient and instantaneous effect coefficient based on the time-series variation trend of the global equivalent damage factor, and construct the initial mapping relationship for the evolution equation of the remaining strength of the tubing based on the cumulative effect coefficient and instantaneous effect coefficient. Substitute the global equivalent damage factor at the current moment into the initial mapping relationship for iterative correction to obtain the evolution equation of the remaining strength of the tubing. Specifically, this includes: firstly, collecting global equivalent damage factor data in continuous time series, with a collection interval of 1 hour, and the collection duration covering at least 3 complete production cycles, with a single cycle duration of not less than 7 days; performing moving average smoothing on the collected time series data to eliminate random noise and abnormal fluctuation interference; and decomposing the change of the global equivalent damage factor into a long-term cumulative effect component and a short-term instantaneous effect component through piecewise trend fitting. The quantities are as follows: the long-term cumulative impact component is defined as the trend of change over 30 consecutive days with a change rate ≤ 0.001 / day, corresponding to the linear change portion of the tubing string under long-term production, corrosion, and fatigue. The short-term instantaneous impact component is defined as the change portion within 24 hours with a change rate > 0.01 / hour, corresponding to the nonlinear instantaneous change portion caused by changes in operating conditions such as production system switching, pressure shock, and sudden flow changes. The cumulative effect coefficient is calculated using the long-term cumulative impact component, i.e., cumulative effect coefficient = long-term cumulative impact component ÷ cumulative production time. The instantaneous effect coefficient is calculated using the short-term instantaneous impact component, i.e., instantaneous effect coefficient = short-term instantaneous impact component ÷ operating condition fluctuation amplitude, where the operating condition fluctuation amplitude is the maximum value of the difference between the casing pressure and the oil pressure at the current moment and the previous moment.

[0045] Using the cumulative effect coefficient (denoted as α) and the instantaneous effect coefficient (denoted as β) as the core modeling parameters, and combining the time variable (denoted as t) and the global equivalent damage factor at the current moment (denoted as D), an initial mapping relationship for the evolution equation of the remaining strength of the tubing string is constructed. The formula for the initial mapping relationship is: R(t,D)=R0×(1-α×t-β×D), where R(t,D) is the remaining strength of the tubing string at time t, R0 is the initial strength of the tubing string, which is taken from the material performance parameters of the tubing string at the factory and the well completion structure design parameters, and is a known fixed value; α is the cumulative effect coefficient, β is the instantaneous effect coefficient, t is the cumulative production time, and D is the global equivalent damage factor at the current moment. The actual remaining strength of the tubing string is measured using downhole testing instruments such as a wellbore caliper and electromagnetic thickness gauge. This value is then converted to an equivalent remaining strength value based on the tubing material strength standard. Specifically, the equivalent remaining strength value is calculated by combining the tensile strength and yield strength parameters specified in the tubing material strength standard: Equivalent Remaining Strength = Initial Strength of Tubing String × (Current Effective Wall Thickness of Tubing String ÷ Initial Wall Thickness of Tubing String). The predicted value is directly calculated from the current mapping relationship. The current global equivalent damage factor is substituted into the initial mapping relationship, and the iteration count is set to be no less than 10. Each iteration... The absolute deviation between the predicted value and the measured residual strength is calculated by the algebraic average. The deviation calculation formula is: Deviation = |Predicted residual strength - Actual measured residual strength|. When the deviation is greater than 0.01, the cumulative effect coefficient α and the instantaneous effect coefficient β are simultaneously corrected according to the deviation ratio. The correction method is to decrease the coefficient by a step size of 0.001 when the deviation increases and increase the coefficient by a step size of 0.001 when the deviation decreases, thus gradually reducing the prediction error. When the deviation after iterative calculation is less than the preset allowable deviation of 0.01, the iteration terminates, and finally the residual strength evolution equation of the tubing string that is adapted to the actual working conditions of the wellbore and meets the fitting accuracy standard is obtained.

[0046] Step 501: Input the globally equivalent damage factor calculated at the current moment into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment; compare the wall thickness loss rate and residual strength coefficient at the current moment with the safety thresholds for the wall thickness loss rate and residual strength coefficient, respectively, to determine the deviation of the wall thickness loss rate and the residual strength coefficient at the current moment. Specifically, this includes: using the globally equivalent damage factor (D) calculated at the current moment as the only input parameter, substituting it into the residual strength evolution equation of the tubing obtained in step 500 for solving, and directly calculating the wall thickness loss rate (denoted as η) and residual strength coefficient (denoted as K) of the tubing at the current moment, where the calculation formulas for the wall thickness loss rate and residual strength coefficient are respectively, η=1- K= In the formula, η is the current wall thickness loss rate of the tubing, K is the current remaining strength coefficient of the tubing, and R(t,D) is the current remaining strength of the tubing. The initial strength of the tubing string is defined. Preset safety thresholds for wall thickness loss rate and residual strength coefficient, conforming to downhole tubing string engineering safety standards. The wall thickness loss rate safety threshold is set to 0.15 (15%), and the residual strength coefficient safety threshold is set to 0.85 (85%). Both thresholds are determined based on the tubing string material's tolerance limit, well completion safety requirements, and field engineering experience, serving as fixed safety judgment standards. The difference between the current wall thickness loss rate and the wall thickness loss rate safety threshold is calculated to obtain the deviation degree of the wall thickness loss rate; the difference between the current residual strength coefficient and the residual strength coefficient safety threshold is calculated to obtain the deviation degree of the residual strength coefficient. A positive deviation degree for the wall thickness loss rate indicates that the current wall thickness loss rate exceeds the safety threshold, and the larger the value, the more severe the deviation from the safety standard. A positive deviation degree for the residual strength coefficient indicates that the current residual strength coefficient is below the safety threshold, and the larger the value, the more severe the deviation from the safety standard. Negative values ​​for both indicate that the current parameters are within the safe range.

[0047] Step 502: Based on the deviation of the wall thickness loss rate and the deviation of the residual strength coefficient, and combined with the current casing pressure time series value, oil pressure time series value, and instantaneous production volume time series value, calculate the production system parameters to be adjusted and the adjustment amount. The production system parameters include production differential pressure, pumping unit strokes, gas lift injection volume, or electric submersible pump frequency. Based on the adjustment amount, determine the intelligent decision information containing control instructions, specifically including: establishing a production system parameter optimization calculation model, that is, clarifying the optimization objective, with minimizing the deviation of the wall thickness loss rate and the deviation of the residual strength coefficient as the core optimization objective, constructing an objective function, the objective function being to minimize (the deviation of the wall thickness loss rate + the deviation of the residual strength coefficient), ensuring that the safety risk of the tubing structure is minimized after adjustment; combined with the actual wellbore production control parameters, take the production differential pressure, pumping unit strokes, gas lift injection volume, and electric submersible pump frequency as core optimization variables, clarify the value range of each variable (production differential pressure 0.5 to 3.0 MPa, pumping unit strokes 3 to 10 times / minute, gas lift injection volume 5 to 20 m³ / min). 3 / h, ESP frequency 30 to 60Hz), to avoid variables exceeding the equipment's operating limits; set constraints, using the current collected casing pressure time series value, oil pressure time series value, and instantaneous fluid production time series value as operating condition constraints, and clearly define the constraint range (casing pressure 0.5 to 2.5MPa, oil pressure 1.0 to 3.0MPa, instantaneous fluid production 10 to 50m³ / h). 3 / h), ensuring that the adjusted production regime matches the current downhole conditions and does not affect production continuity; combining the linear characteristics of the optimization objectives, variables, and constraints, a linear programming algorithm is used as the model solution algorithm. Through iterative algorithm solution, the adjustment amount of the production regime parameters that meet the constraints and achieve the optimization objectives is obtained. The specific adjustment amount calculation formula is as follows: Production pressure difference adjustment amount = k1 × wall thickness loss rate deviation + k2 × residual strength coefficient deviation; Pumping unit stroke adjustment amount = k3 × wall thickness loss rate deviation + k4 × residual strength coefficient deviation; Gas lift injection amount adjustment amount = k5 × wall thickness loss rate deviation + k6 × residual strength coefficient deviation; Electric submersible pump frequency adjustment amount = k7 × wall thickness loss rate deviation + k8 × residual strength coefficient deviation; k1 to k8 are adjustment coefficients, set according to field production experience, the degree of tubing damage impact, and parameter control sensitivity, with k1=0. 2. Because the impact of production pressure differential on wall thickness loss is moderate, it needs to be adjusted appropriately to balance damage control and production efficiency; k2=0.3, because the impact of production pressure differential on residual strength is more significant, the weight needs to be increased to ensure strength safety; k3=0.1, because the impact of pumping unit strokes on tubing damage is relatively mild, the adjustment range should not be too large; k4=0.2, because the impact of pumping unit strokes on residual strength is slightly higher than that on wall thickness loss, the weight needs to be increased appropriately; k5=0.15, because the impact of gas lift injection volume on tubing damage is moderate, and gas lift efficiency needs to be considered, the adjustment range needs to be moderate; k6=0.25, because fluctuations in gas lift injection volume can easily lead to pressure shocks, which have a significant impact on residual strength, the weight needs to be increased; k7=0.1, because the impact of ESP frequency on tubing damage is minimal, only slight adjustment is needed; k8=0.2, because changes in ESP frequency will indirectly affect the stress on the tubing, which has a certain impact on residual strength, the adjustment needs to be moderate.

[0048] Based on the calculated adjustment amounts of various production system parameters, corresponding control instructions are generated. Each control instruction specifies the target parameters, such as production differential pressure, pumping unit strokes, adjustment range (i.e., the calculated adjustment amount), and execution timing. The adjustment instruction is executed within 10 minutes of being issued, and the execution standard is verified after execution. Relevant parameters are collected 30 minutes after execution to verify the adjustment effect. All control instructions are sorted by priority, with higher priority for the greater the deviation. This information is integrated to form complete intelligent decision-making information, which is directly used to guide the automated control and safety management of the wellbore production system. This ensures that the adjusted production system can effectively reduce tubing damage and guarantee tubing safety and production stability.

[0049] This embodiment constructs a residual strength evolution equation based on the temporal variation of the global equivalent damage factor, enabling accurate prediction of tubing structural damage from static evaluation to dynamic evolution. It can predict the strength decay pattern of the tubing in advance, improving the foresight and reliability of safety warnings. By quantifying the structural safety risk through the deviation between the wall thickness loss rate and the residual strength coefficient, and combining real-time production condition data to optimize production system parameters, a dynamic balance between structural safety and production efficiency is achieved, avoiding overproduction that exacerbates tubing damage. A full-process intelligent decision-making system is formed, from damage assessment and strength prediction to parameter control, which can automatically output control commands and adjustments without manual intervention, improving the intelligence level and response speed of wellbore production management. By comprehensively considering the cumulative and instantaneous effects of damage, and balancing the long-term structural safety and short-term production stability of the tubing, the system effectively slows down the rate of tubing damage development, extends the service life of the tubing, and reduces downhole maintenance costs and safety risks.

[0050] like Figure 2 As shown, embodiments of the present invention also provide a real-time intelligent analysis system for oil and gas well production data, including: The module is used to collect and preprocess multi-source time-series data from oil and gas wells to obtain a standard dataset. From the standard dataset, the real-time difference rate of change of casing pressure and oil pressure, the short-time fluctuation amplitude of instantaneous production, the centroid offset of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference of wellhead and bottom hole temperature are extracted to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions. The completion structure parameters and tubing string geometric parameters of the oil and gas well are obtained to construct a three-dimensional digital model of the wellbore. The analysis module is used to perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; and to construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore using the first, second, third, and fourth nucleation base points as vertices. The partitioning module is used to discretize the spatial tetrahedral domain into several grid cells and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. The calculation module is used to calculate the local coupling damage index of each grid cell, and to perform a volume-weighted integral based on the local coupling damage indices of all grid cells to obtain the global equivalent damage factor. The adjustment module is used to input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision information containing control instructions.

[0051] It should be noted that this system is a system corresponding to the above method. All implementation methods in the above method embodiments are applicable to this embodiment and can achieve the same technical effect.

[0052] Embodiments of the present invention also provide a computer-readable storage medium storing instructions that, when executed on a computer, cause the computer to perform the method described above. All implementations in the above method embodiments are applicable to this embodiment and can achieve the same technical effects.

[0053] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A real-time intelligent analysis method for oil and gas well production data, characterized in that, The method includes: Step 1: Collect and preprocess multi-source time-series data of oil and gas wells to obtain a standard dataset. Extract the real-time difference rate of change of casing pressure and oil pressure, the short-term fluctuation amplitude of instantaneous production, the centroid offset of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference of wellhead and bottom hole temperature from the standard dataset to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions. Obtain the completion structure parameters and tubing string geometric parameters of the oil and gas well to construct a three-dimensional digital model of the wellbore. Step 2: Perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; using the first, second, third, and fourth nucleation base points as vertices, construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore; Step 3: Discretize the spatial tetrahedral domain into several grid cells, and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. Step 4: Calculate the local coupling damage index of each mesh element, and perform a volume-weighted integral based on the local coupling damage indices of all mesh elements to obtain the global equivalent damage factor. Step 5: Input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision-making information containing control instructions.

2. The real-time intelligent analysis method for oil and gas well production data according to claim 1, characterized in that, The multi-source time-series data of the oil and gas wells includes: casing pressure time-series value, oil pressure time-series value, instantaneous production time-series value, wellhead temperature time-series value, bottom hole temperature time-series value, and tubing string axial vibration time-series signal.

3. The real-time intelligent analysis method for oil and gas well production data according to claim 2, characterized in that, The well completion structural parameters and tubing string geometric parameters include: wellhead flange coordinates, casing coupling position coordinates, perforation section top depth coordinates, packer installation depth coordinates, casing outer diameter, casing inner diameter, tubing outer diameter, tubing inner diameter, and the elastic modulus and Poisson's ratio of the tubing string material.

4. The real-time intelligent analysis method for oil and gas well production data according to claim 3, characterized in that, Step 2 includes: Local extreme points are detected in the real-time difference rate of change of casing pressure and oil pressure in the multidimensional feature matrix. The extreme point features corresponding to the zero crossing of the first derivative are extracted. The first nucleation base point is determined according to the spatial position of the extreme point features in the three-dimensional digital model of the wellbore. The first nucleation base point is located at the center point of the wellhead flange. A joint mutation point analysis was performed on the multidimensional feature matrix to extract the synchronous anomaly segments of the real-time difference change rate sequence of casing pressure and oil pressure and the short-term fluctuation amplitude sequence of instantaneous production volume. The second nucleation base point was determined based on the spatial location of the starting time of the synchronous anomaly segment in the three-dimensional digital model of the wellbore. The second nucleation base point is located at the center point of the top of the downhole perforated section. Cumulative offset analysis is performed on the centroid offset sequence of the axial vibration signal spectrum of the tubing string in the multidimensional feature matrix. The cumulative offset extreme points are extracted. The third nucleation base point is determined according to the spatial position of the time corresponding to the cumulative offset extreme point in the three-dimensional digital model of the wellbore. The third nucleation base point is located at the midpoint of the casing coupling. Inflection point detection is performed on the real-time gradient difference sequence of wellhead and bottomhole temperatures in the multidimensional feature matrix. Gradient change inflection points are extracted. The fourth nucleation base point is determined based on the spatial location of the time corresponding to the gradient change inflection point in the three-dimensional digital model of the wellbore. The fourth nucleation base point is located at the center point of the packer position. Using the first, second, third, and fourth nucleation base points as vertices, a spatial tetrahedral relational domain is constructed in the three-dimensional digital model of the wellbore.

5. The real-time intelligent analysis method for oil and gas well production data according to claim 4, characterized in that, Step 3 includes: The spatial tetrahedral associative domain is discretized into multiple grid cells, and the node coordinates and cell volume of each grid cell are recorded. Based on a standard dataset and a three-dimensional digital model of the wellbore, a spatial interpolation method is used to calculate the real-time physical field parameters at the center point of each grid cell. The real-time physical field parameters include temperature, pressure, and strain. The real-time physical field parameters at the center point of all grid cells are mapped to the corresponding grid cells, and the physical field distribution matrix in the spatial tetrahedral associated domain is constructed based on the mapping results.

6. The real-time intelligent analysis method for oil and gas well production data according to claim 5, characterized in that, The spatial tetrahedral associative domain is discretized into multiple mesh elements, and the node coordinates and element volume of each mesh element are recorded, including: The twelve edges of the spatial tetrahedral associated domain are extracted as boundary curves. Each boundary curve is sampled at equal intervals to obtain a set of sampling points. The curvature value at each sampling point is calculated. The local maxima of curvature are detected from the curvature values ​​as the inflection points of the boundary curves. Each boundary curve is divided into multiple curve segments based on the inflection point, and the length of each curve segment is calculated. Based on the length of the curve segment and the spatial position of the inflection point, a parametric spline surface covering the four triangular faces of the spatial tetrahedral associated domain is constructed with the boundary curve as a constraint, so that the boundary of the corresponding spline surface fits the boundary curve. The parametric spline surface is recursively subdivided to generate a boundary surface mesh composed of multiple surface patches. The normal vector of each surface patch on the boundary surface mesh is calculated, and the surface curvature at each mesh vertex on the boundary surface mesh is estimated. Based on the estimated surface curvature distribution characteristics, regions with drastic surface changes are identified. According to the positional distribution of these regions on the boundary surface mesh, tetrahedral mesh elements are generated in the interior space of the spatial tetrahedral associated domain, so that the generated positions of the tetrahedral mesh elements correspond to the positions of the regions with drastic surface changes. The spatial coordinates of the four vertices of each tetrahedral mesh element in the 3D digital model of the well shaft are determined and used as the node coordinates of the corresponding mesh element. The element volume of each tetrahedral mesh element is calculated based on the spatial coordinates of the four vertices.

7. The real-time intelligent analysis method for oil and gas well production data according to claim 6, characterized in that, Step 4 includes: Temperature, pressure, and strain values ​​of all grid cells within the spatial tetrahedral correlation domain are extracted from the physical field distribution matrix. The mean and standard deviation of the temperature, pressure, and strain values ​​are calculated respectively. The temperature variation coefficient, pressure variation coefficient, and strain variation coefficient are determined based on the ratio of the standard deviation to the mean of the temperature value, the pressure variation coefficient, and the strain variation coefficient. The temperature variation coefficient, pressure variation coefficient, and strain variation coefficient are then normalized to obtain the temperature weighting coefficient, pressure weighting coefficient, and strain weighting coefficient. The average temperature value of all grid cells within the spatial tetrahedral domain is used as the temperature reference value; the average pressure value of all grid cells within the spatial tetrahedral domain is used as the pressure reference value; and the average strain value of all grid cells within the spatial tetrahedral domain is used as the strain reference value. Based on the temperature weighting coefficient, pressure weighting coefficient, and strain weighting coefficient, as well as the temperature reference value, pressure reference value, and strain reference value, the local coupling damage index of each grid cell is calculated. The local coupling damage index characterizes the coupled influence of temperature, pressure, and strain on the damage of the tubing at the corresponding grid cell. The local coupling damage index of all grid cells in the spatial tetrahedral association domain is multiplied by the volume of the corresponding grid cell to obtain the damage contribution volume value of each grid cell. The damage contribution volume values ​​of all grid cells are summed and then compared with the total volume of the spatial tetrahedral association domain to obtain the global equivalent damage factor.

8. The real-time intelligent analysis method for oil and gas well production data according to claim 7, characterized in that, Step 5 includes: The cumulative effect coefficient and the instantaneous effect coefficient are determined based on the temporal variation trend of the global equivalent damage factor. The initial mapping relationship of the residual strength evolution equation of the tubing is constructed based on the cumulative effect coefficient and the instantaneous effect coefficient. The global equivalent damage factor at the current moment is substituted into the initial mapping relationship for iterative correction to obtain the residual strength evolution equation of the tubing. Input the global equivalent damage factor calculated at the current moment into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment; compare the wall thickness loss rate and residual strength coefficient at the current moment with the wall thickness loss rate safety threshold and the residual strength coefficient safety threshold respectively to determine the deviation of the wall thickness loss rate and the residual strength coefficient at the current moment. Based on the deviation of the wall thickness loss rate and the deviation of the residual strength coefficient, combined with the current casing pressure time series value, oil pressure time series value and instantaneous production volume time series value, the production system parameters that need to be adjusted and the adjustment amount are calculated. The production system parameters include production pressure difference, pumping unit strokes, gas lift injection volume or electric submersible pump frequency; based on the adjustment amount, intelligent decision information containing control instructions is determined.

9. A real-time intelligent analysis system for oil and gas well production data, wherein the system implements the method as described in any one of claims 1 to 8, characterized in that, include: The module is used to collect and preprocess multi-source time-series data from oil and gas wells to obtain a standard dataset. From the standard dataset, the real-time difference rate of change of casing pressure and oil pressure, the short-time fluctuation amplitude of instantaneous production, the centroid offset of the axial vibration signal spectrum of the tubing string, and the real-time gradient difference of wellhead and bottom hole temperature are extracted to construct a multi-dimensional feature matrix characterizing the dynamic changes of downhole operating conditions. The completion structure parameters and tubing string geometric parameters of the oil and gas well are obtained to construct a three-dimensional digital model of the wellbore. The analysis module is used to perform local extremum point detection, joint mutation point analysis, cumulative offset analysis, and inflection point detection on the multidimensional feature matrix to determine the first, second, third, and fourth nucleation base points; and to construct a spatial tetrahedral association domain in the three-dimensional digital model of the wellbore using the first, second, third, and fourth nucleation base points as vertices. The partitioning module is used to discretize the spatial tetrahedral domain into several grid cells and map the real-time physical field parameters at the center point of each grid cell to all grid cells to construct the physical field distribution matrix within the spatial tetrahedral domain. The calculation module is used to calculate the local coupling damage index of each grid cell, and to perform a volume-weighted integral based on the local coupling damage indices of all grid cells to obtain the global equivalent damage factor. The adjustment module is used to input the global equivalent damage factor into the residual strength evolution equation of the tubing to obtain the wall thickness loss rate and residual strength coefficient at the current moment, so as to dynamically adjust the production system parameters and obtain intelligent decision information containing control instructions.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 8.

Images