A multivariate time series causal connection identification method, system and storage medium

Abstract

The application discloses a multivariate time series causal connection identification method and system and a storage medium; a processor calculates one-step prediction error by using a prediction regression window up to t-1, independently estimates a candidate edge evaluation full model and a reduced model by using a causal analysis window containing t, calculates a conditional prediction error ratio index and a conditional prediction residual difference signal. Set ensemble empirical mode decomposition is performed on the candidate cause variable signal and the conditional prediction residual difference signal, and original multi-scale causal scores are obtained by pairing according to intrinsic mode function orders; a candidate directed graph is formed by sequentially performing double-threshold joint screening, indirect path batch screening and bidirectional edge direction advantage judgment. The system switches identification states, an effective state and a re-identification state according to a rolling average error stability criterion and a change candidate graph two-out-of-three voting, and outputs a directed adjacency matrix only in the effective state, which is used to generate an industrial fault propagation path or an electroencephalogram connection topological feature.

Description

Technical Field

[0001] This invention relates to the fields of multivariate time series processing, data-driven system identification, and conditional causal inference, and particularly to a method, system, and storage medium for identifying causal connections in multivariate time series signals. In industrial applications, it is used to generate fault propagation paths between industrial process variables, and in EEG applications, it is used to generate directed effective connection networks and EEG connection topology features between EEG channels or brain regions. Background Technology

[0002] Continuous industrial processes, energy systems, and complex electromechanical equipment typically deploy sensors for flow, pressure, level, temperature, current, and vibration. Upstream disturbances, equipment failures, or changes in control parameters propagate along the process chain, causing multiple measuring points to change approximately simultaneously. Establishing causal relationships solely based on correlations or two-variable predictions can easily misjudge indirect correlations propagated from intermediate variables as direct effects, and may also continue to output outdated causal graphs after changes in operating conditions.

[0003] The linear multivariate Granger method determines the directional predictive effect by comparing the residuals of the full model and the reduced model after removing candidate variables, but the linear lag model is difficult to fully fit the nonlinear coupling of state dependence; the comparison of residual variance itself does not provide a multi-scale signal object corresponding to the candidate causal variables on an hourly basis.

[0004] Existing methods for analyzing causal relationships in multivariate time-series signals typically focus on single-scale predictive relationships, identification of static or staged causal structures, or rely on specific models and pre-defined data conditions for connectivity analysis. In situations involving nonlinear coupling, multi-scale dynamic changes, and structural switching, problems may arise such as indirect effects being identified as direct effects, difficulties in determining the direction of bidirectional relationships, and the continued use of invalid causal connection results after system state changes. Furthermore, some existing methods lack multi-scale analysis objects that correspond to the interactions of candidate variables on an hourly basis, and online control mechanisms for the validity of causal connection results.

[0005] Therefore, it is necessary to establish a unified technology chain that is independent of specific industrial variables or EEG channel names and can be repeatedly implemented to complete the online identification of direct directional predictive effects between multivariate time-series signals. This would enable the identification of direct directional predictive effects under complex nonlinear and multi-scale conditions, suppress indirect paths and directional unresolved relationships, and further allow for online judgment and control of the effectiveness of causal connection results when the system structure changes. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method, system and storage medium for identifying causal connections between multivariable time series signals. It utilizes the comparison of nonlinear full model and reduced model, conditional prediction residual difference signal, ensemble empirical mode decomposition, multi-scale graph decision and online state control to identify the direct directional predictive effects between multivariable time series signals online.

[0007] The objective of this invention is achieved through the following technical solution:

[0008] The processor acquires M synchronously observed variables through a multi-channel signal acquisition device 100 and a data acquisition and synchronization interface 110. The multi-channel signal acquisition device can be an industrial process sensor array or an EEG signal acquisition device. The observation vector at discrete time t is represented as:

[0009]

[0010] Where M is the number of variables, This represents the numerical sample of the i-th variable at time t. The window management module 120 maintains a length of... Predictive regression window and length is Causal analysis window The prediction regression window strictly ends at time t-1 and is used to predict the current sample X(t) that has not yet participated in the model estimation; the causal analysis window includes the current sample X(t) and is used to re-estimate the candidate edges to evaluate the model at the candidate graph update time.

[0011]

[0012]

[0013] and All values ​​are preset positive integers, which can be the same or different, and there is no mandatory size relationship between them; the length of each window should not be less than the minimum number of samples required to complete the corresponding lag model estimation, residual calculation, and signal decomposition. When the number of samples is less than the fixed length, the reached samples form an expanded window; after the fixed length is reached, a sliding window is used. (Freezing example) , However, this set of values ​​does not constitute a limit on the window length.

[0014] The preprocessing includes performing timestamp alignment according to a unified clock, marking missing samples, marking or clipping obvious out-of-boundary value or artifact samples, and recording variable scales within the prediction regression window and the causal analysis window respectively. Applicable data are multivariate time series data with a unified sampling time base, preset finite lag order, and sufficient effective variations in corresponding windows. The frozen implementation manner does not use the end segment of the complete sequence or future samples to determine the convergence band.

[0015] Specifically, the online state prediction model and the candidate edge evaluation model are as follows:

[0016] For any ordered variable pair , is referred to as a candidate cause variable, and is referred to as a to-be-explained outcome variable. This designation represents the role of the variables in the current round of test; before completing subsequent graph decision, it is not pre-determined that the two have formed a causal relationship. The same process variable can be used as a candidate cause variable or a to-be-explained outcome variable in different ordered variable pairs respectively.

[0017] Given the maximum lag order p and the highest polynomial order Q, an original nonlinear lag feature vector without intercept is defined. This vector is composed of the 1st to p-th order lag values of each variable, the 2nd to Q-th order powers of the lag values of each variable, and pairwise products of different variables at the same lag order; the intercept of each model is set as an independent constant term:

[0018]

