A method for early warning of ice blockage risk of natural gas hydrates based on a gas-solid equilibrium thermodynamic model

Abstract

This invention discloses a method for early warning of natural gas hydrate ice blockage risk based on a gas-solid equilibrium thermodynamic model, belonging to the field of oil and gas production safety technology. The method includes: acquiring natural gas component data and pressure and temperature time-series data after throttling from a production data acquisition and monitoring system; preprocessing these data to obtain synchronized pressure-temperature sequence pairs; substituting the component data and pressure sequence into the gas-solid equilibrium constant equation and solving for the stability condition to obtain the hydrate formation temperature sequence; calculating the temperature safety margin and comparing it with a preset classification threshold to obtain the risk classification result and risk score; further identifying and blocking drainage periods, performing continuous boundary violation confirmation, and outputting the final alarm level and risk score. This invention balances prediction accuracy and real-time performance, achieving graded alarms and resistance to false alarms, and is applicable to natural gas throttling and transportation processes.

Description

Technical Field

[0001] This invention relates to the field of oil and gas production safety technology, specifically to a method for early warning of natural gas hydrate ice blockage risk based on a gas-solid equilibrium thermodynamic model. Background Technology

[0002] Natural gas hydrates are solid crystals with an ice-like appearance and a cage-like internal structure, formed by guest molecules such as methane, ethane, propane, and hydrogen sulfide in natural gas and water molecules under high pressure and low temperature thermodynamic conditions. They are also known as cage-like hydrates. During natural gas extraction and transportation, throttling and disturbance points such as throttling valves, needle valves, elbows, and orifice flow meters experience drastic temperature drops due to the Joule-Thomson effect. Once the fluid temperature drops below the hydrate formation temperature under the given pressure conditions, hydrates will nucleate, grow, and gradually block pipelines and valves. This can lead to reduced gas transmission capacity and increased pressure drop, or even freezing and pressure buildup after production shutdowns, causing safety incidents. Therefore, online calculation and graded early warning of hydrate formation risks in throttling and transportation processes are key technical requirements for safe winter production and unmanned operation of oil and gas fields. Summary of the Invention

[0003] To address the problems existing in the prior art, this invention provides a method for early warning of ice blockage risk of natural gas hydrates based on a gas-solid equilibrium thermodynamic model.

[0004] To achieve the above objectives, the present invention provides the following technical solution: a method for early warning of ice blockage risk of natural gas hydrates based on a gas-solid equilibrium thermodynamic model, comprising the following steps:

[0005] S1. Acquire natural gas component data, pressure time-series data after throttling, and temperature time-series data after throttling; perform unit normalization, timestamp alignment, missing value imputation, and outlier removal on the pressure time-series data and the temperature time-series data to obtain synchronized pressure-temperature sequence pairs.

[0006] S2. Based on the natural gas component data and the pressure sequence in the synchronized pressure-temperature sequence pair, the hydrate formation temperature at each time point is calculated using a gas-solid equilibrium thermodynamic model to obtain the hydrate formation temperature sequence.

[0007] S3. Based on the temperature sequence in the hydrate formation temperature sequence and the synchronous pressure-temperature sequence pair, calculate the temperature safety margin at each time point to obtain a temperature safety margin sequence; compare each temperature safety margin in the temperature safety margin sequence with a preset classification threshold to obtain a risk classification result, and map the risk classification result to a risk score.

[0008] S4. Identify the drainage period and mask the data points within the drainage period to obtain the masked risk classification result; perform continuous boundary verification on the masked risk classification result, and generate and output the final alarm level and the risk score after verification.

[0009] The steps of this invention are connected end to end with data as the clue: the synchronous pressure-temperature sequence obtained in step S1 is both the input for the hydrate formation temperature calculation in step S2 and the source of the actual temperature data in step S3; the hydrate formation temperature sequence obtained in step S2 is used in step S3 to calculate the temperature safety margin; the risk classification result obtained in step S3 is used in step S4 as the object of shielding and confirmation, and the risk score is corrected in the confirmation process and output along with the final alarm level as a quantitative basis for alarm sorting and handling priority.

[0010] The beneficial effects of this invention are:

[0011] This invention uses a gas-solid equilibrium thermodynamic model as the core criterion and sets up a trainable correction network to absorb systematic biases. The prediction accuracy of hydrate formation temperature is significantly better than that of empirical chart methods and empirical correlation methods, while retaining the extrapolation capability and physical interpretability of the mechanistic model. It provides both analytical fast path for pure methane and numerical root-finding path for multi-component working conditions, taking into account the real-time performance of online high-frequency calculations and the universality of general working conditions. It also provides deterministic handling rules for abnormal situations such as root-finding failure and no roots in the interval, ensuring the continuity of online services.

[0012] This invention uses temperature safety margin as a unified standard to achieve a three-level risk classification of normal, warning, and severe, which is mapped to a quantitative risk score. The alarm level directly corresponds to the intensity of treatments such as alcohol injection and heating, supporting graded treatment. The liquid discharge period shielding mechanism eliminates the false temperature drop caused by separator liquid discharge operations, and the tail continuous cross-boundary confirmation mechanism suppresses false alarms caused by single-point shaking. Together, they significantly reduce the false alarm rate and do not require the introduction of smoothing filters that would delay risk exposure. Multi-point collaborative judgment and quantitative evaluation of treatment effects form a complete closed loop of early warning-treatment-evaluation, providing a directly implementable technical means for ice blockage risk management in unmanned stations. Attached Figure Description

[0013] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0014] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0015] Please see Figure 1 This invention provides a method for early warning of ice blockage risk of natural gas hydrate based on a gas-solid equilibrium thermodynamic model. The method includes steps S1 to S4. The entire method can be deployed on the edge computing device of an oil and gas field station or on a production data analysis server. It pulls recent data from the production data acquisition and monitoring system at a set period (e.g., every 5 minutes) and performs a complete calculation.

[0016] S1. Acquire natural gas component data, pressure time-series data after throttling, and temperature time-series data after throttling; perform unit normalization, timestamp alignment, missing value imputation, and outlier removal on the pressure time-series data and the temperature time-series data to obtain synchronized pressure-temperature sequence pairs.

[0017] The input data for this step comes from two sources. One is the Supervisory Control and Data Acquisition (SCADA) system for oil and gas field production. This system acquires the time-series pressure and temperature data after throttling at the target monitoring points within a set time range. The target monitoring points are, for example, pressure transmitter measuring points and temperature transmitter measuring points installed in the downstream pipe section of the throttling valve.

