VOCs emission early warning method and system based on dcs timing feature mapping

Abstract

The application discloses a VOCs emission early warning method and system based on DCS time sequence feature mapping, belongs to the field of industrial process control and intelligent monitoring, and comprises the following steps: extracting DCS multi-dimensional operation parameters to form an original time sequence data set; calculating the input-output accumulated difference in a time window to obtain material energy residual time sequence; calculating the valve instruction and feedback difference to obtain instruction execution deviation time sequence; performing fast Fourier transform on the deviation time sequence to extract local energy proportion; calculating correlation coefficients to construct a parameter linkage correlation coefficient matrix; respectively solving node differences to obtain three change rate time sequences; and outputting an over-standard early warning when the three change rates synchronously exceed dynamic threshold values. The application adopts a dynamic causal chain diagnosis strategy of tracing from macro-process imbalance to micro-control fault and verifying in a closed loop by using parameter linkage relationship, and can realize high-precision, traceable root source early warning of VOCs emission risk.

Description

Technical Field

[0001] This invention relates to the field of industrial process control and intelligent monitoring, and in particular to a VOCs emission early warning method and system based on DCS time-series feature mapping. Background Technology

[0002] Industrial production processes, such as chemical, pharmaceutical, and metallurgical manufacturing, are often accompanied by the generation and emission of volatile organic compounds (VOCs). To meet increasingly stringent environmental protection requirements, it is essential to monitor and provide early warnings of VOC emission risks from industrial equipment in real time to prevent excessive emissions from harming the ecological environment and human health. Distributed control systems (DCS) are widely used in such industrial production processes, controlling and recording massive amounts of underlying equipment operating parameters in real time.

[0003] In related technologies, Chinese invention patent application CN120848441A discloses a dynamic prediction and energy-saving control system and medium for VOCs emissions of catalytic combustion equipment in a foundry workshop. The system includes a high-frequency VOCs concentration monitoring unit that uses TDLAS and FTIR spectrometers and a data fusion algorithm to obtain high-precision concentration data, where the second derivative signal from the TDLAS analyzer is extracted separately as a key feature; a data preprocessing and feature engineering module that processes multi-source data; a dynamic emission prediction module with a built-in attention mechanism temporal convolutional network (ATCN) model that achieves ultra-early prediction of emission concentration by fusing the second derivative signal with historical data; and a multi-objective optimization control module that uses a model predictive control (MPC) framework with the goal of minimizing overall operating costs. The dynamic virtual carbon penalty coefficient β can be dynamically adjusted according to the workshop production plan, and rolling optimization generates advanced control commands to adjust equipment power and fan frequency.

[0004] Regarding the aforementioned technologies, the prediction process heavily relies on additional dedicated high-precision optical gas monitoring instruments, which are highly susceptible to interference in harsh industrial environments such as high temperatures and dust, and the hardware investment and maintenance costs are high. At the same time, the prediction logic relies on deep learning black-box models, which not only require a large amount of computing power for long-term model training and online incremental updates, but also lack clear physical and industrial mechanism support, making it difficult to directly utilize the common operating parameters already available in industrial sites to achieve low-cost, highly reliable VOCs early warning. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a VOCs emission early warning method and system based on DCS time-series feature mapping. It employs a dynamic causal chain diagnostic strategy that traces macroscopic process imbalances back to microscopic control faults and utilizes parameter linkage relationships for closed-loop verification. This approach enables high-precision, traceable early warning of VOCs emission risks.

[0006] The above objectives can be achieved through the following approach: A VOCs emission early warning method based on DCS time-series feature mapping includes extracting feed flow rate time-series data, discharge flow rate time-series data, heating power time-series data, regulating valve control command time-series data, regulating valve execution feedback time-series data, and in-vessel pressure time-series data from the DCS to form a six-dimensional original time-series dataset; using the process operation cycle as the material calculation window, the feed flow rate time-series data and the heating power time-series data are accumulated window by window, and the accumulated values ​​of the discharge flow rate time-series data and the in-vessel pressure time-series data are subtracted to obtain the material energy residual time-series; the difference between the regulating valve control command time-series data and the regulating valve execution feedback time-series data is calculated point by point to obtain the command execution deviation time-series; using the process response time as the spectrum calculation window, the... The instruction execution deviation time series is subjected to Fast Fourier Transform to extract the peak frequency and the local energy ratio at the peak frequency; through the spectrum calculation window, the one-step lag autocorrelation coefficient and Pearson correlation coefficient are calculated on the six-dimensional original time series dataset to obtain the parameter linkage correlation coefficient matrix; the differences between adjacent time series nodes are calculated on the material energy residual time series, the local energy ratio, and the parameter linkage correlation coefficient matrix to obtain the residual growth rate time series, the spectrum energy change rate time series, and the linkage mutation rate time series in sequence; when the current values ​​of the residual growth rate time series, the spectrum energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectrum calculation window, the VOCs exceedance warning time series stamp and the current values ​​of the three change rates are output.

[0007] Optionally, the six-dimensional original time-series dataset is constructed by: calling the historical database interface of the DCS, parsing the process tag to separate and read the feed flow rate time-series data, discharge flow rate time-series data, heating power time-series data, regulating valve control command time-series data, regulating valve execution feedback time-series data, and in-vessel pressure time-series data to obtain DCS time-series data; performing outlier detection on the DCS time-series data, removing outlier data points and filling missing intervals using spline interpolation to generate a continuous time-series sequence; extracting the highest sampling frequency from the continuous time-series sequence as a reference time axis, performing timestamp alignment and resampling calculations on the remaining continuous time-series sequences, and performing channel dimension splicing to generate a six-dimensional original time-series dataset.

[0008] Optionally, obtaining the material energy residual time series includes: extracting the process operation cycle to construct a material calculation window, extracting the feed flow rate time series data, the heating power time series data, the discharge flow rate time series data, and the in-vessel pressure time series data; performing extreme value normalization mapping on the feed flow rate time series data and the heating power time series data to output a dimensionless input sequence, and performing extreme value normalization mapping on the discharge flow rate time series data and the in-vessel pressure time series data to output a dimensionless output sequence; performing window-by-window integral accumulation calculation on the dimensionless input sequence to generate an input accumulation value, performing window-by-window integral accumulation calculation on the dimensionless output sequence to generate an output accumulation value, and performing difference calculation between the input accumulation value and the output accumulation value to obtain the material energy residual time series.

[0009] Optionally, obtaining the command execution deviation timing includes: extracting the timing data of the control valve control command and the timing data of the control valve execution feedback, performing cross-correlation function calculation to extract the physical response hysteresis time, and using the physical response hysteresis time to perform time axis translation on the timing data of the control valve control command to generate aligned command timing data; performing numerical subtraction calculation on the aligned command timing data and the timing data of the control valve execution feedback at each sampling point to obtain the command execution deviation timing.

[0010] Optionally, the extraction of the peak frequency and the local energy ratio at the peak frequency includes: extracting the process response time to construct a spectrum calculation window; performing a fast Fourier transform on the instruction execution deviation timing to generate an amplitude-frequency feature sequence; extracting the maximum amplitude coordinate frequency of the non-zero frequency band as the peak frequency through the amplitude-frequency feature sequence; locating the effective oscillation frequency band of the amplitude-frequency feature sequence at the peak frequency; calculating the local energy integral of the effective oscillation frequency band and the global energy integral of the amplitude-frequency feature sequence respectively; and performing a division operation to obtain the local energy ratio.

[0011] Optionally, the step of performing a division operation to obtain the local energy percentage includes: extracting the peak amplitude at the peak frequency using the amplitude-frequency feature sequence, calculating the left and right boundary frequency coordinates at half the peak amplitude, defining the frequency band of the left and right boundary frequency coordinates as a full-width half-maximum (WHM) band as an effective oscillation band; performing a numerical definite integral operation on the square of the amplitude within the WHM band to generate a local energy integral, performing a numerical definite integral operation on the square of the frequency amplitude of the amplitude-frequency feature sequence to generate a global energy integral, and calculating the quotient to obtain the local energy percentage.

