A metro SoS anomaly detection method based on physical-data collaborative diffusion model

Abstract

The application discloses a subway SoS anomaly detection method based on a physical-data collaborative diffusion model. The method comprises the following steps: collecting operation data of vehicle, signal, traction power supply and passenger flow subsystems, and obtaining multi-modal time window samples through time alignment and feature processing; constructing a physical constraint set according to the operation mechanism of each subsystem, and jointly training a diffusion model based on data fidelity loss and physical consistency loss; performing deterministic single-step reconstruction on the to-be-detected sample to obtain a reconstructed anomaly score, and simultaneously calculating a physical violation score; standardizing the two types of scores and performing scene adaptive fusion, and determining whether the subway SoS is abnormal according to a dynamic threshold; and when an anomaly occurs, the fault subsystem, equipment and violation rules are located through physical residual decomposition, reconstruction error feature decomposition and causal knowledge graph. The application improves the real-time performance, accuracy and interpretability of cross-subsystem anomaly detection.

Description

Technical Field

[0001] This invention belongs to the field of subway operation and maintenance technology, specifically involving a subway SoS anomaly detection method based on a physical-data collaborative diffusion model. Background Technology

[0002] Subway operation is a typical System of Systems (SoS), encompassing multiple subsystems such as rolling stock, signaling, power supply, and passenger flow. These subsystems are deeply coupled, and even a minor disturbance in any component can trigger a chain reaction through inter-system interfaces. Currently, anomaly detection technologies in this field mainly include single-system physical rule methods based on fixed thresholds and purely data-driven machine learning methods. However, both methods have limitations and struggle to achieve collaborative perception and intelligent diagnosis across the entire subway system.

[0003] Single-system physical rule-based methods using fixed thresholds are typically implemented by various professional maintenance teams, who set fixed alarm thresholds or simple logic rules in their respective system's programmable logic controllers or monitoring servers, based on design specifications and operational experience. Since subsystems such as vehicles, signals, power supply, and passenger flow are generally monitored separately by different professionals, the detection rules for each subsystem primarily judge its own operating parameters, making it difficult to correlate monitoring information from different subsystems. When anomalies are triggered by changes in the coupling relationships between multiple subsystems, and the individual parameters of each subsystem have not yet exceeded their corresponding thresholds, such anomalies are easily missed.

[0004] Setting fixed thresholds presents a trade-off between sensitivity and false alarm rate. Relaxing the threshold to reduce false alarms results in insufficient sensitivity to early anomalies such as slow sensor drift and gradual wear of mechanical components, easily leading to missed alarms. Conversely, tightening the threshold to capture early anomalies can cause frequent alarm triggers under normal but fluctuating operating conditions, such as peak hours, severe weather, or schedule adjustments, generating numerous invalid alarms. Furthermore, alarm information generated by the fixed threshold method typically only indicates that a parameter exceeds the preset range, failing to explain the cause of the anomaly, the associated subsystems, and the scope of impact. Complex faults still require manual backtracking and correlation analysis by multiple technical personnel.

[0005] To compensate for the shortcomings of physical rule-based methods, data-driven methods, represented by deep learning, have been introduced into the field of subway operation and maintenance in recent years. Typical methods include long short-term memory autoencoders and single-class support vector machines.

[0006] Long Short-Term Memory (LSTM) autoencoders are typically trained using historical normal data and determine the presence of anomalies based on the error between the input data and the reconstructed data. This method primarily learns the statistical correlations in the training data, making it difficult to directly represent the physical causal relationships during equipment operation. For infrequent but physically normal operating conditions, the model may misclassify them as anomalies; when unidentified fault modes are mixed into the training data, the model may also learn these fault modes as normal modes. Due to the lack of explicit physical constraints within the model, its reconstruction results may not satisfy the dynamics, safety logic, or energy conservation relationships in actual operation.

[0007] Single-class support vector machines typically construct decision boundaries for normal samples in the feature space and classify data outside these boundaries as anomalies. The detection performance of this method largely depends on manually designed features. For subway SoS (System-on-Site) operations, which contain various types of data such as continuous physical quantities, discrete signal states, and passenger flow statistics, and exhibit strong temporal and dynamic coupling, manually designed features are insufficient to cover all operational modes. When the operating schedule, new vehicle models, seasons, or equipment loads change, the distribution of normal data may drift, and the original decision boundaries may easily become invalid.

[0008] Existing purely data-driven methods generally suffer from insufficient physical consistency and limited anomaly interpretation capabilities. The models typically lack explicit constraints such as train dynamics, signal interlocking logic, traction power supply energy conservation, and passenger flow load consistency, leading to anomaly detection results that may not align with actual system operation patterns. Furthermore, increased reconstruction errors or decision function out-of-bounds errors usually only indicate that the current sample deviates from the normal data distribution; they cannot directly identify the corresponding subsystem, equipment, monitoring variables, or the disrupted physical or safety rules.

[0009] The prior art CN121980476A discloses a multivariate temporal anomaly detection method based on stage-aware transfer diffusion across working conditions. It extracts spatiotemporal features through temporal local normalization, dynamic graph construction, temporal convolutional network and graph attention network, and uses graph variational autoencoder, stage encoder and stage conditional diffusion model for anomaly detection to adapt to the changes in data distribution between different working conditions.

[0010] However, this existing technology mainly relies on the statistical characteristics of sensors to construct dynamic graph relationships and uses latent variables, stage probabilities, normal prototypes, and staged thresholds to complete anomaly judgments. Its detection basis still mainly comes from the statistical correlation between data. For heterogeneous subsystems such as subway vehicles, signaling, traction power supply, and passenger flow, it is difficult to directly determine whether the monitoring data violates specific operational rules such as train force balance, signal interlocking, movement authorization protection, traction power supply energy conservation, or load consistency.

[0011] Furthermore, CN121980476A primarily obtains latent variable reconstruction results through diffusion backsampling and combines reconstruction error, reconstruction uncertainty, and phased thresholds for anomaly detection. This approach focuses on data distribution adaptation across operating conditions, making it difficult to simultaneously utilize data distribution deviation information and physical state deviation information for mutual verification. After outputting anomaly results, it is also difficult to further determine the subsystem, equipment, characteristic variables, and disrupted physical rules involved in the anomaly, and its ability to locate the root cause and analyze the correlation path of anomalies across subsystems remains limited.

[0012] In summary, existing technologies still have the following technical problems: a lack of unified representation and cross-subsystem correlation mechanisms for multi-source heterogeneous data applicable to metro SoS; anomaly detection mainly relies on fixed thresholds or data statistical patterns, making it difficult to guarantee the physical reliability of anomaly judgment results; anomaly detection results under different operating scenarios are easily affected by operating condition switching, load fluctuations, and environmental changes; anomaly results usually remain at the level of parameter out-of-bounds, reconstruction errors, or anomaly scores, making it difficult to automatically locate the subsystem, equipment, and characteristic variables where anomalies occur, and also difficult to determine the disrupted physical rules and anomaly correlations. Summary of the Invention

[0013] The purpose of this invention is to overcome the shortcomings of the existing technology and provide a subway SoS anomaly detection method based on the physical-data co-diffusion model.

[0014] The objective of this invention can be achieved through the following technical solutions: This invention provides a subway SoS anomaly detection method based on a physical-data co-diffusion model, comprising the following steps: Operational data from multiple subsystems in the subway SoS are collected, and time alignment and feature processing are performed on the operational data to obtain multimodal time window samples. A set of physical constraints is constructed based on the operating mechanisms of the multiple subsystems. The set of physical constraints is used to characterize the physical consistency relationships between the multiple subsystems themselves and among the multiple subsystems. Using the multimodal time window samples under normal operating conditions as training samples, the diffusion model is jointly trained based on the data fidelity loss of the training samples and the physical consistency loss determined according to the physical constraint set, to obtain the physical-data co-diffusion model. The multimodal time window sample to be detected is input into the physical-data co-diffusion model for reconstruction. Based on the difference between the multimodal time window sample to be detected and the corresponding reconstructed sample, the reconstruction anomaly score is determined. Physical consistency detection is performed on the multimodal time window samples to be detected based on the physical constraint set to determine the physical violation score; By combining the reconstructed anomaly score and the physical violation score, a comprehensive anomaly score is obtained. Based on the comparison result between the comprehensive anomaly score and the anomaly judgment threshold, it is determined whether the subway SoS has an anomaly.

[0015] Furthermore, the process of collecting operational data from multiple subsystems within the metro SoS, performing time alignment and feature processing on the operational data to obtain multimodal time window samples, specifically includes: Operational data with sampling timestamps are collected from the vehicle subsystem, signaling subsystem, traction power supply system, and passenger flow subsystem, respectively. The running data is timestamped and missing values ​​are processed according to a unified timeline to obtain aligned running data corresponding to each sampling time. According to the data mode to which the running data belongs, the continuous variables in the aligned running data are standardized, and the discrete states in the aligned running data are one-hot encoded and scaled. The processed running data of each subsystem corresponding to the same sampling time are spliced ​​together to obtain the multimodal feature vector corresponding to that sampling time. The multimodal feature vectors of a preset number of consecutive samples are stacked in chronological order to obtain the multimodal time window sample.