[0018] Each time series data point is represented by a binary sequence consisting of a sampling time and a measured value: pressure time series data is denoted as D. P = {(τ k , p k Temperature time series data is denoted as D. T = {(τ k ′, t k )}, where τ is the sampling time (millisecond-level timestamp, the same below), p is the pressure measurement value, t is the temperature measurement value, and k is the data point sequence number. The sampling period is determined according to the on-site acquisition configuration, for example, 10 seconds or 20 seconds; the number of historical cache points obtained in a single round of calculation is, for example, 10080 points (approximately 2.8 days based on 20-second sampling), of which the most recent window for judgment is the most recent 360 points (approximately 2 hours). In addition, when liquid level location codes are configured, the liquid level time series data D of that measurement point is acquired synchronously. L = {(τ k ", l k (corresponding to the liquid level of the separator or metering tank), its purpose is explained in step S4.

[0019] The second is natural gas component analysis data, obtained from gas sample analysis or chromatographic analysis reports. This data is provided as a string format: "Component Name: Mole Percentage," such as "CH4:96.2, C2H6:1.5, C3H8:0.5, CO2:0.8, N2:1.0," where each value is a mole percentage (%). Component data can be updated periodically during periods of stable gas supply. A summary of the data sources, content, and formats involved in this step is shown in Table 1.

[0020] Table 1. Input Data Sources and Formats

[0021] Pressure time series data after throttling Production data acquisition and monitoring system Pressure transmitter measuring point downstream of throttle valve MPa (or kPa / bar / Pa), time series, sampling period 10~20s Temperature time series data after throttling Production data acquisition and monitoring system Temperature transmitter measuring point downstream of throttle valve ℃ (or K), time series, sampling period 10~20s Liquid level time series data (optional) Production data acquisition and monitoring system Separator / metering tank level measuring point % range, timing sequence, for drainage shielding only. Natural gas composition data Gas sample analysis or chromatographic analysis Laboratory report string Mole percentage (%), updated periodically Traffic data (optional) Production data acquisition and monitoring system Flow meter measuring points m³ / d, used only to correct network characteristics

[0022] The above component strings are split and parsed into mole percentages (v) of each component according to the delimiter. i (i iterates through components such as methane, ethane, propane, isobutane, n-butane, carbon dioxide, hydrogen sulfide, and nitrogen). Due to potential rounding errors in laboratory reports, the sum of the components may not strictly equal 100. Furthermore, the calculation only involves known components from the coefficient table. Therefore, the known components used in the calculation are normalized to mole fractions using the following formula:

[0023] y i = v i / Σ j v j ;

[0024] Where j iterates through all known components involved in the calculation, satisfying Σ i y i = 1.

[0025] The normalized component set G = {(i, y} is obtained. i This set will be used as input to the gas-solid equilibrium constant equation in step S2. Unknown components (components not present in the coefficient table) are not included in the normalization denominator to avoid distorting the calculation results.

[0026] Thermodynamic models require temperature to be measured in Kelvin (K) and pressure in megapascals (MPa), but field-reported values ​​may be measured in degrees Celsius (°C), kilopascals (kPa), bars (bar), or pascals (Pa). Therefore, it is necessary to standardize the dimensions first.

[0027] Temperature conversion: If the unit is set to Celsius, convert it to Kelvin using T = t + 273.15; if the unit is set to Kelvin, it remains unchanged; if no unit is set, the median is used for automatic determination: calculate the median med(t) of the temperature series, and if med(t) > 150, determine that the series is already in Kelvin (the temperature of hydrates cannot exceed 150℃), otherwise determine it to Celsius and add 273.15 for conversion.

[0028] Pressure conversion: P(MPa) = P(kPa) / 1000 = P(bar) / 10 = P(Pa) / 10 6 The different reporting units will be unified under the term "megapascal".

[0029] We obtain a uniform pressure and temperature sequence, denoted as P′ = {(τ k , p k )} (MPa) and T′ = {(τ) k ′,t k )}(K).

[0030] Because pressure and temperature measurement points report independently, the sampling time sets of the two sequences are not completely consistent, and directly pairing them by row number will result in misalignment. Furthermore, it is difficult to achieve strict millisecond-level synchronization between the reporting clocks of the two independent measurement points. If only strict equality matching is performed, a large number of time points will not be able to be paired, and temperature values ​​will be set to missing values. Therefore, a tolerance-tolerant nearest neighbor alignment is performed based on the timestamp of the pressure sequence for each τ. k Perform the following three levels of matching in sequence:

[0031] Exact match: Finding a match with τ in a temperature series k The values ​​must be taken at the same time; if a match is found, the pairing is performed directly.

[0032] Tolerance-based nearest neighbor matching: When a match is not found, if |τ′ -τ k | ≤ δ temperature time selected with τ k Pair the samples with the smallest time difference. The tolerance δ is half a sampling period (e.g., δ = 5 seconds when sampling for 10 seconds and δ = 10 seconds when sampling for 20 seconds). The upper limit of half a period ensures that mismatched pairings across sampling periods will not occur.

[0033] No match found: If there is no temperature moment within the tolerance window, set the temperature value of that moment as missing and perform missing value imputation;

[0034] In cases where multiple records appear at the same timestamp in any sequence (e.g., duplicate reporting caused by communication retransmission), the last record with the same timestamp is retained.

[0035] After processing, coarsely aligned sequence pairs are obtained:

[0036] F = {(τ k , p k , t k )}, of which part t k This could be a missing value.

[0037] Missing temperature values ​​in coarsely aligned sequence pairs are filled using linear interpolation:

[0038] Let the missing time τm The most recent valid times before and after are τ a τ b (corresponding value t) a t b ),but:

[0039] t̂′ m = t a + (t b -t a ) × (τ m -τ a ) / (τ b -τ a );

[0040] For missing data at both ends of the sequence where it is impossible to obtain points from both sides, the nearest valid value from the forward or backward direction is used for filling. When the number of consecutive missing points exceeds a preset missing limit (e.g., 2 points), the data is considered broken, and no imputation is performed, but the missing data is retained and marked to avoid introducing false data through long-segment interpolation. Isolated missing data in stress sequences are handled in the same way. After processing, the imputed sequence pair F′ is obtained.

[0041] Then, four types of elimination are performed in sequence:

[0042] Non-positive pressure elimination: Thermodynamic equations require pressure to be positive, pressure p k Points with a value ≤ 0 (usually bad points or points that have stopped transmitting and have been set to zero) are directly discarded;

[0043] Physical range out-of-range values ​​are eliminated: Points that exceed the preset physical range (e.g., pressure exceeding 0~30MPa, temperature exceeding 223.15~353.15K, i.e. -50~80℃; the range is expressed in Kelvin and is consistent with the normalized dimensions) are considered as transmission interference or instrument abnormalities and are eliminated.

[0044] Constant data segment removal: Data segments with no more than a set threshold (e.g., 1, meaning the entire window contains only a single value) indicate that the measuring point is invalid or has not been refreshed for a long time, and are not meaningful for judgment, so they are removed.

[0045] Data volume gating: Count the remaining valid points N. When N reaches the preset data volume threshold (default 10 points; to ensure the significance of hierarchical statistics, it is recommended to fully configure 360 ​​points in the recent window), allow the process to proceed to the next step. Otherwise, return to the data shortage status and wait for the next round of data accumulation without making a risk judgment.

[0046] Finally, we obtained synchronized pressure-temperature sequence pairs (P) with timestamps one-to-one. sync , T sync ), where P sync ={p sync , k (MPa), Tsync = {t sync , k} (K), both share the same timestamp set {τ k This sequence pair will serve as input data for both step S2 (taking the pressure sequence) and step S3 (taking the temperature sequence).

[0047] S2. Based on the natural gas component data and the pressure sequence in the synchronized pressure-temperature sequence pair, the hydrate formation temperature at each time point is calculated using a gas-solid equilibrium thermodynamic model to obtain the hydrate formation temperature sequence.

[0048] The gas-solid equilibrium thermodynamic model used in this step is based on a two-step mechanism of hydrate formation: First, water molecules react with guest molecules through a quasi-chemical reaction to form a basic hydrate with a fixed stoichiometric ratio; second, smaller gas molecules are further adsorbed into the cavities of the basic hydrate, stabilizing the hydrate structure. At phase equilibrium, the distribution of components between the gas phase and the hydrate phase follows a gas-solid equilibrium constant.

[0049] The hydrate formation temperature is defined as the temperature at which, under given pressure P and composition conditions, the system is exactly in the critical state of hydrate formation, denoted as T. h When the actual temperature is below T h At this point, the system enters the hydrate formation region and possesses the thermodynamic conditions for hydrate formation.

[0050] Construct the gas-solid equilibrium constant equation, and for each component i in natural gas, establish the relationship between its gas-solid equilibrium constant and temperature and pressure:

[0051] K i (T, P) = exp(A i + B / T) × P C ;

[0052] Among them, K i Let be the gas-solid equilibrium constant (dimensionless) of the i-th component, T be the thermodynamic temperature (K), P be the pressure (MPa), and A be the solid-gas equilibrium constant (dimensionless). i is the characteristic coefficient of the i-th component, and its value is taken from Table 2 according to the type of hydrate cage structure (structure type I or structure type II) in which the component participates in the formation; B is the temperature sensitivity coefficient shared by all components (Arrhenius-like form, B=-10000.0 in the example), the sign of which makes the equilibrium constant increase when the temperature rises and the system move away from the hydrate region; C is the pressure index shared by all components (C=-0.791 in the example), the sign of which makes the equilibrium constant decrease when the pressure rises and the system more likely to enter the hydrate region.

[0053] Table 2. Example values ​​of characteristic coefficients of gas-solid equilibrium constant (Structure Type I)

[0054] Methane CH4 36.73 Isobutane iC4 31.43 Ethane C2H6 34.83 n-Butane nC4 32.82 Propane C3H8 32.12 Carbon dioxide (CO2) 35.53 Hydrogen sulfide (H2S) 33.73 Nitrogen N2 37.83

[0055] Characteristic coefficient A i The rule for the value of A is: i The smaller the value, the smaller the equilibrium constant of the component, and the easier it is for it to enter the hydrate phase, i.e., hydrate promoters (such as propane and hydrogen sulfide); A i The larger the value, the more likely it is to remain in the gas phase, i.e., a hydrate inhibitor (such as nitrogen). Methane, as the main component of natural gas, is used as the baseline value, and its coefficient can be calibrated by least squares fitting using data on pure methane generation conditions from literature (e.g., 1 MPa corresponds to 273 K, 5 MPa to 282 K, 10 MPa to 287 K, and 20 MPa to 291 K); the coefficients of other components are determined by scaling their relative strength to methane's hydrate promoting / inhibiting effects.

[0056] The above coefficients are approximate engineering values, applicable to operating conditions such as 1~30MPa and -20~+20℃. Specific values ​​can be recalibrated based on the selected thermodynamic fundamental data (such as experimental data or calculation results from commercial software), and this invention does not limit them.

[0057] When the sum of the mole fractions of components other than methane in the normalized component set is less than a set threshold (e.g., 10), -6 When K is in the pure methane condition, the stability condition is determined to be analytically so: Let K CH4 = 1, that is:

[0058] exp(A CH4 + B / T) × P C = 1;

[0059] Taking the natural logarithm of both sides, we get:

[0060] A CH4 + B / T + C × ln P = 0;

[0061] Solving for the temperature yields the hydrate formation temperature at that time point:

[0062] T h = -B / (A CH4 + C × ln P);

[0063] This analytical path requires no iteration, and the calculation time at a single point is on the order of microseconds, making it suitable for online high-frequency computation. As a numerical example: Let A... CH4 = 36.73, B = -10000.0, C = -0.791, when P = 5MPa, ln P = 1.609, substituting these values ​​gives T. h= 10000 / (36.73 -0.791×1.609) = 10000 / 35.46 ≈ 282.0 K (approximately 8.9℃), which is consistent with the literature report of the formation temperature of pure methane at 5MPa, which is 281~283K, verifying the correctness of the coefficient calibration.

[0064] When natural gas is a multi-component mixture, the stability conditions for hydrate formation are:

[0065] Σ i (y i / K i (T h , P)) = 1;

[0066] This involves searching for the temperature at which the sum of the ratios of the mole fractions of each component to their gas-solid equilibrium constant equals 1. To facilitate numerical solutions, a residual function is constructed:

[0067] f(T) = Σ i (y i / K i (T, P)) -1;

[0068] The residual function has a clear physical meaning and monotonicity: a residual greater than zero indicates that the sum of the components is too large at the current temperature, and the system is in the hydrate formation region (temperature is too low); a residual less than zero indicates that the system is outside the formation region (temperature is too high); a residual equal to zero corresponds to the formation temperature.

[0069] Since B < 0, we know that as T increases, exp(A) i + B / T) increase, K i Increase, y i / K i Since the temperature decreases, f(T) decreases monotonically with temperature within the defined interval, and the equation has a unique real root, which guarantees the convergence determinism of numerical root finding.

[0070] The numerical solution employs a hybrid root-finding method (Brent's method) combining inverse quadratic interpolation, the secant method, and the bisection method. This method alternates between inverse quadratic interpolation and the secant method to accelerate convergence, provided the function values ​​at the interval endpoints have opposite signs. When the iteration step size does not meet the safety condition, it degenerates into the bisection method, achieving both superlinear convergence speed and global convergence reliability, without requiring derivative information. The search interval is a preset temperature range [T]. lo , T hi For example, [200, 350] K (corresponding to approximately -73℃ to +77℃, covering all operating conditions of natural gas gathering and transportation); the convergence condition takes the absolute tolerance x. tol = 0.01K and relative tolerance r tol = 10 -6 .

[0071] As a complete example, the following solution is provided for natural gas with a mole fraction of methane 0.962, ethane 0.015, propane 0.005, carbon dioxide 0.008, and nitrogen 0.010, at P = 6.5 MPa:

[0072] Initial interval endpoint residual f(200) = +4.062×10 6 > 0, f(350) = -0.998 < 0, the two have opposite signs to confirm that the root is in the interval; the current interval endpoint, current estimated point, corresponding residual and the sub-method actually used in each subsequent iteration are shown in Table 3.

[0073] Table 3 Iterative trajectory of the Brent method for solving the hydrate formation temperature (P=6.5MPa)

[0074] 0 200.000 350.000 200.000 +4.062e6 -0.998 initial interval 1 350.000 349.995 200.000 -0.998 -0.998 Secant method 2 349.995 274.998 200.000 -0.998 +3.861 Binary search method (step size protection) 3 349.995 334.590 274.998 -0.998 -0.925 Secant method 4 334.590 304.794 274.998 -0.925 -0.861 Binary search method (step size protection) 5 304.794 289.896 274.998 -0.861 -0.250 Binary search method (step size protection) 6 289.896 284.748 274.998 -0.250 +0.399 inverse quadratic interpolation 7 289.896 287.914 284.748 -0.250 -0.049 Secant method 8 287.914 287.484 284.748 -0.049 +0.0018 inverse quadratic interpolation 9 287.484 287.499 287.914 +0.0018 -4.7e-5 Secant method 10 287.499 287.494 287.484 -4.7e-5 +5.6e-4 Secant method 11 287.494 287.499 287.494 +5.6e-4 -4.7e-5 Satisfy tolerance, convergence

[0075] As shown in Table 3, the first two iterations resulted in excessively small secant step sizes due to the massive residuals at the low-temperature end (on the order of 10^6), triggering step size protection and reverting to the bisection method to rapidly narrow the interval. From the sixth iteration onwards, inverse quadratic interpolation became effective, and from the eighth to the eleventh iteration, the residuals rapidly decreased from the order of 10^-2 to the order of 10^-5, exhibiting superlinear convergence. Finally, after 11 iterations, it converged to T. h =287.4990 K (14.35℃), at this time f(T) h ) = -4.7×10 -5 , Σ i y i / K i = The deviation between 0.999953 and 1 is less than one ten-thousandth, which verifies the accuracy of the root finding.

[0076] If the residuals at both ends of the interval have the same sign, it indicates that the generated temperature is outside the preset interval. In this case, a conservative approach should be taken: when f(T) lo ) > 0 and f(T) hi When ) > 0, it indicates that the temperature is still within the generation zone up to the upper limit of the interval, but the actual generation temperature is higher than the upper limit, so we return to the upper limit T. hi As the formation temperature (at which the temperature safety margin is extremely negative, conservatively judged as high risk, and will not be missed); when both residuals are negative, it indicates that hydrates will not form even at the lower bound temperature, and the temperature returns to the lower bound T. lo (At this point, the safety margin is extremely positive, and the algorithm is deemed safe.) This rule ensures that the algorithm has deterministic output under any input conditions and will not interrupt online service due to abnormal values.

[0077] For P syncThe pressure at each time point is selected for root finding according to the component conditions. For time points where root finding fails (e.g., component abnormality or pressure exceeding the applicable range), the missing value is returned and removed from the sequence. This is necessary because the result of comparing the missing value with any threshold is always false. If it is not removed, it will be misjudged as normal, resulting in missed reports.

[0078] After processing, a hydrate formation temperature sequence that corresponds one-to-one with the synchronous pressure sequence in time was obtained:

[0079] T h = {T h ,k}(K), where k is consistent with the set of timestamps of the synchronization sequence.

[0080] As another verification example: Solving for the same composition (methane 0.962, ethane 0.015, propane 0.005, carbon dioxide 0.008, nitrogen 0.010) at P = 5 MPa, we obtain T. h = 285.7955K (12.65℃), which is about 3.8K higher than the 282.0K of pure methane under the same pressure. This is consistent with the physical law that "heavy components increase the hydrate formation temperature", which once again verifies the rationality of the model coefficients.

[0081] To facilitate review and understanding of the roles of each component under stable conditions, Table 4 further provides the values ​​for P = 6.5 MPa and T. h = Mole fraction of each component at 287.4990 K, gas-solid equilibrium constant and its summation term Σ with respect to stability conditions i y i / K i . contributions.

[0082] Table 4. Contributions of each component to the summation term of the stability condition (P = 6.5 MPa, T) h = 287.4990K)