[0012] Optionally, obtaining the parameter linkage correlation coefficient matrix includes: extracting independent dimension sequences using the six-dimensional original time series dataset through the spectrum calculation window; calculating the inner product of a single independent dimension sequence and a single-step delayed sequence to generate a lag one-step autocorrelation coefficient; calculating the covariance normalization value of each pair of independent dimension sequences to generate a Pearson correlation coefficient; mapping the lag one-step autocorrelation coefficient sequentially to the main diagonal coordinate position of the blank matrix; symmetrically mapping the Pearson correlation coefficient to the off-diagonal coordinate position of the blank matrix; and combining them to generate the parameter linkage correlation coefficient matrix.

[0013] Optionally, obtaining the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series sequentially includes: extracting the current node value and the previous node value based on the material energy residual time series and the local energy ratio, performing arithmetic subtraction to output the residual growth rate time series and the spectral energy change rate time series respectively; extracting the current node matrix and the previous node matrix based on the parameter linkage correlation coefficient matrix, calculating the absolute value of the difference between elements at the same row and column coordinate positions item by item, and summing all the absolute values ​​of the difference to obtain the linkage mutation rate time series.

[0014] Optionally, the output of the VOCs exceeding the standard warning time stamp and the current values ​​of the three change rates include: extracting the historical sequence data of the residual growth rate time series, the spectral energy change rate time series, and the linked mutation rate time series in the spectrum calculation window; calculating the arithmetic mean and sample standard deviation of the historical sequence data respectively and performing addition to generate independent dynamic thresholds; obtaining the current values ​​of the three change rates and comparing them with the independent dynamic thresholds respectively; when all three comparison results are determined to be greater than, extracting the current control system time to generate the VOCs exceeding the standard warning time stamp and outputting it synchronously in conjunction with the current values ​​of the three change rates.

[0015] Based on the same inventive concept, this invention also provides a VOCs emission early warning system based on DCS time-series feature mapping. The system includes: a DCS time-series data acquisition module, used to extract feed flow rate time-series data, discharge flow rate time-series data, heating power time-series data, regulating valve control command time-series data, regulating valve execution feedback time-series data, and in-vessel pressure time-series data from the DCS, forming a six-dimensional original time-series dataset; a material energy residual calculation module, used to calculate the feed flow rate time-series data and the heating power time-series data window by window, subtracting the window-by-window accumulation of the discharge flow rate time-series data and the in-vessel pressure time-series data, to obtain the material energy residual time-series; an instruction execution deviation extraction module, used to calculate the difference between the regulating valve control command time-series data and the regulating valve execution feedback time-series data point by point, to obtain the instruction execution deviation time-series; and a spectrum peak energy ratio extraction module, used to calculate the difference between the process operation cycle and the process operation cycle, using the process operation cycle as the material calculation window. The process response time serves as the spectrum calculation window. A Fast Fourier Transform is performed on the instruction execution deviation time series to extract the peak frequency and the local energy proportion at that peak frequency. A parameter linkage correlation coefficient matrix construction module is used to calculate the one-step lag autocorrelation coefficient and Pearson correlation coefficient on the six-dimensional original time series dataset through the spectrum calculation window, obtaining the parameter linkage correlation coefficient matrix. A triple change rate difference calculation module is used to calculate the difference between adjacent time series nodes for the material energy residual time series, the local energy proportion, and the parameter linkage correlation coefficient matrix, sequentially obtaining the residual growth rate time series, the spectrum energy change rate time series, and the linkage mutation rate time series. A triple evidence joint early warning output module is used to output a VOCs exceeding warning time series stamp and the current values ​​of the three change rates when the current values ​​of the residual growth rate time series, the spectrum energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectrum calculation window.

[0016] Compared with the prior art, the present invention has the following advantages: 1. By capturing abnormal signals at the fundamental level of the material energy balance process, focusing on the dynamic performance of the control execution link for in-depth diagnosis, and finally verifying the closed loop through the degree of instability of the internal linkage, this multi-layer evidence chain strong correlation judgment mechanism ensures the high accuracy and reliability of the early warning and effectively avoids false alarms caused by fluctuations in a single indicator. 2. By dynamically mining early signs such as deterioration of control performance and imbalance of materials and energy hidden in DCS data, sufficient reaction time can be provided before VOCs exceed the emission standard, enabling operators to take preventive intervention measures, thereby transforming the traditional passive response environmental management into proactive preventive process safety control. 3. By intelligently segmenting the production process and establishing dynamic, adaptive baselines and thresholds for different process stages, the system can accurately identify real anomalies under specific operating conditions. This dynamic analysis mechanism enables the early warning model to adapt to dynamic changes such as normal switching and load fluctuations in the production process, ensuring the stability and reliability of early warning performance in complex and ever-changing real industrial environments.

[0017] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures pointed out in the description, claims and drawings. Attached Figure Description

[0018] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a flowchart illustrating the VOCs emission early warning method based on DCS time-series feature mapping according to an embodiment of the present invention.

[0020] Figure 2 This is a thermal distribution diagram of the parameter linkage correlation coefficient matrix in an embodiment of the present invention.

[0021] Figure 3 This is a radar distribution diagram of the triple rate of change warning triggering mechanism according to an embodiment of the present invention.

[0022] Figure 4 This is a schematic diagram of the structure of the VOCs emission early warning system based on DCS time-series feature mapping according to an embodiment of the present invention. Detailed Implementation

[0023] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, 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.

[0024] Reference Figure 1One embodiment of the present invention proposes a VOCs emission early warning method based on DCS time-series feature mapping. It adopts a dynamic causal chain diagnosis strategy that traces from macroscopic process imbalance to microscopic control faults and uses parameter linkage relationship for closed-loop verification. This can achieve high-precision and traceable early warning of VOCs emission risks.

[0025] The method described in this embodiment specifically includes: The feed flow rate time series data, discharge flow rate time series data, heating power time series data, regulating valve control command time series data, regulating valve execution feedback time series data, and in-vessel pressure time series data are extracted from the DCS to form a six-dimensional original time series dataset. Using the process operation cycle as the material calculation window, the time series data of the feed flow rate and the time series data of the heating power are accumulated window by window, and the accumulated time series data of the discharge flow rate and the time series data of the pressure inside the vessel are subtracted to obtain the material energy residual time series. The difference between the timing data of the control command and the timing data of the control valve execution feedback is calculated for each sampling point to obtain the command execution deviation timing. Using the process response time as the spectrum calculation window, a fast Fourier transform is performed on the instruction execution deviation timing to extract the peak frequency and the local energy ratio at the peak frequency; Through the spectrum calculation window, the one-step lag autocorrelation coefficient and Pearson correlation coefficient are calculated for the six-dimensional original time series dataset to obtain the parameter linkage correlation coefficient matrix. The differences between adjacent time series nodes are calculated for the material energy residual time series, the local energy ratio and the parameter linkage correlation coefficient matrix, respectively, to obtain the residual growth rate time series, the spectrum energy change rate time series and the linkage mutation rate time series in sequence. When the current values ​​of the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectral calculation window, the VOCs exceedance warning time series stamp and the current values ​​of the three change rates are output.

[0026] Optionally, the six-dimensional original time-series dataset comprises: The historical database interface of the DCS is called to parse the process tag to separate and read the feed flow rate time series data, discharge flow rate time series data, heating power time series data, regulating valve control command time series data, regulating valve execution feedback time series data, and in-vessel pressure time series data to obtain DCS time series data. To transform the dispersed and inconsistent low-level equipment data into a standardized format, the query service of the DCS distributed control system's historical database is invoked through an industry-standard open platform communication interface. Based on the process tag list in the mapping file, the feedback messages from each sensor node and control node are parsed and separated in batches. From these, the raw measurement values ​​and corresponding timestamps of feed flow rate timing data, discharge flow rate timing data, heating power timing data, regulating valve control command timing data, regulating valve execution feedback timing data, and in-vessel pressure timing data are extracted to form the initial DCS timing data.

[0027] For example, in a reactor process section for resin synthesis, data is extracted based on tag numbers. For instance, the feed flow meter with tag FI-101 is read to obtain feed flow rate timing data, and the temperature controller with tag TIC-201 is read to obtain heating power timing data. At the same time, the valve opening command sequence and the actual opening sequence returned by the field valve positioner are read. This raw data is extracted and stored in a local cache array.