[0016] Furthermore, the construction of the physical constraint set based on the operating mechanisms of the multiple subsystems specifically includes: Based on the mechanical relationship between the train's traction, braking, and running resistance, longitudinal dynamic constraints of the train are constructed. Based on the legal transfer relationships of subway signal states and the security relationship between train speed and movement authorization protection range, construct signal security logic constraints; Based on the energy conservation relationship between the injected power of the traction power supply network and the net power consumption of online trains, traction power flow constraints are constructed. Based on the consistency relationship between the train load estimated from passenger flow data and the train load deduced from onboard air spring pressure data, load consistency constraints are constructed. The train longitudinal dynamics constraints, the signal safety logic constraints, the traction power flow constraints, and the load consistency constraints are respectively converted into differentiable residual functions to form the physical constraint set. For segmented relationships containing maximum values, minimum values, or state transitions, corresponding smoothing processing is adopted during the joint training phase to ensure that the gradient can be propagated. Each residual function tends to zero when the corresponding physical consistency relationship is satisfied and increases when the corresponding physical consistency relationship is violated.

[0017] Furthermore, the construction of the train longitudinal dynamic constraints and the signal safety logic constraints specifically includes: The train acceleration is determined based on the train speed at adjacent sampling times, the basic running resistance of the train is determined based on the train speed and preset running resistance parameters, and the longitudinal dynamic residual of the train is determined based on the force balance relationship between the train traction force, train braking force, train mass, the train acceleration, the basic running resistance of the train and the track gradient resistance. The signal state at each sampling time is encoded into a one-hot vector. Based on the preset legal signal state transition matrix and the signal state at the current sampling time, the predicted signal state distribution at the next sampling time is determined. Based on the difference between the predicted signal state distribution and the actual signal state at the next sampling time, the interlocking transfer residual is determined. Based on the train's current position, the endpoint of the movement authorization, the preset emergency braking deceleration, the system response travel distance, and the permissible speed on the line, the limit safe speed corresponding to the train's current position is determined, and the movement authorization protection residual is determined based on the degree to which the train's actual speed exceeds the limit safe speed. The interlocking transfer residual and the mobile authorization protection residual are weighted and combined to obtain the signal security logic residual.

[0018] Furthermore, the determination of the traction power flow constraint, the load consistency constraint, and the physical violation score specifically includes: The total injected power on the feeder side into the traction power supply network is obtained, and the total net power consumption of the traction power supply network is determined based on the traction or regenerative braking power of each online train, the power of other loads, the power loss of the line and converter, and the net charging power of the energy storage device. The traction power flow residual is determined based on the difference between the total injected power on the feeder side and the total net power consumption. The load of the first train is determined based on the passenger flow statistics of the automatic fare collection system, and the load of the second train is determined based on the air spring pressure. The load consistency residual is determined based on the difference between the load of the first train and the load of the second train. The residuals in the physical constraint set are dimensionless by using the mean and standard deviation of the normal historical residuals respectively. Then, the dimensionless residuals are processed by Huber loss. The processing results are weighted according to the safety criticality level of each physical constraint to obtain the physical violation energy corresponding to the current multimodal time window sample. The Huber threshold is a threshold that is preset on the dimensionless scale or calibrated using the validation set. According to a preset forgetting factor, the physical violation energy corresponding to multiple consecutive multimodal time window samples is subjected to an exponential moving average to obtain the physical violation score.

[0019] Furthermore, the diffusion model includes a forward diffusion module, a time-step encoding module, a denoising network, and a clean sample estimation module; The denoising network adopts a one-dimensional temporal U-Net structure, including an encoder, a temporal self-attention bottleneck layer, a decoder, and skip connections connecting network layers of the same scale in the encoder and the decoder. The encoder and the decoder each include multiple one-dimensional convolutional residual blocks. The temporal step encoding module includes a sinusoidal encoding layer and a feature-level linear modulation layer. The feature-level linear modulation layer is connected to each of the one-dimensional convolutional residual blocks. The training samples are input into the forward diffusion module, and Gaussian noise is added to the training samples at selected diffusion time steps according to the preset noise schedule to obtain noisy samples. The diffusion time step is input into the sinusoidal coding layer to obtain the time step code, and the feature-level linear modulation layer modulates the features of each one-dimensional convolutional residual block according to the time step code; The noisy sample is input into the encoder for step-by-step downsampling and feature extraction. The features output by the encoder are input into the temporal self-attention bottleneck layer to extract the correlation features of the multimodal time window sample in the temporal dimension. Then, the decoder performs step-by-step upsampling and fuses the features of the corresponding scale in the encoder through the skip connection to output the predicted noise corresponding to the noisy sample. The clean sample estimation module determines the clean sample estimate corresponding to the training sample based on the noisy sample, the noise scheduling coefficient corresponding to the diffusion time step, and the prediction noise. The prediction noise is used to determine the data fidelity loss, and the clean sample estimate is used to determine the physical consistency loss based on the physical constraint set.

[0020] Furthermore, the joint training of the diffusion model based on the data fidelity loss of the training samples and the physical consistency loss determined according to the physical constraint set specifically includes: The data fidelity loss is determined based on the difference between the Gaussian noise added to the training samples by the forward diffusion module and the predicted noise output by the denoising network. The data fidelity loss is expressed as: in, This indicates the loss of data fidelity; This represents the expectation operation; This represents the multimodal time window sample used as the training sample; This indicates the addition of Gaussian noise to the training samples; Indicates the selected diffusion time step; This represents the noisy sample obtained at the diffusion time step; The parameter is The denoising network outputs predicted noise based on the noisy sample and the diffusion time step; Represents the square of the L2 norm; Based on the physical constraints in the physical constraint set, physical residuals are calculated for the clean sample estimates. Huber loss is applied to each physical residual, and the results are weighted according to the corresponding constraint weights to determine the physical consistency loss. The physical consistency loss is expressed as follows: in, This represents the loss of physical consistency; This indicates the number of physical constraints in the physical constraint set; Indicates the first The constraint weights corresponding to each physical constraint; This indicates the estimated value based on the clean sample. The calculated first One physical residual; Indicates the threshold is Huber loss function; Indicates the first The Huber loss threshold corresponding to each physical residual; The total loss for joint training is determined based on the data fidelity loss and the physical consistency loss, and the weight of the physical consistency loss is adjusted according to the training progress and the noise level corresponding to the diffusion time step. The total loss and the weight of the physical consistency loss are expressed as follows: in, This represents the total loss; Indicates the first The training rounds and the first The weights of the physical consistency loss at each diffusion time step; The underlying weights representing the loss of physical consistency; This represents the warm-up coefficient, which increases linearly with the number of training rounds. The preset noise scheduling indicates that the preset noise scheduling is in the first... The cumulative noise retention factor corresponding to each diffusion time step; During backpropagation based on the total loss, the physical rule parameters in the physical constraint set are fixed, the gradient of the physical rule parameters is not updated, and the gradient of the physical consistency loss is backpropagated to the denoising network through the clean sample estimation module and the calculation path of the prediction noise to update the parameters of the denoising network.

[0021] Furthermore, the step of inputting the multimodal time window sample to be detected into the physical-data co-diffusion model for reconstruction, and determining the reconstruction anomaly score based on the difference between the multimodal time window sample to be detected and the corresponding reconstructed sample, specifically includes: According to the preset noise evaluation step, the multimodal time window sample to be detected is noise-added according to the preset reproducible noise generation rule to obtain the noisy time window sample to be detected. The noisy time window sample to be detected is input into the denoising network to obtain the predicted noise corresponding to the noisy time window sample to be detected. Based on the noise prediction parameterization relationship, the reconstructed sample is generated through a single network forward computation based on the noisy time window sample to be detected and the predicted noise. The reconstruction anomaly score is determined based on the mean squared error between the multimodal time window sample to be detected and the reconstructed sample; While determining the reconstruction anomaly score, the physical consistency of the multimodal time window sample to be detected is performed using the differentiable residual functions in the physical constraint set, and the physical violation score is determined based on the calculation results of each differentiable residual function.

[0022] Furthermore, the process of fusing the reconstructed anomaly score and the physical violation score to obtain a comprehensive anomaly score, and determining whether the subway SoS has anomalies based on the comparison result of the comprehensive anomaly score and the anomaly judgment threshold, specifically includes: A sliding window for reconstructed abnormal scores and a sliding window for physical violations, containing only historical scores under normal operating conditions, are established respectively. The corresponding score sliding windows are initialized using the reconstructed abnormal scores and physical violation scores corresponding to historical multimodal time window samples that have been confirmed to be under normal operating conditions. Based on the mean and standard deviation of each historical score within the corresponding score sliding window, the current reconstructed abnormal score and the physical violation score are standardized to obtain the data abnormality deviation and physical abnormality deviation: The fusion weight is determined based on the current operational scenario, and the data anomaly deviation and the physical anomaly deviation are fused according to the fusion weight to obtain the comprehensive anomaly score: in, This represents the current reconstruction anomaly score; and These represent the mean and standard deviation of the sliding window of scores corresponding to the reconstructed outlier scores, respectively. This represents the current physical violation score; and These represent the mean and standard deviation of the sliding window of scores corresponding to physical violation scores, respectively. This indicates the degree of abnormal deviation of the data; Indicates the degree of deviation of the physical anomaly; This represents the overall anomaly score; Indicates the fusion weights; It is a small constant that is numerically stable; The current operating scenarios include off-peak operating scenarios, peak operating scenarios, operating scenarios during severe weather, and low-density operating scenarios at night; Establish a sliding window for comprehensive anomaly scores for different operating scenarios, and determine the preset quantile of the historical comprehensive anomaly scores in the sliding window corresponding to the current operating scenario as the anomaly judgment threshold. When the overall anomaly score is greater than the anomaly determination threshold, it is determined that the subway SoS has an anomaly; when the overall anomaly score is not greater than the anomaly determination threshold, it is determined that the subway SoS does not have an anomaly, and the current reconstruction anomaly score, physical violation score and overall anomaly score are updated to the reconstruction anomaly score sliding window, the physical violation score sliding window and the overall anomaly score sliding window corresponding to the current operation scenario, respectively.