[0019] where D is the number of original non-intercept features, q ranges from 1 to M, ranges from 1 to p, r ranges from 2 to Q, and a and b satisfy 1≤a<b≤M. Each non-intercept feature records the variable set it contains; accordingly, when deleting the candidate cause variable, its linear term, univariate polynomial term, and pairwise product terms containing this variable can be deleted at the same time. The frozen implementation manner does not include interaction terms of three or more elements, nor does it include cross-lag order product terms.

[0020] The online state prediction model is only used to evaluate the current prediction error and control the online state. At each sampling time t, the processor uses the prediction regression window strictly cut off at t-1 to construct a design matrix. Using the sample mean and sample standard deviation of each non-intercept feature in this window as normalization parameters, the original feature vector φ is converted into a prediction window normalized feature vector , and the online state prediction model containing all variable features is independently estimated; the current sample that has not participated in this model estimation has a one-step prediction value of:

[0021]

[0022] The superscript (on) indicates online, i.e., online status prediction; and For online state prediction models targeting outcome variables Independently estimated intercept and non-intercept coefficient vectors, This is the current time-time feature vector standardized using the predictive regression window statistics. It is only used for the instantaneous prediction error and rolling error statistics described later, and is not used as a candidate edge in the conditional prediction residual difference signal to evaluate the full model prediction value; the superscript (on) is independent of p, which represents the maximum lag order.

[0023] like Figure 3 As shown, the candidate edge evaluation full model and the candidate edge evaluation reduction model only utilize the causal analysis window containing the current sample X(t) when updating the candidate graph or checking the internal graph. Re-evaluate. Let The time-aligned index set is obtained by subtracting the initial samples required to form the first p-order lag features from this window; the causal window standardized feature vector is obtained by using the sample mean and sample standard deviation of each non-intercept feature in this index set as standardization parameters. Unexplained outcome variable Candidate edge evaluation full model Retaining all D non-intercept features, their in-window fitted values ​​and residuals are as follows:

[0024]

[0025]

[0026] The superscript (f) indicates full; and These represent the estimated intercept and D-dimensional non-intercept coefficient vector of the candidate edge evaluation full model, respectively. The fitted value at index τ within the same causal analysis window. This corresponds to the full model residuals. This model is used to generate the in-window residuals required for candidate edge comparison, and is not used to predict the current samples outside the prediction window.

[0027] For candidate cause variables By removing all components containing j from the variable set in the causal window standardized eigenvector, we obtain... ,in Less than D. Candidate edge evaluation reduction model The design matrix was independently re-estimated using the same causal analysis window, the same time-aligned index, and the deleted design matrix, without using candidate edges to evaluate the coefficients of any retained features in the full model; the in-window fitted values ​​and residuals were as follows:

[0028]

[0029]

[0030] Here, the superscript (r,j) indicates reduced-by-j, that is, for the candidate causal variable. The resulting reduced model; and The intercept and the independent estimates of the reduced model are respectively 3D non-intercept coefficient vector This corresponds to the reduced model residuals. Both the full model residuals and the reduced model residuals are in the same... The calculation is performed on top and aligned time-by-time; the former only requires calculation on each result variable. The latter estimates once, for each candidate causal variable. Remove the features and re-estimate them.

[0031] The online state prediction model, the full candidate edge evaluation model, and the candidate edge evaluation reduction model all employ ridge regression with unpenalized intercept. For any model m, let... The column of all-1 intercepts is combined with the normalized non-intercept eigenma matrix of the model to form the augmented design matrix. , corresponding result variable The target vector is composed of samples at the effective index of the model. And the intercept and non-intercept coefficients are combined into an augmented coefficient vector. Then the estimation formula can be uniformly written as:

[0032]

[0033] Where, if model m has One valid regression sample and If there are non-intercept features, then λ is a non-negative ridge parameter. Its first diagonal element is 0 to ensure the intercept is not penalized. The online model m=on uses a predictive regression window and does not include X(t); the full model m=f and the reduced model m=(r,j) use the same causal analysis window containing X(t) and the same time-aligned index, the only difference being that the reduced model removes... Relevant features are re-evaluated independently.

[0034] Conditional prediction error ratio index and conditional prediction residual difference signal

[0035] Candidate causal variables Pointing to the variable of the result to be explained The conditional prediction error ratio is defined as the logarithm of the ratio of the variance of the residuals in the candidate edge evaluation reduced model to the variance of the residuals in the candidate edge evaluation full model after hysteresis alignment of the current causal analysis window.

[0036]

[0037] in, To prevent positive numbers with a denominator of zero. It is the dimensionless logarithm ratio of the residual variance: a positive value indicates deletion. The variance of the residuals in the reduced model increases after considering relevant features; a larger value indicates... right The more important the prediction step, the better. This indicator is calculated directly from the residual variance ratio, and is not mutual information, transitive entropy, or other Shannon information; this statement is only used to define its mathematical properties, and this distinction will not be repeated later.

[0038] To transform the candidate variable removal operation into a signal object that can be decomposed over time at multiple scales, a conditional prediction residual difference signal is constructed:

[0039]

[0040] Conditional prediction of residual difference signal Indicates deletion After considering the relevant features, the reduced model adds a residual component at time τ compared to the full model. This signal has the same characteristics as... With the same physical dimensions, they can be directly used as time-sequential multi-scale decomposition objects.

[0041] Empirical Mode Decomposition of Two Object Sets

[0042] like Figure 4 As shown, the candidate causal variable signals are aligned by lag within the causal analysis window. Conditional prediction residual difference signal As two decomposed objects processed separately. To ensure repeatability and avoid the two objects directly reusing the same generated noise array, the processor uses a deterministic base seed. Candidate edges And the signal type identifier 1 or 2, respectively generating the first pseudo-random noise seed and the second pseudo-random noise seed:

[0043]

[0044] Where H is a deterministic hash mapping, number 1 represents the candidate causal variable signal, and number 2 represents the conditional prediction residual difference signal. The freezing software combines the basic seed, candidate edge identifier, and signal type identifier into a hash input, and constructs a 32-bit unsigned integer seed using the first 4 bytes of the SHA-256 digest. The same basic seed, candidate edges, and signal type produce the same decomposition results; the two signal types respectively initialize the corresponding ensemble empirical mode decomposition process and generate reproducible auxiliary noise sequences.