[0028] Outlier detection is performed on the DCS time series data to remove outlier data points and spline interpolation is used to fill in missing intervals to generate a continuous time series. Due to electromagnetic interference or momentary sensor malfunctions, DCS time-series data often contains noise points with abrupt numerical jumps. Outlier detection based on a dynamic statistical window is performed on the data sequences of each dimension. The detection principle is to calculate the dynamic boundary using the statistical characteristics of the data within the sliding window. When the data value of a sampling point exceeds this dynamic boundary, it is identified as an outlier and removed. The formula for calculating the upper limit of the dynamic boundary is: , in, This represents the upper limit value of the dynamic boundary, used to define the maximum allowable value for normal fluctuations; It represents the arithmetic mean of all sampled values ​​within the current sliding time window, calculated by data sliding and truncation, and reflects the baseline level of the current working condition; The sample standard deviation represents the values ​​of all sampling points within the current sliding time window, reflecting the degree of fluctuation in the current operating conditions. The confidence multiplier coefficient is 3, directly set to this value based on the statistical principle of normal distribution, to cover the vast majority of normal process fluctuation ranges without requiring manual intervention. For the missing intervals left after removing outlier data points, a cubic spline interpolation algorithm is used to construct piecewise smooth polynomial curves using valid data points before and after the missing intervals for numerical extrapolation and filling, generating a continuous time series without data breaks and conforming to the physical laws of smooth process evolution.

[0029] The highest sampling frequency is extracted from the continuous time series as the reference time axis. Timestamp alignment and resampling calculations are performed on the remaining continuous time series. Channel dimension splicing is then performed to generate a six-dimensional original time series dataset.

[0030] Different types of industrial instrument hardware have different inherent sampling periods. The system iterates through the generated continuous time series of the six channels, identifying and extracting the time series with the smallest sampling interval (i.e., the highest sampling frequency), and uses the timestamp set of this sequence as the global reference time axis. For the remaining continuous time series with sampling frequencies lower than this highest sampling frequency, timestamp alignment and linear resampling are performed using this reference time axis to ensure that at each time node of the reference time axis, all six physical variables have corresponding valid values. Finally, the six timestamp-aligned one-dimensional sequences are concatenated along the column channel direction to construct a structurally consistent two-dimensional matrix, i.e., the six-dimensional original time series dataset.

[0031] For example, assume that the sampling frequency of the valve execution feedback time series data is the highest, 10 times per second, with a time interval of 100 milliseconds; while the sampling frequency of the feed flow rate time series data is once per second, with a time interval of 1000 milliseconds. Use the 100-millisecond interval as the baseline time axis. For the feed flow rate time series data, between every two original values ​​at 1000-millisecond intervals, generate nine interpolated values ​​at equal intervals using a linear interpolation algorithm, increasing its frequency to be completely consistent with the baseline time axis. After resampling all variables, these six sequences are concatenated along the channel dimension, ultimately generating a six-dimensional original time series dataset containing six columns of data, with the number of rows matching the number of nodes on the high-frequency time axis, providing a well-organized data source for matrix operations.

[0032] Optionally, obtaining the material energy residual time series includes: The process operation cycle is extracted to construct a material calculation window, and the feed flow rate time series data, the heating power time series data, the discharge flow rate time series data, and the pressure time series data inside the vessel are extracted. First, the average time taken for a complete production batch of the target process equipment to run continuously is read from the underlying production scheduling database. This average time is defined as the process operation cycle, and a material calculation window is constructed based on the length of this process operation cycle. Then, using this material calculation window as the time span, a sliding truncation operation is performed on the time axis. From the original time-series dataset containing six dimensions, the feed flow rate time-series data and heating power time-series data belonging to the input end, and the discharge flow rate time-series data and in-vessel pressure time-series data belonging to the output end and internal status feedback, are selectively separated and extracted.

[0033] For example, in a specific chemical resin synthesis reaction section, the average time for a complete reaction batch is two hours. This two-hour period is extracted as the process operation cycle, and a material calculation window with a width of two hours is constructed. As production time progresses, time-series data records of feed flow rate, heating power, discharge flow rate, and reactor pressure are continuously extracted from the high-frequency dataset of the DCS at fixed two-hour intervals, providing time-strictly aligned data slices for energy and mass conservation calculations.

[0034] The feed flow rate time series data and the heating power time series data are subjected to extreme value normalization mapping to output a dimensionless input sequence, and the discharge flow rate time series data and the in-vessel pressure time series data are subjected to extreme value normalization mapping to output a dimensionless output sequence. Because flow velocity, power, and pressure have completely different physical dimensions, they cannot be directly added or subtracted mathematically. Therefore, extreme value normalization mapping is performed on the four types of time-series data, proportionally mapping all the original measurement values ​​with engineering physical units to a pure numerical range of 0 to 1. The extreme value normalization mapping formula used in the normalization calculation unit is: , in, This represents the dimensionless value output after the extreme value normalization mapping; This represents the actual measured value of the time-series data read at the current sampling point; This represents the minimum observed extreme value of the corresponding physical quantity in the historical steady-state operation record; This represents the maximum observed extreme value of the corresponding physical quantity in the historical steady-state operation record. Both the minimum and maximum values ​​are objectively obtained by performing an extreme value retrieval algorithm on all data from the DCS's historical normal, leak-free operation, ensuring data-driven operation without manual intervention. The feed flow rate time-series data, after formula transformation, is combined with the heating power time-series data to output a dimensionless input sequence; similarly, the transformed discharge flow rate time-series data is combined with the in-vessel pressure time-series data to output a dimensionless output sequence.

[0035] For example, suppose the minimum observed extreme value of the feed flow rate in historical normal records is 10 cubic meters per hour, and the maximum observed extreme value is 50 cubic meters per hour. When the actual measured value of the feed flow rate read at the current sampling point is 20 cubic meters per hour, substituting it into the formula yields a dimensionless value of 0.25. Through this mathematical mapping, the feed flow rate, heating power, discharge flow rate, and pressure inside the vessel, which originally had different physical units, are all converted into relative intensity indicators between 0 and 1, completely eliminating the physical dimension barrier.

[0036] A window-by-window integration and accumulation calculation is performed on the dimensionless input sequence to generate an input accumulation value. A window-by-window integration and accumulation calculation is performed on the dimensionless output sequence to generate an output accumulation value. The difference between the input accumulation value and the output accumulation value is calculated to obtain the material energy residual time series.

[0037] To quantify the macroscopic material and energy imbalance within the entire material calculation window, a discrete integral summation is performed on the dimensionless sequence over time. First, the values ​​of all sampling points in the dimensionless input sequence are accumulated to obtain the input accumulated value, representing the relative total input intensity within the calculation window. Then, the values ​​of all sampling points in the dimensionless output sequence are accumulated to obtain the output accumulated value, representing the relative total consumption and state output intensity within the calculation window. Finally, the input accumulated value is subtracted from the output accumulated value. As the material calculation window slides continuously, this difference calculation will continuously output a series of results that change over time, thus obtaining a continuous material-energy residual time series. The underlying mathematical formula for this calculation logic is: , in, This represents the residual material energy value calculated at the current sliding moment t. This represents performing a summation operation on all discrete sampling points within the current material calculation window; This represents the dimensionless value of the feed velocity in the dimensionless input sequence. This represents the dimensionless value of the heating power in the dimensionless input sequence. This represents the dimensionless value of the discharge velocity in the dimensionless output sequence. This represents the dimensionless value of the pressure inside the vessel in the dimensionless output sequence. Physically, if the relative total input of a closed industrial environment consistently exceeds the combined characteristics of the relative total output and the state pressure, the resulting positive residual indicates that the material and energy have not been completely converted into the target product, but are highly likely to be undergoing unexpected physical dissipation in invisible gaseous forms such as volatile organic compounds.