[0023] Furthermore, after determining that the subway SoS is abnormal, the following steps are also included: Based on the physical residuals corresponding to each physical constraint in the physical constraint set, determine the contribution of each physical constraint to the physical violation score, and identify the physical constraint with the highest contribution as the violated physical constraint. The reconstruction error between the multimodal time window sample to be detected and the reconstructed sample is decomposed along the feature dimension to determine the contribution of each feature variable to the reconstruction anomaly score. One or more feature variables with the highest contribution are identified as anomalous feature variables, and the subsystem and monitoring equipment corresponding to the anomalous feature variables are identified. The violated physical constraints and the abnormal feature variables are used as abnormal evidence and input into a preset causal knowledge graph. The causal knowledge graph uses subsystems, equipment, physical constraints and observations in the subway SoS as nodes, and the relationship between each node is represented by preset association edges. Based on the anomalous evidence, a correlation path matching the violated physical constraint and the anomalous feature variable is retrieved in the causal knowledge graph. Based on the correlation path, the faulty subsystem, faulty equipment, and anomalous causal chain corresponding to the subway SoS anomaly are determined, and anomaly diagnostic information containing the violated physical constraint, the anomalous feature variable, and the anomalous causal chain is generated.

[0024] Compared with the prior art, the present invention has the following advantages: (1) Existing subway anomaly detection typically processes data separately according to specialties such as vehicle, signal, traction power supply, and passenger flow. The data from different subsystems differ in sampling time, data type, and dimension, which can easily lead to the fragmentation of cross-subsystem operational relationships. This invention performs time alignment and feature processing on the operational data of each subsystem, concatenates data from the same sampling time into a multimodal feature vector, and stacks multiple consecutive multimodal feature vectors in chronological order into a multimodal time window sample. This allows the horizontal coupling relationship between different subsystems and the vertical evolution relationship of system state to be preserved simultaneously in the same data primitive, reducing false correlations caused by asynchronous sampling and dimensional differences, and improving the ability to identify coupling anomalies where individual parameters have not exceeded limits but cross-subsystem collaborative relationships have changed.

[0025] (2) Existing data-driven anomaly detection models mainly learn statistical correlations in historical data. They are prone to false alarms for statistically rare but physically normal operating states, while they may miss anomalies that have not changed much in value but have disrupted the system's operating rules. This invention constructs physical constraints such as train longitudinal dynamics, signal safety logic, traction power flow, and load consistency based on the operating mechanisms of multiple subsystems. It also uses the physical consistency loss and data fidelity loss together for diffusion model training, so that the model is constrained by the actual operating rules while learning the joint distribution of normal data.

[0026] Meanwhile, this invention adjusts the weights of the physical consistency loss according to the training process and the level of diffused noise, and performs gradient isolation only on updates to physical rule parameters and statistics, ensuring that the physical consistency loss is backpropagated to the denoising network via clean sample estimates. This reduces the interference of unreliable clean sample estimates on training during high-noise phases, while simultaneously making the model output in low-noise phases approach the physical feasible region. A comparative study of the same structure in Example 2 shows that, compared to M4, M5 using the aforementioned combined technical features improves precision by 0.0269, F1 score by 0.0091, accuracy by 0.0020, and specificity by 0.0028, while reducing the average false positive rate by approximately 70.37%, thus providing quantitative evidence for the comprehensive technical effect produced by the combined technical features.

[0027] (3) Existing time-series models struggle to simultaneously extract short-term local changes, coupling relationships between multiple data channels, and state-response relationships spanning longer periods. This invention employs a one-dimensional temporal U-Net as the denoising network. It extracts local features at different time scales through one-dimensional convolutional residual blocks, extracts long-range correlation features through a time-dimensional self-attention bottleneck layer, and preserves low-level details through skip connections. Simultaneously, it utilizes time-step encoding and feature linear modulation to adjust each network layer. This approach can take into account transient fluctuations, cross-channel collaboration, and time delay relationships, reduce the loss of weak anomalous features during the step-by-step downsampling process, and improve the modeling capability for complex time-series operation modes of multiple subway subsystems.

[0028] (4) Existing diffusion models typically require multiple reverse denoising operations, resulting in significant computational overhead, and random sampling may lead to highly fluctuating reconstruction results for the same operating state. This invention utilizes a fixed global seed and time window number to generate reproducible evaluation noise, and employs a deterministic single-step denoising process with noise prediction parameterization to generate reconstruction samples. Compared to multi-step reverse sampling, this approach reduces the number of online iterations; the same evaluation noise is used when inputting repeatedly within the same time window, thus enabling the acquisition of reproducible reconstruction anomaly scores.

[0029] (5) Existing technologies typically rely solely on reconstruction error or physical thresholds for anomaly detection. Using only reconstruction error is susceptible to interference from changes in normal operating conditions, while using only physical rules may miss combined anomalies that haven't yet significantly disrupted individual physical relationships. Furthermore, the two types of scores have different dimensions and numerical ranges, making direct fusion difficult. This invention determines reconstruction anomaly scores and physical violation scores in parallel for the same multimodal time window sample to be detected. These scores are standardized using historical scores from normal operating conditions, and then fused based on the current operating scenario. This allows data distribution deviation information and physical state deviation information to be mutually verified on a unified scale: when data distribution changes due to factors such as passenger flow and weather, but the system still conforms to operating rules, false alarms caused by statistical fluctuations can be reduced; when the joint relationship of multiple variables is abnormal but the individual physical residual is not yet significant, the reconstruction anomaly score can provide supplementary judgment criteria.

[0030] Furthermore, this invention sets corresponding fusion weights for different operating scenarios such as off-peak, peak, severe weather, and low-density nighttime operation, and maintains the statistical baseline of the comprehensive anomaly score and the dynamic anomaly judgment threshold respectively, so that the anomaly judgment standard can adapt to the changes in the normal fluctuation range under different scenarios, thereby taking into account both the false alarm suppression capability under high fluctuation conditions and the detection sensitivity of weak anomalies under stable conditions.

[0031] (6) To address the issues of instantaneous noise easily causing false alarms and simple averaging easily weakening asymptotic anomalies, this invention employs Huber loss to process each physical residual and performs exponential moving average on the physical violation energy corresponding to continuous multimodal time window samples. This can suppress the excessive influence of occasional outliers on the detection results, while accumulating continuously occurring small deviations, thereby improving the early identification capability of slowly changing anomalies such as sensor drift, braking force attenuation, and power supply performance degradation.

[0032] (7) Existing anomaly detection systems typically only output anomaly labels or scores, failing to explain the source of the anomaly or the basis for alarms. This invention, after determining that anomalies exist in the subway's SoS (Solar System), decomposes the physical violation energy corresponding to each physical constraint to identify the violated physical or safety rules; simultaneously, it decomposes the error between the multimodal time window samples to be detected and the reconstructed samples along the feature dimension to identify anomaly feature variables and their corresponding subsystems and monitoring equipment. Both types of anomaly evidence are then input into a causal knowledge graph to retrieve matching association paths. This allows for the combination of which monitoring data deviate from normal conditions and which system operating rules are disrupted, outputting faulty subsystems, faulty equipment, anomaly variables, and anomaly causal chains. Attached Figure Description

[0033] Figure 1 This is a flowchart of the subway SoS anomaly detection method according to an embodiment of the present invention; Figure 2 This is a flowchart illustrating the training process of the physical-data co-diffusion model according to an embodiment of the present invention. Figure 3 This is a flowchart illustrating the online anomaly scoring and judgment process according to an embodiment of the present invention. Figure 4 This is a flowchart of the two-layer decomposition and anomaly tracing process according to an embodiment of the present invention; Figure 5 This is a schematic diagram of the hardware deployment structure according to an embodiment of the present invention; Figure 6 This is a schematic diagram illustrating the fixed-point detection performance of various comparison methods in embodiments of the present invention; Figure 7 This is a comparison chart of the average false positive rates of various comparison methods in the embodiments of the present invention across five time-series test blocks. Figure 8 This is a matrix diagram of the M5 classification results in an embodiment of the present invention; Figure 9 This is a comparison chart of the F1 scores of various comparison methods in this embodiment of the invention on five time-sequence test blocks. Detailed Implementation