[0045] Ensemble empirical mode decomposition is performed on the candidate causal variable signals and the conditional prediction residual difference signals, respectively. For each object, white noise generated by its corresponding seed is added in each ensemble trial, and the ensemble mean is calculated for components of the same order. The two sets of intrinsic mode functions are arranged according to the order of the decomposition outputs; if the number of modes are respectively... and Then keep A pair of modes of the same order. The freezing implementation only pairs modes by order, without performing center frequency or energy centroid matching.

[0046] Multiscale phase coherence, mean modal variance and variance weighting

[0047] The intrinsic mode function of the m-th candidate causal variable signal and the intrinsic mode function of the conditional prediction residual difference signal Perform Hilbert transform to obtain the instantaneous phase. and The modal phase coherence value is:

[0048]

[0049] Where N is the number of window samples participating in the coherent computation, and abs represents the complex modulus function. Located in a closed interval To simultaneously reflect the amplitude changes of both the candidate causal variable mode and the residual difference mode, the arithmetic mean of the two-sided variances is first calculated:

[0050]

[0051] The sum of all mean modal variances is ,in Index for modal summation. This applies when the sum of the variances is greater than the positive constant. When the average modal variance is normalized, the sum of the variances is not greater than 100%. At this time, all weights and aggregate scores are set to zero, instead of using uniform weights:

[0052]

[0053]

[0054] When the sum of the above variances is not greater than At that time, all and . The candidate edge's original multi-scale causal score lies within a closed interval. The raw scores are used for subsequent path comparisons and direction comparisons across the variables to be explained.

[0055] VII. Dual threshold joint screening, indirect path batch screening, and bidirectional edge direction dominance determination

[0056] For the same unexplained outcome variable The original multi-scale causal scores of each candidate causal variable are normalized into incoming edge weights:

[0057]

[0058] Here, q iterates through all candidate cause variables except i. Indicate candidate edges The original multiscale causal fraction. When the denominator is not greater than season Normalized incoming edge weights Located in a closed interval And it is only used for comparing candidate sources within the same explained outcome variable. When Not less than the incoming edge threshold and Not less than the error ratio threshold At that time, Retained to the initial candidate directed graph This process is called dual-threshold joint screening. The freezing implementation uses a "not less than" comparison, and both the code and formulas include equal signs as boundaries.

[0059] The fixed candidate directed graph obtained after dual-threshold joint screening Complete all indirect path determinations above. Let This represents nodes from j to i, with lengths from 2 to M-1, where nodes are unique and do not contain any direct edges to be compared. The set of all indirect directed paths. For any path in this set... Path bottleneck score Defined as the minimum of the original multi-scale causal scores of all candidate edges of the path:

[0060]

[0061] Candidate edges Maximum indirect explanation score Defined as the maximum value among all qualified indirect paths:

[0062]

[0063] If at least one of the aforementioned indirect paths exists and Then the direct candidate edges Add to the set to be deleted, where The indirect path suppression coefficient is positive. The processor operates on a fixed candidate directed graph. After determining all candidate edges and all indirect paths, all candidate edges in the set to be deleted are removed at once; this process is called batch indirect path removal, and the deletion result does not depend on the candidate cause variable, the variable to be explained, or the traversal order of the paths. For the three-variable case, the indirect paths degenerate into... .

[0064] After completing the batch filtering of indirect paths, for those that still exist... and The variables are used to determine the dominance of the bidirectional edge direction, where The directional dominance coefficient is greater than 1. If Only retain ;like Only retain Otherwise, mark the variable pair as having undecided directions and delete both directions. This yields the directed adjacency matrix:

[0065]

[0066] in, This indicates that candidate causal variables are identified under the current valid state. Treating the explanatory outcome variable Directed predictive function; This indicates that the signal failed the dual-threshold joint screening, was deleted through batch screening of indirect paths, was marked as direction undecided after bidirectional edge direction dominance determination, or the direct edge does not exist. The direct directed prediction effect is the connection relationship determined by the conditional prediction and multi-scale graph decision rules of this invention; in EEG applications, when it is referred to as directed effective connectivity, the calculation result is not directly equivalent to anatomical connectivity.

[0067] Rolling average error stability criterion and change candidates Figure 2 Take three votes

[0068] At each sampling time t, the processor uses a prediction regression window that is strictly limited to t-1. Establish an online state prediction model to obtain one-step predictions for current samples that have not yet been included in this model estimation. The standard deviation of each variable within the predicted regression window is used as the basis for the prediction. and positive scale lower limit Normalized prediction error, instantaneous multivariate normalized error is:

[0069]

[0070] The rolling average error of length L is:

[0071]

[0072] Let the most recent H rolling average errors form a set:

[0073]

[0074] The rolling error statistic is considered to satisfy the rolling average error stability criterion and has accumulated a complete causal analysis window if and only if the following three conditions are met simultaneously:

[0075]

[0076] in, As the starting point for current identification, This is the length of the causal analysis window; during the initial identification period, This serves as the initial identification starting point; during the re-identification process, Candidates for change Figure 2 The moment of change is determined by three votes. The third condition indicates that... The cumulative number of samples has reached Therefore, the entire causal analysis window used for causal identification consists of data from the current identification stage. The freeze implementation does not estimate the so-called final stable band from the end of the complete sequence.

[0077] like Figure 5 As shown, when the system first uses the three states, they are defined as the identification state (IDENTIFY), the valid state (VALID), and the re-identification state (REIDENTIFY). The system enters the valid state when all three conditions are met in the identification state or the re-identification state.

[0078] In the valid state (VALID), the processor computes the directed adjacency matrix once upon entering this state, and then performs an internal candidate graph update every U samples; the previous valid directed adjacency matrix is ​​maintained between adjacent graph update times. This maintenance behavior is not equivalent to replacing the valid graph with the changed candidate graph.

[0079] when The error alarm threshold is not less than D consecutive time intervals. When this occurs, a prediction error alarm is generated and an internal candidate graph check is immediately triggered, but the prediction error alarm itself does not directly declare the valid causal graph invalid. If the internal candidate graph is different from the currently valid one... Figure 1 If the change is consistent, then the candidate records for change are cleared and the record remains valid.