[0038] For example, within a two-hour material calculation window, the dimensionless values ​​of all sampling points within the window are summed. The calculated cumulative input value corresponding to the feed and heating power is 1500.5, and the cumulative output value corresponding to the discharge and in-vessel pressure is 1490.0. A differential calculation is performed, subtracting 1490.0 from 1500.5 to obtain the material energy residual value at the corresponding moment in the window, which is 10.5. As time progresses, the window slides forward according to the sampling step size, continuously repeating the above addition and subtraction calculations, ultimately outputting a material energy residual time series that continuously evolves with the process progress, providing the first macroscopic data evidence for leak early warning.

[0039] Optionally, obtaining the instruction execution deviation timing includes: Extract the timing data of the control command of the regulating valve and the timing data of the execution feedback of the regulating valve, perform cross-correlation function calculation to extract the physical response hysteresis time, and use the physical response hysteresis time to perform time axis translation on the timing data of the control command of the regulating valve to generate aligned command timing data; The aim is to objectively quantify the dynamic difference between control commands and actual actions, eliminating the inherent time delay caused by mechanical inertia. Timing data of control commands and feedback from the control valve are extracted from the underlying database. The control command timing data represents the target valve position opening sequence output by the controller, while the feedback timing data represents the actual opening sequence returned by the valve position sensor. Since there is an inherent time delay in the mechanical structure from receiving the electrical signal to pneumatic execution and stabilization, direct comparison will result in physical misalignment. Discrete cross-correlation function calculations are performed on these two sets of sequences. By sliding along the time axis to calculate the point-by-point product sum of the two sets of sequences, the time offset corresponding to the maximum peak value of the product sum is found. This time offset corresponding to the peak value of the product sum is extracted as the physical response hysteresis time. After obtaining this delay time, a time axis translation operation is performed on the control command timing data, delaying its timestamps by the length of one physical response hysteresis time, thereby generating aligned command timing data. The mathematical formula for this translation operation is: , in, This represents the specific value of the alignment instruction timing data generated at the current time node t; This represents the original value of the timing data of the control command of the regulating valve before the historical shift occurred. This represents the physical response hysteresis time objectively extracted through the peak value calculation of the cross-correlation function; Represents the current time point; This represents the historical time point obtained by subtracting the physical response lag time from the current time point. The physical meaning of the formula is that it compares the control commands issued in the past with the actual mechanical feedback in the same time dimension. Its value is determined entirely by the cross-correlation peak calculation of the historical data of the actual valve action in the field, avoiding the errors caused by artificially setting fixed delay parameters.

[0040] For example, suppose in an industrial operation log, the controller issues a control command to adjust the valve opening to 50% at the 10th second, while the actual valve opening at 50% is not shown in the feedback data until the 12th second. Performing a cross-correlation function calculation on these two continuous sequences lasting several minutes reveals that the waveform overlap is highest when the command sequence is shifted forward by 2 seconds. Therefore, the physical response hysteresis time is extracted to be 2 seconds. Subsequently, the value at the 10th second of the command sequence is shifted and aligned to the 12th second time node, generating aligned command timing data. Through this timeline shift, the aligned command value at the 12th second is 50%, achieving synchronization with the actual 50% value feedback at the 12th second on the timeline.

[0041] The timing data of the alignment command is subtracted from the timing data of the control valve execution feedback at each sampling point to obtain the command execution deviation timing.

[0042] After eliminating physical inertia delay, the internal processing unit performs a basic arithmetic subtraction calculation on a sample-by-sample basis between the alignment instruction timing data and the control valve execution feedback timing data for data nodes at the same timestamp. This subtraction operation directly eliminates the normal mechanical motion lag, and the remaining pure numerical difference is the instruction execution deviation timing. The mathematical formula for this calculation process is: , in, This represents the numerical value of the instruction execution deviation timing calculated and output at the current time node t; This represents the value of the alignment instruction timing data read at the current time node t; This represents the numerical value of the feedback timing data of the control valve execution actually collected at the current time node t. The physical meaning of the formula lies in measuring the actual frictional resistance and mechanical wear of the actuator. If this difference deviates from zero for a long time or exhibits high-frequency oscillation, it directly indicates that there are substantial mechanical faults inside the pneumatic valve, such as valve stem jamming or seal damage, thus eliminating the interference of normal process delays.

[0043] For example, after the translation and alignment operation, at the sampling point of 15 seconds, the alignment command timing data value at that moment is read as 60%. Simultaneously, the actual valve execution feedback timing data value transmitted back by the field sensor at 15 seconds is read as 58%. Substituting these values ​​into the difference formula, the difference between 60% and 58% is calculated, yielding an instantaneous difference value of 2% at 15 seconds. As the time axis progresses through each sampling point, this subtraction operation is continuously performed. Finally, the instantaneous difference values ​​obtained from all sampling points are arranged in chronological order, outputting a complete and continuous command execution deviation timing sequence. This provides a high-frequency waveform data source reflecting the characteristics of microscopic mechanical faults for the Fast Fourier Transform.

[0044] Optionally, the extracted peak frequency and the proportion of local energy at the peak frequency include: The process response time is extracted to construct a spectrum calculation window. A fast Fourier transform is performed on the instruction execution deviation timing to generate an amplitude-frequency feature sequence. The maximum amplitude coordinate frequency of the non-zero frequency band is extracted as the peak frequency through the amplitude-frequency feature sequence. First, the process response time is extracted by reading historical open-loop step test records stored in the underlying database. The process response time represents the time span required from a step change in the manipulated variable to the controlled variable reaching a new steady state. A spectrum calculation window is constructed based on the length of this process response time, and using this window as a fixed data length boundary, instruction execution deviation timing within this span is truncately extracted on the time axis. Subsequently, the underlying digital signal processing chip performs a Fast Fourier Transform (FFT) algorithm on the extracted instruction execution deviation timing. The FFT algorithm calculates the amplitude-frequency characteristic sequence by orthogonally projecting discrete time-domain data points onto frequency-domain basis functions, which is used to map the magnitude of the deviation signal at each discrete frequency point. To remove the DC component representing the static steady-state error, the amplitude-frequency characteristic sequence is traversed, skipping the zero-frequency point, and searching only within the non-zero frequency band for the abscissa frequency corresponding to the highest amplitude. This abscissa frequency is extracted as the peak frequency to objectively characterize the most dominant mechanical oscillation or closed-loop instability frequency in the current control loop.

[0045] For example, for the pressure control loop inside a typical chemical reactor, by analyzing historical step test data, the process response time of this process object is extracted to be 600 seconds, and a 600-second spectrum calculation window is constructed based on this. The instruction execution deviation timing data within the most recent 600 seconds is extracted and subjected to a Fast Fourier Transform. The amplitude-frequency characteristic sequence output after the transform shows that the DC component amplitude at 0 Hz is 2.0, the amplitude at 0.05 Hz is 0.5, and the oscillation amplitude at 0.1 Hz reaches a maximum peak value of 5.0. The 0 Hz DC data is automatically filtered out, and the coordinates corresponding to the maximum amplitude value of 5.0 are locked within the non-zero frequency band, thus accurately extracting 0.1 Hz as the peak frequency.

[0046] Locate the effective oscillation frequency band of the amplitude-frequency characteristic sequence at the peak frequency, calculate the local energy integral of the effective oscillation frequency band and the global energy integral of the amplitude-frequency characteristic sequence respectively, and perform a division operation to obtain the local energy ratio.

[0047] After identifying the dominant oscillation frequency, the degree of dominance of this oscillation in the overall control deviation energy is quantified. In the amplitude-frequency characteristic sequence, a numerical traversal is performed to the left and right of the peak frequency in both low-frequency and high-frequency directions to find the two frequency boundary points corresponding to the amplitude decreasing to half of the peak amplitude. The range enclosed by these two frequency boundary points is defined as the effective oscillation frequency band. Subsequently, the local energy integral within the effective oscillation frequency band and the global energy integral covering all non-zero effective calculated frequencies are calculated. Finally, a division operation is performed, dividing the calculated local energy integral by the global energy integral to obtain the local energy percentage in decimal form. The mathematical logic of this calculation process is expressed by the following definite integral calculation formula: , in, It represents the local energy percentage of the final output and is an objectively quantified dimensionless ratio. The representative amplitude-frequency characteristic sequence at any frequency point The corresponding actual amplitude at that location; This represents the lower boundary frequency of the effective oscillation band extracted using the half-amplitude method; The upper boundary frequency representing the effective oscillation band; This represents performing continuous integration or discrete summation operations within the effective oscillation frequency band, which is related to... The calculated numerator is the local energy integral, which physically represents the local mechanical energy contained in that specific oscillation waveform. This represents the smallest non-zero frequency point in the amplitude-frequency characteristic sequence after excluding the zero frequency point. This represents the upper limit frequency of the Nyquist limit supported by the Fast Fourier Transform; This represents the continuous integration or discrete summation operation performed over the entire non-zero full-frequency range. The denominator term calculated is the global energy integral, which physically represents the total abnormal fluctuation energy generated during valve actuation. This ratio is obtained through division. A higher value indicates that the deviation behavior of the actuator is more concentrated at a single specific frequency, which can be used to accurately determine the instability of the control system or the resonance state of the equipment.