[0083] Methane CH4 0.962 36.73 1.5947 0.603258 60.3% Ethane C2H6 0.015 34.83 0.23851 0.062890 6.3% Propane C3H8 0.005 32.12 0.015870 0.315062 31.5% Carbon dioxide (CO2) 0.008 35.53 0.48031 0.016656 1.7% Nitrogen N2 0.010 37.83 4.7907 0.002087 0.2% sum 1.000 — — 0.999953 100%

[0084] As shown in Table 4, propane, with a molar fraction of only 0.5%, contributes as much as 31.5%, exceeding the sum of all other components except methane. This confirms the mechanism by which propane acts as a strong hydrate promoter—which is the physical basis for the requirement that natural gas containing heavy components must retain a multi-component root-finding path in addition to the pure methane analysis path. The sum of the contributions of all components is 0.999953, which deviates from the theoretical value of 1 under stability conditions by less than one ten-thousandth.

[0085] The thermodynamic model is based on idealized assumptions and engineering approximations, and does not include the influence of impurities, flow rates, and heat transfer conditions of specific gas sources. Therefore, there may be systematic deviations between the calculated values ​​and the actual measured conditions. To absorb this deviation, in one optional implementation, a correction network constructed using a multilayer sensing mechanism is set up to correct the residuals of the mechanistic calculation results.