[0080] When an internal candidate graph differs from the currently valid graph, the changed candidate graph is recorded. The freeze implementation method continuously saves the three most recent changed candidate graphs and executes the changed candidate graph. Figure 2Three-way voting: For any candidate edge, if the edge appears at least twice in the three candidate change graphs, it is included in the voting graph. If the voting graph is still different from the current valid graph, the state enters the REIDENTIFY state, and the earliest time of the candidate change graph participating in the voting is recorded as the change evidence time. The same voting criterion is used for periodic internal chart checks and internal checks triggered by prediction error alarms.

[0081] During the REIDENTIFY phase, one more step of prediction and model identification is performed, but the changed candidate map is not considered as a valid output. The cumulative number of samples since the change evidence time reaches [number missing]. And once the rolling average error stability criterion is satisfied again, the state is restored to a valid state (VALID), and this moment is recorded as _____. If the physical switching time Given from the truth experiment, the detection delay and recovery delay are respectively:

[0082]

[0083] In real-world online deployments, the physical switchover time is usually unknown, and therefore the output data packets are recorded. , Status and reason codes; reported in verification experiments with known switching times. and .

[0084] IX. System and Output Interfaces

[0085] like Figure 2 As shown, the online causal connection identification system for multivariable time-series signals includes a multi-channel signal acquisition device 100, a data acquisition and synchronization interface 110, a processor, a memory, and a result output interface 180. The processor executes a window management module 120, a state prediction and full model and reduced model identification module 130, a conditional prediction calculation module 140, an ensemble empirical mode decomposition and multi-scale coherent aggregation module 150, a graph decision module 160, and a convergence determination and state control module 170.

[0086] The memory stores the algorithm version, deterministic foundation seed, window length, model order, threshold, current state, previous valid directed adjacency matrix, and state transition log. The result output interface 180 outputs a result data packet containing the sequence number, timestamp, algorithm version, state enumeration, valid flag, and reason code. In the IDENTIFY and REIDENTIFY states, the valid flag is 0, the result output interface 180 does not provide a change candidate graph, and prohibits the corresponding downstream analysis modules from using the change candidate graph. In the VALID state, the valid flag is 1, and the result output interface 180 provides a valid directed adjacency matrix in the result data packet. When the multi-channel signal acquisition device 100 is an industrial process sensor group, the result output interface 180 provides the valid directed adjacency matrix and valid flag to the industrial fault propagation monitoring module; the industrial fault propagation monitoring module determines the directed connections corresponding to elements with a value of 1 in the matrix as direct fault propagation edges between industrial process variables, and generates industrial fault propagation paths according to the direction of the directed edges and the connection relationship between the first and last nodes. When the multi-channel signal acquisition device 100 is an EEG signal acquisition device, the result output interface 180 provides the effective directed adjacency matrix and effective markers to the EEG intention state analysis module. The EEG intention state analysis module determines the effective directed adjacency matrix as a directed effective connection network between EEG channels or brain regions, and generates EEG connection topology features based on this network as input data for intention state analysis. This invention can also be implemented as a program stored in a storage medium. When the processor executes this program, it completes data window updates, online state prediction, candidate edge evaluation model estimation, conditional predictor calculation, dual-object EEMD, multi-scale aggregation, deterministic graph decision, state control, and effective result output.

[0087] The beneficial effects of this invention are: compared with the scheme that only uses linear residual variance comparison, this invention improves the nonlinear one-step prediction capability through polynomial lag terms and paired coupling terms, and ensures that the residual comparison corresponds to the candidate variable deletion operation through the reduction model of independent re-estimation.

[0088] Compared with the scheme of separately decomposing the full model residuals and the reduced model residuals, the present invention directly decomposes the conditional prediction residual difference signal formed by subtracting the full model residuals from the reduced model residuals, so that the multi-scale analysis object focuses on the time-by-time residual part added after deleting the relevant features of the candidate causal variables.

[0089] By initializing the empirical mode decomposition process of the set for candidate causal variable signals and conditional prediction residual difference signals based on deterministic foundation seeds, candidate edge identifiers, and signal type identifiers, and generating reproducible auxiliary noise sequences, it is possible to avoid the two decomposition objects directly reusing the same generated noise array, while preserving experimental reproducibility.

[0090] By performing batch filtering of indirect paths on a fixed candidate directed graph, namely, determining indirect paths with no duplicate nodes, calculating path bottleneck scores and deleting the set to be deleted at once, and performing bidirectional edge direction dominance determination, the graph decision results do not depend on the traversal order and can handle indirect interpretation paths with a length of more than two hops when there are more than three variables.

[0091] By strictly using only historical data, the rolling average error stability criterion, prediction error alarm, and change candidates are established. Figure 2 The dual-evidence approach with three-vote voting, the complete causal window constraint, and the invalid state output gating can pause the changed candidate graph after structural changes and restore the output of the valid directed graph after re-identification. The same core processing chain can process industrial process variables and multi-channel EEG signals respectively without changing the candidate edge evaluation, dual-object EEMD, and state control rules. It can also convert the directed adjacency matrix output only in the valid state into industrial fault propagation paths or EEG connection topology features, thereby avoiding downstream modules from performing fault propagation analysis or intention state analysis based on candidate graphs that have not yet stabilized or have changed. Attached Figure Description

[0092] Figure 1 This is a flowchart of the online causal connection identification method for multivariable time-series signals of the present invention, and shows the downstream processing branches for generating industrial fault propagation paths and generating EEG connection topology features under effective conditions.

[0093] Figure 2 This is a block diagram of the online causal connection identification system for multivariable time-series signals of the present invention, and shows the data connection relationship between the result output interface 180 and the external industrial fault propagation monitoring module and the EEG intention state analysis module, respectively.

[0094] Figure 3 This diagram illustrates the construction of the causal analysis window, the candidate edge evaluation full model, the candidate edge evaluation reduced model, the conditional prediction error ratio index, and the conditional prediction residual difference signal.