[0048] For example, in a scenario involving an internal pressure control loop, the highest amplitude measured at a peak frequency of 0.1 Hz is 5.0. Searching to both sides, the amplitude at 0.08 Hz and 0.12 Hz is found to decay to exactly 2.5, thus precisely defining the 0.08 Hz to 0.12 Hz range as the effective oscillation band. Discretely summing the squares of the amplitudes within this effective oscillation band yields a local energy integral of 150.0. Simultaneously, summing the squares of the amplitudes across the entire non-zero frequency band from 0.01 Hz to the Nyquist threshold frequency of 0.5 Hz yields a global energy integral of 200.0. Performing a division operation, the result is... The final output local energy ratio is 0.75. This high value directly reveals that 75% of the error energy in the loop is caused by a single 0.1 Hz oscillation, providing solid data evidence for troubleshooting mechanical jamming or control parameter divergence.

[0049] Optionally, the step of performing a division operation to obtain the local energy percentage includes: The peak amplitude is extracted at the peak frequency by the amplitude-frequency feature sequence, and the left and right boundary frequency point coordinates at half the peak amplitude are calculated. The frequency band of the left and right boundary frequency point coordinates is defined as the half-width full-width frequency band as an effective oscillation frequency band. First, the amplitude-frequency characteristic sequence is read, and the vertical coordinate value corresponding to the peak frequency is located. This vertical coordinate value is then extracted as the peak amplitude. Next, half of this peak amplitude is calculated using a basic division operator, and this half value is used as the baseline threshold. Centered on the peak frequency, numerical comparisons are performed on discrete frequency points in both the low-frequency and high-frequency directions to locate the two critical frequency points where the curve descends and intersects the baseline threshold. The horizontal coordinates corresponding to these two critical frequency points are extracted as the left and right boundary frequency coordinates, respectively. The continuous frequency span between the left and right boundary frequency coordinates is directly defined as the full-width half-maximum (HWHM) bandwidth, and this HWHM bandwidth is established as the effective oscillation band, thus completely eliminating the black-box operation of manually specifying a fixed bandwidth.

[0050] For example, in a control loop abnormal fluctuation analysis scenario, the peak amplitude extracted at the peak frequency of 0.5 Hz is 10.0. Half of 10.0 is calculated to obtain the baseline threshold of 5.0. Searching from 0.5 Hz to the left in the low-frequency region, it is found that the amplitude decays to 5.0 at 0.4 Hz, and 0.4 Hz is extracted as the left boundary frequency point coordinate; searching to the right in the high-frequency region, it is found that the amplitude decays to 5.0 at 0.6 Hz, and 0.6 Hz is extracted as the right boundary frequency point coordinate. The frequency band from 0.4 Hz to 0.6 Hz is defined as the full-width half-maximum frequency band, which is used as the effective oscillation frequency band objectively intercepted for integral calculation.

[0051] A local energy integral is generated by performing a numerical definite integral operation on the square of the amplitude within the full width at half maximum (WHM) frequency band, and a global energy integral is generated by performing a numerical definite integral operation on the square of the frequency point amplitude of the amplitude-frequency characteristic sequence. The quotient is then calculated to obtain the local energy percentage.

[0052] After establishing objective integration boundaries, the final quantization calculation is performed on the energy characteristics. In discrete signal processing, the physical energy carried by a signal is strictly proportional to the square of its vibration amplitude. First, for all discrete frequency points within the full width at half maximum (FWHM) bandwidth, the square of the amplitude corresponding to each frequency point is calculated, and all squared values ​​are discretely summed. This discrete summation operation is equivalent to the numerical definite integral operation of a continuous function in discrete calculus, and the accumulated result directly generates the local energy integral. Similarly, for the entire effective frequency range from the starting frequency excluding the DC zero frequency point to the Nyquist limit frequency in the amplitude-frequency characteristic sequence, the same discrete summation operation is performed on the squares of the amplitudes of all frequency points to generate the global energy integral. Finally, the local energy integral is divided by the global energy integral to obtain a pure decimal quotient representing the concentration of oscillations, thus obtaining the local energy proportion. The underlying core mathematical formula of this calculation process is: , in, It represents the local energy percentage of the final output, and is a dimensionless pure numerical value. Its physical meaning is the proportion of energy occupied by a specific oscillation in the overall control valve fluctuation. The representative amplitude-frequency characteristic sequence in the frequency variable The corresponding actual objective amplitude; The square of the amplitude is used to characterize the anomalous physical energy carried at that specific frequency. This represents the coordinates of the left boundary frequency point obtained from the retrieval. This represents the coordinates of the right boundary frequency point obtained in the previous step; This represents the summation operation performed for each discrete frequency step within the effective oscillation frequency band, and its overall molecular calculation result represents the local energy integral. Represents the coordinates of the smallest effective discrete frequency point greater than zero, used to completely filter out static, non-fluctuating DC offset signals; This means that the Fast Fourier Transform is limited by the sampling theorem, which determines the coordinates of the highest effective frequency point that can be analyzed. This represents the summation operation performed within the global non-zero fluctuation frequency band, and the overall calculation result of its denominator represents the global energy integral. All coordinate limits and amplitude variables are directly derived from the high-frequency objective physical sampling of the underlying sensors, completely eliminating the need to rely on expert experience to set weights or guess the filtering frequency band.

[0053] For example, suppose the defined full-width half-maximum (FWHM) band range of 0.4 Hz to 0.6 Hz contains five discrete frequency points. The amplitudes of these five frequency points are extracted sequentially, multiplied by themselves, squared, and summed to calculate a local energy integral of 80.0. Next, the amplitudes of all 100 non-zero discrete frequency points from 0.01 Hz to the highest frequency of 1.0 Hz are squared and summed to calculate a global energy integral of 200.0. The internal arithmetic logic unit executes a basic division instruction, dividing 80.0 by 200.0, ultimately yielding a local energy percentage of 0.4. This value clearly quantifies that up to 40% of the total deviation energy during the device's operating cycle is highly concentrated in a narrowband oscillation at a specific frequency of 0.5 Hz, providing an irrefutable data basis for early warning of mechanical resonance.

[0054] Optionally, the obtained parameter linkage correlation coefficient matrix includes: Through the spectrum calculation window, independent dimension sequences are extracted using the six-dimensional original time series dataset. The inner product of a single independent dimension sequence and a single-step delayed sequence is calculated to generate a one-step lag autocorrelation coefficient. The covariance normalization value of each pair of independent dimension sequences is calculated to generate a Pearson correlation coefficient. Using a defined spectral calculation window, the original six-dimensional time-series dataset is decomposed along the column vector direction to extract six single independent dimension sequences. For indices requiring characterization of single-parameter motion inertia, each single independent dimension sequence is extracted, and its time axis is shifted backward by one underlying hardware sampling step to generate a single-step delay sequence. Subsequently, the values ​​of the original sequence and the single-step delay sequence at the same time are multiplied term by term, summed to calculate the inner product, and then divided by the sequence variance to complete numerical normalization, generating a one-step lag autocorrelation coefficient. This coefficient objectively quantifies the memory effect of the physical parameter itself. For indices requiring characterization of the synchronous linkage effect of different parameters, pairwise independent dimension sequences are randomly selected from the six independent dimension sequences, and the covariance between them is calculated. This covariance is divided by the product of the standard deviations of the two sequences to obtain the normalized covariance value, thus generating the Pearson correlation coefficient. The mathematical formula for generating the Pearson correlation coefficient is: , in, The Pearson correlation coefficient, representing the generated output, is a dimensionless pure numerical value ranging from negative one to positive one, characterizing the linear coupling strength of synchronous fluctuations between two different process physical parameters. Representing the The independent dimension sequence and the first The covariance of each independent dimension sequence is calculated by the expected value of the product of the difference between the discrete values ​​of the two independent dimension sequences and their arithmetic mean within the bottom sampling window. Representing the The sample standard deviation of each independent dimension sequence within the current spectrum calculation window; Representing the The standard deviation of each independent dimension sequence within the current spectrum calculation window. All variables in the formula are directly derived from the original six-dimensional time-series dataset acquired objectively at high frequencies. The entire calculation process is entirely driven by on-site data, without any subjective weighting based on expert experience.