[0034] 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, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0035] Example 1: like Figure 1 As shown, this embodiment provides a subway SoS anomaly detection method based on a physical-data collaborative diffusion model. This method considers the vehicle subsystem, signaling subsystem, traction power supply system, and passenger flow subsystem of a subway line as interconnected system components. It forms multimodal time window samples using a unified time reference, and the same multimodal time window samples are processed in parallel by a physical constraint calculation branch and a diffusion model reconstruction branch. Then, the two anomaly scores are fused to complete anomaly determination. When an anomaly is determined, the physical violation energy and reconstruction error are further decomposed to determine the physical constraints, feature variables, subsystems, and monitoring equipment involved in the anomaly. Thus, a complete processing chain is formed from data acquisition, physical consistency calculation, diffusion reconstruction, score fusion to anomaly localization.

[0036] In this embodiment, the detection system for the execution method includes a multi-source real-time data acquisition module, a multi-modal preprocessing and alignment module, a physical rule base and constraint calculation module, a diffusion model anomaly scoring module, a comprehensive anomaly score fusion and threshold determination module, and an anomaly location and interpretation module. The multi-source real-time data acquisition module is connected to the train control network, the signal system interface, the traction power supply monitoring interface, the automatic fare collection system interface, and the on-board suspension monitoring interface, respectively. The physical rule base and constraint calculation module and the diffusion model anomaly scoring module receive the same multi-modal time window samples to avoid inconsistencies in the time range or monitoring objects corresponding to the two scores.

[0037] The method in this embodiment includes the following steps.

[0038] S1. Collect operational data from multiple subsystems in the subway SoS, perform time alignment and feature processing on the operational data, and obtain multimodal time window samples.

[0039] In one specific implementation, the operational data of the vehicle subsystem includes train identification, sampling timestamp, train position, train speed, traction force, braking force, train mass, and operating conditions; the operational data of the signaling subsystem includes signal equipment identification, signal status, route status, movement authorization endpoint, and temporary speed limit; the operational data of the traction power supply system includes feeder voltage, feeder current, substation injected power, and the traction or regenerative braking status of online trains; the operational data of the passenger flow subsystem includes the number of passengers entering the station, the number of passengers exiting the station, the number of passengers transferring, and the estimated train passenger capacity. Onboard air spring pressure, as load cross-validation data, is obtained from the vehicle subsystem or an independent onboard monitoring interface. Video surveillance data is not a required input; in another implementation, numerical features such as platform passenger density and door area congestion are first extracted from the video frames, and then the extracted numerical features are used as part of the passenger flow modality, instead of directly inputting the original video frames into the subsequent physical constraint module.

[0040] To maintain a fixed feature dimension, detection areas are pre-divided and a maximum number of online trains is set. For multiple online trains within the same detection area, they are sorted according to their spatial position along the track and written into their corresponding train slots; empty slots are filled with preset default values ​​and a train presence flag is set. Another approach is to use a single train as the detection object, associating the train's associated forward signals, power supply sections, and passenger flow information with that train's feature vector. The former method preserves the collaborative relationships between multiple trains within a section, while the latter method occupies fewer feature dimensions, making it suitable for onboard edge deployment.

[0041] All data sources use a unified clock for time synchronization. Preferably, the scheduling center's time server or a precise time protocol clock is used as the unified time reference, and the original sampling timestamps of each data source are mapped to an interval of [missing information]. A unified time axis is used. For continuous variables with sampling frequencies higher than the unified frequency, the mean, median, or last value can be taken within the corresponding time interval; for continuous variables with lower sampling frequencies, linear interpolation can be used within the allowed time interval; for discrete variables such as signal states, the most recent previous value is retained. If the duration of missing data for a continuous variable exceeds a preset missing duration, or if the proportion of invalid features within the same time window exceeds a preset proportion, then the samples in that time window are marked as unscored samples and are not written into the normal historical score sliding window. This avoids large-scale interpolated data being mistakenly identified as real normal data.

[0042] In one example, uniform sampling interval Take 1 second; when it is necessary to detect transient anomalies in traction or braking, The sampling interval can range from 100ms to 500ms. The selection of a uniform sampling interval should ensure that the main time scales of train dynamics, signal state changes, and power supply response can be covered by multimodal time window samples.

[0043] Regarding the first Each continuous feature is used to calculate the mean from training data that has been confirmed to be under normal operating conditions. and standard deviation And standardize it according to the following formula: in, Indicates the sampling time The A continuous feature, This represents the corresponding standardized feature. This represents a stability constant to prevent the denominator from reaching zero. For discrete states such as signal state, route state, and traction / regenerative braking state, one-hot encoding is used, and the encoded values ​​are scaled to a range similar to that of continuous features. The mean, standard deviation, one-hot encoding order, and feature index are stored together in a feature mapping table to ensure that the same data definitions are used in the training phase, online detection phase, and anomaly feature localization phase.

[0044] Continuous The multimodal feature vectors from each sampling time point are stacked in chronological order to obtain: in, Represents multimodal time window samples, Indicates the length of the time window. Preferably, Multimodal time-window samples preserve the lateral coupling relationships between different subsystems at the same sampling time in the feature dimension, and retain continuous time in the temporal dimension. The state evolution relationship at each sampling time.

[0045] S2. Construct a set of physical constraints based on the operating mechanisms of multiple subsystems, and determine the physical residuals, physical violation energy, and physical violation scores based on the set of physical constraints.

[0046] The physical constraint set consists of multiple physical residual functions The residuals are structured as follows: when the inputs satisfy the corresponding physical consistency relationship, the residuals approach zero; when the correspondence is violated, the residuals increase. Since the diffusion model uses standardized data, while dynamics and energy conservation calculations require actual dimensions, continuous features are de-standardized according to the feature mapping table before physical constraint calculations, and the soft-coded or one-hot coded discrete states are converted into state probabilities. During the training phase, continuous probabilities are used for discrete state channels, and Softplus-type smoothing is applied to segmented operations such as maximum and minimum values; during the online detection phase, the collected one-hot states and corresponding hard-decision forms are used. To ensure that force, power, mass, and probability residuals are weighted on a uniform scale, the original residuals of each class are dimensionless using statistics from the normal training data.

[0047] in, Let represent the j-th type of original physical residual; the mean and standard deviation of this type of normal training residual are used for centering and scaling, respectively, with a stability constant added to the denominator; the quantity before the equal sign represents the dimensionless residual. The dynamic, power supply, and load standardized residuals in the following violation energy formulas are all calculated according to this formula; the interlocking transfer residuals can also be dimensionless according to their normal historical statistics.

[0048] Huber loss is used to suppress occasional outliers in each physical residual. The threshold is... The Huber loss function is expressed as: in, This represents the physical residual to be processed. This represents the Huber loss threshold. Compared to directly using squared loss, this function maintains sensitivity to early biases when residuals are small, and limits their excessive influence on physical violation scores when residuals are caused by communication glitches or isolated outliers.

[0049] For the longitudinal dynamic constraints of the train, let the first... The train speed at each sampling time is The train acceleration is expressed as: The basic running resistance of the train is determined using the Davis equation: in, , and These are the operating resistance parameters determined based on vehicle type, train formation, and track testing. Track gradient resistance can be expressed as... ,in, Indicates the first Train mass at each sampling time, Represents gravitational acceleration. This indicates the track gradient angle corresponding to the train's current position. It is obtained by adding the basic running resistance and the gradient resistance. Considering the slewing mass coefficient Subsequently, the longitudinal dynamic residuals of the single-step train Represented as: in, Indicates the train's traction force. This indicates the train's braking force. If the equivalent mass is already given in the vehicle model parameters, then the equivalent mass is used directly instead. The dynamic violation energy corresponding to a multimodal time window sample is expressed as: in, The Huber loss threshold, representing the dynamic residual, is preferably set to 0.7 times the standard deviation of the normal historical dynamic residual. This constraint can detect states that disrupt the force balance, such as traction response hysteresis, braking force decay, mechanical braking jamming, and significant errors in train mass estimation.

[0050] For signal security logic constraints, the first The signal state at each sampling time is represented as follows: dimensional state vector Generate a valid state transition matrix based on the route interlocking table. The predicted signal state probability at the next sampling time is expressed as: The interlocking transfer residual is expressed as: in, The first digit represents the actual signal state vector at the next sampling time. One portion, The first digit represents the probability of the predicted signal state. Each component has a valid transition matrix, which can be directly generated from the existing interlocking table or stored separately for different stations. When the operation diagram or interlocking table is updated, only the corresponding matrix and version number are updated, without changing the subsequent diffusion model interface.

[0051] Movement authorization protection constraints are used to determine whether the train speed meets the braking protection requirements between the current position and the movement authorization endpoint. Let the movement authorization endpoint be... The train's current location is The preset emergency braking deceleration is The system's response distance is The maximum safe speed corresponding to the current position is expressed as: Mobile authorization protection residuals are expressed as: In a preferred embodiment, for Then use the threshold as The Huber loss is processed, and The standard deviation of the normal historical protection residual is taken as 0.5 times. The interlocking transfer residual and the mobile authorization protection residual are aggregated within the time window and then combined using preset weights to form the signal security logic violation energy. Interlocking transfer is used to identify illegal transitions in discrete states, while movement authorization protection is used to identify abnormal relationships in continuous speed protection. The combination of the two can simultaneously cover discrete security logic and continuous protection curves.