[0095] Figure 4 This is a schematic diagram of dual-object EEMD using separately generated and reproducible auxiliary noise sequences, paired by IMF order, and multi-scale coherent aggregation.

[0096] Figure 5 This diagram illustrates the transitions between the IDENTIFY state, the VALID state, and the REIDENTIFY state.

[0097] Figure reference numerals: 100-Multi-channel signal acquisition device; 110-Data acquisition and synchronization interface; 120-Window management module; 130-State prediction and full model and reduced model identification module; 140-Conditional prediction quantity calculation module; 150-Ensemble empirical mode decomposition and multi-scale coherent aggregation module; 160-Graph decision module; 170-Convergence determination and state control module; 180-Result output interface. Detailed Implementation

[0098] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings, but the scope of protection of the present invention is not limited to the following description.

[0099] This invention utilizes a prediction regression window up to the current time to complete online state prediction; it independently estimates candidate edges using a causal analysis window including the current time to evaluate the full and reduced models and generate time-aligned residuals; it forms a conditional prediction residual difference signal with a time index from the residual difference; it obtains multi-scale components using separately generated and reproducible EEMD random streams; it suppresses indirect paths and undecided edges using explicit, ergonomically independent graph decision rules; and it controls when the directed connection graph can be provided as a valid result to the industrial fault propagation monitoring module or the EEG intention state analysis module through an online state machine that strictly uses only historical data, and uses it to generate industrial fault propagation paths or EEG connection topology features, respectively. The overall flow of this invention's multivariate time-series signal online causal connection identification method is as follows: Figure 1 As shown in the figure, the structural block diagram of the online causal connection identification system for multivariable time-series signals of the present invention is as follows. Figure 2 As shown, the construction diagrams of the causal analysis window, the full candidate edge evaluation model, the reduced candidate edge evaluation model, the conditional prediction error ratio index, and the conditional prediction residual difference signal are as follows. Figure 3 As shown, a schematic diagram of dual-object EEMD, pairing according to IMF order, and multi-scale coherent aggregation using separately generated and reproducible auxiliary noise sequences is illustrated below. Figure 4 As shown in the diagram, the transitions between the IDENTIFY state, VALID state, and REIDENTIFY state are illustrated below. Figure 5 As shown.

[0100] The following examples illustrate the implementable process of the present invention. Unless otherwise specified, the frozen software version mapped and re-verified item by item in this application is mRCAD-EEMD-FROZEN-1.1.0, and the Python version used for verification is 3.12.13. Fixed parameters include: p=2, Q=2, λ=10 -6 The EEMD ensemble has 5 trials, a noise width of 0.2, and retains a maximum of 3 IMFs, providing a deterministic foundational seed. Inbound threshold Error ratio than threshold Indirect path suppression coefficient directional dominance coefficient .

[0101] Online parameters include: minimum number of training samples for prediction (80), and prediction regression window. Causal analysis window Rolling length L=20, stable holding length H=20, convergence threshold Stable bandwidth Error alarm threshold The alarm maintains a length D=3, continuously saves 3 candidate graphs of change and performs a 2-out-of-3 voting process edge by edge, the graph update cycle U=10, and the scale lower limit. .

[0102] Example 1: Three-variable stochastic nonlinear propagation chain:

[0103] The following stochastic nonlinear system is used:

[0104]

[0105]

[0106]

[0107] The initial state is [0.15, -0.10, 0.20]; , and The samples are independent and follow a Gaussian distribution with a mean of 0 and a standard deviation of 0.10. The first 100 transitional samples are discarded, and the subsequent 300 samples are used. The true direct edges are... and , Only by mediating variables Forming an indirect path.

[0108] Ten replicates were performed using the baseline data seeds 240519 to 240528. The mean precision, recall, F1 score, and accuracy of the frozen method were all 1.000, and the mean NRMSE for regression was 0.698. The mean precision of the linear multivariate Granger baseline was 0.333, the mean recall was 1.000, the F1 score was 0.500, and the mean accuracy was 0.333. NRMSE represents the root mean square error of the one-step fit, normalized to the standard deviation of each channel sample.

[0109] Example 2: Coupled Logistic Propagation Chain:

[0110] The following coupled Logistic system is adopted:

[0111]

[0112]

[0113]