[0055] For example, suppose that within the current spectrum calculation window, independent dimension sequences representing the feed flow rate and independent dimension sequences representing the vessel pressure are extracted. The inner product of the feed flow rate sequence and its own sequence delayed by one time step is calculated, yielding an autocorrelation coefficient of 0.85. This indicates that the current feed flow rate possesses strong smoothness and inertia, with no significant drastic fluctuations. Simultaneously, the covariance between the feed flow rate sequence and the vessel pressure sequence is calculated to be 4.0. The standard deviations of the feed flow rate and vessel pressure are calculated to be 2.0 and 2.5 respectively. Substituting these values ​​into the above formula, the Pearson correlation coefficient is obtained. This high value directly reveals a very strong and positive physical correlation between the feed flow rate and the pressure inside the vessel.

[0056] The lag one-step autocorrelation coefficients are sequentially mapped to the main diagonal coordinates of the blank matrix, and the Pearson correlation coefficients are symmetrically mapped to the off-diagonal coordinates of the blank matrix, thus generating a parameter linkage correlation coefficient matrix.

[0057] After obtaining all correlation coefficient quantification results, a six-row, six-column two-dimensional array structure is allocated and initialized in computer memory, defined as a blank matrix. Following a pre-determined fixed index arrangement of the six physical variables, the six calculated one-step lag autocorrelation coefficients are sequentially written into the main diagonal coordinates of this blank matrix. The main diagonal coordinates are the coordinates where the row and column indices in the two-dimensional array are exactly equal. Next, the Pearson correlation coefficients between all different variable pairs are interleaved and filled into the off-diagonal coordinates of the blank matrix according to the row and column indices of the variables involved in the calculation, and a symmetric mapping operation is enforced. The off-diagonal coordinates are the coordinates where the row and column indices are not equal. Through the symmetric mapping logic, the value in the first row, second column is guaranteed to be strictly equal to the value in the second row, first column. After filling all thirty-six coordinate positions, the two-dimensional array structure is packaged to generate a parameter-linked correlation coefficient matrix. This matrix presents a strictly symmetrical array at any time slice, forming a digital fingerprint of the subtle global coupling variations in the process system.

[0058] For example, a 6x6 initial blank matrix is ​​created in memory. For the first dimension variable representing the feed flow rate, its calculated one-step autocorrelation coefficient of 0.85 is directly written to the [1,1] main diagonal coordinate position with row index 1 and column index 1. Subsequently, for the Pearson correlation coefficient of 0.8 calculated between the feed flow rate and the pressure inside the vessel, it is simultaneously written to the [1,6] off-diagonal coordinate position in the first row and sixth column, and the [6,1] off-diagonal coordinate position in the sixth row and first column. After thirty-six such absolutely objective row and column addressing and assignment instructions, a complete six-dimensional parameter linkage correlation coefficient matrix is ​​assembled and output. This matrix organically weaves six isolated fluid and thermodynamic parameters into a mathematical grid, preparing a multidimensional data carrier for calculating the overall mutation rate of the spatial topology. Figure 2 As shown, the Pearson correlation and one-step autocorrelation between six key process parameters, such as feed flow rate and pressure inside the reactor, within a specific spectrum calculation window are intuitively quantified by the grayscale depth, which objectively reflects the tightness of synchronous linkage between multidimensional parameters.

[0059] Optionally, obtaining the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series sequentially includes: Based on the material energy residual time series and the local energy ratio, extract the current node value and the previous node value, and perform arithmetic subtraction to output the residual growth rate time series and the spectrum energy change rate time series respectively. For a one-dimensional scalar time series, at each set discrete calculation step node, the latest generated value of the current node in the material energy residual time series and the value of the immediately preceding node on the time axis are read synchronously. Simultaneously, the latest generated value of the current node in the local energy percentage and the value of the immediately preceding node are read. Then, using the subtraction instructions of the internal arithmetic logic unit, the value of the current node is subtracted from the value of the preceding node, performing basic arithmetic subtraction operations. This mathematical action of directly calculating the first-order forward difference strips away the absolute amplitude of the data, retaining only the relative slope of the numerical change. The difference results of the material energy residual time series are arranged in chronological order, outputting the residual growth rate time series; similarly, the difference results of the local energy percentage are arranged in chronological order, outputting the spectral energy change rate time series.

[0060] For example, when the process reaches the calculation node of the tenth minute, the current node value of the material energy residual time series is extracted as 10.5, and the previous node value of the ninth minute is extracted as 10.0. The arithmetic subtraction is performed to obtain a difference of 0.5, and this 0.5 is recorded as the residual growth rate at that moment. Similarly, for the local energy proportion, the current node value of the tenth minute is extracted as 0.75, and the previous node value of the ninth minute is 0.60. The arithmetic subtraction is performed to obtain a difference of 0.15, and this 0.15 is recorded as the spectral energy change rate at that moment.

[0061] Based on the parameter linkage correlation coefficient matrix, extract the current node matrix and the previous node matrix, calculate the absolute value of the difference between the elements at the same row and column coordinate positions, and sum all the absolute values ​​of the difference to obtain the linkage mutation rate time series.

[0062] At the current time point, the most recently assembled parameter linkage correlation coefficient matrix in memory is extracted as the current node matrix, and the parameter linkage correlation coefficient matrix of the immediately preceding time point stored in the historical cache is extracted as the previous node matrix. To capture relationship mutations in any direction, for these two 6x6 two-dimensional arrays, a double-loop traversal algorithm is used to extract the corresponding coordinate positions of elements with identical row and column indices, perform subtraction, and take the absolute value of the difference. Finally, the calculated absolute values ​​of the thirty-six differences are globally summed to obtain a single scalar value, namely the linkage mutation rate at the current time, which is then appended to form a continuous linkage mutation rate time series. The mathematical formula for this core matrix dimensionality reduction calculation process is: , in, The linkage mutation rate, representing the calculated output, is a dimensionless scalar value. Its physical meaning is the total distorted dynamic energy of the topological relationship of the six core parameters of the entire industrial process. The number of dimensions representing process parameters is fixed at 6 in this algorithm structure, representing the number of physical variables in the six-dimensional original time series dataset; and These represent the row index and column index in the two-dimensional matrix, respectively, with values ​​ranging from 1 to 6. Represents the node located at the current node in the matrix. Line number The objective numerical values ​​of the elements in the column; Represents the position of the node in the previous node matrix. Line number The objective numerical values ​​of the elements in the column; This represents the result of taking the absolute value of the difference between the numerical values ​​of elements at the same coordinate position. By designing an operator that takes the absolute value, whether two parameters suddenly change from positive to negative correlation, or from uncorrelated to strongly correlated, all absolute changes will be equivalently accumulated into a positive anomalous indicator. This completely avoids the mathematical loophole of positive and negative differences canceling each other out in matrix summation, and perfectly reduces the high-dimensional sixth-order square matrix, which is difficult to judge intuitively, into a one-dimensional time series feature.