[0052] Regarding traction power flow constraints, in the first At each sampling time, the total electrical power injected into the feeder lateral traction power supply network is calculated based on the voltage and current of each feeder. The total net power consumption is determined based on the traction or regenerative braking conditions of each online train. Among them, traction power consumption, other loads, losses, and energy storage charging power are taken as positive values, while regenerative braking feedback and energy storage discharging power are taken as negative values. The single-step traction power flow residual is expressed as: The corresponding traction power flow violation energy is represented as: in, The optimal value is 0.8 times the standard deviation of the normal historical traction power supply residual. For lines with multiple traction substations, the power supply residual is calculated separately for each power supply section, and then the illegal energy of each power supply section is weighted and summed to avoid incorrect power cancellation between different power supply sections.

[0053] For load consistency constraints, the load of the first train is determined based on passenger flow statistics from the automatic fare collection system and the train operation plan. For example, the empty weight of the train is added to the product of the estimated number of passengers and the average passenger weight. The load capacity of the second train is determined based on the onboard air spring pressure and the vehicle calibration curve. The single-step load consistency residual is expressed as: The corresponding load consistency violation energy is represented as: in, The optimal value is 0.8 times the standard deviation of the historical load consistency residual. If AFC data can only be updated by station section, the first train load can be updated when the train departs and maintained within adjacent station sections; air spring pressure is continuously updated along a uniform time axis. This process can identify weighing sensor anomalies or passenger flow statistics failures using load estimation results from two independent sources.

[0054] The physical violation energy is obtained by weighting the violation energy corresponding to each physical constraint in the current multimodal time window sample according to the safety criticality level: in, Indicates the window sequence number. Indicates the number of physical constraints. Indicates the first The weight of each physical constraint, Indicates the first The first time window The violation energy corresponds to each physical constraint. Weights can be pre-set and normalized according to the safety criticality level, ensuring the sum of all weights is 1. An exponential moving average is applied to the physical violation energy over consecutive time windows to obtain the physical violation score. in, The forgetting factor, preferably 0.9, is used. The first effective window initializes the physical violation score with its physical violation energy. Huber loss suppresses individual outliers, and exponential moving average accumulates continuous small deviations. The combination of these two methods can improve sensitivity to asymptotic anomalies while avoiding instantaneous false alarms.

[0055] In another specific implementation, when some lines do not yet have complete AFC or air spring pressure interfaces, a physical constraint set is constructed using train longitudinal dynamics constraints, signal safety logic constraints, and traction power flow constraints, and the weights of the missing constraints are renormalized. This approach is compatible with phased deployment; adding load consistency constraints after the data interface is complete can further enhance the ability to detect cross-system data inconsistencies.

[0056] S3. Using multimodal time window samples under normal operating conditions as training samples, the diffusion model is jointly trained based on data fidelity loss and physical consistency loss to obtain a physical-data co-diffusion model.

[0057] like Figure 2 As shown, the diffusion model includes a forward diffusion module, a time-step encoding module, a denoising network, and a clean sample estimation module. Training samples are denoised by the forward diffusion module and then input into the denoising network, which outputs predicted noise. One training branch determines the data fidelity loss based on the real Gaussian noise and the predicted noise. Another training branch obtains clean sample estimates from the clean sample estimation module and inputs them into the physical constraint set to determine the physical consistency loss. Finally, the denoising network parameters are updated based on the total loss from joint training.

[0058] In one specific implementation, the forward diffusion module uses the total number of diffusion steps. A fixed Markov chain. Let the first... The noise variance of each diffusion time step is , , Then the training sample is in the th The noisy sample corresponding to each diffusion time step is represented as follows: Cosine noise modulation is preferred, and the cumulative noise retention factor is expressed as: in, This is the cosine scheduling offset, which can be set to 0.008. During training, it ranges from 1 to... Uniform extraction diffusion time step This enables the denoising network to learn the recovery relationship under different noise levels.

[0059] The denoising network employs a one-dimensional temporal U-Net structure, including an encoder, a temporal self-attention bottleneck layer, a decoder, and skip connections connecting network layers of the same scale in the encoder and decoder. Both the encoder and decoder comprise multiple one-dimensional convolutional residual blocks. Each one-dimensional convolutional residual block sequentially includes a one-dimensional convolutional layer, a normalization layer, a non-linear activation layer, and another one-dimensional convolutional layer, with the input features added to the convolutional output through residual connections. Preferably, the convolutional kernel length is 3, the non-linear activation function is SiLU, and the normalization layer uses group normalization. Using one-dimensional convolution allows for local feature extraction along the temporal dimension, with each multimodal feature serving as channel input, thus enabling simultaneous perception of adjacent temporal state changes and cross-modal channel combinations.

[0060] In a more complete server implementation, the encoder comprises four scales with channel numbers of 64, 128, 256, and 512, respectively. Each scale includes two one-dimensional convolutional residual blocks and is downsampled using one-dimensional convolution with a stride of 2. The decoder employs the corresponding four scales, recovering temporal resolution through upsampling and one-dimensional convolution. The temporal self-attention bottleneck layer uses a multi-head self-attention structure with eight attention heads to extract long-range correlations between signal openness, train traction response, and power supply load changes along the time axis. Skip connections concatenate or add features from the encoder at the corresponding scale to the decoder features to preserve weak transient information.

[0061] The time-step encoding module consists of a sinusoidal encoding layer and a feature linear modulation layer. The sinusoidal encoding layer encodes the diffusion time steps into time-step vectors, which are then used by a multilayer perceptron to generate the scaling and offset for each one-dimensional convolutional residual block. For the intermediate features of the residual block... The characteristic linear modulation is expressed as: in, and These represent the scaling and offset generated by the time step encoding, respectively. This represents element-wise multiplication. Therefore, denoising tasks at different levels of diffused noise can share the main network structure, while the features of each layer are adjusted by time-step encoding.

[0062] Noisy samples and diffusion time step Input parameters are The denoising network is used to obtain the predicted noise. The clean sample estimation module determines the clean sample estimate based on the noisy sample, the cumulative noise retention coefficient, and the predicted noise: Data fidelity loss is expressed as: in, This represents the expectation operation. This represents the square of the L2 norm. This loss directly constrains the predicted noise of the denoising network to be close to the Gaussian noise actually added by the forward diffusion module.

[0063] The clean sample estimates are input into the physical constraint set, and the physical residuals corresponding to each physical constraint are calculated. These residuals are then weighted according to the constraint weights. The physical consistency loss is expressed as: in, This represents the first value calculated based on the clean sample estimate. One physical residual, This represents the corresponding Huber loss threshold. For physical constraints that generate multiple single-step residuals within a time window, the violation energy of the constraint is first obtained according to the time aggregation method of the corresponding constraint, and then written into the above weighted sum.

[0064] The weights of the total loss from joint training and the physical consistency loss are expressed as follows: in, Indicates training rounds, The basic weight representing the physical consistency loss is preferably set to 0.005; This represents the warm-up coefficient, which increases linearly from 0 to 1 over the first 50 training rounds. This represents the noise level adjustment factor. In the high-noise stage, the estimation error of clean samples is relatively large, and the adjustment factor makes the weight of the physical consistency loss relatively low; in the low-noise stage, the estimation of clean samples is closer to the original samples, and the effect of the physical consistency loss is enhanced, thereby avoiding unreliable high-noise estimation from imposing too strong constraints on training.

[0065] Vehicle parameters, route parameters, interlocking matrices, braking parameters, and sensor calibration parameters from the physical constraint set are stored as non-training parameters. During backpropagation, gradient isolation is performed on these physical rule parameters, while maintaining... The computational path between the target noise and the target noise is continuous, making The gradient can be fed back to the denoising network via the clean sample estimation module. This gradient path is represented as: During implementation, it was also used Update the normal historical standard deviation, Huber threshold, or statistical monitoring quantity of the physical residuals; Not used as a substitute for those participating in backpropagation This approach retains the gradient-free statistical update method from the attached implementation plan while avoiding the physical consistency loss that would prevent training the denoising network if the clean sample estimates were truncated as a whole.

[0066] In one training example, manually reviewed multimodal time-window samples under normal operating conditions are used, with 70% used as the training set, 15% as the validation set, and 15% as the test set. The same time-window samples are not reused across sets; the training segment only retains manually reviewed normal operating condition samples. The batch size is 64, the optimizer is AdamW, and the initial learning rate is... The system is trained for 200 epochs, and early stopping is performed based on the total loss of the validation set. For edge devices with limited computing resources, a one-dimensional temporal U-Net with three scales is used, with the number of channels being 64, 128, and 256 respectively, and the number of attention heads is reduced to 4. The server implementation has strong modeling capabilities for long-range associations, while the lightweight implementation has low inference resource consumption.

[0067] S4. Input the multimodal time window samples to be detected into the physical-data co-diffusion model for reconstruction, and determine whether there are anomalies in the subway SoS based on the reconstruction anomaly score and the physical violation score.

[0068] like Figure 3 As shown, the multimodal time window samples to be detected simultaneously enter the DDIM single-step reconstruction branch and the physical consistency calculation branch. The two branches output reconstruction anomaly scores and physical violation scores, respectively. These scores are then standardized through their respective normal historical score sliding windows to obtain data anomaly deviation and physical anomaly deviation. Finally, they are fused and thresholded based on the current operational scenario.