[0086] For each time point, construct a 6-dimensional feature vector:

[0087] x k = [p k , t k , y CH4 , y C2 , ΔT k , q k ];

[0088] Where, p k The pressure (MPa) at that time point, t k The actual temperature (K) at that time point, y CH4 y C2 The mole fractions of methane and ethane are respectively (ethane is used as a heavy component to characterize the influence of other components), ΔT k = t k - T h,k This is the initial value of the temperature safety margin at that time point (its introduction allows the network to learn the "margin-deviation" relationship rather than simply memorizing temperature and pressure), q k To normalize the flow rate, calculate using the following formula:

[0089] q k = min(max(Q k / Q max , 0), 1), Q max = 10000 m³ / d;

[0090] Among them, Q k For traffic at that point in time, q is used when traffic data is missing. k Take 0.

[0091] Because the dimensions of the original features mentioned above differ significantly (pressure is approximately 0-30 MPa, temperature is approximately 200-350 K, and mole fraction and normalized flow rate are 0-1), directly inputting them into a network using hyperbolic tangent as the activation function would cause the activation values ​​to saturate around ±1 due to excessively large input amplitudes in the intermediate layers, resulting in vanishing gradients and difficulty in training convergence. Therefore, before inputting the features into the network, each feature is standardized using the following formula:

[0092] x̂′ k,d = (x k,d - μ d ) / σd , d = 1, 2, …, 6;

[0093] Where, x k,d Let μ be the d-th dimension of the original feature of the k-th sample. d With σ d These are the mean and standard deviation of this feature on the training set, respectively (when σ d = 0 is taken as 1 to avoid division by zero), and both are saved together with the correction network; the same μ is used in the online inference stage as in the training stage. d σ d This ensures the consistency of feature transformation.

[0094] After standardization, each feature dimension has approximately zero mean and unit variance, falling within the input-sensitive region of the hyperbolic tangent function. In subsequent steps, the standardized feature vector will be denoted as x̂′. k .

[0095] The network is modified to a multilayer perceptron, with the following example structure: Input layer (6-dimensional) — Fully connected layer (32 neurons, hyperbolic tangent activation) — Fully connected layer (16 neurons, hyperbolic tangent activation) — Output layer (1 neuron, linear output z). The layer-by-layer forward propagation process is as follows:

[0096] First hidden layer:

[0097] h (1) = tanh(W (1) x̂ + b (1) );

[0098] Among them, W (1) ∈ R 32×6 The first layer weight matrix is ​​32 rows and 6 columns, with a total of 192 weights, where x̂ is the standardized eigenvector and b... (1) ∈ R 32 h is the first layer bias vector (32 biases). (1) ∈ R 32 This is the output of the first hidden layer.

[0099] Second hidden layer:

[0100] h (2) = tanh(W (2) h (1) + b (2) );

[0101] Among them, W (2) ∈ R 16×32 (512 weights), b (2) ∈ R 16 (16 biases), h (2)∈ R 16 .

[0102] Output layer:

[0103] z = w (3)T h (2) + b (3) ;

[0104] Among them, w (3) ∈ R 16 (16 weights), b (3) This is a scalar bias.

[0105] The total number of parameters to be trained in the entire network is 192 + 32 + 512 + 16 + 16 + 1 = 769. The hyperbolic tangent activation function is:

[0106] tanh(z) = (e z - e -z ) / (e z + e -z );

[0107] The output layer applies a constraint transformation to the linear output to obtain the correction amount:

[0108] r = 3 × tanh(z);

[0109] Since tanh(z) ∈ (-1, 1), the correction amount r is naturally constrained to the range of -3℃ to +3℃. This structural constraint ensures that the data-driven correction amount does not deviate too far from the mechanistic value: even if the network extrapolates from outside the training distribution, its output is limited to a safe bound, preventing the neural network correction from undermining the reliability of the mechanistic model.

[0110] The corrected generation temperature is:

[0111] T′ h,k = T h,k + r k。

[0112] Based on the field-measured hydrate formation condition samples (obtained from field sampling experiments or operational condition backtracking annotations), residual labels are constructed:

[0113] r k label = T h,k obs - T h,k mech ;

[0114] Among them, T h,k obs To measure the actual generation temperature, T h,kmech These are values ​​calculated using the thermodynamic model; samples with labels exceeding [-3, +3]℃ are truncated to the boundary (consistent with the network output domain).

[0115] The training objective is to minimize the mean squared error loss.

[0116] L = (1 / N) × Σ k (r k label -r k pred ) 2 ;

[0117] The training uses an adaptive moment estimation optimizer (Adam), whose parameter update rule is as follows:

[0118] m t = β1·m t-1 + (1-β1)·g t ;

[0119] v t = β2·v t-1 + (1-β2)·g t 2 ;

[0120] m̂′ t = m t / (1 - β1 t ), v̂′ t = v t / (1 - β2 t );

[0121] θ t = θ t-1 - η·m̂′ t / (√v̂′ t + ε);

[0122] Among them, g t Let m be the gradient of the current mini-batch. t v t For the exponential moving average of the first and second moments, m̂ t v̂ t It is the estimate after bias correction (t is the number of iterations, and the correction term is 1 / (1-β)). t (Compensation for the systematic underestimation caused by zero initialization in the early stages of training) β1 = 0.9, β2 = 0.999, ε = 10 -8 The learning rate η = 0.001.

[0123] The gradient of the loss function with respect to the linear output of the output layer is explicitly given by the chain rule:

[0124] ∂L / ∂z k = -(2 / N)·(r k label -r k pred )·3·(1 - tanh 2 (z k ));

[0125] Among them, tanh 2 (z k The term comes from the derivative of the constraint transformation r = 3·tanh(z), dr / dz = 3·(1 - tanh). 2 (z)); The gradients of the weights and biases of each hidden layer are calculated layer by layer by backpropagation, that is, they are obtained by multiplying the error term of the next layer with the input of the current layer.

[0126] The training samples are randomly divided into a training set and a validation set in an 8:2 ratio (with a fixed random seed, e.g., 42). The validation set is used only for early stopping detection and learning rate decay monitoring, and does not participate in parameter updates. To enhance training stability, the following measures are also adopted:

[0127] Weight initialization: The weights of each fully connected layer are uniformly distributed U(-√(6 / (n)). in +n out )), √(6 / (n in +n out Random initialization (n) in n out (These are the input and output dimensions of this layer, respectively). The bias is initialized to 0 so that the variance of the forward output of each layer remains approximately unchanged, avoiding saturation or decay upon initial activation.

[0128] Gradient clipping: When the gradient norm exceeds the set value, it is scaled by g ← g·min(1, γ / ‖g‖), where γ = 1.0, to prevent gradient explosion caused by individual abnormal samples;

[0129] Learning rate adaptive decay: If the validation loss does not decrease for 5 consecutive rounds, the learning rate is multiplied by 0.5, with a lower decay limit of 10. -6 ;

[0130] Early stopping and restoration of optimal weights: Training is terminated when the validation loss has not been updated to the optimal value for 10 consecutive rounds, and the historical optimal weights are restored to prevent overfitting;

[0131] Mini-batch training: Batch size 256, maximum number of training rounds 500, and training is reproducible with a fixed random seed.

[0132] The corrected hydrate formation temperature sequence T′ was obtained after processing. h , replacing T h This step participates in the calculation of step S3. In implementations where training labels are not obtained, this step can be disabled, and the mechanistic calculation value can be used directly. In this case, the method can still run completely (the correction amount is always 0), which demonstrates the good deployability and scalability of this solution.

[0133] S3. Based on the temperature sequence in the hydrate formation temperature sequence and the synchronous pressure-temperature sequence pair, calculate the temperature safety margin at each time point to obtain a temperature safety margin sequence; compare each temperature safety margin in the temperature safety margin sequence with a preset classification threshold to obtain a risk classification result, and map the risk classification result to a risk score.

[0134] Read the hydrate formation temperature sequence and the synchronization temperature sequence, and calculate for each time point k:

[0135] M k = t sync,k -T h,k ;

[0136] Since both are on the Kelvin temperature scale, the difference is equal to the difference in degrees Celsius. k The unit is denoted by °C. Temperature safety margin (the opposite of supercooling) characterizes the degree to which the actual operating temperature exceeds the hydrate formation temperature: M k The larger the value, the farther away from the generation conditions and the safer it is; M k ≤ 0 indicates that the actual temperature has fallen below the formation temperature, entering the hydrate formation region. The temperature safety margin sequence M = {M} is obtained by iterating through each time point. k}

[0137] Preset grading thresholds include warning threshold M warn (Default 5℃) and critical threshold M crit (Default 0℃), and should satisfy M crit ≤ M warn When configuration is reversed (M) warn ≤ M crit When the temperature reaches 0℃, it automatically reverts to the safe default value (5℃ / 0℃) to prevent false negatives caused by configuration errors. The threshold can be adjusted according to the season and operating conditions; for example, in winter when the water content is high, it can be set to M. warn = 5℃, M crit = -10℃.

[0138] Temperature safety margin M at the end of the recent window (current time) cur Perform the following judgment:

[0139] When M cur > Mwarn When the condition is far from the conditions for hydrate formation, it is considered normal: the operating condition is far from the conditions for hydrate formation and no intervention is required.

[0140] When M crit < M cur ≤ M warn When the actual temperature is close to the production temperature, it is determined to be a warning: the actual temperature is close to the production temperature, and it is necessary to intensify monitoring and prepare for alcohol injection and heating.

[0141] When M cur ≤ M crit When the temperature reaches a critical level, it is determined that the actual temperature has entered the hydrate formation zone, and immediate measures such as alcohol injection, heating, or adjustment of operating conditions are required.

[0142] To provide a quantitative benchmark for alarm prioritization and response assessment, the current temperature safety margin is mapped to a risk score ranging from 0 to 1 using a piecewise linear function:

[0143] s = 1, when M cur ≤ M crit ;

[0144] s = (M warn - M cur ) / (M warn - M crit When M crit < M cur < M warn ;

[0145] s = 0, when M cur ≥ M warn ;

[0146] The closer the score is to 1, the higher the risk. When the risk score reaches the preset concern score (e.g., 0.3), it is initially judged as abnormal and enters the blocking and confirmation stage in step S4.

[0147] The margin at the current moment only reflects the instantaneous state. If a depth exceedance occurred within the window and the current situation happens to recover, the instantaneous score will underestimate the persistent risk. Therefore, we introduce the recent window minimum margin M. min = min k M k Bonus rules:

[0148] If M min ≤ M crit If , then s = max(s, 0.9);

[0149] If M min ≤ M warnThen s = max(s, 0.7 × (M) warn - M min ) / (M warn - M crit ));

[0150] Finally, the score is truncated to the interval [0, 1], and the risk classification result (normal / warning / critical) and risk score s are obtained after processing. Both are used as the processing objects in step S4. At the same time, the window statistics mean margin and minimum margin are output for reference in alarm details.

[0151] S40. Identify the drainage period and mask the data points within the drainage period to obtain the masked risk classification result; perform continuous boundary verification on the masked risk classification result, and generate and output the final alarm level after verification.

[0152] During separator draining (venting and draining) operations, a rapid rise in liquid level and a rapid drop in temperature occur simultaneously. This temperature drop stems from the throttling and evaporative cooling effects of the draining operation, rather than from hydrate formation. Without shielding, this will generate regular false alarms. The required liquid level timing data for identification comes from the synchronized liquid level sequence L obtained in step S1 and preprocessed in the same manner. sync = {l k (This step is skipped if no liquid level is configured, and the masking process will not be performed).

[0153] To ensure that the judgment thresholds have a unified semantics across different sampling periods, the rate of change is first normalized based on a ten-minute baseline:

[0154] Δt̄′ = (τ N -τ1) / (N - 1), k 10 = Δt̄′ / 600;

[0155] Among them, τ1, τ N The recent window's first and last times, N is the number of points within the window, and k is the number of points within the window. 10 The normalization factor is the ten-minute normalization factor (e.g., k when sampling for 20 seconds). 10 = 1 / 30). Based on this, calculate the pointwise rate of change:

[0156] L′ k = (l k - l k-1 ) / k 10 , T′ k = (t k -t k-1 ) / k 10 ;

[0157] The criteria for determining the drainage period are that both conditions must be met simultaneously:

[0158] L′ k > λ L And T′ k < -λ T , λ L = 5.0 (% range / 10min), λ T = 2.0 (℃ / 10min);

[0159] That is, the liquid level rises at a rate exceeding 5% of the range every ten minutes, while the temperature drops at a rate exceeding 2°C every ten minutes.

[0160] Considering that single-point communication spikes in the liquid level channel may momentarily meet the above conditions, to avoid missed reporting of actual ice blockage periods due to false blocking, in a preferred embodiment, the corresponding time point is marked as a drainage period only when m consecutive sampling points (default m = 2) simultaneously meet the above conditions; if only a single point meets the conditions but subsequent sampling points do not, it is judged as a data spike and is not blocked. Data points within the drainage period are not included in alarm determination; when all data points in the recent window are in the drainage period, the system directly returns to the normal level and notes "in the drainage stage, ice blockage risk is not assessed at this time"; when liquid level data is missing, blocking is skipped and the system is treated as if there is no drainage.

[0161] Instantaneous data jitter may cause single-point exceedances. In order to suppress such noise without introducing smoothing filters (which would delay the exposure of real risks), continuous exceedance statistics are performed on the tail of recent windows.

[0162] First, construct the out-of-bounds marker sequence:

[0163] e k = 1, if M k ≤ M warn And k is not in the drainage period; otherwise e k = 0;

[0164] Then, scanning backwards from the end of the window, the number of out-of-bounds points with consecutive 1s at the end is counted:

[0165] C = max{ c : e N-j = 1, j = 0, 1, …, c-1};

[0166] When C is less than the preset confirmation threshold N min (Default 2, meaning at least two consecutive points exceed the limit) When the result is downgraded to non-abnormal, the score is reduced to no more than half of the preset attention score (e.g., s ≤ 0.15); when C ≥ N min The alarm is confirmed to be active in time.

[0167] After confirmation, the final alarm level (normal / warning / critical) and alarm details are generated and output. The risk score, after processing and correction, is output along with the final alarm level as one of the fields in the alarm details (the "score" field in Table 5), allowing the upper-level alarm platform to sort alarms from multiple locations and determine their priority based on the score. The names, meanings, value formats, and examples of each field in the alarm details are shown in Table 5. The above results are also output in structured data format for consumption by the upper-level alarm platform and the work order system.

[0168] Table 5. Field Definitions and Examples for Alarm Output

[0169] alarm_level Final Alarm Level Enumeration of the three states: normal, warning, and critical. critical margin_current Current temperature safety margin Real number, unit: °C -9.35 hydrate_temp hydrate formation temperature Real number, unit °C (converted from K). 14.35 actual_temp Current actual temperature Real number, unit: °C 5.0 pressure Current operating pressure Real number, unit MPa 6.5 margin_mean / margin_min Recent window mean margin / minimum margin Real number, unit: °C -4.2 / -9.8 score Risk Score real numbers from 0 to 1 1.0 violation_start The moment the boundary was crossed began Timestamp of the start of the continuous boundary crossing segment at the tail 2026-01-15 08:23:40 duration Duration of crossing the boundary The difference between the current time and the time when the boundary violation began, in minutes. 36 confirm_count Tail consecutive cross-boundary points positive integer 108 drainage_masked Number of drainage shielding points in this window non-negative integers 0 advice Recommendation document for handling Suggested string matching the level Immediately add methanol / ethylene glycol inhibitor

[0170] The start time of the boundary violation is determined by the start time of the processed tail segment, and the duration of the boundary violation is the difference between the current time and the start time. The mapping relationship between the final alarm level and the handling recommendation is shown in Table 6.

[0171] Table 6. Risk Level and Response Measures Mapping

[0172] normal The margin is greater than the warning threshold and the confirmation fails. Maintain regular monitoring cycle No intervention required Warning The margin falls within (critical threshold, warning threshold) and is confirmed to pass. Increase monitoring frequency; verify inhibitor reserves and heating equipment availability in advance; prepare for alcohol injection. Handled within 24 hours critical The margin is not greater than the critical threshold and the confirmation is passed. Immediately add methanol / ethylene glycol inhibitors, or adjust the throttling pressure drop and increase the temperature after throttling. Immediate action

[0173] For multiple monitoring points within the same pipeline or station, steps S1 to S4 are executed separately, sorted by severity: normal < warning < critical. The highest level among all points is taken as the overall alarm conclusion, and the margin of each point is listed in the alarm details to locate high-risk nodes. After the alarm triggers the response measures, the changes in margin are continuously tracked: the effectiveness of the measures is quantitatively evaluated by the difference in margin between the time the measures were initiated and the current time (the recovery rate). When the recovery rate exceeds a set value (e.g., 3°C), the effectiveness of the measures is automatically confirmed. Normal monitoring is resumed after the margin recovers to a level greater than the warning threshold, forming a closed loop of "early warning - response - assessment".

[0174] Thus, the complete data chain for a single calculation is as follows: The raw time-series data from the production data acquisition and monitoring system and the analytical component data are processed in step S1 to obtain synchronized pressure-temperature sequences and normalized component sets; step S2 uses the pressure sequence and component set to calculate the hydrate formation temperature sequence; an optional step calculates corrections based on these sequences along with temperature, pressure, component, and flow characteristics to obtain the corrected formation temperature; step S3 uses the formed temperature sequence and the temperature sequence to calculate the temperature safety margin sequence and complete grading and scoring; step S4 uses the temperature safety margin sequence and the synchronized liquid level sequence to complete drainage shielding and continuous boundary crossing confirmation, and outputs the final alarm. A summary of the input sources, processing operations, and output data and their destinations for each step of the above data chain is shown in Table 7. As can be seen from Table 7, each step has a source for inputs and a destination for outputs, and there is no isolated data in the entire text.

[0175] Table 7. Summary of Input-Processing-Output Data Links for Each Step

[0176] S1 Data Acquisition and Preprocessing Pressure / temperature time series data (production data acquisition and monitoring system), component data (laboratory report), liquid level time series data (optional) Component analysis normalization, unit normalization, timestamp alignment, linear interpolation, four-category elimination, and data volume gating Synchronized pressure-temperature sequence pairs (P_sync, T_sync), normalized component set G, and synchronized liquid level sequence L_sync (for use by S20 / S30 / S40). S2 generation temperature calculation P_sync and G (from S10) Construct the gas-solid equilibrium constant equation; for pure methane, use the analytical approach; for multi-component methane, use the Brent method for numerical root determination. Hydrate formation temperature sequence T_h (for use by S30 / S25) Correction (optional) T_h, P_sync, T_sync, G, Traffic data (from S10 / S20) Constructing 6-dimensional feature vectors; forward propagation using a multilayer perceptron; ±3℃ constraint on output correction. The corrected temperature sequence T′_h is generated (replacing T_h for use by S30). S3 Margin Calculation and Classification T_h (or T′_h) and T_sync (from S10 / S20 / S25) Point-by-point calculation of temperature safety margin; threshold verification; three-level classification; piecewise linear scoring; minimum margin addition. Risk classification results, risk score s, and window statistics (for use by S40). S4 shielding and confirmation Margin sequences M and L_sync (from S30 / S10) Drainage period identification and shielding; boundary crossing marker construction; statistics and confirmation of continuous tail boundary crossings; scoring correction. Final alarm level, risk score, and alarm details (output to the alarm platform and ticket system)

[0177] In this implementation, all parameters can be adjusted through the configuration file. Table 8 provides the meaning and example values ​​of each parameter.

[0178] Table 8. Examples of Main Parameter Configurations

[0179] Recent window points Number of data points used in the classification and confirmation process (approximately 2 hours based on 20-second sampling). 360 points Historical cache points Maximum amount of historical data to be retrieved in a single round of calculation 10080 points Data volume threshold If the number of valid points is lower than this value, the data is considered insufficient. 10 o'clock Upper limit of consecutive missing numbers No imputation is performed if the number of consecutive missing points exceeds this limit. 2 o'clock Timestamp alignment tolerance The time tolerance for nearest neighbor pairing based on temperature is half a sampling period. 5 s (10-second sampling) Constant determination threshold The window is considered a constant segment if the number of different values ​​within it does not exceed this value. 1 Temperature sensitivity coefficient B Common temperature coefficient of gas-solid equilibrium constant equation -10000.0 Stress Index C Common pressure coefficient of gas-solid equilibrium constant equation -0.791 Root search interval Temperature search range for numerical root finding (Kelvin) [200, 350] Root absolute tolerance Absolute temperature tolerance of iterative convergence 0.01 K Root relative tolerance Relative tolerance of iterative convergence 1e-6 Warning threshold Upper limit of temperature safety margin warning judgment 5 ℃ Critical threshold Upper limit of temperature safety margin for severe judgment 0 ℃ Follow rating This score indicates an initial abnormality. 0.3 Minimum margin addition lower limit Lower bound of scoring when the minimum margin of the window is exceeded 0.9 Confirmation threshold Minimum number of points required for consecutive tail crossings to be confirmed 2 o'clock Liquid level change threshold Threshold for the rate of liquid level rise in determining drainage 5% range / 10min Temperature drop threshold Temperature drop rate threshold for drainage determination 2 ℃ / 10min Drainage confirmation of consecutive points The number of sampling points required to continuously meet the conditions during the drainage period. 2 o'clock Effective recovery rate of disposal Confirmation of the effectiveness of the response measures and the resulting margin recovery threshold 3 ℃ Hidden layer structure Correct the configuration of neurons in the hidden layer of the network. 32、16 Correction quantity constraint boundary Correct the upper limit of the absolute value of the network output. 3 ℃ Training learning rate Correcting the initial learning rate during network training 0.001 Training batch size / maximum number of rounds Adjust the mini-batch and epoch limits for network training. 256 / 500

[0180] The above values ​​are illustrative examples. Those skilled in the art can adjust the above parameters based on the characteristics of the gas source components, sampling cycle, process conditions, and risk preferences, without departing from the concept of this invention.

[0181] Example 1

[0182] Throttling valve Joule-Thomson effect ice blockage early warning. Under winter operating conditions, after natural gas from a gas well is throttled by a throttle valve, the pressure is 6.5 MPa and the temperature is 5℃. The molar fraction of the natural gas composition is methane 0.962, ethane 0.015, propane 0.005, carbon dioxide 0.008, and nitrogen 0.010.

[0183] After executing step S1, the synchronization sequence pairs and the normalized component set are obtained;

[0184] The multi-component root-finding path of step S2 (the complete iterative trajectory is shown in Table 3) was executed. After 11 iterations, the hydrate formation temperature T_h under this pressure was found to be 287.4990K, which is about 14.35℃ (the contribution of each component is shown in Table 4).

[0185] Step S3 calculates the temperature safety margin M = 5 - 14.35 = -9.35℃, which is not greater than the critical threshold of 0℃. Therefore, it is judged as serious and the score is 1.

[0186] Step S4 counts the number of consecutive out-of-bounds points at the tail end, confirming the alarm is established. According to Table 6, a severe alarm and handling suggestion are output: "Immediately add methanol / ethylene glycol inhibitor, or adjust the throttling pressure drop." After adding methanol as suggested on-site, the temperature rises after throttling, and the margin increases to above the warning threshold, thus clearing the alarm.

[0187] Example 2

[0188] Rapid calculation under pure methane operating conditions. The composition of a certain treated dry gas is close to pure methane (the sum of non-methane components is less than 10⁻⁶), and the pressure is 4.0 MPa.

[0189] Step S2 determines to use the analytical fast path: T_h = 10000 / (36.73 - 0.791×ln4.0) = 10000 / (36.73 - 1.0966) ≈ 280.64K (approximately 7.5℃). The actual temperature is 12℃, with a margin of approximately 4.5℃, falling between the critical threshold of 0℃ and the warning threshold of 5℃. Therefore, a warning is issued, and it is recommended to increase monitoring and prepare inhibitors. This path's single-point calculation time is in the microsecond range, significantly reducing online computational overhead compared to the numerical root-finding path (approximately 11 iterations). For stations operating under pure methane conditions at all points, the overall calculation time is reduced by approximately an order of magnitude.

[0190] Example 3

[0191] Multi-point collaborative monitoring of pipelines in winter. In a certain winter operating condition with an ambient temperature of -15℃, steps S1 to S4 are executed for 5 monitoring points of the same pipeline, and warning thresholds of 5℃ and critical thresholds of -10℃ are configured according to the winter operating conditions.

[0192] Results: The margin at the point downstream of the throttle valve was -12.5℃, not exceeding the critical threshold, and was judged as severe; the margin at a certain low point was -8.2℃, between the two thresholds, and was judged as a warning; the margins at the remaining points were greater than the warning threshold and were judged as normal. Based on the multi-point coordination rule, the highest severity was taken as the overall conclusion, identifying the point downstream of the throttle valve as a high-risk node and prioritizing alcohol injection, achieving precise allocation of disposal resources; the margins at each point were output in the same table for easy comparison and analysis by operators.

[0193] Example 4

[0194] Discharge period shielding. A separator performs a periodic discharge operation. During the discharge period, the liquid level rises at a rate of approximately 8% of the capacity every 10 minutes, while the temperature drops at a rate of approximately 3°C every 10 minutes. If not shielded, the rapid temperature drop during the discharge period will cause the margin to momentarily fall into the over-limit zone and trigger a false alarm.

[0195] After executing step S4, if the liquid level change rate of 8% of the range / 10min is greater than the threshold of 5% of the range / 10min, and the temperature change rate of -3℃ / 10min is less than -2℃ / 10min, the corresponding time points are marked as drainage periods and masked. For the remaining data points, continuous out-of-bounds errors at the tail end are counted, but the confirmation threshold is not reached, and no alarm is generated. After drainage is completed and the liquid level and temperature return to normal, the system automatically removes the mask. Statistics show that after enabling the drainage mask, false alarms of ice blockage caused by drainage operations at this station were completely eliminated.

[0196] Example 5

[0197] Correction for network bias absorption. A gas source contains small amounts of nitrogen and carbon dioxide. The thermodynamic model calculations show a systematic deviation of approximately +1.8℃ from the measured generation conditions backtracked from the field (the model is too high).

[0198] Residual labels were constructed from collected field samples, and a correction network was trained (32 and 16 hidden layers, learning rate 0.001, batch size 256, and 10 early stop rounds). After training, the mean square error of the validation set converged to below 0.4. During online operation, the network output a correction of approximately -1.6℃ for the current operating condition (within the ±3℃ constraint). After correction, the deviation between the generated temperature and the measured conditions decreased to within 0.5℃, and the classification judgment and scoring were subsequently corrected, significantly improving the accuracy of the early warning.

[0199] Example 6

[0200] Quantitative evaluation of the effectiveness of alcohol injection measures. After a critical alarm is triggered at a certain location, alcohol injection is initiated at a rate of 50 L / h. At the time of the alarm, the margin was -6.5°C; after 30 minutes, the margin recovered to -2.3°C, an increase of 4.2°C. The system automatically confirms the effectiveness of the measures according to the treatment evaluation logic in step S4, and continues to track the margin until it recovers to above the warning threshold, after which normal monitoring resumes. All data from the entire process is accumulated into a treatment knowledge base, providing quantitative reference for the selection of measures for similar alarms in the future.

[0201] Test case

[0202] To verify the effectiveness of this invention, a test set covering various gas components (methane molar fraction 85%~99%, containing ethane, propane, carbon dioxide, and nitrogen), pressures of 1~30 MPa, and temperatures of -20~+30℃ was constructed. Ablation tests, comparative tests, threshold sensitivity tests, and computation time tests were conducted. Evaluation indicators included: mean absolute error (MAE) and root mean square error (RMSE, in °C) of hydrate formation temperature, accuracy of risk classification (ACC), false alarm rate (FAR), and false negative rate (MDR).

[0203] Ablation test: Based on the complete method, individual steps were removed separately, and the contribution of each step was examined. The results are shown in Table 9.

[0204] Table 9 Ablation Test Results

[0205] Remove correction network 2.45 3.12 0.823 0.112 0.089 Remove drainage shield 1.02 1.35 0.801 0.178 0.041 Remove consecutive out-of-bounds confirmations 1.01 1.33 0.795 0.163 0.039 Remove minimum margin bonus 1.05 1.41 0.856 0.089 0.078 Complete method of the present invention 0.89 1.23 0.934 0.045 0.038

[0206] As shown in Table 9, the generated temperature error increased significantly after removing the correction network, indicating that the correction network contributed the most to absorbing systematic biases. After removing the drain shield or continuous boundary crossing confirmation, the false alarm rate increased from 0.045 to 0.178 and 0.163, respectively, indicating that the two anti-false alarm mechanisms were significant and independent of each other. After removing the minimum margin addition, the false alarm rate increased from 0.038 to 0.078, indicating that the mechanism effectively compensated for the underestimation of continuous risk by the instantaneous scoring.

[0207] Comparative experiment: The results are shown in Table 10, compared with existing representative methods.

[0208] Table 10 Comparison results with existing methods

[0209] Experienced diagram method 4.56 5.89 0.712 0.198 0.123 Empirical correlation method 3.78 4.67 0.756 0.156 0.098 Thermodynamic model (uncorrected) 2.45 3.12 0.823 0.112 0.089 Complete method of the present invention 0.89 1.23 0.934 0.045 0.038

[0210] As shown in Table 10, the average absolute error of the generation temperature of the complete method of the present invention decreased from 4.56℃ in the empirical chart method and 2.45℃ in the uncorrected thermodynamic model to 0.89℃, the classification accuracy increased from 0.823 to 0.934, and the false alarm rate decreased from 0.112 to 0.045. All indicators are optimal.

[0211] Threshold sensitivity test: The combination of warning threshold and critical threshold was adjusted on the winter operating condition test subset, and the results are shown in Table 11.

[0212] Table 11 Threshold combination sensitivity test (winter operating condition subset)

[0213] 5 0 0.908 0.052 0.061 5 -10 0.941 0.043 0.035 7 -10 0.928 0.031 0.058 3 -10 0.922 0.087 0.026

[0214] As shown in Table 11, the combination of warning threshold of 5℃ and critical threshold of -10℃ performs best under winter conditions. Although tightening the warning threshold reduces the false alarm rate, it significantly increases the false alarm rate, which verifies the rationality of the default threshold combination and also shows that the threshold mechanism of the present invention can be flexibly configured according to the season.

[0215] Computation time test: For a single round of full calculation of 10080 historical data points, the average single-point calculation time for the pure methane analytical path and the multi-component numerical root-finding path is 3 microseconds and 41 microseconds respectively (including multiple calculations of the residual function). The total calculation time for a single round is less than 0.5 seconds, which meets the real-time requirement of a 5-minute online calculation cycle. The inference time of the correction network is negligible.

[0216] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for early warning of natural gas hydrate ice blockage risk based on a gas-solid equilibrium thermodynamic model, characterized in that, Includes the following steps: Step S1: Obtain natural gas component data, pressure time series data after throttling, and temperature time series data. Perform unit normalization, timestamp alignment, missing value imputation, and outlier removal on the pressure time series data and temperature time series data to obtain synchronized pressure-temperature sequence pairs. Step S2: Based on the pressure sequence in the natural gas composition data and pressure-temperature sequence pair, the hydrate formation temperature at each time point is calculated using a gas-solid equilibrium thermodynamic model to obtain the hydrate formation temperature sequence. Step S3: Based on the temperature sequence in the hydrate formation temperature sequence and the pressure-temperature sequence pair, calculate the temperature safety margin at each time point to obtain the temperature safety margin sequence. Each temperature safety margin in the temperature safety margin sequence is compared with a preset classification threshold to obtain the risk classification result, and the risk classification result is mapped to a risk score. Step S4: Identify the drainage period in the risk classification results and mask the data points within the drainage period to obtain the masked risk classification results; Continuous boundary verification is performed on the risk classification results after shielding. Once the verification is successful, the final alarm level and risk score are generated and output.

2. The method according to claim 1, characterized in that, In step S1, natural gas component data is acquired, including: receiving natural gas component strings recorded in the format "component name: mole percentage", and parsing them to obtain the mole percentage v of each component. i , and press y i = v i / Σ j v j Normalized to the mole fraction of each component, where i and j traverse the known components involved in the calculation.

3. The method according to claim 2, characterized in that, In step S1, a synchronized pressure-temperature sequence pair is obtained, specifically including: The units of temperature time series data are uniformly converted to Kelvin, and the units of pressure time series data are uniformly converted to MPa. Based on the common timestamp of pressure time series data and temperature time series data, pressure values ​​and temperature values ​​at the same time are paired. For duplicate data at the same timestamp, the last one is retained to obtain coarse alignment sequence pairs. Missing values ​​in coarsely aligned sequence pairs are filled with linear interpolation, and when the number of consecutive missing points exceeds the preset missing limit, the missing points are retained and marked; data points with non-positive pressure values, outlier data points that exceed the preset physical range, and constant data points whose values ​​do not change within the window are removed. After the above processing and the number of valid points reaches the preset data volume threshold, pressure-temperature sequence pairs corresponding to timestamps are obtained.

4. The method according to claim 3, characterized in that, Step S2 specifically includes: Construct the gas-solid equilibrium constant equation: K i (T , P) = exp(A i + B / T) × P C ; Among them, K i Let A be the gas-solid equilibrium constant of the i-th component, T be the thermodynamic temperature in Kelvin, and P be the pressure in megapascals; i B is the characteristic coefficient determined by the i-th component based on the hydrate cage structure type, C is the temperature sensitivity coefficient shared by all components, and C is the pressure index shared by all components. When the natural gas composition data indicates that the natural gas is pure methane, substitute the coefficient corresponding to methane into the analytical formula T. h =-B / (A CH4 The hydrate formation temperature T at that time point can be directly obtained by solving for C × ln P. h ; When natural gas composition data indicates that the natural gas is multi-component, the stability condition Σ is solved using the numerical root-finding method within a preset temperature range. i (y i / K i (T h , P))=1, to obtain the hydrate formation temperature at that time point, where y i Let be the mole fraction of the i-th component; By iterating through each time point in the pressure sequence and summarizing them, a hydrate formation temperature sequence is obtained.

5. The method according to claim 4, characterized in that, The stability conditions are solved using numerical root-finding methods within a preset temperature range, including: Construct the residual function f(T)=Σ i (y i / K i (T, P))-1, the residual function decreases monotonically with temperature within a preset temperature range, and a residual greater than zero corresponds to the hydrate formation region, while a residual less than zero corresponds to the outside of the formation region; A hybrid root-finding method combining inverse quadratic interpolation, secant method, and bisection method is used to find roots within a preset temperature range. The convergence conditions include absolute tolerance and relative tolerance. When the residuals at both ends of the preset temperature range have the same sign and are both positive, the hydrate formation temperature is determined to be higher than the upper limit of the preset temperature range, and the upper limit is returned as the hydrate formation temperature at that time point; when the residuals at both ends have the same sign and are both negative, the hydrate formation temperature is determined to be lower than the lower limit of the preset temperature range, and the lower limit is returned. The time point at which the root search fails is not included in the subsequent risk classification.

6. The method according to claim 5, characterized in that, In step S3, before calculating the temperature safety margin at each time point, a correction is made to the hydrate formation temperature sequence: Construct the feature vector x = [P, T, y] CH4 , y C2 , ΔT , q]; Where P is pressure, T is actual temperature, and y is... CH4 With y C2 ΔT represents the mole fractions of methane and ethane, ΔT is the initial value of the temperature safety margin at that time point, and q is the normalized flow rate; The feature vector is input into a pre-trained correction network, which includes an input layer, at least two hidden layers, and an output layer connected in sequence. Each hidden layer uses a hyperbolic tangent activation function, and the output layer outputs a correction amount with r = 3·tanh(z), where z is the linear output of the output layer, constraining the correction amount to between -3℃ and +3℃. The correction network is trained using the difference between the field-measured hydrate formation temperature sample and the value calculated by the thermodynamic model as training labels. The corrected hydrate formation temperature is used in subsequent temperature safety margin calculations.

7. The method according to claim 6, characterized in that, Step S3 specifically includes: Calculate the temperature safety margin M at each time point. t = t t - T h , t , where t t T represents the actual temperature at that time point. h , t This refers to the hydrate formation temperature at that time point; When M t If the value is greater than the warning threshold, it is considered normal. When M t A warning is issued when M is greater than the critical threshold but not greater than the warning threshold. t If the value is not greater than a critical threshold, it is determined to be serious; wherein, the critical threshold is less than or equal to the warning threshold, and when the configured warning threshold is not greater than the critical threshold, it reverts to the preset safe default value; M is mapped by piecewise linear mapping t Mapped to a risk score with values ​​ranging from 0 to 1: when M t The risk score is 1 when M is not greater than the critical threshold. t The risk score is 0 when it is not less than the warning threshold, and is adjusted according to the time interval. t = (M warn -M t ) / (M warn -M crit )calculate; If the minimum temperature safety margin M within the recent window min If the value is not greater than the critical threshold, the risk score will be increased to not less than 0.9; if M min If the risk score is not greater than the warning threshold, then the risk score will be increased to not less than 0.7 × (M) warn -M min ) / (M warn -M crit The calculated value; the risk score is truncated to the 0~1 range.

8. The method according to claim 7, characterized in that, In step S4, the risk classification result after shielding is obtained, specifically including: Acquire liquid level time series data and preprocess it in the same way as the pressure time series data to obtain a synchronized liquid level sequence; Based on the sampling interval Δt̄′ = (τ) between the synchronized liquid level sequence and the synchronized temperature sequence N -τ1) / (N -1), calculate the ten-minute normalization factor k =Δt̄′ / 600, where τ1, τ N The first and last moments of the recent window are represented by N, where N is the number of points within the window. Calculate the rate of change of liquid level L′ t = (L t - L t-1 ) / k and the rate of temperature change T′ t = (T t - T t-1 ) / k; when L′ t Greater than the preset liquid level change threshold and T′ t When the value is less than the negative of the preset temperature drop threshold, the time point is marked as the drainage period; Data points within the drainage period are not included in alarm determination. When all data points within the recent window are within the drainage period, the alarm is directly returned to the normal level. When the liquid level data is missing, the drainage shielding process is skipped.

9. The method according to claim 8, characterized in that, In step S4, continuous boundary crossing confirmation is performed on the risk classification results after shielding, including: For each time point within the recent window, an out-of-bounds marker is constructed. When the temperature safety margin at that time point is not greater than the warning threshold and the time point is not marked as a drainage period, the out-of-bounds marker is set to 1; otherwise, it is set to 0. Scan backwards from the end of the recent window and count the number of out-of-bounds points that are consecutively 1; When the number of out-of-bounds points is less than the preset confirmation threshold, the risk classification result is downgraded to non-abnormal, and the risk score is reduced to no more than half of the preset concern score; otherwise, the risk classification result is maintained, and alarm details including the start time of the out-of-bounds violation, duration, current temperature safety margin and risk score are generated and output together with the final alarm level and the risk score.

10. The method according to claim 9, characterized in that, Also includes: For multiple monitoring points in the same pipeline or the same site, steps S1 to S4 are executed respectively. The final alarm level with the highest severity among all monitoring points is taken as the overall alarm conclusion to locate high-risk nodes. Based on the increase in the temperature safety margin before and after the alarm, the effectiveness of alcohol injection and heating treatment measures is quantitatively evaluated.

Images