[0063] For example, when comparing the parameter linkage correlation coefficient matrix of two adjacent time steps, these thirty-six elements are traversed. Assume that at the coordinate position in the first row and sixth column, the Pearson correlation coefficient recorded by the previous node matrix is ​​0.8, while the correlation value at the current node matrix drops sharply to 0.2 due to the sudden breakdown of the physical coupling between the feed and pressure at the site. The difference between these two is calculated and its absolute value is taken, yielding an absolute difference of 0.6 for this coordinate point. The same absolute difference calculation is then performed on the remaining thirty-five coordinate positions, assuming the total change at the remaining positions is 1.8. Finally, an accumulation and summation calculation is performed, adding 0.6 to 1.8 to obtain 2.4, outputting the linkage mutation rate value at this moment as 2.4. As time progresses, this scalar value is continuously generated, eventually forming a complete linkage mutation rate time series, providing extremely sensitive comprehensive high-dimensional variation evidence for ultimately judging complex faults such as control loop decoupling or material imbalance.

[0064] Optionally, the output VOCs exceedance warning time stamp and the current values ​​of the three change rates include: Extract the historical sequence data of the residual growth rate time series, the spectral energy change rate time series, and the linked mutation rate time series in the spectrum calculation window, calculate the arithmetic mean and sample standard deviation of the historical sequence data respectively, and perform addition to generate independent dynamic thresholds; For the time series of residual growth rate, spectral energy change rate, and linked mutation rate, historical sequence data of these three indicators within the time span of the spectral calculation window are extracted. For each type of historical sequence data, its arithmetic mean and sample standard deviation are independently calculated using underlying statistical functions. Subsequently, for the same indicator, the calculated arithmetic mean and sample standard deviation are added together, and the sum is used to generate an independent dynamic threshold for each rate of change indicator that perfectly matches the fluctuations of the current actual operating conditions. The mathematical formula for this calculation process is expressed as follows: , in, Represented by a specific rate of change indicator The calculated independent dynamic threshold, where These can refer to the residual growth rate, the rate of change of spectral energy, or the linkage mutation rate, respectively. Represents this specific indicator The arithmetic mean within a fixed historical time window, in physical terms, is the basic drift center of the underlying process in the current operating stage, which can automatically and smoothly move with the rise and fall of the process load; This represents the sample standard deviation of a specific indicator within the same historical time window, and its physical meaning is the degree of dispersion of the underlying data around the mean. The physical and statistical basis for adding the two to generate a threshold is that, according to the characteristics of a normal distribution, the mean plus one standard deviation constitutes a one-sided upper limit dynamic envelope covering approximately 84% of the confidence interval. When actual operating data exceeds this envelope, it mathematically proves that the deviation is no longer a normal process background noise disturbance.

[0065] For example, assume the set spectrum calculation window duration is 15 minutes. Extract the time-series data of the residual growth rate over the past 15 minutes, calculate the arithmetic mean of the residual growth over these 15 minutes as 0.20, and the sample standard deviation as 0.05. Perform addition, adding 0.20 and 0.05 to obtain 0.25, thus setting the independent dynamic threshold for the residual growth rate to 0.25. Similarly, independently calculate the independent dynamic threshold for the spectral energy change rate as 0.60, and the independent dynamic threshold for the linked mutation rate as 2.00. These three thresholds are not hard-coded into the program beforehand, but are objectively derived entirely from the actual operating data of the past 15 minutes.

[0066] The current values ​​of the three rates of change are obtained and compared with the independent dynamic thresholds respectively. When all three comparison results are determined to be greater than the threshold, the current control system time is extracted to generate a VOCs exceedance warning time stamp and output synchronously in conjunction with the current values ​​of the three rates of change.

[0067] At the latest time calculation node, the current values ​​of three change rates—residual growth rate, spectral energy change rate, and linkage mutation rate—are synchronously acquired. These three current change rate values ​​are then compared with the three corresponding independent dynamic thresholds generated in the previous step. To completely eliminate false alarms caused by single sensor failures or localized random process disturbances, a strict AND gate architecture is employed. Only when the residual growth rate, spectral energy change rate, and linkage mutation rate all exceed their respective thresholds are the three comparison results considered true. Once these triple concurrent triggering conditions are confirmed, the hardware real-time clock interface of the industrial field control bus is immediately invoked to read the current control system time, accurate to milliseconds, and encapsulated to generate a unique and tamper-proof volatile organic compound (VOC) exceedance warning time stamp. Finally, this time stamp is concatenated with the current values ​​of the three change rates that triggered the alarm, forming a complete diagnostic evidence chain that is synchronously output to the host computer or manufacturing execution system interface.

[0068] For example, at the current time point, the latest values ​​of the three change rates are calculated in real time: residual growth rate is 0.30, spectral energy change rate is 0.80, and linked mutation rate is 3.50. These are then compared with the previously calculated corresponding independent dynamic thresholds (0.25, 0.60, and 2.00). Since 0.30 is greater than 0.25, 0.80 is greater than 0.60, and 3.50 is greater than 2.00, the Boolean comparison results for all three dimensions are true. Therefore, a high-confidence risk of abnormal VOCs emission is immediately determined. The system then reads the underlying RTC chip time as "February 27, 2026, 10:15:30.250 milliseconds", generates a high-precision warning time stamp, and packages the data message "Time: [2026-02-27 10:15:30.250], residual growth rate: 0.30, spectral energy change rate: 0.80, linkage mutation rate: 3.50" and outputs it synchronously to the central control room's large screen. This output, carrying a complete set of underlying traceability parameters, allows on-site engineers to quickly pinpoint whether the problem is a pipeline leak or a valve jamming based on specific exceedances. Figure 3 As shown, the diagram uses an intuitive spatial polygonal shape to compare and demonstrate the synchronous breakthrough states of the residual growth rate, spectral energy change rate, and linkage mutation rate relative to the adaptive mean plus standard deviation dynamic threshold boundary under normal operating conditions and VOCs exceeding the standard conditions, effectively illustrating the high-confidence multidimensional concurrent verification logic.

[0069] Based on the same inventive concept, this invention also provides a VOCs emission early warning system based on DCS time-series feature mapping, such as... Figure 4 As shown, the system includes: The DCS timing data acquisition module is used to extract feeding flow rate timing data, discharging flow rate timing data, heating power timing data, regulating valve control command timing data, regulating valve execution feedback timing data, and in-vessel pressure timing data from the DCS to form a six-dimensional original timing dataset. The material energy residual calculation module is used to calculate the material energy residual time series by adding the feed flow rate time series data and the heating power time series data window by window, and subtracting the accumulated values ​​of the discharge flow rate time series data and the pressure time series data inside the vessel from the window-by-window values ​​of the process operation cycle. The instruction execution deviation extraction module is used to calculate the difference between the timing data of the control instruction of the regulating valve and the timing data of the execution feedback of the regulating valve at each sampling point to obtain the instruction execution deviation timing. The spectrum peak energy ratio extraction module is used to perform a fast Fourier transform on the instruction execution deviation timing with the process response time as the spectrum calculation window, and extract the peak frequency and the local energy ratio at the peak frequency. The parameter linkage correlation coefficient matrix construction module is used to calculate the one-step lag autocorrelation coefficient and Pearson correlation coefficient of the six-dimensional original time series dataset through the spectrum calculation window to obtain the parameter linkage correlation coefficient matrix. The triple rate of change difference calculation module is used to calculate the difference between adjacent time series nodes for the material energy residual time series, the local energy ratio and the parameter linkage correlation coefficient matrix, respectively, and obtain the residual growth rate time series, the spectrum energy change rate time series and the linkage mutation rate time series in sequence. The triple evidence joint early warning output module is used to output a VOCs exceeding warning time stamp and the current values ​​of the three change rates when the current values ​​of the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectrum calculation window.

[0070] It should be noted that the functional division and information interaction between the various modules described above are logical, but in terms of physical implementation, they can be integrated on the same software platform or deployed in a distributed manner. The connections between them represent data flow and control flow, aiming to collaboratively achieve the objectives of this invention. The above descriptions are merely exemplary embodiments of this invention and should not be construed as limiting the scope of protection of this invention.

Claims