[0069] Set the noise evaluation step to the 150th diffusion time step, i.e. To reduce the impact of repeated random sampling on the reconstruction results of the same sample to be detected, evaluation noise with the same dimension as the multimodal time window sample to be detected is obtained according to a preset reproducible evaluation noise generation rule, and noise is added according to the following formula: The noisy time window sample to be detected is input into the denoising network to obtain the predicted noise, and the reconstructed sample is generated by deterministic single-step inversion based on the denoising diffusion implicit model: The reconstructed outlier scores were obtained using the mean squared error: in, This represents the square of the Frobenius norm. While determining the reconstructed anomaly score, samples from the same multimodal time window to be detected are input into the physical constraint set. The current physical violation energy is determined according to the method in S2, and an exponential moving average is performed combining the physical violation scores from the previous time window to obtain the current physical violation score. .

[0070] Score sliding windows containing only historical scores under normal operating conditions are established for both reconstructed abnormal scores and physical violation scores, with a preferred window length of 500. During initial system deployment, initialization is performed using scores corresponding to historical multimodal time windows confirmed to be under normal operating conditions. The current scores for both channels are standardized based on the mean and standard deviation of historical scores within the corresponding sliding windows. in, and Let represent the mean and standard deviation of the reconstructed outlier score sliding window, respectively. and Let represent the mean and standard deviation of the sliding window for physical violation scores, respectively. Indicates the degree of abnormal deviation in the data. This indicates the degree of deviation from physical anomalies. Including a stability constant can handle situations where the standard deviation is close to zero during the system initialization phase or when fractional fluctuations are extremely small.

[0071] The fusion weight is determined based on the current operational scenario, and the overall anomaly score is calculated: in, This indicates the fusion weight. The current operating scenario is determined jointly based on the operating schedule, passenger load factor, weather alarms, and dispatch status. When the current operating scenario is an off-peak operating scenario, When operating during peak hours or in severe weather, When operating in a low-density, nighttime scenario, To avoid frequent scene switching under critical conditions, the scene label should only be updated after the new scene conditions are met in multiple consecutive time windows, or different hysteresis intervals should be set for the entry and exit thresholds.

[0072] The system establishes a sliding window for comprehensive anomaly scores for different operational scenarios, and determines the 99th percentile of the historical comprehensive anomaly scores within the window corresponding to the current operational scenario as the anomaly judgment threshold. .when When an anomaly is detected in the subway's SoS (Solar System), an alarm is triggered; when At that time, it is determined that there are no anomalies in the subway SoS, and the current reconstruction anomaly score, physical violation score, and comprehensive anomaly score are updated to the corresponding sliding windows. Updating the sliding window only with samples deemed normal can reduce the risk of anomalies gradually contaminating the normal statistical baseline.

[0073] In another specific implementation, the fusion weights are calibrated based on verification data within a preset range, but are still stored separately for different operating scenarios; the anomaly detection threshold can also be a quantile between the 98th and 99.9th percentiles. Using fixed scenario weights facilitates deployment and interpretation, while using verification data calibration can adapt to different line noise levels and equipment configurations.

[0074] S5. After determining that there is an anomaly in the subway SoS, perform a two-level decomposition of physical violation energy and reconstruction error, and generate anomaly diagnosis information based on the causal knowledge graph.

[0075] like Figure 4 As shown, the first layer is physical residual decomposition, used to identify the violated physical constraints or safety rules; the second layer is reconstruction error feature decomposition, used to identify feature variables that deviate from the normal distribution, as well as the corresponding subsystems and monitoring equipment. Both types of results serve as independent anomalous evidence input into the causal knowledge graph to generate anomalous causal chains.

[0076] No. The contribution of each physical constraint to the current physical violation energy can be expressed as: The physical constraints are sorted from highest to lowest contribution, and those with the highest contribution, or multiple physical constraints whose contribution exceeds a preset contribution threshold, are identified as the violated physical constraints. Using multiple constraints as the common output method is suitable for situations where traction power supply anomalies simultaneously cause abnormal train dynamic response, and can preserve evidence of cross-system linkage.

[0077] For the The error between the multimodal time window samples to be detected and the reconstructed samples is aggregated along the time dimension using 10 feature variables, resulting in: The corresponding feature contribution is expressed as: Based on the feature mapping table, the feature variables with higher contribution are mapped to the corresponding subsystems, monitoring equipment, and physical quantities. For example, traction force and braking force characteristics correspond to vehicle traction and braking equipment, feeder voltage and feeder current correspond to traction substations or power supply sections, signal status and movement authorization endpoints correspond to signal equipment, and air spring pressure corresponds to vehicle suspension monitoring equipment.

[0078] The causal knowledge graph uses subsystems, devices, physical constraints, and observations as nodes, and establishes paths through associative edges such as devices belonging to subsystems, devices generating observations, physical constraints using observations, and device failures affecting observations. Violated physical constraints and anomalous characteristic variables are used as anomalous evidence. Paths covering both types of evidence are retrieved, and candidate paths are ranked based on path length, node matching degree, and contribution of both types of evidence. The final generated anomaly diagnostic information can include anomaly occurrence time, comprehensive anomaly score, main physical violation rules, main anomalous characteristic variables, candidate faulty subsystems, candidate faulty devices, and anomaly causal chains.

[0079] For example, when the deviation between the feeder injected power and the net power consumed by the online train increases significantly, and the reconstruction error of the contact network voltage characteristics contributes the most, a candidate causal chain can be formed: traction power supply system - contact network or substation - feeder voltage observation - traction power supply energy conservation constraint. When the load estimated by AFC is inconsistent with the load converted by air spring pressure, and the air spring pressure characteristics are abnormal, an anomaly in the on-board pressure sensor or suspension monitoring equipment can be prioritized. Mutual corroboration of physical and data evidence can reduce misjudgments caused by relying solely on a single out-of-limit parameter for positioning.

[0080] S6. Deploy and execute the subway SoS anomaly detection method in the hardware environment.

[0081] In one specific implementation, such as Figure 5 As shown, the input hardware includes train control equipment, trackside sensing equipment, video surveillance equipment, and a dispatching and command center interface; the computing end includes servers with CPUs and GPUs, memory for storing model parameters and historical sliding windows, communication buses, and network buses; the output end includes a metro SoS simulation test platform, an operation and maintenance terminal, and an alarm log storage.

[0082] The CPU can perform time alignment, missing value handling, physical residual calculation, sliding window updates, and causal knowledge graph retrieval, while the GPU can perform training of the one-dimensional temporal U-Net and single-step inference of DDIM. The physical constraint calculation branch and the diffusion model reconstruction branch are associated through the same time window number, and score fusion is only performed after both calculations are completed. If either calculation times out or the input data integrity does not meet the requirements, an output data quality alarm is issued, rather than directly classifying the time window as a system anomaly.

[0083] Model parameters, feature mapping tables, physical rule parameters, legal signal state transition matrices, fusion weights for each operational scenario, sliding window data, and anomaly detection thresholds are all managed through version control. After updating the runtime diagram, vehicle parameters, or interlocking rules, the new version is first replayed and verified on the simulation testing platform before switching to the online version. This approach avoids anomaly scoring distortion caused by incompatibility between the model version and the physical rule version.

[0084] In another implementation method where onboard computing resources are limited, data acquisition, basic preprocessing, and lightweight diffusion model reconstruction can be completed on the onboard unit, while power flow constraints, cross-train information fusion, and causal knowledge graph analysis can be performed on the ground server. The onboard unit and the ground server exchange intermediate results using a unified time window number. Centralized deployment facilitates the comprehensive utilization of data from multiple subsystems, while distributed deployment can reduce the amount of raw high-frequency data transmission and improve continuous detection capabilities under partial communication interruption conditions.

[0085] Example 2: Based on the aforementioned embodiments, this embodiment experimentally verifies the subway SoS anomaly detection method based on the physical-data co-diffusion model to examine the method's performance in multimodal cross-system coupled anomaly detection, false alarm control capability, and stability over continuous time intervals.

[0086] The experiment used actual operation data of a city's subway system to construct a dataset, which was then divided into training, validation, and test sets according to 70%, 15%, and 15% of the time sequence, respectively. The training segment used only manually verified normal operating condition samples; after the time segmentation, four types of cross-system coupled anomalies (A, B, C, and D) were injected into the validation and test segments based on a physical model to prevent anomaly injection information from entering the training process. The multimodal time window length L was set to 64. The unified test set contained 37,333 time window samples, including 33,546 normal samples and 3,787 abnormal samples. To verify temporal stability, the unified test set was divided into 5 consecutive test blocks without changing the sample order; the first 3 test blocks each contained 7,467 samples, and the last 2 test blocks each contained 7,466 samples, totaling the same number of samples as the unified test set.

[0087] For comparison, the isolated forest is denoted as M1, the single-class support vector machine as M2, the long short-term memory autoencoder as M3, the pure diffusion model without physical consistency loss and scene adaptive fusion mechanism as M4, and the physical-data co-diffusion model of this invention as M5. In M5, the basic weight of the physical consistency loss is set to 0.005, a linear warm-up is performed in the first 50 training rounds, and a gradient isolation mechanism is used. The Huber loss threshold corresponding to each physical residual is set according to the aforementioned implementation method. The evaluation noise step t_eval for online detection is set to 150. The standardized sliding window length corresponding to the reconstruction anomaly score and physical violation score is set to 500. The fusion weight is determined according to the operational scenario, and the anomaly judgment threshold is the 99th percentile of the comprehensive anomaly score sliding window under the corresponding operational scenario.