[0114] Where r = 3.8, , The initial state is [0.21, 0.37, 0.49]. , and The trajectories are independent and follow a Gaussian distribution with a mean of 0 and a standard deviation of 0.01. At each step, the trajectories are clipped to... [10 -6 ,1-10 -6 The first 200 transitional samples are discarded, and the subsequent 300 samples are used. The true direct edges are... and .

[0115] The same 10 data seeds and algorithm parameters as in Example 1 were used. The mean precision, recall, F1, and accuracy of the frozen method were all 1.000, and the mean regression NRMSE was 0.033; the mean precision of the linear multivariate Granger baseline was 0.333, the mean recall was 1.000, the F1 was 0.500, and the mean accuracy was 0.333.

[0116] Example 3: A continuous process with four operating conditions under real topological changes:

[0117] To distinguish it from experiments that only change the dynamics but maintain the same causal graph, a four-condition process with 1600 consecutive points was constructed, with each condition containing 400 points. The physical switching time was... The initial state is [0.15, -0.10, 0.20]. Each channel is superimposed with independent Gaussian perturbations with a mean of 0 and a standard deviation of 0.05. The data seed is 240519. The state is not manually reset during switching.

[0118] The first working condition adopts Transmission chain:

[0119]

[0120]

[0121]

[0122] The second working condition uses a cross diagram. and :

[0123]

[0124]

[0125]

[0126] The third operating condition uses a reverse propagation chain. :

[0127]

[0128]

[0129]

[0130] The fourth working condition adopts Transmission chain:

[0131]

[0132]

[0133]

[0134] Initial identification resulted in a valid state (VALID) at t=181. The three structural changes occurred at... Candidates subject to change Figure 2 The detection was carried out by taking three votes, with detection delays of 127, 98, and 87 points respectively; respectively in The recovery latency for the VALID state, calculated from the physical handover time, is 380, 257, and 246 points, respectively, and the re-identification duration from the detection of recovery is 253, 159, and 159 points, respectively.

[0135] The average F1 score for all four operating conditions at the post-recovery graph update time is 1.000, and the F1 score for the last valid graph in each segment is also 1.000. The average value of the entire time trajectory also includes the unavoidable detection interval where the old valid graph is retained after the physical switch and before causal detection; therefore, this average value is not presented as the post-recovery recognition performance.

[0136] Example 4: Deployment of an industrial three-tank liquid level process:

[0137] In an industrial three-tank level process deployment, the industrial process sensor group 100 includes an inlet flow transmitter FT-101, a first tank level transmitter LT-102, and a second tank level transmitter LT-103, which respectively form variable... , and The data acquisition and synchronization interface 110 acquires 16-bit digital data at a sampling period of 1 second, with a timestamp alignment error of no more than 50 milliseconds, and sends the synchronization sample to the processor.

[0138] Window Management Module 120 Maintenance Length The predictive regression window must strictly adhere to the current moment, and its length. The system includes a causal analysis window for the current moment. The processor executes online state prediction, candidate edge evaluation model estimation, two-object EEMD, graph decision, and state control according to the aforementioned implementation parameters. In the VALID state, the result output interface 180 sends the valid directed adjacency matrix, valid flag 1, and cause code to the industrial fault propagation monitoring module. The industrial fault propagation monitoring module determines the directed connections corresponding to the matrix elements with a value of 1 in the valid directed adjacency matrix as direct fault propagation edges between industrial process variables. It then performs a directed traversal according to the direction of the direct fault propagation edges and the connection relationship between the first and last nodes, generating an industrial fault propagation path consisting of at least one direct fault propagation edge. Simultaneously, it outputs the industrial process variables traversed by the path and their propagation direction. In the IDENTIFY or REIDENTIFY state, the result output interface 180 only sends the valid flag 0 and cause code, without sending the change candidate graph. The industrial fault propagation monitoring module does not use the change candidate graph to generate an industrial fault propagation path for the current operating condition.

[0139] In this Windows 11 and Python 3.12.13 verification environment, the average update time of a single causal graph in the time-varying freeze experiment was 0.697 seconds, with a maximum of 1.469 seconds. Therefore, when the actual sampling period is greater than the maximum update time measured by this device, a graph update can be completed before the next sample arrives. For shorter sampling periods, higher variable dimensions, or different hardware, the graph update period should be increased, parallel computing should be used, or computing resources should be increased based on the measured time of the deployed device. The time consumption in this set is only used to illustrate the computational overhead of the freeze embodiment and does not constitute a general guarantee of real-time performance.

[0140] Example 5: Directed Effective Connectivity Analysis of Multichannel EEG:

[0141] In a directed and efficient connectivity analysis deployment for electroencephalography (EEG), a multi-channel signal acquisition device 100 is used for EEG signal acquisition, and all electrode channels employ a unified sampling clock. The data acquisition and synchronization interface 110 performs timestamp alignment, power frequency notch filtering, causal digital filtering, and artifact labeling on continuous EEG data; online processing uses only data arriving at the current time and before, excluding samples arriving after the current time. Each EEG electrode channel can be directly used as an observation variable. Alternatively, multiple brain region temporal signals can be generated first based on preset channel grouping or source space mapping, and then each brain region temporal signal can be used as an observation variable. .

[0142] Window management module 120 sets the prediction regression window length based on the sampling rate and the minimum sample requirement of the lag model. With causal analysis window length Both are preset positive integers and are not limited to 120 and 180. The processor performs the same online state prediction, independent estimation of the full model and reduced model for candidate edge evaluation, construction of conditional prediction residual difference signal, dual-object EEMD, multi-scale coherent aggregation, dual-threshold joint screening, batch screening of indirect paths, determination of bidirectional edge direction dominance and state control on the temporal signals of each EEG channel or brain region as in the aforementioned embodiments.

[0143] In the valid state (VALID), the matrix elements in the valid directed adjacency matrix that have a value of 1 are... This indicates that, under the current window and model conditions, historical information of EEG channel or brain region j has a direct, directed predictive effect on EEG channel or brain region i through joint criteria. The result output interface 180 sends the effective directed adjacency matrix and effective marker 1 to the EEG intention state analysis module. The EEG intention state analysis module determines the effective directed adjacency matrix as a directed effective connection network between EEG channels or brain regions, and calculates at least one of the following based on the directed effective connection network: in-degree, out-degree, edge duration, network path, and connection changes between preset brain regions. The calculation results are combined to form EEG connection topology features, which serve as input data for intention state analysis.

[0144] In the IDENTIFY or REIDENTIFY state, the result output interface 180 only sends the valid marker 0 and the reason code, and does not send the change candidate map; the EEG intention state analysis module does not use the change candidate map to generate EEG connection topology features, nor does it use the change candidate map for intention state analysis. This embodiment outputs the directed predicted connections and their validity markers obtained from EEG signal processing, and does not directly provide disease diagnosis, treatment plans, or anatomical connection conclusions; the training method and classification labels of the intention state classifier can be set separately by the downstream module without changing the core directed connection recognition processing chain of this invention.

Claims

1. A method for identifying causal connections in multivariable time-series signals, characterized in that: Executed by a processor that is communicatively connected to a multi-channel signal acquisition device, including: Obtain synchronized digital time-series samples of M observed variables, and maintain a length of [value missing] at the current evaluation time t. Predictive regression window and length is Causal analysis window Where M is an integer not less than 2. and The value is a preset positive integer; the prediction regression window only contains historical samples before time t, and the causal analysis window contains the current sample at time t and its historical samples before time t; At each sampling time t, an online state prediction model is established using the prediction regression window to obtain the one-step prediction value of each observed variable at time t, and the rolling average error and rolling error statistic are calculated based on the one-step prediction value. For any ordered pair of variables ,Will As a candidate causal variable, As the variable to be explained; at the candidate graph update time, the causal analysis window is used as... Establish a full candidate edge evaluation model incorporating multivariable nonlinear lag characteristics and obtain the full model residuals; delete all candidate edges containing... The nonlinear lag characteristics are used to independently re-estimate the candidate edges to evaluate the reduced model and obtain the residuals of the reduced model; the two types of residuals are aligned according to the same unexplained outcome variable and the same sampling time. The conditional prediction error ratio index is determined by the logarithm of the ratio of the variance of the reduced model residuals to the variance of the full model residuals, and the conditional prediction residual difference signal at each sampling time is formed by subtracting the full model residuals from the reduced model residuals. The candidate causal variable signal and the conditional prediction residual difference signal are subjected to ensemble empirical mode decomposition respectively. The two decompositions use reproducible auxiliary noise sequences generated according to the same deterministic foundation seed, candidate edge identifier and different signal type identifiers respectively. The two sets of intrinsic mode functions are aligned according to the order of the intrinsic mode functions, and mode pairs with the same order number are retained. Perform Hilbert transform on each mode pair and calculate the phase coherence value. Determine the mode weight based on the arithmetic mean of the variances on both sides of the mode pair. Weighted aggregate the phase coherence values ​​of each mode to obtain the original multi-scale causality score of the candidate edge. For the same unexplained outcome variable, the original multi-scale causal scores of each candidate edge are normalized, and a dual-threshold joint screening is performed when the normalized inbound edge weight and the conditional prediction error ratio both reach the corresponding thresholds; indirect path batch screening is performed on the fixed candidate directed graph obtained by the dual-threshold joint screening; bidirectional edge direction dominance is determined for candidate edges in opposite directions to form a candidate directed adjacency matrix. Under the IDENTIFY state, only if the rolling error statistic satisfies the rolling average error stability criterion, and the cumulative number of samples since the current identification starting point reaches the causal analysis window length. When the time is reached, the system enters the effective state VALID; in the effective state, the candidate directed adjacency matrix is ​​updated according to the preset graph update cycle, and the previous effective directed adjacency matrix is ​​maintained between adjacent update times; The three most recent candidate graphs that differ from the current valid graph are continuously saved, and a 2-out-of-3 voting process is performed on each candidate edge to form a voting graph. When the voting graph still differs from the current valid graph, the system enters the REIDENTIFY state and stops outputting candidate graphs. The cumulative number of samples since the change evidence time reaches [a certain value]. After the rolling average error stability criterion is satisfied again, the effective directed adjacency matrix output is restored. When the device is active, one of the following downstream processes is performed depending on the type of the multi-channel signal acquisition device: When the multi-channel signal acquisition device includes an industrial process sensor group, the effective directed adjacency matrix and the effective marker are output to the industrial fault propagation monitoring module. When the effective marker indicates an effective state, the industrial fault propagation monitoring module determines the directed connection between the candidate cause variable and the variable to be explained, corresponding to the matrix element with a value of 1 in the effective directed adjacency matrix, as the direct fault propagation edge between the industrial process variables, and generates an industrial fault propagation path according to the direction of the direct fault propagation edge and the connection relationship between the first and last nodes. When the multi-channel signal acquisition device includes an EEG signal acquisition device, the effective directed adjacency matrix and effective markers are output to the EEG intention state analysis module. When the effective markers indicate an effective state, the EEG intention state analysis module determines the effective directed adjacency matrix as a directed effective connection network between EEG channels or brain regions, and calculates at least one of the following based on the directed effective connection network: in-degree, out-degree, edge duration, network path, and connection changes between preset brain regions, to generate EEG connection topology features, and uses the EEG connection topology features as input data for intention state analysis. In the identification state and the re-identification state, an invalid flag is output to the corresponding industrial fault propagation monitoring module or EEG intention state analysis module, and the industrial fault propagation monitoring module or EEG intention state analysis module is prohibited from using the changed candidate graph for fault propagation analysis or intention state analysis.

2. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: The nonlinear lag features include: linear terms of each variable lagged from the 1st to the pth order, univariate polynomial terms of each variable lagged from the 2nd to the Qth order, and pairwise product terms between different variables of the same lag order, where p and Q are both positive integers and Q is not less than 2; each feature records the set of variables it contains, and the candidate edge evaluation reduction model deletes all features in the variable set that contain candidate causal variables; the features of the candidate edge evaluation full model and the candidate edge evaluation reduction model are standardized with the mean and standard deviation of the corresponding features within the causal analysis window, and are independently estimated by ridge regression with the intercept term unpenalized.

3. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: Candidate causal variables Pointing to the variable of the result to be explained Conditional prediction error ratio index Conditional prediction residual difference signal Determine them respectively using the following formulas: ； ； in, This is the current causal analysis window. From The time-aligned index set obtained after subtracting the initial samples required to form the first p-order hysteresis features. For the variable to be explained The candidate edges are evaluated to assess the residuals of the entire model. To remove candidate cause variables The candidate edge evaluation reduction model residuals are obtained by independently estimating relevant features. It is a positive number; the conditional prediction error ratio index reflects the deletion. The residual variance after relevant features evaluates the change in the full model residual variance relative to the candidate edges, and the conditional prediction residual difference signal has the same characteristics as... Same physical dimensions.

4. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: The deterministic basic seed is Candidate edges First pseudo-random noise seed Second pseudo-random noise seed Generate using the following formula: ； Where H is a deterministic hash map, j and i identify candidate edges, and numbers 1 and 2 identify candidate causal variable signals and conditional prediction residual difference signals, respectively; two seeds initialize the corresponding ensemble empirical mode decomposition processes and generate reproducible auxiliary noise sequences; the two sets of intrinsic mode functions are arranged according to the decomposition output order, when their quantities are respectively and At that time, only the first one is retained. For a pair of modes of the same order, center frequency matching is not performed.

5. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: Phase coherence value of the m-th mode pair Mean modal variance Modal weights Original multi-scale causal scores of candidate edges Determine them respectively using the following formulas: ； ； ； ； Where N is the number of samples participating in coherent calculation after hysteresis alignment in the current causal analysis window. and Let be the m-th eigenmode functions of the candidate causal variable signal and the conditional prediction residual difference signal, respectively. and For the corresponding instantaneous phase, K is the number of aligned mode pairs. For modal summation index, Let abs be a positive number, representing the complex modal function; the sum of all mean modal variances is... When the sum of the variances is not greater than At that time, set all modal weights and the original multiscale causality score to zero; Located in a closed interval .

6. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: Candidate edges Normalized incoming edge weights Determine using the following formula: ； Here, q iterates through all candidate cause variables except i. Candidate edges The original multiscale causal fraction, It is a positive number; when the denominator is not greater than season ,and Located in a closed interval ;when Not less than the incoming edge threshold and Not less than the error ratio threshold When candidate edges are selected, they are retained, thus forming a fixed candidate directed graph after dual threshold joint screening; normalized inbound edge weights are only used for comparing candidate sources for the same unexplained result variable.

7. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: The fixed candidate directed graph obtained after dual threshold joint screening consists of nodes from j to i, with lengths from 2 to M-1, where nodes are unique and do not contain direct edges to be compared. The set of all indirect directed paths; paths Path bottleneck score and candidate edges Maximum indirect explanation score Determine them respectively using the following formulas: ； ； in, To form a path The minimum value among the original multi-scale causal scores of each candidate edge, i.e., the path bottleneck score; The maximum value among all eligible indirect paths' bottleneck scores, i.e., the maximum indirect explanation score; when Existence and Not greater than At that time, Add to the deletion set; after completing the determination of all paths, delete all candidate edges in the deletion set at once, completing the batch screening of indirect paths and making the result independent of the traversal order; then perform the bidirectional edge direction dominance determination: when and When they exist simultaneously, if Not less than Then only retain ,like Not less than Then only retain Otherwise, mark the variable pair as direction undecided and delete both directions, where A positive indirect path suppression coefficient. The directional dominance coefficient is greater than 1.

8. The method for identifying causal connections of multivariable time-series signals according to claim 1, characterized in that: At time t, the predictive regression window is used. The established online state prediction model obtains one-step prediction values. And utilize the scale of each variable within the predicted regression window. Normalize the one-step prediction error; instantaneous multivariate normalized error. and rolling average error of length L Determine using the following formula: ； Let the most recent H rolling average errors form the set of rolling error statistics. The rolling average error stability criterion is determined by the following formula: ； in, The lower limit of the positive scale. As the starting point for current identification, The length of the causal analysis window; during the initial identification period. As the initial identification starting point, during the re-identification period For moments of change as evidence; when The error alarm threshold is not less than D consecutive time intervals. The system immediately performs an internal graph check, but does not declare the causal graph invalid based solely on this check; it performs a two-out-of-three voting on the three most recent candidate graphs that are different from the current valid graph, and includes the same candidate edge in the voting graph if the same candidate edge appears at least twice; when the voting graph is still different from the current valid graph, it enters a re-identification state, and no candidate graph is output during the re-identification period.

9. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: The multi-channel signal acquisition device includes an industrial process sensor group. The observed variables are at least two industrial process variables selected from flow rate, pressure, liquid level, temperature, current, vibration, and control quantity. In the valid state, the industrial fault propagation monitoring module receives the valid directed adjacency matrix and valid markers, determines the directed connections corresponding to the matrix elements with a value of 1 in the valid directed adjacency matrix as direct fault propagation edges between industrial process variables, and performs directed traversal according to the direction of the direct fault propagation edges and the connection relationship between the first and last nodes to generate an industrial fault propagation path consisting of at least one direct fault propagation edge. The industrial fault propagation path includes the industrial process variables traversed by the path and their propagation direction. In the identification state and re-identification state, the industrial fault propagation monitoring module receives invalid markers and prohibits the generation of industrial fault propagation paths using change candidate graphs.

10. The method for identifying causal connections in multivariable time-series signals according to claim 1, characterized in that: The multi-channel signal acquisition device includes an electroencephalogram (EEG) signal acquisition device, and the observed variable is the EEG electrode channel timing signal, or the brain region timing signal formed by channel grouping or source spatial mapping of multiple EEG electrode channels; In the effective state, the EEG intention state analysis module receives the effective directed adjacency matrix and the effective marker, determines the effective directed adjacency matrix as a directed effective connection network between EEG channels or brain regions, and calculates at least one of the following based on the directed effective connection network: in-degree, out-degree, edge duration, network path, and connection changes between preset brain regions, generates EEG connection topology features, and uses the EEG connection topology features as input data for intention state analysis; In the identification state and the re-identification state, the EEG intention state analysis module receives invalid tags and prohibits the generation of EEG connection topology features or the performance of intention state analysis using the changed candidate graph.

11. A causal connection identification system for multivariable time-series signals, characterized in that: The system includes a multi-channel signal acquisition device 100, a data acquisition and synchronization interface 110, a processor, a memory, and a result output interface 180. The processor functionally includes a window management module 120, a state prediction and full model and reduced model identification module 130, a conditional prediction calculation module 140, an ensemble empirical mode decomposition and multi-scale coherent aggregation module 150, a graph decision module 160, and a convergence decision and state control module 170. The memory stores a program, deterministic foundation seeds, model parameters, thresholds, and a state log. When the program is executed by the processor, it causes the processor to implement the method described in any one of claims 1 to 10. The result output interface 180 outputs a result data packet including a timestamp, algorithm version, status, valid flag, and reason code. When the valid flag is 1, the result data packet also includes a valid directed adjacency matrix. When the multi-channel signal acquisition device includes an industrial process sensor group, the valid directed adjacency matrix is ​​provided to the industrial fault propagation monitoring module to generate an industrial fault propagation path. When the multi-channel signal acquisition device includes an EEG signal acquisition device, the valid directed adjacency matrix is ​​provided as a directed valid connection network between EEG channels or brain regions to the EEG intention state analysis module to generate EEG connection topology features. When the valid flag is 0, no change candidate graph is provided, and the corresponding downstream analysis module is prohibited from using the change candidate graph for fault propagation analysis or intention state analysis.

12. A storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method described in any one of claims 1 to 10.

Images