1. A VOCs emission early warning method based on DCS time-series feature mapping, characterized in that, The method includes: The feed flow rate time series data, discharge flow rate time series data, heating power time series data, regulating valve control command time series data, regulating valve execution feedback time series data, and in-vessel pressure time series data are extracted from the DCS to form a six-dimensional original time series dataset. Using the process operation cycle as the material calculation window, the time series data of the feed flow rate and the time series data of the heating power are accumulated window by window, and the accumulated time series data of the discharge flow rate and the time series data of the pressure inside the vessel are subtracted to obtain the material energy residual time series. The difference between the timing data of the control command and the timing data of the control valve execution feedback is calculated for each sampling point to obtain the command execution deviation timing. Using the process response time as the spectrum calculation window, a fast Fourier transform is performed on the instruction execution deviation timing to extract the peak frequency and the local energy ratio at the peak frequency; Through the spectrum calculation window, the one-step lag autocorrelation coefficient and Pearson correlation coefficient are calculated for the six-dimensional original time series dataset to obtain the parameter linkage correlation coefficient matrix. The differences between adjacent time series nodes are calculated for the material energy residual time series, the local energy ratio and the parameter linkage correlation coefficient matrix, respectively, to obtain the residual growth rate time series, the spectrum energy change rate time series and the linkage mutation rate time series in sequence. When the current values ​​of the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectral calculation window, the VOCs exceedance warning time series stamp and the current values ​​of the three change rates are output.

2. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The six-dimensional original time-series dataset includes: The historical database interface of the DCS is called to parse the process tag to separate and read the feed flow rate time series data, discharge flow rate time series data, heating power time series data, regulating valve control command time series data, regulating valve execution feedback time series data, and in-vessel pressure time series data to obtain DCS time series data. Outlier detection is performed on the DCS time series data to remove outlier data points and spline interpolation is used to fill in missing intervals to generate a continuous time series. The highest sampling frequency is extracted from the continuous time series as the reference time axis. Timestamp alignment and resampling calculations are performed on the remaining continuous time series. Channel dimension splicing is then performed to generate a six-dimensional original time series dataset.

3. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The obtained material energy residual time series includes: The process operation cycle is extracted to construct a material calculation window, and the feed flow rate time series data, the heating power time series data, the discharge flow rate time series data, and the pressure time series data inside the vessel are extracted. The feed flow rate time series data and the heating power time series data are subjected to extreme value normalization mapping to output a dimensionless input sequence, and the discharge flow rate time series data and the in-vessel pressure time series data are subjected to extreme value normalization mapping to output a dimensionless output sequence. A window-by-window integration and accumulation calculation is performed on the dimensionless input sequence to generate an input accumulation value. A window-by-window integration and accumulation calculation is performed on the dimensionless output sequence to generate an output accumulation value. The difference between the input accumulation value and the output accumulation value is calculated to obtain the material energy residual time series.

4. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The obtained instruction execution deviation timing includes: Extract the timing data of the control command of the regulating valve and the timing data of the execution feedback of the regulating valve, perform cross-correlation function calculation to extract the physical response hysteresis time, and use the physical response hysteresis time to perform time axis translation on the timing data of the control command of the regulating valve to generate aligned command timing data; The timing data of the alignment command is subtracted from the timing data of the control valve execution feedback at each sampling point to obtain the command execution deviation timing.

5. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The extracted peak frequency and the local energy ratio at the peak frequency include: The process response time is extracted to construct a spectrum calculation window. A fast Fourier transform is performed on the instruction execution deviation timing to generate an amplitude-frequency feature sequence. The maximum amplitude coordinate frequency of the non-zero frequency band is extracted as the peak frequency through the amplitude-frequency feature sequence. Locate the effective oscillation frequency band of the amplitude-frequency characteristic sequence at the peak frequency, calculate the local energy integral of the effective oscillation frequency band and the global energy integral of the amplitude-frequency characteristic sequence respectively, and perform a division operation to obtain the local energy ratio.

6. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 5, characterized in that, The local energy percentage obtained by performing the division operation includes: The peak amplitude is extracted at the peak frequency by the amplitude-frequency feature sequence, and the left and right boundary frequency point coordinates at half the peak amplitude are calculated. The frequency band of the left and right boundary frequency point coordinates is defined as the half-width full-width frequency band as an effective oscillation frequency band. A local energy integral is generated by performing a numerical definite integral operation on the square of the amplitude within the full width at half maximum (WHM) frequency band, and a global energy integral is generated by performing a numerical definite integral operation on the square of the frequency point amplitude of the amplitude-frequency characteristic sequence. The quotient is then calculated to obtain the local energy percentage.

7. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The obtained parameter linkage correlation coefficient matrix includes: Through the spectrum calculation window, independent dimension sequences are extracted using the six-dimensional original time series dataset. The inner product of a single independent dimension sequence and a single-step delayed sequence is calculated to generate a one-step lag autocorrelation coefficient. The covariance normalization value of each pair of independent dimension sequences is calculated to generate a Pearson correlation coefficient. The lag one-step autocorrelation coefficients are sequentially mapped to the main diagonal coordinates of the blank matrix, and the Pearson correlation coefficients are symmetrically mapped to the off-diagonal coordinates of the blank matrix, thus generating a parameter linkage correlation coefficient matrix.

8. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The time series of residual growth rate, spectral energy change rate, and linkage mutation rate obtained sequentially include: Based on the material energy residual time series and the local energy ratio, extract the current node value and the previous node value, and perform arithmetic subtraction to output the residual growth rate time series and the spectrum energy change rate time series respectively. Based on the parameter linkage correlation coefficient matrix, extract the current node matrix and the previous node matrix, calculate the absolute value of the difference between the elements at the same row and column coordinate positions, and sum all the absolute values ​​of the difference to obtain the linkage mutation rate time series.

9. The VOCs emission early warning method based on DCS time-series feature mapping according to claim 1, characterized in that, The output VOCs exceedance warning time stamp and the current values ​​of the three change rates include: Extract the historical sequence data of the residual growth rate time series, the spectral energy change rate time series, and the linked mutation rate time series in the spectrum calculation window, calculate the arithmetic mean and sample standard deviation of the historical sequence data respectively, and perform addition to generate independent dynamic thresholds; The current values ​​of the three rates of change are obtained and compared with the independent dynamic thresholds respectively. When all three comparison results are determined to be greater than the threshold, the current control system time is extracted to generate a VOCs exceedance warning time stamp and output synchronously in conjunction with the current values ​​of the three rates of change.

10. A VOCs emission early warning system based on DCS time-series feature mapping, applied to the VOCs emission early warning method based on DCS time-series feature mapping as described in any one of claims 1-9, characterized in that, The system includes: The DCS timing data acquisition module is used to extract feeding flow rate timing data, discharging flow rate timing data, heating power timing data, regulating valve control command timing data, regulating valve execution feedback timing data, and in-vessel pressure timing data from the DCS to form a six-dimensional original timing dataset. The material energy residual calculation module is used to calculate the material energy residual time series by adding the feed flow rate time series data and the heating power time series data window by window, and subtracting the accumulated values ​​of the discharge flow rate time series data and the pressure time series data inside the vessel from the window-by-window values ​​of the process operation cycle. The instruction execution deviation extraction module is used to calculate the difference between the timing data of the control instruction of the regulating valve and the timing data of the execution feedback of the regulating valve at each sampling point to obtain the instruction execution deviation timing. The spectrum peak energy ratio extraction module is used to perform a fast Fourier transform on the instruction execution deviation timing with the process response time as the spectrum calculation window, and extract the peak frequency and the local energy ratio at the peak frequency. The parameter linkage correlation coefficient matrix construction module is used to calculate the one-step lag autocorrelation coefficient and Pearson correlation coefficient of the six-dimensional original time series dataset through the spectrum calculation window to obtain the parameter linkage correlation coefficient matrix. The triple rate of change difference calculation module is used to calculate the difference between adjacent time series nodes for the material energy residual time series, the local energy ratio and the parameter linkage correlation coefficient matrix, respectively, and obtain the residual growth rate time series, the spectrum energy change rate time series and the linkage mutation rate time series in sequence. The triple evidence joint early warning output module is used to output a VOCs exceeding warning time stamp and the current values ​​of the three change rates when the current values ​​of the residual growth rate time series, the spectral energy change rate time series, and the linkage mutation rate time series all exceed the mean plus standard deviation dynamic threshold of the spectrum calculation window.

Images