[0088] M4 uses the same diffusion model network structure, training data, and evaluation noise step as M5, but does not use physical consistency loss, calculate physical violation scores, or employ dual-score scene adaptive fusion and dynamic thresholding; its decision threshold is determined on the validation set with the objective of maximizing the F1 score. M1 to M3 use the same multimodal time window samples and the same training, validation, and testing time divisions, and configure hyperparameters and determine thresholds on the validation set according to their respective methods. This setup is used to compare the overall performance of the M5 combination scheme with each control model, without interpreting the difference between M4 and M5 as an independent contribution of any single technical feature.

[0089] The experiment was evaluated from three dimensions: First, precision, accuracy, and F1 score were used to evaluate the fixed-point detection capability under a fixed decision threshold; second, the average false positive rate on five time-series test blocks was further compared; third, the decision results of M5 were decomposed into sample counts of normal decisions, false positives, false negatives, and abnormal decisions to verify the results. Figure 7 The source of the proportional indicators. Fourth, considering that the aggregated results cannot individually demonstrate the stability across different time periods, Figure 9 The F1 scores of each model on five consecutive test blocks were evaluated. The experimental results are as follows: Figures 6 to 9 As shown.

[0090] like Figure 6 As shown, M5 achieved the highest precision (0.9881), F1 score (0.9301), accuracy (0.9866), and specificity (0.9988) among the five models on the unified test set; its false positive rate (0.0012) was the lowest among the five models. Therefore, Figure 6 M5 achieved the best results across all listed indicators.

[0091] Further comparison of M4 and M5 with the same structure showed that M5 improved precision by 0.0269, F1 score by 0.0091, accuracy by 0.0020, and specificity by 0.0028 compared to M4; the false positive rate decreased from 0.0040 to 0.0012, a reduction of approximately 70.37%. These results indicate that the synergistic setting of physical consistency constraints and scene-adaptive dynamic thresholds can improve alarm precision, overall detection performance, and normal operating condition recognition capabilities. Since both have the same diffusion network structure, training data, and evaluation noise step, these results support the comprehensive improvement achieved by the combined techniques of physical co-training, dual-score fusion, and scene-adaptive dynamic thresholds added to M5 compared to M4. Because subway online alarms also require reducing false triggers under normal operating conditions, [further details are needed]. Figure 6 Further derivation to Figure 7 The false positive rate was examined separately.

[0092] like Figure 7 As shown, after calculating the average false positive rate for each of the five time-sequenced test blocks, the average false positive rate for M5 was 0.0012, lower than that for M3 (0.0038) and M4 (0.0040). Compared to M4, the average false positive rate of M5 decreased by approximately 70.37%; compared to M3, it decreased by approximately 68.25%. Therefore, under the continuous time segmentation conditions of this embodiment, the suppression effect of M5 on false triggering under normal operating conditions remained stable. Considering that the false positive rate is a proportional indicator, its value still needs to be verified by the original classification count. Therefore, the next step is to... Figure 8 Decompose the judgment result of M5.

[0093] like Figure 8 As shown, out of 33,546 normal samples, M5 correctly identified 33,506 as normal and only 40 as abnormal; out of 3,787 abnormal samples, it correctly identified 3,327 as abnormal. This yields a precision of 0.9881, a specificity of 0.9988, and a false positive rate of 0.0012. Compared to M4, M5 reduced the number of false positive samples from 135 to 40, a reduction of 95, or approximately 70.37%. Among the five models, M5 had the fewest false positive samples, the highest precision, and the highest specificity. However, the aggregate count on the unified test set cannot indicate whether the effect is concentrated in a specific time period; therefore, further analysis is needed... Figure 8 Derivation to Figure 9 Perform time-series stability verification.

[0094] like Figure 9As shown, the F1 scores of M5 on test blocks 1 through 5 are 0.9349, 0.9364, 0.9267, 0.9306, and 0.9227, respectively, representing the highest values ​​among the five models for each test block. The average F1 score across the five test blocks is 0.9303, with a standard deviation of 0.0057, and the difference between the maximum and minimum scores is 0.0137. This demonstrates that after splitting the test set chronologically, M5 maintains relatively small performance fluctuations.

[0095] In summary, based on the above results, M5 is the optimal model among the five models in terms of precision, F1 score, accuracy, specificity, and false positive rate listed in this embodiment. Specifically, M5 achieved the highest precision of 0.9881, the highest F1 score of 0.9301, the highest accuracy of 0.9866, the highest specificity of 0.9988, and the lowest false positive rate of 0.0012. The results on the five time-series test blocks also show that M5 achieved the highest F1 score in each test block. These results support the synergistic effect of physical consistency constraints, data reconstruction scores, and scene-adaptive dynamic thresholds in improving the accuracy, stability, and engineering usability of subway SoS anomaly detection.

[0096] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0097] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any person skilled in the art can easily conceive of various equivalent modifications or substitutions within the technical scope disclosed in the present invention, and these modifications or substitutions should all be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A subway SoS anomaly detection method based on a physical-data co-diffusion model, characterized in that, Includes the following steps: Operational data from multiple subsystems in the subway SoS are collected, and time alignment and feature processing are performed on the operational data to obtain multimodal time window samples. A set of physical constraints is constructed based on the operating mechanisms of the multiple subsystems. The set of physical constraints is used to characterize the physical consistency relationships between the multiple subsystems themselves and among the multiple subsystems. Using the multimodal time window samples under normal operating conditions as training samples, the diffusion model is jointly trained based on the data fidelity loss of the training samples and the physical consistency loss determined according to the physical constraint set, to obtain the physical-data co-diffusion model. The multimodal time window sample to be detected is input into the physical-data co-diffusion model for reconstruction. Based on the difference between the multimodal time window sample to be detected and the corresponding reconstructed sample, the reconstruction anomaly score is determined. Physical consistency detection is performed on the multimodal time window samples to be detected based on the physical constraint set to determine the physical violation score; By combining the reconstructed anomaly score and the physical violation score, a comprehensive anomaly score is obtained. Based on the comparison result between the comprehensive anomaly score and the anomaly judgment threshold, it is determined whether the subway SoS has an anomaly.

2. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, The process of collecting operational data from multiple subsystems within the subway SoS (Solar Operation System), performing time alignment and feature processing on the operational data to obtain multimodal time window samples, specifically includes: Operational data with sampling timestamps are collected from the vehicle subsystem, signaling subsystem, traction power supply system, and passenger flow subsystem, respectively. The running data is timestamped and missing values ​​are processed according to a unified timeline to obtain aligned running data corresponding to each sampling time. According to the data mode to which the running data belongs, the continuous variables in the aligned running data are standardized, and the discrete states in the aligned running data are one-hot encoded and scaled. The processed running data of each subsystem corresponding to the same sampling time are spliced ​​together to obtain the multimodal feature vector corresponding to that sampling time. The multimodal feature vectors of a preset number of consecutive samples are stacked in chronological order to obtain the multimodal time window sample.

3. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, The construction of the physical constraint set based on the operating mechanisms of the multiple subsystems specifically includes: Based on the mechanical relationship between the train's traction, braking, and running resistance, longitudinal dynamic constraints of the train are constructed. Based on the legal transfer relationships of subway signal states and the security relationship between train speed and movement authorization protection range, construct signal security logic constraints; Based on the energy conservation relationship between the injected power of the traction power supply network and the net power consumption of online trains, traction power flow constraints are constructed. Based on the consistency relationship between the train load estimated from passenger flow data and the train load deduced from onboard air spring pressure data, load consistency constraints are constructed. The train longitudinal dynamics constraints, the signal safety logic constraints, the traction power flow constraints, and the load consistency constraints are respectively converted into differentiable residual functions to form the physical constraint set. For segmented relationships containing maximum values, minimum values, or state transitions, corresponding smoothing processing is adopted during the joint training phase to ensure that the gradient can be propagated. Each residual function tends to zero when the corresponding physical consistency relationship is satisfied and increases when the corresponding physical consistency relationship is violated.

4. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 3, characterized in that, The construction of the train longitudinal dynamic constraints and the signal safety logic constraints specifically includes: The train acceleration is determined based on the train speed at adjacent sampling times, the basic running resistance of the train is determined based on the train speed and preset running resistance parameters, and the longitudinal dynamic residual of the train is determined based on the force balance relationship between the train traction force, train braking force, train mass, the train acceleration, the basic running resistance of the train and the track gradient resistance. The signal state at each sampling time is encoded into a one-hot vector. Based on the preset legal signal state transition matrix and the signal state at the current sampling time, the predicted signal state distribution at the next sampling time is determined. Based on the difference between the predicted signal state distribution and the actual signal state at the next sampling time, the interlocking transfer residual is determined. Based on the train's current position, the endpoint of the movement authorization, the preset emergency braking deceleration, the system response travel distance, and the permissible speed on the line, the limit safe speed corresponding to the train's current position is determined, and the movement authorization protection residual is determined based on the degree to which the train's actual speed exceeds the limit safe speed. The interlocking transfer residual and the mobile authorization protection residual are weighted and combined to obtain the signal security logic residual.

5. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 3, characterized in that, The determination of the traction power flow constraints, the load consistency constraints, and the physical violation scores specifically includes: The total injected power on the feeder side into the traction power supply network is obtained, and the total net power consumption of the traction power supply network is determined based on the traction or regenerative braking power of each online train, the power of other loads, the power loss of the line and converter, and the net charging power of the energy storage device. The traction power flow residual is determined based on the difference between the total injected power on the feeder side and the total net power consumption. The load of the first train is determined based on the passenger flow statistics of the automatic fare collection system, and the load of the second train is determined based on the air spring pressure. The load consistency residual is determined based on the difference between the load of the first train and the load of the second train. The residuals in the physical constraint set are dimensionless by using the mean and standard deviation of the normal historical residuals respectively. Then, the dimensionless residuals are processed by Huber loss. The processing results are weighted according to the safety criticality level of each physical constraint to obtain the physical violation energy corresponding to the current multimodal time window sample. The Huber threshold is a threshold that is preset on the dimensionless scale or calibrated using the validation set. According to a preset forgetting factor, the physical violation energy corresponding to multiple consecutive multimodal time window samples is subjected to an exponential moving average to obtain the physical violation score.

6. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, The diffusion model includes a forward diffusion module, a time-step encoding module, a denoising network, and a clean sample estimation module; The denoising network adopts a one-dimensional temporal U-Net structure, including an encoder, a temporal self-attention bottleneck layer, a decoder, and skip connections connecting network layers of the same scale in the encoder and the decoder. The encoder and the decoder each include multiple one-dimensional convolutional residual blocks. The temporal step encoding module includes a sinusoidal encoding layer and a feature-level linear modulation layer. The feature-level linear modulation layer is connected to each of the one-dimensional convolutional residual blocks. The training samples are input into the forward diffusion module, and Gaussian noise is added to the training samples at selected diffusion time steps according to the preset noise schedule to obtain noisy samples. The diffusion time step is input into the sinusoidal coding layer to obtain the time step code, and the feature-level linear modulation layer modulates the features of each one-dimensional convolutional residual block according to the time step code; The noisy sample is input into the encoder for step-by-step downsampling and feature extraction. The features output by the encoder are input into the temporal self-attention bottleneck layer to extract the correlation features of the multimodal time window sample in the temporal dimension. Then, the decoder performs step-by-step upsampling and fuses the features of the corresponding scale in the encoder through the skip connection to output the predicted noise corresponding to the noisy sample. The clean sample estimation module determines the clean sample estimate corresponding to the training sample based on the noisy sample, the noise scheduling coefficient corresponding to the diffusion time step, and the prediction noise. The prediction noise is used to determine the data fidelity loss, and the clean sample estimate is used to determine the physical consistency loss based on the physical constraint set.

7. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 6, characterized in that, The joint training of the diffusion model based on the data fidelity loss of the training samples and the physical consistency loss determined according to the physical constraint set specifically includes: The data fidelity loss is determined based on the difference between the Gaussian noise added to the training samples by the forward diffusion module and the predicted noise output by the denoising network. The data fidelity loss is expressed as: in, This indicates the loss of data fidelity; This represents the expectation operation; This represents the multimodal time window sample used as the training sample; This indicates the addition of Gaussian noise to the training samples; Indicates the selected diffusion time step; This represents the noisy sample obtained at the diffusion time step; The parameter is The denoising network outputs predicted noise based on the noisy sample and the diffusion time step; Represents the square of the L2 norm; Based on the physical constraints in the physical constraint set, physical residuals are calculated for the clean sample estimates. Huber loss is applied to each physical residual, and the results are weighted according to the corresponding constraint weights to determine the physical consistency loss. The physical consistency loss is expressed as follows: in, This represents the loss of physical consistency; This indicates the number of physical constraints in the physical constraint set; Indicates the first The constraint weights corresponding to each physical constraint; This indicates the estimated value based on the clean sample. The calculated first One physical residual; Indicates the threshold is Huber loss function; Indicates the first The Huber loss threshold corresponding to each physical residual; The total loss for joint training is determined based on the data fidelity loss and the physical consistency loss, and the weight of the physical consistency loss is adjusted according to the training progress and the noise level corresponding to the diffusion time step. The total loss and the weight of the physical consistency loss are expressed as follows: in, This represents the total loss; Indicates the first The training rounds and the first The weights of the physical consistency loss at each diffusion time step; The underlying weights representing the loss of physical consistency; This represents the warm-up coefficient, which increases linearly with the number of training rounds. The preset noise scheduling indicates that the preset noise scheduling is in the first... The cumulative noise retention factor corresponding to each diffusion time step; During backpropagation based on the total loss, the physical rule parameters in the physical constraint set are fixed, the gradient of the physical rule parameters is not updated, and the gradient of the physical consistency loss is backpropagated to the denoising network through the clean sample estimation module and the calculation path of the prediction noise to update the parameters of the denoising network.

8. The subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, The process involves inputting the multimodal time-window sample to be detected into the physical-data co-diffusion model for reconstruction, and determining the reconstruction anomaly score based on the difference between the multimodal time-window sample to be detected and the corresponding reconstructed sample. Specifically, this includes: According to the preset noise evaluation step, the multimodal time window sample to be detected is noise-added according to the preset reproducible noise generation rule to obtain the noisy time window sample to be detected. The noisy time window sample to be detected is input into the denoising network to obtain the predicted noise corresponding to the noisy time window sample to be detected. Based on the noise prediction parameterization relationship, the reconstructed sample is generated through a single network forward computation based on the noisy time window sample to be detected and the predicted noise. The reconstruction anomaly score is determined based on the mean squared error between the multimodal time window sample to be detected and the reconstructed sample; While determining the reconstruction anomaly score, the physical consistency of the multimodal time window sample to be detected is performed using the differentiable residual functions in the physical constraint set, and the physical violation score is determined based on the calculation results of each differentiable residual function.

9. A subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, The process of fusing the reconstructed anomaly score and the physical violation score to obtain a comprehensive anomaly score, and determining whether the subway SoS has anomalies based on the comparison result of the comprehensive anomaly score and the anomaly judgment threshold, specifically includes: A sliding window for reconstructed abnormal scores and a sliding window for physical violations, containing only historical scores under normal operating conditions, are established respectively. The corresponding score sliding windows are initialized using the reconstructed abnormal scores and physical violation scores corresponding to historical multimodal time window samples that have been confirmed to be under normal operating conditions. Based on the mean and standard deviation of each historical score within the corresponding score sliding window, the current reconstructed abnormal score and the physical violation score are standardized to obtain the data abnormality deviation and physical abnormality deviation: The fusion weight is determined based on the current operational scenario, and the data anomaly deviation and the physical anomaly deviation are fused according to the fusion weight to obtain the comprehensive anomaly score: in, This represents the current reconstruction anomaly score; and These represent the mean and standard deviation of the sliding window of scores corresponding to the reconstructed outlier scores, respectively. This represents the current physical violation score; and These represent the mean and standard deviation of the sliding window of scores corresponding to physical violation scores, respectively. This indicates the degree of abnormal deviation of the data; Indicates the degree of deviation of the physical anomaly; This represents the overall anomaly score; Indicates the fusion weights; It is a small constant that is numerically stable; The current operating scenarios include off-peak operating scenarios, peak operating scenarios, operating scenarios during severe weather, and low-density operating scenarios at night; Establish a sliding window for comprehensive anomaly scores for different operating scenarios, and determine the preset quantile of the historical comprehensive anomaly scores in the sliding window corresponding to the current operating scenario as the anomaly judgment threshold. When the overall anomaly score is greater than the anomaly determination threshold, it is determined that the subway SoS has an anomaly; when the overall anomaly score is not greater than the anomaly determination threshold, it is determined that the subway SoS does not have an anomaly, and the current reconstruction anomaly score, physical violation score and overall anomaly score are updated to the reconstruction anomaly score sliding window, the physical violation score sliding window and the overall anomaly score sliding window corresponding to the current operation scenario, respectively.

10. A subway SoS anomaly detection method based on a physical-data co-diffusion model according to claim 1, characterized in that, After determining that the subway SoS is abnormal, the following steps are also included: Based on the physical residuals corresponding to each physical constraint in the physical constraint set, determine the contribution of each physical constraint to the physical violation score, and identify the physical constraint with the highest contribution as the violated physical constraint. The reconstruction error between the multimodal time window sample to be detected and the reconstructed sample is decomposed along the feature dimension to determine the contribution of each feature variable to the reconstruction anomaly score. One or more feature variables with the highest contribution are identified as anomalous feature variables, and the subsystem and monitoring equipment corresponding to the anomalous feature variables are identified. The violated physical constraints and the abnormal feature variables are used as abnormal evidence and input into a preset causal knowledge graph. The causal knowledge graph uses subsystems, equipment, physical constraints and observations in the subway SoS as nodes, and the relationship between each node is represented by preset association edges. Based on the anomalous evidence, a correlation path matching the violated physical constraint and the anomalous feature variable is retrieved in the causal knowledge graph. Based on the correlation path, the faulty subsystem, faulty equipment, and anomalous causal chain corresponding to the subway SoS anomaly are determined, and anomaly diagnostic information containing the violated physical constraint, the anomalous feature variable, and the anomalous causal chain is generated.

Images