WO2020090770A1 - 異常検出装置、異常検出方法、およびプログラム - Google Patents

異常検出装置、異常検出方法、およびプログラム Download PDF

Info

Publication number
WO2020090770A1
WO2020090770A1 PCT/JP2019/042252 JP2019042252W WO2020090770A1 WO 2020090770 A1 WO2020090770 A1 WO 2020090770A1 JP 2019042252 W JP2019042252 W JP 2019042252W WO 2020090770 A1 WO2020090770 A1 WO 2020090770A1
Authority
WO
WIPO (PCT)
Prior art keywords
sensor
unit
mode information
abnormality
sensors
Prior art date
Application number
PCT/JP2019/042252
Other languages
English (en)
French (fr)
Inventor
誠司 堤
美樹 平林
Original Assignee
国立研究開発法人宇宙航空研究開発機構
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立研究開発法人宇宙航空研究開発機構 filed Critical 国立研究開発法人宇宙航空研究開発機構
Priority to EP19880340.5A priority Critical patent/EP3876056A4/en
Priority to US17/288,611 priority patent/US11669080B2/en
Priority to JP2020553912A priority patent/JP7340265B2/ja
Priority to CN201980071540.XA priority patent/CN112955839A/zh
Publication of WO2020090770A1 publication Critical patent/WO2020090770A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0259—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterized by the response to fault detection
    • G05B23/0283—Predictive maintenance, e.g. involving the monitoring of a system and, based on the monitoring results, taking decisions on the maintenance schedule of the monitored system; Estimating remaining useful life [RUL]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D21/00—Measuring or testing not otherwise provided for
    • G01D21/02—Measuring two or more variables by means not covered by a single other subclass
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D3/00—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups
    • G01D3/08—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups with provision for safeguarding the apparatus, e.g. against abnormal operation, against breakdown
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0221—Preprocessing measurements, e.g. data collection rate adjustment; Standardization of measurements; Time series or signal analysis, e.g. frequency analysis or wavelets; Trustworthiness of measurements; Indexes therefor; Measurements using easily measured parameters to estimate parameters difficult to measure; Virtual sensor creation; De-noising; Sensor fusion; Unconventional preprocessing inherently present in specific fault detection methods like PCA-based methods
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B23/00—Testing or monitoring of control systems or parts thereof
    • G05B23/02—Electric testing or monitoring
    • G05B23/0205—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults
    • G05B23/0218—Electric testing or monitoring by means of a monitoring system capable of detecting and responding to faults characterised by the fault detection method dealing with either existing or incipient faults
    • G05B23/0224—Process history based detection method, e.g. whereby history implies the availability of large amounts of data
    • G05B23/024—Quantitative history assessment, e.g. mathematical relationships between available data; Functions therefor; Principal component analysis [PCA]; Partial least square [PLS]; Statistical classifiers, e.g. Bayesian networks, linear regression or correlation analysis; Neural networks
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D2218/00—Indexing scheme relating to details of testing or calibration
    • G01D2218/10—Testing of sensors or measuring arrangements

Definitions

  • the present invention relates to an abnormality detection device, an abnormality detection method, and a program.
  • the present application claims priority based on Japanese Patent Application No. 2018-204515 filed in Japan on October 30, 2018, and the content thereof is incorporated herein.
  • the present invention has been made in consideration of such circumstances, and provides an abnormality detection device, an abnormality detection method, and a program capable of estimating a sensor that may be abnormal without setting a threshold value in advance.
  • One of the purposes is to provide.
  • an acquisition unit that acquires a plurality of detection results that each of the plurality of sensors has detected the state of the detection target every predetermined time, and to each of the plurality of detection results acquired by the acquisition unit
  • a plurality of mode information obtained by performing a plurality of mode information acquisition processing for extracting a characteristic according to the nature of the abnormality, of the plurality of sensors, a specific sensor with a possibility of abnormality is estimated.
  • an estimating unit for performing the abnormality detection is performed.
  • a sensor with a possibility of abnormality can be estimated without setting a threshold value in advance.
  • FIG. 6 is a diagram showing an example of the contents of detection result information 300.
  • FIG. It is a figure which shows an example of the content of the correlation information 202. It is a figure which shows an example of the content of the reference information 206. It is a figure which shows an example of a "1st type mode information acquisition process.” It is a figure which shows an example of a "2nd kind mode information acquisition process.” It is a graph which shows an example of the time-dependent change of the sensor measurement value in a normal state. It is a figure which shows an example of the shape of the orbit on the phase plane in a normal state.
  • FIG. 8 is a diagram (dendrogram) showing an example of the complex power cepstrum distance Dc classified by the dissimilarity calculation unit 108.
  • FIG. 6 is a diagram showing an example of a state transition matrix generated by an estimation unit 109. 6 is a flowchart showing a series of processing flow of the abnormality detection device 1. It is a flow chart which shows a series of flows of presumption processing of a specific sensor concerning Step S102. It is a flow chart which shows a series of flows of presumption processing of a specific sensor concerning Step S302.
  • FIG. 1 is a diagram illustrating an example of the configuration of the abnormality detection device 1 of the embodiment.
  • the anomaly detection device 1 detects an anomaly that has occurred in an anomaly detection target device or system.
  • the abnormality detection device 1 includes, for example, a control unit 100 and a storage unit 200.
  • the control unit 100 includes an acquisition unit 102, an analysis unit 104, and a derivation unit 106, for example, when a hardware processor such as a CPU (Central Processing Unit) executes a program (software) stored in the storage unit 200.
  • the dissimilarity calculation unit 108, the estimation unit 109, and the output unit 110 are realized as functional units.
  • circuits such as LSI (Large Scale Integration), ASIC (Application Specific Integrated Circuit), FPGA (Field-Programmable Gate Array), and GPU (Graphics Processing Unit). Part; including circuitry), or may be realized by cooperation of software and hardware.
  • LSI Large Scale Integration
  • ASIC Application Specific Integrated Circuit
  • FPGA Field-Programmable Gate Array
  • GPU Graphics Processing Unit
  • the storage unit 200 may be realized by, for example, a storage device (a storage device including a non-transitory storage medium) such as a HDD (Hard Disk Drive) or a flash memory, and a detachable DVD or CD-ROM. It may be realized by a storage medium (non-transitory storage medium) or may be a storage medium mounted in the drive device. Further, part or all of the storage unit 200 may be an accessible external device of the abnormality detection device 1, such as a NAS (Network Attached Storage) or an external storage server.
  • the storage unit 200 stores, for example, correlation information 202, trajectory shape information 204, and reference information 206. Details of various information will be described later.
  • the acquisition unit 102 is information indicating a plurality of detection results obtained by detecting a state of a detection target by each of a plurality of sensors (not shown) arranged in a device or a system of abnormality detection (hereinafter, detection result information). 300) is acquired.
  • the detection target is, for example, a physical quantity that causes a change in the device or system that is the abnormality detection target
  • the state of the detection target is, for example, the value of the physical quantity.
  • This physical quantity may indicate, for example, at least one of mechanical properties, electromagnetic properties, thermal properties, acoustic properties, chemical properties, etc., or spatial information and time indicated by them. It may be information.
  • This physical quantity includes, for example, pressure, flow rate, valve opening, atmospheric pressure, outside air temperature, water temperature, blood pressure, pulse, and the like.
  • the device or system for which an abnormality is to be detected is, for example, a rocket (for example, a reusable rocket engine system), a transportation system, a power generation device, a life support device, or the like.
  • the abnormality detection target system is a rocket engine system
  • the plurality of sensors are mounted on the rocket and detect a detection target (physical quantity) for knowing the state of the engine system.
  • the device or system that is the target of abnormality detection is not limited to the engine system and power generation device of these rockets, and may be any device or system that can acquire the physical quantity that is the detection target.
  • the acquisition unit 102 may be configured to acquire information necessary for detecting an abnormality or estimating the type and cause of the abnormality as necessary.
  • the acquisition unit 102 may acquire the detection result information 300 from, for example, a plurality of sensors connected to the abnormality detection device 1, and the detection result information 300 can be acquired from a connected sensor capable of transmitting and receiving information via a network.
  • the detection history of the sensor included in the operation history (for example, log information) of the abnormality detection target that has completed the operation may be acquired (extracted) as the detection result information 300.
  • the sensor is an example of a “detection unit”.
  • FIG. 2 is a diagram showing an example of the content of the detection result information 300.
  • the detection result information 300 includes, for example, information that can identify the sensor measuring the physical quantity (sensor ID shown in the figure), the sensor name, and the detection result obtained by the sensor detecting the physical quantity that is the detection target at predetermined time intervals. Information in which and are associated with each other is stored separately for each sensor.
  • the sensor with the sensor ID “0001” and the sensor name “A” is the time series information of the abnormality detection target, and the detection result is the measurement value measured (detected) by the sensor. (For example, temperature, pressure, flow rate, etc.) are shown.
  • the analysis unit 104 analyzes the detection results of the plurality of univariates acquired by the acquisition unit 102 by a matrix decomposition method (for example, principal component analysis) used for multivariate analysis, and thus the abnormal feature Acquires a plurality of mode information according to.
  • a matrix decomposition method for example, principal component analysis
  • the deriving unit 106 visualizes the mode information by deriving an orbit on the phase plane or the phase space based on the mode information acquired by the analyzing unit 104, for example.
  • the dissimilarity calculation unit 108 is dissimilar to the derivation unit 106 that calculates an index (hereinafter, dissimilarity) used when the estimation unit 109 estimates the specific sensor based on the visualization information derived by the derivation unit 106. Details of the degree calculation unit 108 will be described later.
  • the estimation unit 109 may, for example, the trajectory derived by the derivation unit 106 based on at least one of the detection result information 300 acquired from the acquisition unit 102 or the mode information acquired by the analysis unit 104, or a dissimilarity. Based on at least one of the dissimilarities calculated by the degree calculation unit, it is estimated whether the sensor detecting each detection result is a sensor showing an abnormal value (hereinafter, a specific sensor).
  • the abnormalities to be detected include various abnormalities such as changes in periodicity that do not normally occur, mixing of irregular noise, and vibration of measured values.
  • the type of abnormality includes, for example, abnormality of the device or system that is the abnormality detection target (hereinafter, system abnormality) and sensor abnormality (hereinafter, sensor abnormality).
  • the output unit 110 outputs information about the sensor estimated to be the specific sensor by the estimation unit 109.
  • the output unit 110 may display, for example, an image showing the specific sensor on a display device (not shown) connected to the abnormality detection device 1, and information indicating the specific sensor on another device connected via the network. May be output, or information indicating the specific sensor may be stored in the storage unit 200.
  • the output unit 110 includes the detection result acquired by the acquisition unit 102, the mode information acquired by the analysis unit 104, the information indicating the trajectory derived by the derivation unit 106, and the non-similarity calculated by the dissimilarity calculation unit 108. The degree of similarity may be stored in the storage unit 200.
  • the storage unit 200 may be configured separately from the abnormality detection device 1.
  • the analysis unit 104 uses, for example, the sliding window method, which is a known technique, on the detection result acquired by the acquisition unit 102 to measure univariate time series information in the vicinity of measurement values called partial time series. Convert to a set of vectors with as an element. As a result, the analysis unit 104 reconstructs the detection result as multivariate data having the measured values in the vicinity as variables.
  • the analysis unit 104 adds the time structure (lag structure) to the principal component analysis result.
  • the estimation unit 109 uses the analysis unit. Various abnormalities are determined based on each mode information acquired by 104.
  • the processing is performed by the analysis unit 104, and time series information of a single sensor is vectorized using a lag structure to reconfigure into multivariate data, and multiple mode information can be acquired.
  • a mode information acquisition process performed by using such a multivariate analysis method (for example, principal component analysis) is referred to as a “first type mode information acquisition process”.
  • it is a process performed by the analysis unit 104, and time series information of two or more sensors is collected at each time and vectorized to reconfigure into multivariate data, which is then converted into a plurality of independent multivariate data.
  • the process of acquiring mode information performed by using a general multivariate analysis method (for example, principal component analysis) in the time domain, which is executed by treating as, is referred to as “second type mode information acquisition process”.
  • the analysis unit 104 acquires mode information from a different viewpoint from the “first type mode information acquisition process” when performing multivariate analysis by collectively using two or more sensors by the “second type mode information acquisition process”. be able to.
  • the analysis unit 104 performs the “second type mode information acquisition process” using the principal component analysis, and the sliding window method in the “first type mode information acquisition process” using the principal component analysis. Unlike the case where the mode information is extracted by the number of window widths, a plurality of pieces of mode information can be acquired by the number of sensors used for reconstruction of multivariate data.
  • the analysis unit 104 acquires the mode information based on the detection results of a plurality of sensors at once, so the “first type mode information” is acquired for each sensor.
  • the number of sensors included in the abnormality detection target is larger than that in the “acquisition process”, the calculation can be saved.
  • the analysis unit 104 standardizes the detection result before the mode information acquisition process of the detection result of the single sensor shown in the detection result information 300 (hereinafter, referred to as 0th-order mode information) or the detection result of the sensor.
  • the data (hereinafter, standardized 0th order mode information) is set as the target data of the “first type mode information acquisition process”.
  • the analysis unit 104 for example, "0th-order mode information”, “standardized 0th-order mode information”, “first-to-higher-order mode information which is a result of the acquisition of the first-type mode information", and "standardized”
  • the first-to-higher-order mode information that is the result of the first-class mode information acquisition "and the like are the target data of the" second-class mode information acquisition process ". Is.
  • the analysis unit 104 performs the “second type mode information acquisition process” by using the data subjected to the “first type mode information acquisition process” as described above, thereby performing the “second type mode information acquisition process”. It is possible to allow the capture of the lag structure lost in. This is for the following reasons.
  • the lag structure is incorporated in the "first type mode information acquisition processing" data obtained by reconstructing the data and performing the multivariate analysis using the sliding window method.
  • the analysis unit 104 collects the time-series information of two or more sensors for each time and re-configures them into multivariate data. The time correlation of is lost.
  • the analysis unit 104 uses the result of the “first type mode information acquisition process”, in which the lag structure is introduced, as the time series information of the sensor in the “second type mode information acquisition process”, so that the analysis unit 104 loses it when the mode information is acquired. It becomes possible to recover the time correlation that is shown.
  • the analysis unit 104 is executed for the “first type mode information acquisition process” executed for each sensor, or for a combination of two or more sensors, based on the detection result information 300 acquired by the acquisition unit 102. Performs "type 2 mode information acquisition processing" and uses multiple mode information to create various abnormalities such as those for which thresholds are not set or those for which thresholds are difficult to set (changes in periodicity, mixing of irregular noise, etc.) Extract the features of.
  • the method by which the analysis unit 104 acquires the mode information that captures the features of various abnormalities is also referred to as (1) “mode information acquisition method for extracting features of various abnormalities”.
  • the derivation unit 106 visualizes the mode information by a known method (for example, the phase plane method) based on the plurality of mode information acquired by the analysis unit 104.
  • a known method for example, the phase plane method
  • the derivation unit 106 derives the trajectory on the phase plane or the phase space based on the plurality of mode information acquired by the analysis unit 104 and visualizes the mode information.
  • the derivation unit 106 can expand the characteristics of the change in the detection result by the derived trajectory and allow the observer of the trajectory to detect even a slight change without overlooking.
  • the deriving unit 106 allows the observer of the orbit to detect the abnormality, estimate the type and cause of the abnormality, etc. with good visibility while retaining a lot of physical information based on the derived trajectory. You can In the following description, the deriving unit 106 derives a trajectory on a phase plane or a phase space based on a plurality of pieces of mode information acquired by the analysis unit 104. It is also described as "visualization method of mode information".
  • the dissimilarity calculation unit 108 calculates the Euclidean distance, the logarithmic likelihood ratio distance, the complex power cepstrum distance, and the complex power calculated from the trajectory shape derived by the deriving unit 106 on the phase plane based on the complex autoregressive coefficient, for example.
  • the power mel cepstrum distance, dynamic time warping method, or neural network calculate the dissimilarity that digitizes the difference of the trajectory shape from the normal time, and perform clustering based on this numerical value and this dissimilarity. By carrying out, an abnormal case and a normal case are classified.
  • the estimation unit 109 uses this clustering result to estimate the presence / absence of an abnormality regardless of the threshold, and also estimates the type and cause of the abnormality.
  • the method in which the dissimilarity calculation unit 108 digitizes the trajectory shape is referred to as (3-1) “quantification method of trajectory shape generated by mode information”.
  • the dissimilarity calculation unit 108 calculates the dissimilarity between the digitized trajectory shape and the digitized trajectory shapes in a plurality of normal cases, and based on this, clustering is performed to classify abnormal cases and normal cases.
  • the method to be performed is also described as (3-2) “calculation of dissimilarity based on digitized information of trajectory shape and clustering method based on it”.
  • the estimation unit 109 refers to, for example, reference information 206 that is generated in advance, reference information 206 that stores reference information 206 that is helpful in determining an abnormality such as a sensor component replacement history, and the like.
  • a state transition matrix that integrates the analysis results
  • anomalies are detected, failures in correlation between sensors are estimated, and causes of anomalies are determined.
  • a method will be described in which the estimation unit 109 uses a state matrix in which analysis results are integrated to estimate the presence or absence of abnormality and the type and cause of abnormality based on the estimation of breakdown of correlation between sensors. , (4) "Method of integrating the calculation results of the dissimilarity calculating unit by the state transition matrix to perform detailed estimation".
  • FIG. 3 is a diagram showing an example of the content of the correlation information 202.
  • the correlation information 202 is information indicating the correlation between the sensors.
  • the correlation information 202 is, for example, information generated in advance based on the detection result of the sensor at normal time (hereinafter, normal data).
  • normal data the detection result of the sensor at normal time
  • the sensors that are correlated with each other are shown for each sensor.
  • the sensors whose sensor ID is “0001” and whose sensor name is “sensor A” and the detection results are correlated with each other have sensor IDs “0002” to “0006” and whose sensor names are “sensor B”, respectively. “Sensor C”, “Sensor D”, and “Sensor F”.
  • a sensor that is correlated with a certain sensor is a detection target of the same type (for example, “outlet voltage”) that changes according to a change of the detection target (for example, “inlet voltage”) of a certain sensor (for example, “voltage sensor”). Or a sensor (for example, “current sensor”) that detects a different type of detection target (for example, “entrance current”).
  • the correlation information 202 may be information generated in advance, or may be generated using the function of the abnormality detection device 1 by a method described later. In this case, the abnormality detection device 1 may detect the abnormality by generating the correlation information 202 and concurrently.
  • the estimation unit 109 refers to the correlation information 202 and estimates whether the sensor has an abnormal value.
  • the estimation unit 109 determines that “Sensor A” and the like based on the mode information acquired by the analysis unit 104, the trajectory derived by the derivation unit 106, the dissimilarity calculated by the dissimilarity calculation unit 108, and the like. If it is estimated that the combination of the sensors of 1 shows an abnormal value, and that the sensor correlated with the “sensor A” does not show the abnormal value of the other combination, the “sensor A” is given. Since the correlation with the sensor that should have the correlation is broken, "sensor A” is estimated to be "sensor failure". In other words, the estimation unit 109 estimates that “sensor A” is a sensor that has an abnormal value due to a sensor failure.
  • the estimation unit 109 correlates with “Sensor A” based on the mode information acquired by the analysis unit 104, the trajectory derived by the derivation unit 106, the dissimilarity calculated by the dissimilarity calculation unit 108, and the like.
  • the correlation of other related sensors is not broken, it is estimated to be "system abnormality".
  • the processing in the analysis unit 104, the derivation unit 106, the dissimilarity calculation unit 108, the estimation unit 109, and the like to only the sensors having the correlation, the time required for the abnormality detection can be shortened. it can. Details of the processing of the estimation unit 109 will be described later.
  • the estimation unit 109 estimates “system abnormality”
  • the reference information 206 is used in addition to the correlation information 202 to identify the sensor indicating the abnormal value and estimate the cause of the abnormality. Details of the above procedure of the estimation unit 109 will be described later.
  • FIG. 4 is a diagram showing an example of the contents of the reference information 206.
  • the reference information 206 includes a sensor ID that can identify the sensor, a sensor name, and a name of a component that configures the device and the system near the measurement point of the sensor (for example, if the sensor is a flow meter installed in the pipeline, Other than the maintenance history of sensor peripheral parts, such as the pipe material that constitutes the sensor), the date of replacement of the part, the manufacturer (manufacturer) of the part, etc.
  • the information necessary for estimating the cause of the abnormality detection such as the replacement date and the sensor manufacturer is information associated with each other for each sensor. Note that the content of the reference information 206 is an example, and the present invention is not limited to this, and information that is considered to be useful to describe in the reference information 206 can be described.
  • an example of the above-described abnormality estimation process is (1) “method of acquiring mode information for extracting various abnormal features", (2) “visualization method of mode information acquisition using trajectory shape", (3 -1) “Quantification method of trajectory shape generated by mode information acquisition information” (3-2) “Calculation of dissimilarity based on digitized information of trajectory shape and clustering method based on it”, (4) "Non The method of integrating the calculation results of the similarity calculation unit by the state transition matrix to perform detailed estimation ”will be described in detail.
  • the analysis unit 104 applies the sliding window method to the time series information of the single sensor acquired by the acquisition unit 102 to convert the time series information into a set of vectors, It can be reconstructed as multivariate data having nearby time information as a variable. Thereby, the analysis unit 104 can extract anomalies in a plurality of modes by the number of window widths by using the matrix decomposition method applied to the multivariate analysis.
  • the analysis unit 104 can also incorporate the lag structure into the mode information at the same time. The process of incorporating the lag structure into the mode information at the same time is the above-described “first type mode information acquisition process”.
  • FIG. 5 shows that the sliding window method with a window width of 5 is used to convert the time series information of a single sensor into a set of vectors, reconstruct it as multivariate data, and perform mode information acquisition processing by principal component analysis. It is a figure which shows an example of a "1st type mode information acquisition process.”
  • the analysis unit 104 Furthermore, in the data reconstructed by the analysis unit 104, not only the features of abnormalities having different characteristics are extracted in each mode, but also in the above-described “second type mode information acquisition process” When the detection results are collected and reconstructed into multivariate data, by using the information subjected to the "type 1 mode information acquisition process", the reconstructed multivariate data is used by ignoring the time correlation. It is possible to recover the time correlation that is lost as the multivariate analysis is performed.
  • the analysis unit 104 converts the univariate data obtained from the detection result of the single sensor into the multivariate data by introducing the lag structure in the “first type mode information acquisition process”, for example. Then, a multivariate analysis method such as a principal component analysis method, a kernel principal component analysis method, a fuzzy principal component analysis method, a sparse principal component analysis method, a stochastic principal component analysis method, or a robust principal component analysis method is applied. As a result, the analysis unit 104 acquires a plurality of pieces of mode information.
  • a multivariate analysis method such as a principal component analysis method, a kernel principal component analysis method, a fuzzy principal component analysis method, a sparse principal component analysis method, a stochastic principal component analysis method, or a robust principal component analysis method is applied.
  • the analysis unit 104 acquires a plurality of pieces of mode information.
  • the analysis unit 104 uses the sliding window method as in the above-described example in the “first type mode information acquisition process”, a plurality of pieces of mode information corresponding to the window width size m (m is a natural number) are provided. Can be obtained. Since different abnormal characteristics are extracted from each of the plurality of mode information, the abnormality detecting device 1 can perform abnormality detection without overlooking by performing abnormality detection using the plurality of mode information.
  • FIG. 6 is a diagram showing an example of “second type mode information acquisition processing”. Is shown.
  • the “second type mode information acquisition process” means that the analysis unit 104 converts the detection results of two or more sensors at the same time into a set of vectors having elements and re-converts them into an independent data set. By configuring, multi-variate analysis is used to perform feature extraction in multiple modes.
  • FIGS. 6A to 6C are three graphs showing an example of the phase surface created by the detection result before the mode information acquisition process, which is named “0th order mode information” before “sensor A” and “sensor B”. is there. As shown in FIGS. 6A to 6C, the abnormality is not clear in the 0th-order mode information.
  • FIGS. 6A to 6C the abnormality is not clear in the 0th-order mode information.
  • 6D to 6E show scores on the first main axis and the second main axis extracted by the “second type mode information acquisition process” using the principal component analysis (hereinafter, primary mode information and 2 respectively). It is a graph showing the next mode information). As shown in FIGS. 6D to 6E, it can be seen that the abnormal feature is extracted on the second main axis (secondary mode). Similarly in the "first type mode information acquisition process”, different abnormal characteristics can be extracted according to the acquired mode. In the “second type mode information acquisition process”, for example, assuming that one data set includes s sensors, different abnormal features are captured by the number of sensors s (s is a natural number). It is possible to perform feature extraction in multiple modes.
  • the abnormality for which the threshold value is set is a mode suitable for the 0th-order mode to detect the abnormality.
  • the primary mode is a mode suitable for detecting the abnormality.
  • the secondary mode is a mode suitable for detecting anomalies.
  • the anomaly detection apparatus 1 can select a mode or an extraction method to be used in advance according to the nature of the anomaly based on the characteristics of the mode information or the extraction method, without performing an exhaustive analysis. It is possible to quickly perform processing such as abnormality detection and estimation of a sign of abnormality according to the purpose.
  • the reference information 206 may include information indicating a mode suitable for detecting each abnormality, and the analysis unit 104 refers to the information and analyzes the mode suitable for detecting each abnormality. It may be configured to. By performing feature extraction in a plurality of modes in this manner, the anomaly detection apparatus 1 can perform various anomalies (abnormalities for which no threshold is set, as well as abnormalities in which it is difficult to set a threshold such as the inclusion of irregular noise). ) Can be detected without overlooking.
  • the multivariate data used by the analysis unit 104 in the “second type mode information acquisition process” is reconstructed into multivariate data by collecting a plurality of detection results (0th order mode information) acquired by the acquisition unit 102.
  • By performing the "type 2 mode information acquisition process” by performing the "type 2 mode information acquisition process” it is possible to recover the time correlation lost when performing the multivariate analysis assuming the independence of the data.
  • the analysis unit 104 may unify the mode types of all the sensors to be combined when reconstructing into multivariate data using the mode information for which the “first type mode information acquisition process” has been performed. (For example, all first-order modes), using the 0th-order mode information for each sensor, changing the mode and reconstructing into multivariate data, etc. You may.
  • the phase plane method has the effect of magnifying on the phase plane even a slight change that is difficult to distinguish by a normal plot with time on the horizontal axis and measurement values on the vertical axis. Is a well-known method.
  • the derivation unit 106 can perform not only the trajectory using the phase plane method, but also the trajectory generation by the time series plot such as the horizontal axis time and the vertical axis principal component score as shown in FIGS. 6D and 6E. Also in the trajectory based on the time-series plot, the abnormality can be detected through the dissimilarity calculation unit 108 and the estimation unit 109 by applying the detection procedure used for the trajectory using the phase plane method.
  • the generated phase plane is not limited to two dimensions. Although it is difficult to visually recognize in a high dimension, it is possible to perform anomaly detection by applying a process executed on a trajectory on a two-dimensional plane by generating a trajectory on a phase plane space of three dimensions or more.
  • the deriving unit 106 visualizes the result of the “first type mode information acquisition process” by using the phase plane method, the “first type univariate phase plane trajectory generation process” and the “first type multivariate phase”. There are two methods of "plane trajectory generation processing”. Details of these processes will be described later.
  • the deriving unit 106 uses a high-dimensional phase plane trajectory that is difficult to visualize. It may be generated and the subsequent processing may be performed.
  • the target for generating the phase plane trajectory in the derivation unit 106 is not limited to the result of the “first type mode information acquisition process”.
  • the phase plane trajectory can be generated including the 0th-order mode information before acquisition of the mode information.
  • the derivation unit 106 uses the result of performing the “first type mode information acquisition process” for the window width m, it can be used for trajectory generation by including 0th-order mode information in m pieces of mode information.
  • the s modes can be used for trajectory generation by including the 0th-order mode information. ..
  • the derivation unit 106 for example, for a single “sensor A”, the horizontal axis represents the primary mode information of the “sensor A” (score on the first principal component axis), and the vertical axis represents the secondary mode of the “sensor A”.
  • a two-dimensional phase plane is generated by a method of plotting information related to the detection result of the same sensor on both axes like information (score on the second principal component axis). In the following description, this generation method will be referred to as “first kind univariate phase plane trajectory generation processing”.
  • the above-described combination of plot targets is an example, and the present invention is not limited to this.
  • the combination can be a combination of modes in which abnormalities can be better detected. For example, when abnormalities often appear in the first-order mode and the third-order mode, the characteristics of the abnormality can be effectively visualized by plotting them in combination as in the horizontal-axis first-order mode and the vertical-axis third-order mode. Will be possible.
  • the mode suitable for estimating the anomaly differs depending on the nature of the anomaly (change of periodicity, mixing of irregular noise, vibration of measured value, etc.), and therefore the deriving unit 106 generates phase plane trajectories in a plurality of modes,
  • the estimation unit 109 uses the obtained index (for example, the degree of dissimilarity) through a process described below using the generated phase plane trajectory to estimate whether or not an abnormality is recognized.
  • the derivation unit 106 comprehensively combines the modes to generate a phase plane when the nature of the abnormality is unknown, but when the nature of the abnormality to be detected is determined, the derivation unit 106 determines the nature of the abnormality.
  • a phase plane is generated in an appropriate mode according to the above.
  • the derivation unit 106 when it is desired to detect an arrhythmia in an electrocardiogram, the derivation unit 106 generates a phase plane in a mode optimal for detecting an abnormality, depending on the type of arrhythmia. Therefore, when the nature of the abnormality to be detected is predetermined, the derivation unit 106 can shorten the time required for abnormality detection by limiting the mode information used for the generated trajectory to a specific one. it can.
  • the deriving unit 106 performs abnormality detection based on changes in a plurality of sensors.
  • the derivation unit 106 is not limited to the example described above based on the change of the single sensor, and may use the same mode information including the 0th-order mode information for all the axes, and combines the mode information as appropriate to obtain the phase. Generate a plane orbit.
  • the derivation unit 106 limits the trajectory generation to a mode suitable for the feature of the abnormality, thereby shortening the time required to detect the abnormality.
  • the deriving unit 106 plots information obtained from different sensors instead of plotting information obtained from the same sensor on both axes. Then, a two-dimensional phase plane can be generated. In the following description, this generation method will be referred to as “first kind multivariate phase plane trajectory generation processing”.
  • the deriving unit 106 also selects a combination of modes to be plotted in the “first kind multivariate phase plane trajectory generation process” as in the “first kind univariate phase plane trajectory generation process” according to the nature of the abnormality. Can effectively detect abnormalities.
  • first type multivariate phase plane trajectory generation process it is possible to detect an abnormality for which a threshold value is not set and to estimate the cause of the abnormality by detecting an abnormality based on a change between a plurality of sensors.
  • the derivation unit 106 performs the “first kind multivariate phase plane trajectory generation process” will be described in detail.
  • a method for generating a phase plane by the “first kind multivariate phase plane trajectory generation process” is as follows, for example.
  • the derivation unit 106 uses the primary mode information of “sensor A” acquired by the “type 1 mode information acquisition process” by the analysis unit 104 on the horizontal axis and the primary mode information of “sensor B” on the vertical axis. Generate the plotted phase plane.
  • the deriving unit 106 can also generate the phase plane using the 0th-order mode information before the mode information acquisition process. Further, as described above, the combination of mode information is not limited to this.
  • the mode in which the anomaly is extracted differs depending on the nature of the anomaly. Therefore, when generating the phase-plane orbit, by selecting the optimal combination of the mode information including the 0th-order mode information, the difference from the normal time becomes more significant. It clearly appears, and it is possible to efficiently detect an abnormality and estimate the type and cause of the abnormality.
  • the deriving unit 106 uses the first-order mode information of the “sensor A” and the first-order mode information of the “sensor B” extracted in the “first-type mode information acquisition process” to calculate the “first-type multivariate phase”.
  • An example of performing "plane trajectory generation processing" will be described.
  • the derivation unit 106 can exhaustively set all sensor combinations as the derivation target of the trajectory, or derivation of only the sensor group previously selected using the correlation information 202 according to the purpose of abnormality detection. Alternatively, it is possible to limit the combination of the target sensors and set the target as the trajectory derivation target.
  • the derivation unit 106 If the derivation unit 106 generates a trajectory for a sensor group having a positive correlation, the abnormality on the phase plane trajectory can be enlarged and expressed, and the abnormality detection and the abnormality type estimation can be easily performed. Can be made If the derivation unit 106 generates a trajectory for a sensor group having a negative correlation, the derivation of the anomaly can be expanded and expressed by devising a plot such as reversing one of the positive and negative signs. Can be easily detected and the type of abnormality can be easily estimated. Further, in the following description, a case will be described as an example in which the combination of the primary mode information of the “sensor A” and the primary mode information of the “sensor B” among all the combinations is the derivation target of the trajectory.
  • the first-order mode information of the “sensor A” selected as the target of the “first-type multivariate phase plane trajectory generation process” will be referred to as mode 1 result x, and the first-order mode information of the “sensor B” will be described.
  • the mode information is described as mode 1 result y.
  • the derivation unit 106 includes each of the mode 1 results x on the phase plane with the mode 1 result x as the first axis (for example, the horizontal axis) and the mode 1 result y as the second axis (for example, the vertical axis).
  • the detection timings match between the elements ⁇ x 1 , x 2 , ..., X n ⁇ (n is a natural number) and each element ⁇ y 1 , y 2 , ..., Y n ⁇ included in the mode 1 result y.
  • Coordinates indicated by the elements are plotted and connected by straight lines in time series order to generate a trajectory.
  • FIG. 7 is a graph showing an example of a change with time of the sensor in a normal state.
  • the waveform W2 shown in FIG. 7 is a waveform showing a change with time of the detection result of the "sensor A”
  • the waveform W1 is a waveform showing a change with time of the detection result of the "sensor B”.
  • FIG. 8 is a diagram showing an example of the shape of the trajectory on the phase plane in the normal state.
  • the deriving unit 106 uses the mode 1 result x derived based on the detection result of the “sensor A” of the waveform W2 as the horizontal axis, and the mode 1 result y derived based on the detection result of the “sensor B” of the waveform W1.
  • the coordinates indicated by the elements and mode 1 result y components of the mode 1 result x illustratedrated, the coordinates (x 1, y 1), ..., the coordinates (x m, y m)
  • a trajectory (trajectory Orb1 shown in the figure) that connects the points with a straight line in the order of time series is generated.
  • the derivation unit 106 describes an example in which the first-order mode information is used when deriving the orbit, but the vertical axis and the horizontal axis of the phase plane orbit are not limited to this combination.
  • the mode information used by the analysis unit 104 a combination of 0th-order mode and usable higher-order mode information can be freely selected for the vertical axis and the horizontal axis.
  • the estimating unit 109 determines whether or not an abnormal value is detected based on the dissimilarity between the trajectory shape derived by the deriving unit 106 and the trajectory shape information 204 that is information indicating the trajectory shape in a normal state. To estimate. For example, in the trajectory shape information 204 in which the trajectory shape information in the normal state is stored, similar to the detection result information 300, all sensor information acquired from a certain time to a certain time is set as one case, and a case is obtained for each acquired case. A number is added, and the result of the “first kind multivariate phase plane orbit generation process” in normal time is recorded together with the acquisition time information and the like.
  • the trajectory shape information 204 includes information indicating the trajectory shape related to the combination of “sensor A” and “sensor B” as shown in FIG.
  • the phase plane orbit in the normal state derived by using the primary mode information of the “sensor A” and the “sensor B” of the case number 1 (hereinafter referred to as the case No. 1) in the phase plane orbit shape of FIG. If it is a shape, the track shape information 204 includes the case number. 1 to case No.
  • phase plane trajectory information is stored for a plurality of cases up to k (k is a natural number).
  • the estimation unit 109 performs clustering of the trajectory shapes of the plurality of normal states and the trajectory shape of the inspection target by the dissimilarity calculation unit 108 based on the dissimilarity, and the trajectory shape of the inspection target is classified into the normal state. Based on the result of checking whether or not it is possible to detect an abnormality for which a threshold value is not set. This is because it can be determined that they do not belong to any of the normal cases if they do not resemble any of the normal cases.
  • a serial number case is described as an example here, the case taken as a normal case does not need to be a serial number, and an appropriate case can be freely selected from past normal cases. In the following, a case will be described in which a serial number of normal cases is selected.
  • the estimation unit 109 uses the case No. stored in the trajectory shape information 204. 1 to case No. From the trajectory shapes in the normal cases up to k, the trajectory shapes corresponding to the sensor combinations (“sensor A” and “sensor B” in this example) corresponding to the sensor pair of interest are extracted, and the case No. . 1 to case No. By comparing the dissimilarity with the phase plane trajectory shape created by the sensor pair corresponding to the sensor pair of interest in the normal cases up to k, it is estimated whether or not the sensor is measuring an abnormal value. As described above, the abnormality detection device 1 analyzes the abnormality using the general time series graph as shown in FIG.
  • the deriving unit 106 derives a phase plane trajectory whose change is magnified and visualized more than in a general time series graph, which causes a slight abnormality. It is possible to provide information (that is, the orbital shape) that enables detection of the sign of (3) without fail. Then, the estimation unit 109 determines that the phase pair trajectory shape of the sensor pair of interest is the case No.
  • the dissimilarity can be estimated qualitatively by the observer, but here we explain a method of digitizing the orbital shape, automating the determination, and reporting the results. Note that the degree of dissimilarity of the orbit shapes is low not only when the shapes completely match each other, but also when the deviations of the orbits are within a predetermined range or when the observers can regard them as the same. Details of the dissimilarity estimation method will be described later.
  • FIG. 9 is a diagram showing an example of a change with time of the sensor in an abnormal state.
  • a waveform W4 shown in FIG. 9 is a waveform showing a temporal change of the detection result of the "sensor A”
  • a waveform 3 is a waveform showing a temporal change of the detection result of the "sensor B”.
  • “Sensor A” is operating abnormally (that is, is a specific sensor), and as shown by the waveform W4, the cycle of change is faster and the amplitude is also smaller.
  • FIG. 10 is a diagram showing an example of the shape of the trajectory on the phase plane in the abnormal state.
  • the deriving unit 106 derives based on the detection result of the "sensor B" of the waveform W3 with the horizontal axis of the mode 1 result x derived based on the detection result of the "sensor A" of the waveform W4 shown in FIG.
  • the coordinates illustrated, coordinates (x 1 , y 1 ), ..., Coordinates (x n , Y n )
  • the shape of the trajectory Orb2 shown in FIG. 10 is different from the shape of the trajectory Orb1 shown in FIG.
  • the estimation unit 109 It can be estimated that either or both of the “sensor A” and the “sensor B” related to the trajectory Orb2 show an abnormal value. As described above, by making a comprehensive judgment by comparing with a plurality of normal cases, it is possible to detect anomalies that do not depend on the threshold.
  • the estimation unit 109 digitizes the phase plane trajectory shape by the dissimilarity calculation unit 108, and based on the degree of dissimilarity to the numerical values obtained from the phase plane trajectory shapes of a plurality of past normal cases, whether the result of clustering is abnormal or not. Presume.
  • the estimation unit 109 compares the result of multivariate analysis (in this case, the analysis based on the change of “sensor A” and “sensor B”) using the phase mode trajectories of a plurality of modes with a plurality of normal data, By estimating whether or not an abnormality is detected, even if it is difficult to estimate whether or not it is abnormal by analyzing the time-series change of a single sensor because the threshold value is not set or the threshold value cannot be set, it is possible to detect the abnormality. It will be possible. Further, since it becomes possible to take out a large change due to anomaly by acquiring the mode information and further magnify and visualize the change by the phase plane, the estimating unit 109 may detect a slight sign of the abnormality without missing it. It will be possible.
  • the derivation unit 106 can provide detailed analysis materials that enable not only abnormality detection but also failure prediction by analysis using phase mode trajectories of multiple modes.
  • the mode information used in the estimation unit 109 for estimating the abnormality is a combination of 0th-order mode information to available higher-order mode information.
  • the estimation unit 109 may perform standardization processing that aligns the average and the variance before the high-order mode extraction processing (or after the mode information acquisition processing, if necessary), depending on the nature of the abnormality. The detection sensitivity can be increased.
  • the derivation unit 106 generates a phase plane trajectory based on the “first type mode information acquisition processing” result, but the information used to generate the phase plane trajectory is not limited to this.
  • the derivation unit 106 generates a phase plane trajectory using an appropriate mode according to the nature of the anomaly, such as a change in periodicity or the inclusion of irregular noise, because the mode in which the anomaly is extracted differs depending on the nature of the anomaly.
  • the dissimilarity calculation unit 108 needs to calculate the dissimilarity based on the trajectory generated by the derivation unit 106.
  • the dissimilarity calculation unit 108 converts the shape of the trajectory derived by the derivation unit 106 into numerical information by an image extraction method that uses a complex autoregressive coefficient for each trajectory shape. In the following description, a method in which the dissimilarity calculation unit 108 calculates the dissimilarity to the normal state by using the complex power cepstrum distance calculated from the complex autoregressive coefficient will be described.
  • the dissimilarity calculation unit 108 uses the information obtained by digitizing the trajectory shape by the complex autoregressive coefficient, and the information obtained by similarly digitizing the trajectory shapes of a plurality of normal states, between the two complex autoregressive coefficients.
  • the complex power cepstrum distance which is the difference (distance) of, is calculated. Further, the dissimilarity calculation unit 108 performs clustering based on the calculated complex power cepstrum distance to classify a normal case and an abnormal case.
  • the complex power cepstrum distance is an example of “calculation of dissimilarity between phase plane trajectory shapes”, and the method of “calculation of dissimilarity between phase plane trajectory shapes” is not limited to this.
  • the dissimilarity calculation unit 108 can calculate the dissimilarity using a method such as Euclidean distance, log likelihood ratio distance, complex power cepstrum distance, complex power mel cepstrum distance.
  • the dissimilarity calculation unit 108 calculates the details of the method of extracting the shape of the trajectory from the phase plane trajectory derived by the derivation unit 106 in digitizing the phase plane trajectory shape, based on the complex autoregressive coefficient.
  • the complex power cepstrum distance is shown below as an example.
  • the dissimilarity calculation unit 108 converts the trajectory shown on the phase plane derived by the derivation unit 106 into an image of a predetermined size (or a predetermined number of pixels), and converts the converted image into grayscale.
  • the shape (contour line (edge)) of the trajectory is extracted based on the brightness (pixel value) of each pixel in the case.
  • the dissimilarity calculation unit 108 sets the point sequence obtained by tracing the contour line of the extracted trajectory as Expression (1),
  • the m-th order complex autoregressive coefficient is derived when the expression is represented by the expression (2).
  • the m-th order complex autoregressive coefficient is defined as a model that approximates contour points by a linear combination of up to m contour points, as shown in Expression (3).
  • the dissimilarity calculation unit 108 extracts trajectory shape information by the above method for each phase plane trajectory generated by the derivation unit 106, and converts the extracted trajectory shape information into a complex autoregressive coefficient. Further, the dissimilarity calculation unit 108 uses the complex autoregressive coefficient z1 of the converted orbit and the corresponding plurality of normal phase plane orbit information stored in the orbit shape information 204 in the same procedure as above. The autoregressive coefficient is derived, and the complex power cepstrum distance Dc (z1, z2) with the complex autoregressive coefficient z2 of the orbit at the normal time is derived.
  • the complex power cepstrum distance Dc (z1, z2) is represented by the equation (4) as the mean square distance of the logarithm of the spectrum envelope between the complex autoregressive coefficient z1 and the complex autoregressive coefficient z2.
  • the dissimilarity calculation unit 108 calculates the complex power cepstrum distance Dc using the converted information. ..
  • the complex power cepstrum distance is calculated by comprehensively combining two of the complex autoregressive coefficients to be detected and a plurality of normal complex autoregressive coefficients added.
  • the normal phase plane trajectory information created by “Sensor A” and “Sensor B” is the case number. 1 to case No.
  • the dissimilarity calculation unit 108 determines that the complex autoregressive coefficients z1 to z18 derived from the normal phase plane trajectory shape and the “sensor A” to be inspected acquired by the acquisition unit 102.
  • the complex power cepstrum distance Dc is calculated by extracting two each from the total of 19 complex autoregressive coefficients z19 based on the phase plane trajectory created by the "sensor B".
  • the dissimilarity calculation unit 108 classifies the complex power cepstrum distance Dc into an abnormal case and a normal case by hierarchical clustering.
  • FIG. 11 is a diagram showing an example of the complex power cepstrum distance Dc classified by the dissimilarity calculation unit 108.
  • FIG. 11 is a diagram called a dendrogram showing the result of clustering all 1 to 19 cases including the detection target based on the complex power cepstrum distance Dc calculated by the dissimilarity calculation unit 108.
  • the complex power cepstrum distance Dc is shown by the length of the branch of the dendrogram.
  • the estimation unit 109 calculates “case No. 19” through the analysis unit 104, the derivation unit 106, and the dissimilarity calculation unit 108 based on the information obtained from the acquisition unit 102.
  • the complex autoregressive coefficient of the phase-plane orbit shape named Case and No. 1 to case No.
  • the dissimilarity calculation unit 108 calculates the Euclidean distance, the logarithmic likelihood ratio distance, the complex power cepstrum distance, or the complex power mel cepstrum distance based on the complex autoregressive coefficient obtained by converting the trajectory shape into a numerical value. It may be configured to calculate.
  • the Euclidean distance, the log likelihood ratio distance, the complex power cepstrum distance, and the complex power mel cepstrum distance are examples of “differences between complex autoregressive coefficients”.
  • the model used by the dissimilarity calculation unit 108 for digitizing the shape of the trajectory is not limited to the method using the complex autoregressive coefficient.
  • the method of calculating the dissimilarity is not limited to the above.
  • the method that uses the complex autoregressive coefficient is a method that extracts only the outline of a figure, so if more detailed discrimination is required, for example, a non-linear shape of a phase plane orbit such as the dynamic time warping method is used.
  • a method such as a neural network capable of extracting the internal shape can be used. Based on this result, the dissimilarity can be calculated as in the case of using the complex autoregressive coefficient.
  • the estimation unit 109 calculates the difference between the numerical value obtained from the abnormality detection target data and the numerical value obtained from the normal data based on the numerical value obtained by quantifying the trajectory shape calculated by the dissimilarity calculation unit 108 (an example of this). Then, the complex power cepstrum distance Dc) is larger than any difference calculated between certain normal data and other normal data, and as shown in FIG. 11, the abnormality detection target data and the normal data group are different groups. When clustered as, it is determined that the dissimilarity between the two is large, and one or both of the sensors (“sensor A” and “sensor B” in this example) related to the abnormality detection target data indicate an abnormal value. Presumed to be a specific sensor.
  • the mode information acquisition method can expand and extract anomalies as shown in FIG. 6E, by plotting this mode information acquisition information on a phase plane that has the effect of expanding the trajectory, the time is plotted on the horizontal axis, It is possible to detect slight changes that are difficult to distinguish by a normal plot in which measured values are plotted on the vertical axis.
  • the estimation unit 109 can detect the occurrence of an abnormality at the sign stage.
  • what the estimation unit 109 estimates as the specific sensor is a combination of sensors.
  • the estimating unit 109 can detect the abnormality of the single sensor instead of the abnormality of the combination of the sensors by the same method. A method for identifying the sensor indicating the abnormal value when the estimation unit 109 estimates that the sensor combination of the abnormality detection target data indicates the abnormal value will be described later.
  • the estimation unit 109 integrates the calculation results of the dissimilarity calculation unit 108 using the state transition matrix to identify the sensor indicating an abnormality and determine whether the sensor abnormality or the system abnormality. To do.
  • the estimation unit 109 generates a state transition matrix and arranges the results obtained by estimating whether or not the combination of the detection target sensors indicates an abnormal value by using the result of the clustering based on the dissimilarity. It is possible to perform a detailed analysis such as identifying a sensor indicating an abnormal value and estimating the type and cause of the abnormality.
  • the state transition matrix is obtained by writing the abnormality estimation result in association with the sensors as the vertical axis and the horizontal axis, respectively. The method of generating the state transition matrix will be described below.
  • FIG. 12 is a diagram showing an example of a state transition matrix generated by the estimation unit 109.
  • the state transition matrix is vertically symmetrical with respect to the diagonal elements whose elements shown on the ordinate and the abscissa coincide with each other, and in this embodiment, a combination of the ordinate "sensor A” and the reading axis "sensor B". Since the combination of the vertical axis “sensor B” and the horizontal axis “sensor A” is not distinguished, the lower half including the diagonal element may be generated as shown in FIG. As illustrated in FIG.
  • the estimation unit 109 includes, among the elements of the state transition matrix, an element having the same sensor name on the vertical axis and the horizontal axis (for example, the vertical axis “sensor A” and the horizontal axis “sensor A”).
  • the element of “” is an element indicating whether or not each sensor measures an abnormal value.
  • elements having the same sensor name will be described as "sensor name on the horizontal axis-sensor name on the vertical axis", such as “sensor A-sensor A” ... "Sensor F-sensor F”.
  • the element with which the sensor name matches is marked with a circle, and if the sensor indicates an abnormal value that exceeds the threshold, "X" is added to the element.
  • a threshold value for example, using the result of multivariate analysis focusing on the change between two sensors by the method described below. It estimates whether or not each sensor is measuring an abnormal value.
  • a state transition matrix for six sensors “sensor A” to “sensor F” is generated as shown in FIG. It is assumed that the only sensors that have a correlation are “sensor A” and “sensor B”.
  • the estimating unit 109 identifies the sensor indicating an abnormal value, fills the blanks, and performs detailed analysis such as estimation of the type of abnormality such as sensor failure or system abnormality and the cause of the abnormality. ), (2), (3-1), and (3-2) based on the results obtained by the procedure, it is possible to organize the contents estimated by the estimation unit 109 in the state transition matrix as write information. For example, the content of the estimation unit 109 estimating the presence / absence of abnormality from the result of clustering indicates that among the elements of the state transition matrix, the elements whose sensor names shown on the vertical axis and the horizontal axis do not match (for the area AR2 shown in the figure). Element).
  • the estimation unit 109 based on the result of the clustering performed by the dissimilarity calculation unit 108, a group of the phase plane trajectories and the phase plane trajectories in the normal state generated by the pair of the sensor shown on the vertical axis and the sensor shown on the horizontal axis. If the dissimilarity to the phase plane trajectory shape generated by the sensor pair of is low and both are clustered in the same group, it is estimated that the detection result of this sensor pair does not indicate an abnormal value.
  • the estimating unit 109 repeats the above processing until the blanks in the elements whose sensor names do not match on the state transition matrix are filled.
  • identify the sensor that shows an abnormal value and fill in all the blanks of the element where the sensor name of the state transition matrix matches with the following procedure to complete it. ..
  • all elements that combine "sensor B” and other sensors show "x" (different from normal times), and all other combinations show " ⁇ " (no change from normal times). )
  • a method of filling a blank space of an element having a matching sensor name will be described.
  • the abnormality is “sensor B” itself, and it is concluded that the other sensors do not show the abnormal value. It is presumed to be a failure of "sensor B". Based on this estimation, among the elements of the state transition matrix, "x”, “sensor A-sensor” in the “sensor B-sensor B” column, which is an element with the same sensor name on the vertical axis and the horizontal axis. Write “ ⁇ ” in the "A” column.
  • the following example describes the case where the sensor is not faulty. Similar to the above, consider the case where the state transition matrixes for the six sensors “sensor A” to “sensor F” are generated. Similar to the above, when the threshold value is not set and it is difficult to detect an abnormality by the sensor alone, and among the elements of the state transition matrix, the elements whose sensor names shown on the vertical axis and the horizontal axis match ( "Sensor A-Sensor A" to "Sensor F-Sensor F”) are blank. Further, similarly to the above, based on the correlation information 202, the only sensor having the correlation with the “sensor A” is the “sensor B”, and the sensors having the correlation with each other are the “sensor A” and the “sensor B”.
  • the estimation unit 109 writes the results of the processing of (1) to (3) in the state transition matrix.
  • the estimation unit 109 writes the results of the processing of (1) to (3) in the state transition matrix.
  • all the elements of the combination of the four sensors (sensors C to F) except “sensor A” and “sensor B” are “ ⁇ ” and the other elements are “x”.
  • the "sensor A” having a correlation with the "sensor B” also shows an abnormal value, no causal failure is observed, and the "sensor A” and the "sensor B” having a correlation with each other are not detected.
  • the cause of the abnormality can be estimated using the completed state transition matrix by the method described below.
  • the system abnormality is estimated by the above procedure, the cause of the abnormality can be estimated using the completed state transition matrix by the method described below.
  • both “sensor A” and “sensor B” are estimated to have abnormal values.
  • the “sensor B” changes due to the change of the “sensor A”. It is estimated that the cause is around “Sensor A”.
  • the reference information 206 is also examined, and if there is a maintenance history such that the part has not been replaced around "Sensor A" for a long time and the replacement time is approaching, it is highly likely that the deterioration of this part is the cause. Since it is estimated, it can be determined that the inspection of this component should be performed first based on the result of the abnormality detection. Similarly, when there are a plurality of sensors having direct correlation or indirect correlation, the cause of the abnormality can be estimated from the state transition matrix using the correlation information 202, the reference information 206 and the like.
  • FIG. 13 is a flowchart showing a series of flow of processing of the abnormality detection device 1.
  • the acquisition unit 102 acquires the detection result information 300 from the log information of the abnormality detection target (step S100).
  • the estimation unit 109 specifies an abnormal value based on the analysis result of the detection result information 300 acquired by the acquisition unit 102 via the analysis unit 104, the derivation unit 106, and the dissimilarity calculation unit 108.
  • the sensor is estimated (step S102). Details of the process of step S102 will be described later.
  • the output unit 110 outputs information regarding the sensor estimated to be the specific sensor by the estimation unit 109 (step S104).
  • the mode information acquisition process for extracting the features of various abnormalities is performed, the mode information acquisition information is visualized using the trajectory shape, and the visualized trajectory shape is quantified.
  • a process of calculating a dissimilarity based on the above, and generating a state transition matrix using the result of clustering based on the calculated dissimilarity (2) In the order of detailed determination using the state transition matrix, the process of step S102 An example of a method of executing will be described.
  • FIG. 14 is a flowchart showing a series of flow of the estimation process of the specific sensor according to step S102.
  • the analysis unit 104 extracts an abnormal feature from the time series information obtained from each sensor by the “first type mode information acquisition process”, and the derivation unit 106 performs the “mode information acquisition process”.
  • a two-dimensional phase plane trajectory is generated by selecting an appropriate combination from a plurality of modes including a zero-order mode, and the trajectory shape is digitized by the dissimilarity calculation unit 108.
  • the estimation unit 109 generates a state transition matrix using the result of clustering based on the dissimilarity calculated as above.
  • FIG. 13 a case where a threshold value is not set and it is difficult to estimate whether or not the sensor shows an abnormal value from the time series information of a single sensor will be described.
  • the threshold value is not set, the abnormality estimation of the diagonal element in which the elements indicated by the vertical axis and the horizontal axis match is performed after the completion of the state transition matrix of the portion below the diagonal element.
  • the state transition matrix is vertically symmetrical with respect to the diagonal elements where the elements shown on the vertical axis and the horizontal axis match, and in this embodiment, "sensor A-sensor B" and “sensor B” are used. Since “B-sensor A” is not distinguished, only the lower half including the diagonal elements need be generated as shown in FIG. Further, as shown in the above embodiment, assuming that the abnormality estimation is performed on the measurement values of the six sensors from “sensor A” to “sensor F” as in FIG. The following will be described. For example, in order to perform the “first type mode information acquisition process”, the analysis unit 104 converts the detection result information 300 acquired by the acquisition unit 102 into a time vector for each sensor using the sliding window method (step S300). ..
  • the analysis unit 104 considers that the data sets reconstructed into multivariate data by this conversion are independent of each other, and uses the multivariate analysis method using the matrix decomposition method for each sensor to identify various abnormal features. It is extracted as a plurality of mode information (step S302).
  • the derivation unit 106 uses the multiple mode information for each sensor acquired by the analysis unit 104 to select a mode to be used on each of the vertical axis and the horizontal axis to derive the phase plane trajectory, and the vertical axis, And two abscissa sensors are combined to derive a phase plane trajectory representing a change between the two sensors on a two-dimensional phase plane (step S304).
  • the phase plane trajectory derived by the deriving unit 106 has, for example, the mode 1 result x of the “sensor A” as the first axis (eg, the horizontal axis) and the mode 1 result y of the “sensor B” as the second axis (eg, the vertical axis). Axis) and the phase plane.
  • the deriving unit 106 outputs the elements ⁇ x 1 , x 2 , ..., X n ⁇ included in the mode 1 result x and the elements ⁇ y 1 , y 2 , ..., Y n ⁇ included in the mode 1 result y.
  • the trajectory is derived by connecting the coordinates (coordinates (x 1 , y 1 ), ..., Coordinates (x n , y n )) indicated by the elements whose detection timing coincides with each other by a straight line in the order of detection time.
  • the dissimilarity calculation unit 108 digitizes the derived phase plane trajectory shape by a quantitative method such as a complex autoregressive coefficient (step S306).
  • the dissimilarity calculation unit 108 is a dissimilarity between the quantified phase surface trajectory shape and the plurality of quantified normal trajectory shapes shown in the trajectory shape information 204 (for example, complex power cepstrum distance Dc). Is calculated, the measurement results are clustered according to the size of the distance from the normal information, and the clustering result is visualized using a visualization method such as the dendrogram shown in FIG. 11 (step S308).
  • the estimating unit 109 determines that the degree of dissimilarity to the normal case is large, and either of the sensor pair or When both are estimated to measure abnormal values, and when they are classified into the same group as the group created by the normal case, it is considered that the degree of dissimilarity to the normal case is small and neither of the sensor pairs has measured the abnormal value. Estimate (step S310).
  • the state transition matrix is a matrix in which the target sensor pair is associated with the vertical axis and the horizontal axis and the abnormality estimation result is written.
  • the corresponding part of the state transition matrix “X” is entered in (elements of “sensor A” on the horizontal axis and “sensor B” on the vertical axis, columns of “sensor A ⁇ sensor B”) (step S312).
  • "o" is entered in the relevant part of the state transition matrix (step S314).
  • the estimating unit 109 determines whether or not the state transition matrix generated in steps S312 and S314 is completed (step S316).
  • the estimation unit 109 repeats the processing of steps S304 to S316 until the state transition matrix is completed.
  • the person who performs the abnormality detection selects a plurality of modes to be used for creating the state transition matrix in advance when performing the abnormality determination, and the estimation unit 109 uses the “sensor” shown in FIG.
  • a method of generating a state transition matrix created by "Sensor F" from "A" will be described. In this case, one state transition matrix is generated for each selected mode.
  • the analysis unit 104 determines whether the state transition matrix is completed in all the selected modes (step S318).
  • the estimation unit 109 repeats the processing of steps S304 to S318 until the state transition matrix is completed in all required modes.
  • the state transition matrix is completed when the entry below the diagonal elements is completed. It is determined that the matrix is completed once, and in the next step S320, the state transition matrix is used to perform abnormality estimation of the diagonal elements having the same sensor names on the vertical axis and the horizontal axis.
  • the threshold is not set, and it is difficult to determine from the time series information of a single sensor whether or not it shows an abnormal value, and the pair of sensors whose sensor names shown on the vertical axis and the horizontal axis match. A case where the corner element is blank will be described. If the abnormality to be detected is determined in advance and the optimum mode to be checked is known in advance, the state transition matrix should be prepared only in the required mode.
  • the estimation unit 109 estimates the element having the same sensor name on the vertical axis and the horizontal axis (for example, the element of “sensor A-sensor A” defined above) as a threshold value. If is set, it is estimated that it is an abnormal value (“ ⁇ ”) if it exceeds the threshold value, and it is not an abnormal value (“ ⁇ ”) if it does not exceed the threshold value, and the result is displayed in the corresponding part of the state transition matrix. Should be entered.
  • the estimation unit 109 does not match the sensor names in the state transition matrix (for example, “sensor A -Sensor B "," Sensor A-Sensor C “, etc.) are all filled, and then the sensor showing an abnormal value is specified together with the following comprehensive judgment of sensor failure or system abnormality.
  • the estimation unit 109 uses the state estimation matrix in which the elements whose sensor names do not match on the vertical axis and the horizontal axis have been entered, and uses the state estimation matrix to store information about the correlation between the sensors stored in the correlation information 202 and the sensor. Based on reference information 206 that stores information such as the maintenance history of components and their peripheral parts, etc., the identification of the sensor indicating an abnormal value (abnormality estimation of the element whose sensor names shown on the vertical and horizontal axes match) ), A sensor failure or a system abnormality is determined, and if a system failure is determined, the cause is estimated (step S320).
  • FIG. 15 is a flowchart showing a series of flow of the estimation process of the specific sensor according to step S320.
  • the correlation information 202 the case where it is known that the sensors having the correlation among them are only the “sensor A” and the “sensor B” will be described as an example. In this example, the case where a plurality of sensor failures do not occur will be described.
  • the estimation unit 109 acquires the state transition matrix information in which the elements whose sensor names do not match have been entered (step S400), and checks whether or not there is a sensor pair with “x” (step S402). If there is no sensor pair with "x”, the estimation unit 109 estimates that there is no abnormality (step S428), and proceeds to the determination of whether or not the state transition matrix of all modes selected in advance has been verified ( Step S420). The procedure after step S420 will be described together with the following steps. If there is a sensor pair with "x”, select one of them (for example, the sensor pair of "sensor A” and “sensor B” from the elements of the vertical axis "sensor A” and the horizontal axis "sensor B”) ( Step S404).
  • One of the selected sensor pairs (for example, “sensor B”) is set as the estimation target sensor (hereinafter referred to as the first sensor) for estimating whether or not it is the specific sensor (step S406).
  • the first sensor is set as the estimation target sensor (hereinafter referred to as the first sensor) for estimating whether or not it is the specific sensor (step S406).
  • the first sensor is set as the estimation target sensor (hereinafter referred to as the first sensor) for estimating whether or not it is the specific sensor.
  • the estimation unit 109 checks whether or not the sensor pair that does not include the first sensor also has “x” (step S408). If the sensor pair that does not include the first sensor has “x”, the correlation information 202 is displayed. On the basis of the above, the sensor having a correlation with the first sensor (“sensor A” in this example) is checked, and all the columns of the sensor pair formed by this sensor with the sensor other than the first sensor are all on the state transition matrix just in case. It is checked whether or not it is “ ⁇ ” (step S410). If all are “ ⁇ ”, the correlation (or the causal relationship) that should exist is lost, and therefore the estimation unit 109 determines that the first sensor is a “sensor failure” (step S412).
  • the estimation unit 109 puts the determination result in “X” in the corresponding column of the state transition matrix (when the first sensor is “sensor B”, the column of “sensor B-sensor B” in the state transition matrix). Is added to indicate that the sensor is faulty (step S414).
  • a different code for example, “#”, may be used to distinguish from a system abnormality around the sensor and indicate that the sensor has a failure. If all of them are not “ ⁇ ” (this should not be the case because there is only “x” in the sensor pair including the first sensor in S408), the result is NO in S408, and therefore the flow merges into step S424. And perform a detailed analysis.
  • step S414 it is confirmed whether both sensors of the sensor pair selected in S404 have been verified as the first sensors. If both sensors have not been verified as the first sensors, the process returns to step S406 for verification. to continue. If the verification has been completed for both sensors, the process returns to S404, and another sensor pair marked with "x" is selected to continue the verification. If all of this is also completed, the process proceeds to step S420, it is determined whether the verification is completed in all the modes selected in advance, and if not completed, the process returns to S400.
  • step S422 the output unit 110 outputs the completed state transition matrix and the determination result.
  • the single sensor failure has been described as an example, but the estimation unit 109 performs the estimation in the same procedure even when a plurality of sensors have failed.
  • the output unit 110 may output all the information related to the sensor estimated to be the specific sensor, or may output the information related to some of the sensors estimated to be the specific sensor. May be.
  • Step S408 If “x” is attached to the sensor pair other than the sensor pair including the first sensor in step S408, it is checked based on the correlation information 202 whether or not a sensor having a correlation with the first sensor is included therein. (Step S424). In the column of the sensor pair formed by the sensor having a correlation with the first sensor and the sensor other than the first sensor, if the “x” is added, the breakdown of the correlation between the sensors is not recognized, and thus the sensor is not detected.
  • step S414 the column of the first sensor of the diagonal element of the state transition matrix (if the first sensor is "sensor A”, the vertical axis is A “x” is added to the "sensor A-sensor A” column in which the sensor names on the horizontal axis match, and a system failure is recorded instead of a sensor failure, and the processing from step S416 is performed in the same manner as described above. Run.
  • the cause of sensor failure is checked with the actual sensor for the sensor that is estimated to be sensor failure during inspection work, and sensor parts may be used as necessary. And the sensor itself are dealt with.
  • system abnormality the state transition matrix is used to identify the sensor that has detected the abnormal value, and the correlation information 202 between the sensors and the devices and systems around the sensor that indicate the abnormal value.
  • the cause of abnormality can be comprehensively estimated based on the reference information 206 in which the maintenance history and the like are stored. For example, if the system is diagnosed as abnormal and the pressure sensor and the flow meter for measuring internal leakage show abnormal values, and it is found from the reference information 206 that the replacement deadline for the seal parts is approaching, the seal wears out. Since internal leakage due to suspicion is suspected, the system can be efficiently maintained by proceeding with the inspection work to identify the faulty part from the inspection of the seal part.
  • step S102 of estimating the specific sensor (1) “mode information acquisition for extracting features of various abnormalities” in the analysis unit 104 (2) “mode information acquisition information using trajectory shape” Visualization ”(3-1)“ Quantification of trajectory shape generated by mode information acquisition information ”(3-2)“ Calculation of dissimilarity based on digitized information of trajectory shape and clustering based on it ”(4)
  • the flow chart in the case of sequentially performing the “detailed estimation processing in which the calculation result of the dissimilarity calculation unit is integrated by the state transition matrix and executed” has been described, but the embodiment of the present invention is not limited to this.
  • step S102 among the procedures (1) to (4), an appropriate procedure can be combined and the flowchart can be executed according to the purpose of abnormality detection.
  • the method described here as an example was a method suitable for analyzing recorded information offline, for example, if an emergency stop due to quick abnormality detection online is the objective, the target failure mode By limiting the number of sensor pairs to be analyzed, it is possible to shorten the time required for abnormality detection. If the purpose is to detect anomalies and estimate their causes for efficient maintenance, or to detect the signs of anomalies if detection that competes for 0.01 seconds is not required, it is comprehensive even if it takes some time. It is desirable to use the method described above. Further, before the mode information acquisition process, a process other than the mode information acquisition process such as standardization for aligning the average and variance of the detection results may be added, or such a process may be performed after the mode information acquisition process.
  • the mode information acquisition process also includes the “first type mode information acquisition process”, the “second type mode information acquisition process”, and a combination thereof, and any mode can be used for deriving the phase plane trajectory. It also states that the trajectory can be derived by combining arbitrary information depending on the information selected, and that the dissimilarity can be calculated by similar processing from the trajectory created by normal time series information that does not use the phase plane.
  • FIG. 11 shows the result of determination of the dissimilarity of the phase plane trajectories formed by “Sensor A” and “Sensor B” among a plurality of normal cases.
  • the phase-plane orbital shapes are generated for the combinations of all sensors in the normal case, and the correlation is calculated from the physical relational expression.
  • the correlation matrix between the sensors can be generated using the abnormality detection device.
  • the correlation table can be generated in parallel with the abnormality determination, but may be executed in advance using normal data and recorded in the correlation information 202.
  • the data obtained by converting the time series data for each sliding window into a set of vectors which is named “type 1 mode information acquisition processing”, is a multivariate having time information in the vicinity as a variable.
  • type 1 mode information acquisition processing is a multivariate having time information in the vicinity as a variable.
  • various types of anomalies including "type 2 mode information acquisition processing" have been described.
  • the method of extracting mode information is not limited to principal component analysis.
  • the derivation unit 106 is a matrix factorization method capable of acquiring mode information similar to principal component analysis such as factor analysis, sparse principal component analysis, fuzzy principal component analysis, kernel principal component analysis, stochastic principal component analysis, and robust principal component analysis.
  • An abnormal feature extraction may be performed by applying a method capable of extracting features in a plurality of modes by multivariate analysis using. Further, as described in the embodiments, the result of "type 1 mode information acquisition" is used for the known method of principal component analysis and factor analysis in the time domain named "type 2 mode information acquisition". Since the lag structure can be introduced, it is also effective in detecting anomalies when the measured values of the sensors interact with each other. By selecting a mode suitable for the nature of the abnormality to be detected, such as the frequency change and the abnormality that occurs with a lag, overlooking and false detection can be reduced, and accurate abnormality detection can be performed.
  • a lag structure is introduced to reconstruct multivariate data, and a series of processes for applying principal component analysis produces an effect like a noise filter.
  • the analysis unit 104, the derivation unit 106, the dissimilarity calculation unit 108, and the estimation unit 109 estimate the specific sensor according to the procedure of (1) to (4). Is not limited to this.
  • the analysis unit 104, the derivation unit 106, the dissimilarity calculation unit 108, and the estimation unit 109 acquire information such as the detection result information 300, they measure an abnormal value having a high dissimilarity with the detection result at the normal time.
  • the specific sensor may be estimated based on a learning model learned to output the existing sensor as the specific sensor.
  • the estimation unit 109 inputs the detection result indicated by the detection result information 300 to the learning model and estimates the sensor obtained as the output as the specific sensor.
  • DESCRIPTION OF SYMBOLS 1 ... Anomaly detection device, 100 ... Control part, 102 ... Acquisition part, 104 ... Analysis part, 106 ... Derivation part, 108 ... Dissimilarity degree calculation part, 109 ... Estimating part, 110 ... Output part, 200 ... Storage part, 202 ... Correlation information, 204 ... Orbit shape information, 206 ... Reference information, 300 ... Detection result information, Dc ... Complex power cepstrum distance, Orb1 ... Orbit, Orb2 ... Orbit, x, y ... Mode 1 result, z1, z2 ... Complex Autoregressive coefficient, z1, z2 ... Complex power cepstrum distance

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Testing And Monitoring For Control Systems (AREA)
  • Testing Or Calibration Of Command Recording Devices (AREA)
  • Information Retrieval, Db Structures And Fs Structures Therefor (AREA)

Abstract

異常検出装置は、複数のセンサの各々が検出対象の状態を所定時間毎に検出した複数の検出結果を取得する取得部と、前記取得部によって取得された複数の前記検出結果から得られた異常の性質に応じた多様な特徴抽出を可能にする複数モード情報取得処理データに基づいて、前記複数のセンサのうち、特定センサを推定する推定部と、を備える。

Description

異常検出装置、異常検出方法、およびプログラム
　本発明は、異常検出装置、異常検出方法、およびプログラムに関する。
　本願は、２０１８年１０月３０日に、日本に出願された特願２０１８－２０４５１５号に基づき優先権を主張し、その内容をここに援用する。
　従来、各種システム（原子炉、航空機のエンジン、鉄道網、生命維持装置等）において、異常がもたらすダメージを最小限に抑えるために、システムに設けられたセンサの検出結果を取得し、取得した検出結果が示す値が、予め設定した閾値を超えた場合に異常として検出する手法がある（例えば、特許文献１）。
特開２０１７－１５６２０１号公報
　しかしながら、従来の技術では、運転条件や環境の変化に対応して閾値が変化することなどから、異常を検出することが難しい場合があった。また、閾値の設定が困難である場合には、異常の発生を予兆の段階で判定することが難しかった。
　本発明は、このような事情を考慮してなされたものであり、閾値を予め設定することなく、異常の可能性があるセンサを推定することができる異常検出装置、異常検出方法、およびプログラムを提供することを目的の一つとする。
　本発明の一態様は、複数のセンサの各々が検出対象の状態を所定時間毎に検出した複数の検出結果を取得する取得部と、前記取得部によって取得された複数の前記検出結果のそれぞれに対して、異常の性質に応じた特徴を抽出する複数のモード情報取得処理を行って得られた複数のモード情報に基づいて、前記複数のセンサのうち、異常の可能性がある特定センサを推定する推定部と、を備える異常検出装置である。
　本発明によれば、閾値を予め設定することなく、異常の可能性があるセンサを推定することができる。
実施形態の異常検出装置１の構成の一例を示す図である。 検出結果情報３００の内容の一例を示す図である。 相関関係情報２０２の内容の一例を示す図である。 参考情報２０６の内容の一例を示す図である。 「第１種モード情報取得処理」の一例を示す図である。 「第２種モード情報取得処理」の一例を示す図である。 通常状態におけるセンサ計測値の経時変化の一例を示すグラフである。 通常状態における位相平面上の軌道の形状の一例を示す図である。 異常状態におけるセンサ計測値の経時変化の一例を示すグラフである。 異常状態における位相平面上の軌道の形状の一例を示す図である。 非類似度算出部１０８によって分類された複素パワーケプストラム距離Ｄｃの一例を示す図（デンドログラム）である。 推定部１０９によって生成された状態遷移マトリクスの一例を示す図である。 異常検出装置１の処理の一連の流れを示すフローチャートである。 ステップＳ１０２に係る特定センサの推定処理の一連の流れを示すフローチャートである。 ステップＳ３０２に係る特定センサの推定処理の一連の流れを示すフローチャートである。
＜実施形態＞
　以下、図を参照して本発明の実施形態について説明する。
［全体構成］
　図１は、実施形態の異常検出装置１の構成の一例を示す図である。異常検出装置１は、異常検出対象の装置やシステムにおいて生じた異常を検出する。異常検出装置１は、例えば、制御部１００と、記憶部２００とを備える。制御部１００は、例えば、ＣＰＵ（Central Processing Unit）等のハードウェアプロセッサが記憶部２００に記憶されるプログラム（ソフトウェア）を実行することにより、取得部１０２と、解析部１０４と、導出部１０６と、非類似度算出部１０８と、推定部１０９と、出力部１１０とを機能部として実現する。また、これらの構成要素のうち一部又は全部は、ＬＳＩ（Large Scale Integration）やＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field-Programmable Gate Array）、ＧＰＵ（Graphics Processing Unit）等のハードウェア（回路部；circuitryを含む）によって実現されてもよいし、ソフトウェアとハードウェアの協働によって実現されてもよい。
　記憶部２００は、例えば、ＨＤＤ（Hard Disk Drive）やフラッシュメモリなどの記憶装置（非一過性の記憶媒体を備える記憶装置）により実現されてもよく、ＤＶＤやＣＤ－ＲＯＭなどの着脱可能な記憶媒体（非一過性の記憶媒体）により実現されてもよく、ドライブ装置に装着される記憶媒体であってもよい。また、記憶部２００の一部又は全部は、ＮＡＳ（Network Attached Storage）や外部のストレージサーバ等、異常検出装置１のアクセス可能な外部装置であってもよい。記憶部２００には、例えば、相関関係情報２０２と、軌道形状情報２０４、参考情報２０６とが記憶される。各種情報の詳細については、後述する。
　取得部１０２は、異常検出対象の装置やシステムに配置された複数のセンサ（不図示）の各々が検出対象の状態を所定時間毎に検出した複数の検出結果を示す情報（以下、検出結果情報３００）を取得する。検出対象は、例えば、異常検出対象の装置やシステムにおいて、変化が生じる物理量であり、検出対象の状態は、例えば、物理量の値である。この物理量は、例えば、機械的性質、電磁気的性質、熱的性質、音響的性質、化学的性質等のうち、少なくともいずれかを示すものであってもよく、或いはそれらで示される空間情報や時間情報であってもよい。この物理量は、例えば、圧力、流量、弁開度、大気圧、外気温、水温、血圧、脈拍等を含む。異常検出対象の装置やシステムは、例えば、ロケット（例えば、再使用ロケットのエンジンシステム）、交通システム、発電装置、生命維持装置等である。異常検出対象のシステムがロケットのエンジンシステムである場合、複数のセンサは、ロケットに搭載され、エンジンシステムの状態を知るための検出対象（物理量）を検出する。
　なお、異常検出対象の装置やシステムは、これらロケットのエンジンシステムや発電装置等に限定されず、検出対象である物理量が取得可能な装置やシステムであればよい。また、取得部１０２は、検出結果情報３００の他、必要に応じて、異常検出、又は異常の種類や原因の推定に必要な情報などを取得する構成であってもよい。
　取得部１０２は、例えば、異常検出装置１に接続される複数のセンサから検出結果情報３００を取得してもよく、ネットワークを介して情報の送受信が可能な接続されたセンサから検出結果情報３００を取得してもよく、動作を完了した異常検出対象の動作履歴（例えば、ログ情報）に含まれるセンサの検出履歴を、検出結果情報３００として取得（抽出）してもよい。以降の説明では、取得部１０２が、データロガーに記憶されるログ情報から検出結果情報３００を取得する場合について説明する。センサは、「検出部」の一例である。
　図２は、検出結果情報３００の内容の一例を示す図である。検出結果情報３００には、例えば、物理量を計測しているセンサを識別可能な情報（図示するセンサＩＤ）と、センサ名と、当該センサが所定時間毎に検出対象である物理量を検出した検出結果とが互いに対応付けられた情報が、センサ毎に区別して格納される。図２において、センサＩＤが「０００１」であり、センサ名が「Ａ」のセンサは、検出対象が異常検出対象の時系列情報であり、検出結果には、センサが測定（検出）した計測値（例えば、温度、圧力、流量等）が示される。
　図１に戻り、解析部１０４は、取得部１０２によって取得された複数の単変量の検出結果を、多変量解析に用いる行列分解法（例えば、主成分分析）によって解析することにより、異常の特徴に応じた複数のモード情報を取得する。
　導出部１０６は、例えば、解析部１０４によって取得されたモード情報に基づいて、位相平面上、又は位相空間上における軌道を導出することにより、モード情報を可視化する。非類似度算出部１０８は、導出部１０６が導出した可視化情報に基づいて、推定部１０９が特定センサを推定する際に用いる指標（以下、非類似度）を算出する導出部１０６と、非類似度算出部１０８の詳細については、後述する。
　推定部１０９は、例えば、取得部１０２から取得された検出結果情報３００、又は解析部１０４によって取得されたモード情報のうち、少なくとも一方に基づいて、導出部１０６によって導出された軌道、又は非類似度算出部によって算出された非類似度の少なくとも一方に基づいて、各検出結果を検出したセンサが、異常な値を示しているセンサ（以下、特定センサ）であるか否かを推定する。検出する異常には、例えば、通常時には発生しない周期性の変化、不定期ノイズの混入、計測値の振動等の多様な異常が含まれる。また、異常の種類には、例えば、異常検出対象の装置やシステムの異常（以下、システム異常）と、センサの異常（以下、センサ異常）とが含まれる。
　出力部１１０は、推定部１０９によって特定センサであると推定されたセンサに係る情報を出力する。出力部１１０は、例えば、異常検出装置１に接続される表示装置（不図示）に特定センサを示す画像を表示してもよく、ネットワークを介して接続される他の装置に特定センサを示す情報を出力してもよく、記憶部２００に特定センサを示す情報を記憶させてもよい。なお、出力部１１０は、取得部１０２によって取得された検出結果、解析部１０４によって取得されたモード情報、導出部１０６によって導出された軌道を示す情報、非類似度算出部１０８によって算出された非類似度等を記憶部２００に記憶させてもよい。また、記憶部２００は、異常検出装置１と別体に構成されていてもよい。
［解析部１０４の処理］
　以下、解析部１０４の処理の詳細について説明する。まず、解析部１０４は、例えば、取得部１０２によって取得された検出結果に、公知の技術である滑走窓法などを利用して単変量の時系列情報を部分時系列とよばれる近傍の計測値を要素とするベクトルの集まりに変換する。これにより、解析部１０４は、検出結果を近傍の計測値を変数としてもつ多変量データとして再構成する。滑走窓法により再構成された検出結果に対して、複数のモード情報取得が可能な行列分解法を用いて多変量解析を行う場合、解析部１０４は、主成分分析結果に時間構造（ラグ構造）を取り込むことができるだけでなく、且つ窓幅の数だけ、検出結果から複数のモード情報を取得することができる。滑走窓法の詳細については、後述する。このようにして得られた複数のモード情報には、それぞれ異なる性質の異常（周期性の変化、不定期ノイズの混入、計測値の振動等）が示されており、推定部１０９は、解析部１０４によって取得された各モード情報に基づいて、多様な異常を判定する。
　以降の説明において、解析部１０４が行う処理であり、単一のセンサの時系列情報を、ラグ構造を利用してベクトル化することによって多変量データに再構成し、複数のモード情報取得が可能な多変量解析法（例えば、主成分分析）を利用して行うモード情報の取得処理を、「第１種モード情報取得処理」と記載する。これに対し、解析部１０４が行う処理であり、２つ以上のセンサの時系列情報を時間毎に集めてベクトル化することによって多変量データに再構成し、これを複数の独立した多変量データとして扱うことにより実行される一般的な時間領域での多変量解析法（例えば、主成分分析）を利用して行うモード情報の取得処理を「第２種モード情報取得処理」と記載する。
　解析部１０４は、「第２種モード情報取得処理」によって、２つ以上のセンサをまとめて多変量解析を行う場合、「第１種モード情報取得処理」とは異なる視点によるモード情報を取得することができる。また、解析部１０４は、主成分分析を用いた「第２種モード情報取得処理」を行う場合、主成分分析を用いた「第１種モード情報取得処理」において滑走窓法を用いた場合において窓幅の数だけモード情報が抽出される場合とは異なり、多変量データの再構成に用いたセンサの数だけ複数のモード情報を取得することができる。また、解析部１０４は、「第２種モード情報取得処理」において、一度に複数のセンサの検出結果に基づいてモード情報を取得するため、センサ毎にモード情報を取得する「第１種モード情報取得処理」に比べて、異常検出対象が備えるセンサ数が多い場合には、計算を省力化することができる。
　解析部１０４は、例えば、検出結果情報３００に示される単一センサの検出結果のモード情報取得処理前の検出結果（以下、０次モード情報と呼ぶ）や、当該センサの検出結果に標準化を行ったデータ（以下、標準化０次モード情報）を、「第１種モード情報取得処理」の対象データとする。また、解析部１０４は、例えば、「０次モード情報」、「標準化０次モード情報」、「第１種モード情報取得の結果である１次から高次までのモード情報」、「標準化された第１種モード情報取得の結果である１次から高次までのモード情報」等を、「第２種モード情報取得処理」の対象データとする。
である。
　また、解析部１０４は、前記のように「第１種モード情報取得処理」を行ったデータを用いて「第２種モード情報取得処理」を行うことにより、「第２種モード情報取得処理」で失われるラグ構造の取り込みを可能にすることができる。これは以下のような理由による。滑走窓法を利用して、データを再構成し、多変量解析を行った「第１種モード情報取得処理」データには、ラグ構造が取り込まれる。一方、「第２種モード情報取得処理」の場合、解析部１０４は、２つ以上のセンサの時系列情報を時間毎に集めてベクトル化することによって多変量データに再構成するため、データ間の時間相関は失われる。解析部１０４は、ラグ構造が導入された「第１種モード情報取得処理」の結果をセンサの時系列情報として、「第２種モード情報取得処理」に用いることでモード情報取得の際に失われる時間相関を回復することが可能になる。
　解析部１０４は、取得部１０２によって取得された検出結果情報３００に基づいて、センサ毎に実行される「第１種モード情報取得処理」、或いは２つ以上のセンサの組合せに対して実行される「第２種モード情報取得処理」を行い、複数のモード情報により、閾値の設定されていない異常や閾値の設定が困難な異常（周期性の変化や不定期ノイズの混入など）など多様な異常の特徴を抽出する。以降の説明において、解析部１０４が、多様な異常の特徴を捉えたモード情報を取得する手法を、（１）「多様な異常の特徴を抽出するためのモード情報取得法」とも記載する。
［導出部１０６の処理］
　以下、導出部１０６の処理の詳細について説明する。導出部１０６は、解析部１０４によって取得された複数のモード情報に基づいて、公知の手法（例えば、位相面法）によりモード情報を可視化する。以降は、導出部１０６は、解析部１０４によって取得された複数のモード情報に基づいて、位相平面、又は位相空間上に軌道を導出し、モード情報を可視化する場合について説明する。これにより、導出部１０６は、導出した軌道により検出結果の変化の特徴を拡大し、軌道の観察者に対してわずかな変化でも見逃しなく検知させることができる。更には、導出部１０６は、導出した軌道により、軌道の観察者に対して多くの物理的な情報を保持したまま視認性良く、異常の検知、異常の種類及び原因の推定などを行わせることができる。また、以降の説明において、導出部１０６が、解析部１０４によって取得された複数のモード情報に基づいて、位相平面、又は位相空間上に軌道を導出する手法を、（２）「軌道形状を利用したモード情報の可視化法」とも記載する。
［非類似度算出部１０８の処理］
　また、非類似度算出部１０８は、例えば、導出部１０６によって位相平面に導出された軌道形状を、複素自己回帰係数に基づいて計算したユークリッド距離、対数尤度比距離、複素パワーケプストラム距離、複素パワーメルケプストラム距離、動的時間伸縮法、或いはニューラルネットワークなどを用いて、通常時との軌道形状の差を数値化する非類似度を算出し、この数値とこの非類似度に基づいてクラスタリングを行うことにより、異常事例と通常事例を分類する。推定部１０９は、このクラスタリング結果を用いて、閾値によらずに異常の有無を推定するとともに、異常の種類や原因の推定などを行う。以降の説明において、非類似度算出部１０８が軌道形状を数値化する手法を、（３－１）「モード情報により生成される軌道形状の定量法」と記載する。同じく非類似度算出部１０８が数値化した軌道形状を複数の通常事例における数値化した軌道形状との非類似度を計算し、これを基にクラスタリングを実行して異常事例と通常事例の分類を行う手法を、（３－２）「軌道形状の数値化情報に基づく非類似度の算出とこれに基づくクラスタリング法」とも記載する。
　推定部１０９は、例えば、事前に生成されたセンサ間の相関関係情報２０２やセンサ部品の交換履歴など異常を判定する際の参考になる情報が格納された参考情報２０６などを参考にして、例えば、解析結果を統合した状態遷移マトリクスを利用（生成）し、異常の検知やセンサ間の相関関係の破綻の推定、異常原因の判定などを行う。以降の説明において、推定部１０９が解析結果を統合した状態マトリクスを利用して、異常の有無の推定やセンサ間の相関関係の破綻の推定に基づく異常の種類や原因の推定などを行う手法を、（４）「非類似度算出部の計算結果を状態遷移マトリクスにより統合して詳細な推定を行う方法」とも記載する。
　図３は、相関関係情報２０２の内容の一例を示す図である。相関関係情報２０２は、センサ間の相関関係を示した情報である。相関関係情報２０２は、例えば、通常時のセンサの検出結果（以下、通常データ）を基づいて、予め生成された情報である。図３には、互いに相関関係にあるセンサが、センサ毎に示されている。センサＩＤが「０００１」であり、センサ名が「センサＡ」であるセンサと検出結果が相関関係にあるセンサは、センサＩＤが「０００２」～「０００６」であり、センサ名がそれぞれ「センサＢ」「センサＣ」「センサＤ」「センサＦ」のセンサである。あるセンサと相関関係にあるセンサは、あるセンサ（例えば、「電圧センサ」）の検出対象（例えば、「入口電圧」）の変化に応じて変化する同種の検出対象（例えば、「出口電圧」）を検出するセンサまたは異種の検出対象（例えば、「入口電流」など）を検出するセンサ（例えば、「電流センサ」）である。
　なお、一方の検出結果が、他方の検出結果の変化に応じて変化するものであっても、設置場所や設置環境に起因して検出結果に相関関係がないと推定されるセンサは、相関関係情報２０２に含まれていなくてもよい。また、相関関係情報２０２は、予め生成される情報であってもよく、後述する方法によって、異常検出装置１の機能を用いて、生成されるものであってもよい。この場合、異常検出装置１は、相関関係情報２０２の生成と同時並行によって異常の検出を行ってもよい。
　推定部１０９は、相関関係情報２０２を参照し、センサが異常値であるか推定する。推定部１０９は、例えば、解析部１０４により取得されたモード情報、導出部１０６により導出された軌道、非類似度算出部１０８により算出された非類似度等に基づいて、「センサＡ」と他のセンサの組合せが異常値を示していると推定し、且つ「センサＡ」と相関関係にあるセンサが他の組合せでは異常値を示してしていないと推定する場合は、「センサＡ」と相関関係にあるべきセンサとの相関関係が破綻しており、「センサＡ」は「センサ故障」であると推定する。換言すると、推定部１０９は、「センサＡ」をセンサ故障により異常値を示しているセンサであると推定する。一方、推定部１０９は、解析部１０４により取得されたモード情報、導出部１０６により導出された軌道、非類似度算出部１０８により算出された非類似度等に基づいて、「センサＡ」と相関関係にある他のセンサの相関関係が破綻していない場合は、「システム異常」と推定する。このとき、例えば、解析部１０４、導出部１０６、非類似度算出部１０８、推定部１０９等における処理を相関関係のあるセンサだけに絞って行うことで、異常検出に要する時間を短縮することができる。推定部１０９の処理の詳細は後述する。
　そして、推定部１０９で、「システム異常」と推定された場合は、相関関係情報２０２の他、参考情報２０６を利用して、異常値を示しているセンサの特定と異常の原因を推定する。推定部１０９の上記手順の詳細は後述する。
　図４は、参考情報２０６の内容の一例を示す図である。参考情報２０６は、センサを識別可能なセンサＩＤと、センサ名と、センサの計測点近傍の装置およびシステムを構成する部品名（例えば、センサが管路に設置された流量計ならば、管路を構成する管材など）、その部品の交換日、部品を製造（又は、販売）する業者（メーカー）等の、センサ周辺部品の整備履歴の他、センサ自身を構成する部品の部品名、センサの交換日、センサ製造者など異常検出の原因の推定に用いるために必要な情報などがセンサ毎に互いに対応付けられた情報である。なお、参考情報２０６の内容は、一例であって、これに限られず、参考情報２０６に記載しておくことが有用であると考えられる情報を記しておくことができる。
　以下、上述の異常推定プロセスの一例を（１）「多様な異常の特徴を抽出するためのモード情報取得法」、（２）「軌道形状を利用したモード情報取得情報の可視化法」、（３－１）「モード情報取得情報により生成される軌道形状の定量法」（３－２）「軌道形状の数値化情報に基づく非類似度の算出とこれに基づくクラスタリング法」、（４）「非類似度算出部の計算結果を状態遷移マトリクスにより統合して詳細な推定を行う方法」に分けて、その詳細を説明する。
［（１）多様な異常の特徴を抽出するためのモード情報取得法について］
　以下、解析部１０４が実行するモード情報取得処理の内容について説明する。上述したように、解析部１０４は、取得部１０２によって取得された単一センサの時系列情報に対して、例えば、滑走窓法を適用して時系列情報をベクトルの集まりに変換することにより、近傍の時間情報を変数としてもつ多変量データとして再構成することができる。これにより、解析部１０４は、多変量解析に適用する行列分解法を利用して窓幅の数だけ、複数モードの異常抽出が可能になる。また、解析部１０４は、同時にモード情報にラグ構造を取り込むこともできる。同時にモード情報にラグ構造を取り込む処理は、上述した「第１種モード情報取得処理」である。
　図５は、窓幅５の滑走窓法を用いて、単一センサの時系列情報をベクトルの集まりに変換し、多変量データとして再構成し、主成分分析によりモード情報取得処理を行う「第１種モード情報取得処理」の一例を示す図である。解析部１０４によって再構成されたデータは、多変量解析法を用いて多変量解析を実行することにより、異なる特徴を抽出したモード情報を得ることが可能になる。さらに、解析部１０４によって再構成されたデータは、各モードにはそれぞれ異なる性質をもった異常の特徴が抽出されるだけでなく、上述の「第２種モード情報取得処理」において複数のセンサの検出結果を集めて多変量データに再構成する際に、「第１種モード情報取得処理」が行われた情報を使うことにより、時間相関を無視して再構成された多変量データを用いて多変量分析が行われることに伴い失われる時間相関を回復することができる。
　以上述べたように、解析部１０４は、「第１種モード情報取得処理」において、例えば、ラグ構造を導入することにより、単一センサの検出結果から得られる単変量データを多変量データに変換して、主成分分析法、カーネル主成分分析法、ファジィ主成分分析法、スパース主成分分析法、確率的主成分分析法、ロバスト主成分分析法のような多変量解析手法を適用する。これにより、解析部１０４は、ことで複数のモード情報を取得する。解析部１０４が、「第１種モード情報取得処理」において上述の例のように滑走窓法を用いた場合には、窓幅のサイズｍ（ｍは、自然数。）だけ、複数のモード情報を取得することができる。複数のモード情報のそれぞれは、異なる異常の特徴が抽出されるため、異常検出装置１は、複数のモード情報を用いて異常検出を行うことで、見逃しのない異常検出が可能になる。
　図６は、「第２種モード情報取得処理」の一例を示す図である。について示したものである。「第２種モード情報取得処理」とは、上述のように解析部１０４が、同時刻の２つ以上のセンサの検出結果を要素にもつベクトルの集まりに変換して、独立なデータセットに再構成することで、多変量解析を利用し、複数モードの特徴抽出を行うものである。図６（Ａ）～（Ｃ）は、「センサＡ」と「センサＢ」の先に０次モード情報と名付けたモード情報取得処理前の検出結果で作る位相面の一例を示す３つのグラフである。図６（Ａ）～（Ｃ）に示される通り、０次モード情報では、異常が明瞭ではない。図６（Ｄ）～（Ｅ）は、主成分分析を利用した「第２種モード情報取得処理」により抽出した第１主軸上と第２主軸上のスコア（以後、それぞれ１次モード情報、２次モード情報ともよぶ）を示すグラフである。図６（Ｄ）～（Ｅ）に示される通り、第２主軸上（２次モード）では、異常の特徴が抽出されていることがわかる。「第１種モード情報取得処理」においても同様に、取得されたモードに応じて異なる異常の特徴を抽出することができる。「第２種モード情報取得処理」では、例えば、一つのデータセットにｓ個のセンサが含まれているとすると、センサの数ｓ（ｓは、自然数。）だけ、異なる異常の特徴をとらえた複数モードの特徴抽出を行うことができる。
　なお、解析部１０４によって取得された複数のモード情報の中で、どのような異常を検出するのにどのモード情報が適しているかは、対象とする装置やシステムのセンサ情報を用いて、事前に調べておくことが可能である。例えば、閾値が設定された異常は、０次モードが異常を検出するのに適したモードである。また、絶対値が大きく変化するような異常は、１次モードが異常の検出に適したモードである。局所的な変化の大きい異常は、２次モードが異常の検出に適したモードである。また、解析部１０４が行う「第２種モード情報取得処理」では、センサ間の相関関係が反映されるため、センサ毎にモード情報取得を行う「第１種モード情報取得処理」では抽出できない異常を検出することが可能になる。異常検出装置１は、これらのモード情報や抽出法がもつ特徴に基づいて、事前に異常の性質に応じて使用するモードや抽出方法を選定しておけば、網羅的な解析をすることなく、目的に応じて異常検出や異常の予兆の推定などの処理を迅速に行うことができる。
　また、参考情報２０６には、各異常の検出に適したモードを示す情報が含まれていてもよく、解析部１０４は、当該情報を参照し、各異常を検出する際に適したモードを解析する構成であってもよい。異常検出装置１は、このようにして、複数モードで特徴抽出を行うことで、多様な異常（閾値が設定されていない異常ばかりでなく不定期ノイズの混入のような閾値の設定が困難な異常）を、見逃しなく検出することが可能になる。
　上述したように、解析部１０４が「第２種モード情報取得処理」に用いる多変量データは、取得部１０２によって取得された検出結果（０次モード情報）を複数集めて多変量データに再構成したものでもよいが、「第１種モード情報取得処理」を行った１次以上のモード情報を複数集めて多変量データに再構成したものでもよい。解析部１０４によって「第１種モード情報取得処理」が行われたデータには、ラグ構造が導入されているため、解析部１０４は、「第１種モード情報取得処理」を行った情報を用いて「第２種モード情報取得処理」を行うことで、「第２種モード情報取得処理」において、データ間の独立を仮定した多変量解析を行う際に失われる時間相関を回復することできるといった利点がある。
　なお、解析部１０４は、「第１種モード情報取得処理」が行われたモード情報を用いて多変量データに再構成する際には、組み合わせるすべてのセンサのモードの種類を統一してもよいし（例えば、すべて１次モード）、センサ毎に０次モード情報も含めて、モードを変えて多変量データに再構成するなど、任意の方法で「第１種モード情報取得処理」情報を利用してもよい。
［（２）軌道形状を利用したモード情報の可視化法について］
　以下、軌道形状を利用したモード情報の可視化法によって、導出部１０６がモード情報取得情報を可視化するための処理の内容について説明する。位相面法は、時間を横軸、計測値を縦軸にとった通常のプロットでは見分けることが困難なわずかな変化でも位相面上に拡大する効果があるため、感度よく状態変化を可視化する手法としてよく知られた手法である。導出部１０６は、位相面法を利用した軌道だけでなく、図６Ｄ、Ｅに示したような横軸時間、縦軸主成分スコアのような時系列プロットによる軌道生成を行うこともできる。時系列プロットによる軌道でも、位相面法を利用した軌道に用いる検知手順を応用して、非類似度算出部１０８、推定部１０９を経て、異常を検知することができる。
　以降の説明では、可視化対象として「第１種モード情報取得処理」の結果を利用して位相面法により、２次元平面上の軌道情報として異常を可視化する場合について説明するが、導出部１０６が生成する位相面は２次元に限定されるものではない。高次元では視認が難しくなるが、３次元以上の位相面空間上に軌道を生成し、２次元平面上の軌道に対して実行される処理を応用して異常検出を行うことは可能である。導出部１０６が「第１種モード情報取得処理」の結果を、位相面法を利用して可視化する場合には、「第１種単変量位相面軌道生成処理」と「第１種多変量位相面軌道生成処理」の２つの方法がある。これらの処理の詳細については、後述する。可視化された結果を最終的に人間が確認する場合には、２次元平面で扱ったほうが、結果を評価しやすい場合が多いが、導出部１０６は、可視化が困難な高次元の位相面軌道を生成し、以降の処理を行ってもよい。
　また、導出部１０６において、位相面軌道を生成する対象は、「第１種モード情報取得処理」の結果に限られない。例えばモード情報取得を行う前の０次モード情報も含めて位相面軌道を生成することができる。導出部１０６は、窓幅ｍの「第１種モード情報取得処理」を行った結果を利用した場合は、ｍ個のモード情報に０次モード情報を含めて軌道生成に利用することができ、ｓ個のセンサの検出結果を組み合わせて「第２種モード情報取得処理」を行った結果を利用した場合には、ｓ個のモードに０次モード情報を含めて軌道生成に利用することができる。
　以下、導出部１０６が、「第１種モード情報取得処理」の結果を利用して２次元位相面を生成する方法の詳細について述べる。導出部１０６は、例えば、単一の「センサＡ」について、横軸に「センサＡ」の１次モード情報（第１主成分軸上のスコア）、縦軸に「センサＡ」の２次モード情報（第２主成分軸上のスコア）のように、両軸に同一センサの検出結果に係る情報をプロットする方法によって、２次元位相平面を生成する。以降の説明において、この生成方法を、「第１種単変量位相面軌道生成処理」と記載する。
　なお、上述したプロット対象の組合せは、一例であってこれに限られない。例えば、横軸に「センサＡ」の２次モード情報、縦軸に「センサＡ」の３次モード情報（第３主成分軸上のスコア）をプロットすることもできる。組み合わせは、異常がよりよく検知できるモードを組み合わせることができる。例えば、１次モードと３次モードに異常がよく表れている場合に、横軸１次モード、縦軸３次モードのように組み合わせてプロットすることで、異常の特徴を効果的に可視化することが可能になる。異常の推定に適切なモードは、異常の性質（周期性の変化、不定期ノイズの混入、計測値の振動など）によって異なるので、導出部１０６は、複数のモードで位相面軌道を生成し、推定部１０９は、生成した位相面軌道を用いた後述のプロセスを経て、得られた指標（例えば、非類似度）を用い、異常が認められるか否かを推定する。導出部１０６は、異常の性質がわからない場合には、網羅的にモードを組み合わせて、位相面を生成することになるが、検知したい異常の性質が決まっている場合には、その異常の性質に応じた適切なモードで位相面を生成する。例えば、心電図において、不整脈を検出したい場合、導出部１０６は、不整脈の種類に応じて、異常の検出に最適なモードで位相面を生成する。したがって、導出部１０６は、検出したい異常の性質が予め決定されている場合には、生成する軌道に使うモード情報を特定のものに限定するなどして、異常検出に要する時間を短縮することができる。
　なお、後述において詳細を説明する「第１種多変量位相面軌道生成処理」において、導出部１０６は、複数センサの変化により、異常検出を行うが、「第１種単変量位相面軌道生成処理」において、導出部１０６は、単一センサの変化に基づき、上述した例に限定されず、すべての軸に０次モード情報を含む同じモード情報を用いてもよく、適宜モード情報を組み合わせて位相面軌道を生成する。導出部１０６は、軌道生成を異常の特徴に適したモードに限定することで、異常の検知に要する時間が短縮される。
　また、導出部１０６は、「第１種モード情報取得処理」の結果を利用する場合には、両軸に同じセンサから得られた情報をプロットする代わりに、異なるセンサから得られた情報をプロットして、２次元位相平面を生成することもできる。以降の説明において、この生成方法を、「第１種多変量位相面軌道生成処理」と記載する。導出部１０６は、「第１種多変量位相面軌道生成処理」においても、「第１種単変量位相面軌道生成処理」と同様にプロットするモードの組み合わせを、異常の性質に応じて選ぶことで効果的な異常検出を行うことができる。「第１種多変量位相面軌道生成処理」では、複数のセンサ間の変化に基づいた異常検出により、閾値が設定されていない異常の検知や、異常の原因の推定が可能になる。次に、導出部１０６が、「第１種多変量位相面軌道生成処理」を行う場合について詳しく説明する。
　導出部１０６において、「第１種多変量位相面軌道生成処理」により位相面を生成する方法は例えば以下のようになる。以降の説明では、説明の便宜上、導出部１０６が「センサＡ」の１次モード情報と「センサＢ」の１次モード情報の組合せで位相面軌道を生成する例について説明する。例えば、導出部１０６は、解析部１０４で「第１種モード情報取得処理」により取得された「センサＡ」の１次モード情報を横軸、「センサＢ」の１次モード情報を縦軸にプロットした位相面を生成する。
　導出部１０６は、モード情報取得処理前の０次モード情報を用いて位相面を生成することも可能である。また、上述したように、モード情報の組合せは、これに限られない。異常の性質に応じて異常が抽出されるモードが異なるため、位相面軌道を生成する際に、０次モード情報も含めた最適なモード情報の組合せを選ぶことにより、通常時との違いがより明瞭に表れ、異常検出や異常の種類や原因の推定を効率よく行うことができる。
　以下では、導出部１０６が、「第１種モード情報取得処理」で抽出した「センサＡ」の１次モード情報と「センサＢ」の１次モード情報を用いて、「第１種多変量位相面軌道生成処理」を行う例を示す。導出部１０６は、網羅的にすべてのセンサの組合せについて軌道の導出対象とすることも、異常検出の目的に応じて、事前に相関関係情報２０２を用いて選定したセンサグループのみを軌道の導出対象とすることも、対象とするセンサの組合せを限定し、軌道の導出対象とすることも可能である。導出部１０６は、正の相関関係にあるセンサ群を対象に軌道を生成すれば、位相面軌道上の異常が拡大して表現することができ、異常の検知や異常の種類の推定が容易にさせることができる。導出部１０６は、負の相関関係にあるセンサ群を対象に軌道を生成すれば、一方の正負の符号を逆転させるなどプロットを工夫することにより、異常を拡大して表現することができ、異常の検知や異常の種類の推定が容易にさせることができる。また、以降の説明では、すべての組合せのうち、「センサＡ」の１次モード情報と「センサＢ」の１次モード情報の組合せを、軌道の導出対象とする場合を一例に説明する。また、以降の説明において、「第１種多変量位相面軌道生成処理」の対象として選定した「センサＡ」の１次モード情報を、モード１結果ｘと記載し、「センサＢ」の１次モード情報を、モード１結果ｙと記載する。
　導出部１０６は、モード１結果ｘを第１軸（例えば、横軸）とし、モード１結果ｙを第２軸（例えば、縦軸）とした位相平面上に、モード１結果ｘに含まれる各要素｛ｘ１、ｘ２、…、ｘｎ｝（ｎは自然数）と、モード１結果ｙに含まれる各要素｛ｙ１、ｙ２、…、ｙｎ｝との中で検出タイミングが合致する要素によって示される座標（図示する、座標（ｘ１，ｙ１）、…、座標（ｘｎ，ｙｎ））をプロットし、時系列順に、直線によって結んで軌道を生成する。
　図７は、通常状態におけるセンサの経時変化の一例を示すグラフである。図７に示される波形Ｗ２は、「センサＡ」の検出結果の経時変化を示す波形であり、波形Ｗ１は、「センサＢ」の検出結果の経時変化を示す波形である。
　図８は、通常状態における位相平面上の軌道の形状の一例を示す図である。導出部１０６は、波形Ｗ２の「センサＡ」の検出結果に基づいて導出されたモード１結果ｘを横軸とし、波形Ｗ１の「センサＢ」の検出結果に基づいて導出されたモード１結果ｙを縦軸とした位相平面において、モード１結果ｘの要素とモード１結果ｙの要素によって示される座標（図示する、座標（ｘ１，ｙ１）、…、座標（ｘｍ，ｙｍ））を時系列順に直線によって結ぶ軌道（図示する軌道Ｏｒｂ１）を生成する。
　なお、上述では、導出部１０６が、軌道の導出に際して１次モード情報を利用した例について説明しているが、位相面軌道の縦軸と横軸は、この組合せに限定されない。先に述べたように、解析部１０４が用いるモード情報は、０次モードから使用可能な高次モード情報まで自由に組み合わせたものを縦軸と横軸に選ぶことができる。
　推定部１０９は、導出部１０６によって導出された軌道の形状と、通常状態における軌道の形状を示す情報である軌道形状情報２０４との非類似度に基づいて、異常値が検出されているか否かを推定する。例えば、通常状態の軌道形状情報が格納された軌道形状情報２０４には、検出結果情報３００と同様に、ある時刻からある時刻までに取得された全センサ情報を１事例として、取得事例ごとに事例番号が付され、取得時刻情報などとともに、通常時の「第１種多変量位相面軌道生成処理」の結果が記録されている。その中には、例えば、通常時の「センサＡ」と「センサＢ」の検出結果を基に０次モードを含む複数のモードで導出した軌道の形状を示す情報が含まれる。具体的には、軌道形状情報２０４には、図８に示すような「センサＡ」と「センサＢ」の組合せに係る軌道の形状を示す情報が含まれる。図８の位相面軌道形状を事例番号１（以降事例Ｎｏ．１のように記述する）の「センサＡ」と「センサＢ」の１次モード情報を用いて導出された通常状態の位相面軌道形状とすると、軌道形状情報２０４には事例Ｎｏ．１～事例Ｎｏ．ｋ（ｋは、自然数。）までの複数の事例を対象にこのような位相面軌道情報が格納されている。推定部１０９は、非類似度算出部１０８による、複数の通常状態の軌道形状と検査対象の軌道形状を非類似度に基づいてクラスタリングを行って、検査対象の軌道形状が通常状態に分類されるかどうかを調べた結果に基づいて、閾値が設定されていない異常の検知が可能になる。複数の通常事例のどれとも似ていなければ、通常事例とは違う分類に属すると判定することができるからである。ここでは、連番事例を例に説明しているが、通常事例としてとりあげる事例は、連番である必要はなく、過去の通常事例から適切なものを自由に選ぶことができる。以下では、連番の通常事例を選んだ場合について説明する。
　非類似度算出部１０８による非類似度の算出法とクラスタリング手法の詳細については後述するが、例えば、推定部１０９は、軌道形状情報２０４に格納されている事例Ｎｏ．１から事例Ｎｏ．ｋまでの通常事例における軌道の形状の中から、注目するセンサ対に対応するセンサの組合せ（この一例では、「センサＡ」と「センサＢ」）に対応する軌道の形状を取り出して、事例Ｎｏ．１から事例Ｎｏ．ｋまでの通常事例における注目するセンサ対に対応するセンサ対がつくる位相面軌道形状との非類似度を比較することにより、センサが異常値を計測しているか否かを推定する。このように、異常検出装置１は、「モード情報取得処理」による異常抽出と「位相面法」を併用することにより、図７に示すような一般的な時系列グラフを用いて異常を解析するよりも、図６Ｅに示したように異常の特徴を抽出し、さらに図１０に示すようにその変化の特徴を位相面軌道上に拡大して検出できるという利点がある。さらに、導出部１０６は、２つのセンサが正の相関関係にある場合は、一般的な時系列グラフよりも、その変化が拡大して可視化される位相面軌道を導出することとなり、わずかな異常の兆候も見逃しなく検知することができるような情報（つまり、軌道形状）を提供することができる。そして、推定部１０９は、注目しているセンサ対の位相面軌道形状が通常時の「センサＡ」と「センサＢ」がつくる事例Ｎｏ．１から事例Ｎｏ．ｋの対応するどの位相面軌道形状と比較しても非類似度が低い場合（すなわちよく似ている場合）、検査対象の「センサＡ」と「センサＢ」は異常値を示していないと推定し、非類似度が高い（異なる形状である）場合、検査対象の「センサＡ」と「センサＢ」は異常値を示していると推定する。推定部１０９は、複数の通常事例と比較することで、事前に閾値を設けることなく、通常事例と異常事例の非類似度による分類が可能になる。非類似度による分類法の詳細は後述する。
　非類似度の推定は、観測者によって定性的に判定することもできるが、ここでは、軌道形状を数値化し、判定を自動化して結果を報告する手法について説明する。なお、軌道の形状の非類似度が低いとは、形状同士が完全一致する場合だけでなく、軌道のずれが所定範囲内にある場合や、観測者がみて同一とみなせる場合も含まれる。非類似度の推定法の詳細は後述する。
　図９は、異常状態におけるセンサの経時変化の一例を示す図である。図９に示される波形Ｗ４は、「センサＡ」の検出結果の経時変化を示す波形であり、波形３は、「センサＢ」の検出結果の経時変化を示す波形である。この一例において、「センサＡ」は異常な動作をしており（つまり、特定センサであり）、波形Ｗ４によって示される通り、変化の周期が速くなり、振幅も小さくなっていることがわかる。
　図１０は、異常状態における位相平面上の軌道の形状の一例を示す図である。導出部１０６は、図９に示した波形Ｗ４の「センサＡ」の検出結果に基づいて導出されたモード１結果ｘを横軸とし、波形Ｗ３の「センサＢ」の検出結果に基づいて導出されたモード１結果ｙを縦軸とした位相平面において、モード１結果ｘの要素とモード１結果ｙの要素によって示される座標（図示する、座標（ｘ１，ｙ１）、…、座標（ｘｎ，ｙｎ））を直線によって時系列順に結ぶ軌道（図示する軌道Ｏｒｂ２）を生成する。
　図１０に示される軌道Ｏｒｂ２の形状は、図８の軌道Ｏｒｂ１の形状とは異なる形状である。軌道形状情報２０４に格納されている他の複数の通常事例の「センサＡ」と「センサＢ」の位相面軌道形状と比較して、どの通常事例とも似ていないとされた時には、推定部１０９は、軌道Ｏｒｂ２に係る「センサＡ」と「センサＢ」のどちらか或いは両方が異常値を示していると推定することができる。上述のように複数の通常事例と比較して総合的な判断をすることで、閾値によらない異常検出が可能になる。推定部１０９は、非類似度算出部１０８による位相面軌道形状を数値化し、複数の過去の通常事例の位相面軌道形状から得られる数値との非類似度に基づき、クラスタリングの結果から異常か否かを推定する。
　推定部１０９が、複数モードの位相面軌道を利用した多変量解析（この場合、「センサＡ」及び「センサＢ」の変化に基づく解析）結果を複数の通常時のデータと比較することにより、異常か否かを推定することで、閾値が未設定、あるいは閾値が設定できないなどの理由で、単一センサの時系列変化の解析では異常か否かの推定が困難な場合でも異常の検知が可能になる。さらに、モード情報取得により異常による変化を大きく取り出し、位相面によりその変化をさらに拡大して可視化することができるようになるため、推定部１０９は、微かな異常の兆候も見逃しなく検知することも可能になる。導出部１０６は、複数モードの位相面軌道を利用した解析により、異常検出のみならず、故障の予知も可能にする詳細な解析材料を提供することができる。
　なお、解析部１０４、及び導出部１０６でのモード情報の取り扱い方法と同様に、推定部１０９で、異常の推定に使うモード情報は、０次モード情報から利用可能な高次モード情報までを組み合わせて網羅的な探索もできるが、検知したい異常の性質に応じた最適モードを予め調べておき、その情報を利用して、特定のモードで得られたモード情報のみについてこれらの処理を実行することができる。また、推定部１０９は、異常の性質によっては、高次モードの抽出処理の前（あるいはそれが必要であれば、モード情報取得処理の後でも）平均と分散をそろえる標準化処理を行うことにより、検出感度を上げることができる。
　以上、導出部１０６が、「第１種モード情報取得処理」結果に基づいて、位相面軌道を生成する場合について説明したが、位相面軌道の生成に用いる情報は、これに限られない。導出部１０６は、異常の性質によって異常が抽出されるモードが異なるため、周期性の変化や不定期ノイズの混入などといった異常の性質に応じて、適切なモードを用いて位相面軌道を生成し、非類似度算出部１０８は、導出部１０６によって生成された軌道に基づいて、非類似度を算出する必要がある。
［（３－１）モード情報取得情報により生成される軌道形状の定量法、（３－２）軌道形状の数値化情報に基づく非類似度の算出とこれに基づくクラスタリング法について］
　以下、モード情報取得情報により生成される軌道形状の定量法と軌道形状の数値化情報に基づく非類似度の算出とこれに基づくクラスタリング法の例として、非類似度算出部１０８が算出する複素自己回帰係数によって位相面軌道形状を数値化し、通常時の位相面軌道情報を記憶した軌道形状情報２０４から算出される通常時の位相面軌道形状の数値と検出対象のデータを用いて生成した位相面軌道から算出した数値を用いて非類似度を算出する処理の内容について説明する。軌道形状情報２０４には、予め以下に述べる方法で軌道形状を数値化した結果を位相面軌道情報と合わせて記録しておき、この数値情報を読み出して非類似度の算出を行ってもよい。これによって位相面軌道形状の数値化にかかる処理時間を短縮することができる。一例として、非類似度算出部１０８は、導出部１０６によって導出された軌道の形状を軌道形状毎に複素自己回帰係数を利用するといった画像抽出手法で数値情報に変換する。以降の説明では、非類似度算出部１０８が、複素自己回帰係数から計算された複素パワーケプストラム距離を用いて、通常状態との非類似度を算出する方法を説明する。非類似度算出部１０８は、複素自己回帰係数により軌道形状を数値化した情報を用いて、複数の通常状態の軌道形状を同様に数値化した情報との間で、２つの複素自己回帰係数間の差分（距離）である複素パワーケプストラム距離を算出する。さらに、非類似度算出部１０８は、算出された複素パワーケプストラム距離に基づいて、クラスタリングを行って、通常事例と異常事例を分類する。複素パワーケプストラム距離は、「位相面軌道形状間の非類似度計算」の一例であり、「位相面軌道形状間の非類似度計算」の手法は、これに限られない。非類似度算出部１０８は、ユークリッド距離、対数尤度比距離、複素パワーケプストラム距離、複素パワーメルケプストラム距離などの手法を用いて非類似度の算出を行うことができる。
　非類似度算出部１０８が、位相面軌道形状の数値化にあたって、導出部１０６によって導出された位相面軌道から、軌道の形状を抽出する手法の詳細を、複素自己回帰係数をもとに計算する複素パワーケプストラム距離を例に以下に示す。例えば、非類似度算出部１０８は、導出部１０６によって導出された位相平面上に示される軌道を所定の大きさ（又は、所定の画素数）の画像に変換し、変換した画像をグレースケールにした場合の各画素の明るさ（画素値）に基づいて、軌道の形状（輪郭線（エッジ））を抽出する。そして、複素自己回帰係数を用いて軌道形状を数値化する場合には、非類似度算出部１０８は、抽出した軌道の輪郭線を追跡して得られる点列を式（１）とし、その複素表現を式（２）とする場合の、ｍ次の複素自己回帰係数を導出する。ｍ次の複素自己回帰係数は、式（３）に示す通り、輪郭点をｍ個前までの輪郭点の線形結合で近似するモデルとして定義される。
Figure JPOXMLDOC01-appb-M000001
Figure JPOXMLDOC01-appb-M000002
Figure JPOXMLDOC01-appb-M000003
　非類似度算出部１０８は、導出部１０６によって生成された位相面軌道毎に、上記の方法で軌道形状情報を抽出し、抽出した軌道の形状情報を複素自己回帰係数に変換する。さらに、非類似度算出部１０８は、変換した軌道の複素自己回帰係数ｚ１と、軌道形状情報２０４に格納された該当する複数の通常時の位相面軌道情報の中から上記と同様の手順で複素自己回帰係数を導出し、通常時の軌道の複素自己回帰係数ｚ２との複素パワーケプストラム距離Ｄｃ（ｚ１，ｚ２）を導出する。複素パワーケプストラム距離Ｄｃ（ｚ１，ｚ２）は、複素自己回帰係数ｚ１と、複素自己回帰係数ｚ２とのスペクトル包絡の対数の平均二乗距離として、式（４）によって示される。なお、軌道形状情報２０４に示される軌道が、予め複素自己回帰係数に変換されている場合には、非類似度算出部１０８は、当該変換された情報を用いて複素パワーケプストラム距離Ｄｃを算出する。
Figure JPOXMLDOC01-appb-M000004
　複素パワーケプストラム距離は、検出対象の複素自己回帰係数に複数の通常時の複素自己回帰係数を加えた中から、網羅的に２つずつを組合せて計算される。例えば、「センサＡ」と「センサＢ」がつくる通常時の位相面軌道情報が事例Ｎｏ．１から事例Ｎｏ．１８まで合計１８個ある場合、非類似度算出部１０８は、この通常時の位相面軌道形状から導出される複素自己回帰係数ｚ１～ｚ１８と取得部１０２で取得した検査対象の「センサＡ」と「センサＢ」がつくる位相面軌道に基づく複素自己回帰係数ｚ１９の計１９個の複素自己回帰係数から、２つずつ取り出して、複素パワーケプストラム距離Ｄｃをそれぞれ算出する。
　非類似度算出部１０８は、この複素パワーケプストラム距離Ｄｃを階層型クラスタリングによって異常事例と通常事例に分類する。図１１は、非類似度算出部１０８によって分類された複素パワーケプストラム距離Ｄｃの一例を示す図である。図１１は、非類似度算出部１０８によって算出された複素パワーケプストラム距離Ｄｃを基に検出対象を含む全１～１９事例のクラスタリングを行った結果を示すデンドログラムとよばれる図である。図１１において、複素パワーケプストラム距離Ｄｃは、樹形図の枝の長さによって示される。一般に、複素パワーケプストラム距離Ｄｃの値が大きいほど、２つの複素自己回帰係数間の違いが大きくなり（つまり、非類似度が大きくなり）、複素パワーケプストラム距離の値が小さいほど、２つの複素自己回帰係数間の違いが小さくなる（つまり、非類似度が小さくよく似た軌道である）。例えば、図１１に示す一例において、推定部１０９は、取得部１０２から得られた情報を基に、解析部１０４、導出部１０６、非類似度算出部１０８を経て算出した「事例Ｎｏ．１９」と名付けられた位相面軌道形状の複素自己回帰係数と事例Ｎｏ．１～事例Ｎｏ．１８までのどの通常時の事例の複素自己回帰係数と比べても複素パワーケプストラム距離Ｄｃが大きいことから、「事例Ｎｏ．１９（検出対象）」は通常時との非類似度が大きく通常時とは別のグループに分類されたと判断し、「事例Ｎｏ．１９」が対象とした「センサＡ」と「センサＢ」の組合せのうち、２つのセンサのどちらかまたは両方が、異常値を示している特定センサとして推定することができる。
　なお、上述では、非類似度算出部１０８は、軌道の形状を数値に変換した複素自己回帰係数に基づいて、ユークリッド距離、対数尤度比距離、複素パワーケプストラム距離、或いは複素パワーメルケプストラム距離を算出する構成であってもよい。ユークリッド距離、対数尤度比距離、複素パワーケプストラム距離、及び複素パワーメルケプストラム距離は、「複素自己回帰係数間の差分」の一例である。
　また、上述で、非類似度算出部１０８が、軌道の形状の数値化に利用するモデルは、複素自己回帰係数を用いた手法に限定されない。また非類似度を算出する手法も上記に限られない。複素自己回帰係数を用いた手法は、図形の外形のみを抽出する手法であるため、より詳細な判別が必要である場合には、例えば、動的時間伸縮法のような位相面軌道形状の非類似度をラグ構造も考慮して算出することができる方法の他、内部形状抽出することが可能なニューラルネットワークなどの方法を用いることもできる。この結果に基づいて、複素自己回帰係数を用いた場合と同様に、非類似度を算出することができる。
　上述のように、推定部１０９は、非類似度算出部１０８によって算出された軌道形状を定量した数値に基づき、異常検出対象データから得られる数値と通常データから得られる数値との差分（この一例では、複素パワーケプストラム距離Ｄｃ）がある通常データと他の通常データ間で算出されるどの差分と比しても大きく、図１１に示すように、異常検出対象データと通常データ群が別のグループとしてクラスタリングされる場合、両者の非類似度が大きいと判断し、異常検出対象データに係るセンサ（この一例では「センサＡ」と「センサＢ」）のどちらかまたは両方が異常値を示している特定センサであると推定する。
　また、モード情報取得法は図６Ｅに示すように異常を拡大抽出することができるので、軌道を拡大する効果のある位相面上にこのモード情報取得情報をプロットすることで、時間を横軸、計測値を縦軸にとった通常のプロットでは見分けることが困難なわずかな変化を検出することが可能になる。これによって、推定部１０９は、異常の発生を予兆の段階で検知することができるようになる。この例で、推定部１０９が、特定センサとして推定するのは、センサの組合せである。導出部１０６が、単一のセンサの軌道を導出した場合は、推定部１０９は、同様の方法でセンサの組合せの異常ではなく、単一センサの異常を検出することも可能である。推定部１０９が異常検出対象データのセンサの組合せが異常値を示していると推定した場合、異常値を示しているセンサを特定する方法については後述する。
［（４）非類似度算出部の計算結果を状態遷移マトリクスにより統合して詳細な推定を行う方法について］
　以下、推定部１０９において、非類似度算出部１０８の計算結果を状態遷移マトリクスを用いて統合することによって、異常を示しているセンサを特定し、センサ異常かシステム異常かを判定する方法について説明する。推定部１０９は、非類似度に基づくクラスタリングの結果を用いて、検出対象のセンサの組み合わせが、異常値を示しているか否かを推定した結果を、状態遷移マトリクスを生成して整理することにより、異常値を示しているセンサを特定し、異常の種類や原因の推定をするなどの詳細な分析を行うことができる。例えば、図１２に示すように、状態遷移マトリクスは、センサをそれぞれ縦軸、及び横軸として対応付けて、異常推定結果を書き込んだものである。以下に状態遷移マトリクスの生成法について説明する。
　図１２は、推定部１０９によって生成された状態遷移マトリクスの一例を示す図である。状態遷移マトリクスは、縦軸と横軸とに示される要素が一致する対角要素をはさんで上下対称であり、本実施例では縦軸「センサＡ」と読軸「センサＢ」の組合せと縦軸「センサＢ」と横軸「センサＡ」の組合せは区別しないので、図１２に示されているように、対角要素を含む下半分を生成すればよい。図１２に示す通り、推定部１０９は、状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素（例えば、縦軸「センサＡ」、横軸「センサＡ」の要素（図示する領域ＡＲ１の要素））は個々のセンサが異常値を計測しているか否かを示す要素である。以後このセンサ名が一致する要素をそれぞれ「横軸のセンサ名－縦軸のセンサ名」として、「センサＡ－センサＡ」…「センサＦ－センサＦ」のように記述する。閾値が設定されている場合は、閾値により異常値が計測されているか否かを簡単に判定することができる。例えば、当該センサが閾値を超える異常値を示していない場合にはこの当該センサ名が一致する要素に「〇」を付し、当該センサが閾値を超えた異常値を示している場合には当該要素に「×」を付す。閾値を用いて個々のセンサが異常値を示しているか否かを推定することができない場合は、例えば、以下に述べる方法で、２つのセンサ間の変化に注目した多変量解析の結果を利用して個々のセンサが異常値を計測しているか否かの推定を行う。以下の説明では、図１２のように「センサＡ」～「センサＦ」までの６個のセンサに関する状態遷移マトリクスを生成する場合を考える。相関関係があるセンサは「センサＡ」と「センサＢ」のみであるとする。またここでは、閾値が設定されていないため、単一センサの挙動のみによって異常の推定を行うことが困難なことから、縦軸と横軸とに示されるセンサ名が一致する要素（「センサＡ－センサＡ」…「センサＦ－センサＦ」）の欄が空欄である場合を例に説明を行う。
　推定部１０９では、異常値を示しているセンサを特定して上記の空欄を埋め、センサ故障かシステム異常かといった異常の種類や異常の原因の推定などの詳細な解析を行うため、上記（１）、（２）、（３－１）、（３－２）の手順で得られた結果から、推定部１０９が推定した内容を、状態遷移マトリクスに書き込み情報を整理することができる。例えば、クラスタリングの結果から、推定部１０９が異常の有無を推定した内容は、状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致しない要素（図示する領域ＡＲ２の要素）に書き込む。以後このセンサ名が一致しない要素を「横軸のセンサ名－縦軸のセンサ名」として、「センサＡ－センサＢ」…「センサＡ－センサＦ」のように記述する。推定部１０９は、非類似度算出部１０８によるクラスタリングの結果から、縦軸に示されるセンサと横軸に示されるセンサとの対で生成された位相面軌道形状と通常時の位相面軌道のグループの当該センサ対で生成された位相面軌道形状との非類似度が低く、両者が同じグループにクラスタリングされた場合には、このセンサ対の検出結果は、異常値を示していないと推定して、対応する状態遷移マトリクスの要素に「〇」を付し、通常時の位相面軌道のグループとの非類似度が高く、異常事例として通常事例のグループとは異なる事例に分類された場合は、このセンサ対の検出結果は異常値を示していると推定して対応する状態遷移マトリクスの要素に「×」を付す。
　推定部１０９は、上述の処理を状態遷移マトリクス上のセンサ名が一致しない要素における空欄が埋まるまで繰り返す。次に、この中に「×」の要素があれば、異常値を示しているセンサを特定し、状態遷移マトリックスのセンサ名が一致する要素の部分の空欄を以下の手順ですべて埋めて完成させる。ここでは、センサ名が一致しない要素において、「センサＢ」と他センサを組み合わせた要素がすべて「×」（通常時と異なる）を示し、他の組み合わせがすべて「〇」（通常時と変わらない）となった場合に、センサ名が一致する要素の空欄を埋める方法について説明する。この時点では、縦軸と横軸とに示されるセンサ名が一致する要素「センサＡ－センサＡ」から「センサＦ－センサＦ」の欄）はすべて空欄であるとする。また相関関係情報２０２に基づいて、「センサＡ」と相関関係があるセンサを調べ、これが「センサＢ」のみであり、なおかつ互いに相関関係にあるセンサは「センサＡ」と「センサＢ」のセンサ対以外にないことがわかっているとする。この時、「センサＢ」と相関関係にあるとわかっている「センサＡ」が「センサＢ」以外の他のセンサとつくる組合せの要素に「×」がないので、推定部１０９は、システム異常の場合に表れるセンサ間の相関関係に基づく異常値が検出されていない状況であると判断し、異常は、「センサＢ」自身であり、他センサは異常値を示していないと結論し、「センサＢ」の故障であると推定する。この推定に基づいて、状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素である「センサＢ－センサＢ」の欄に「×」「センサＡ－センサＡ」の欄に「〇」を書き込む。さらに「センサＢ」を除く５個のセンサ（センサＡ、センサＣ～Ｆ）の組合せがすべて「〇」であることから、「センサＣ」～「センサＦ」も異常値を示していないと推定し、「センサＣ－センサＣ」～「センサＦ－センサＦ」の要素に「〇」を書き込み、状態遷移マトリクスを完成させる。この時、システム異常と区別して、センサ故障であることを示すために「センサＢ－センサＢ」の欄には、「〇」、及び「×」とは異なる符号、（例えば、「＃」）などを用いてもよい。この状態遷移マトリクスの結果に基づいて、センサ故障の場合は、センサの交換や修理など必要な整備を行う。センサ故障でない場合は、センサ名とセンサ名が一致する要素の空欄は、以下の手順で埋める。
　センサ故障ではない場合について以下の例で説明する。上述と同様に、「センサＡ」～「センサＦ」までの６個のセンサに関する状態遷移マトリクスを生成する場合を考える。上述と同様、閾値が設定されておらず、センサ単独では異常の検知が困難な場合であり、状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素（「センサＡ－センサＡ」～「センサＦ－センサＦ」まで）が空欄となっているとする。また上述と同様に、相関関係情報２０２に基づいて、「センサＡ」と相関関係があるセンサが「センサＢ」のみであり、なおかつ互いに相関関係があるセンサはこの「センサＡ」と「センサＢ」のセンサ対以外にないことがわかっているとする。このとき上述と同様に、推定部１０９は、（１）～（３）までの処理の結果を状態遷移マトリクスに書き込んでいく。ここでは、「センサＡ」と「センサＢ」を除く４個のセンサ（センサＣ～Ｆ）の組合せの要素が、すべて「〇」でそれ以外は「×」がついている場合について説明する。この場合「センサＢ」と相関関係のある「センサＡ」も異常値を示していることから、因果関係の破綻は認められず、互いに相関関係にある「センサＡ」と「センサＢ」のまわりでシステム異常が発生していると推定し、状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素の「センサＡ－センサＡ」と「センサＢ－センサＢ」の欄に「×」を付す。さらに、「センサＡ」と「センサＢ」を除く４個のセンサ（センサＣ～Ｆ）の組合せの欄（例えば、「センサＣ－センサＤ」、「センサＣ－センサＦ」など）はすべて「〇」となっていることを確認し、「状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素の「センサＣ－センサＣ」「センサＤ－センサＤ」「センサＥ－センサＥ」「センサＦ－センサＦ」の要素に「〇」を付し、状態遷移マトリクスを完成させる。次に完成した状態遷移マトリクスを利用して、異常の原因を推定してそれに対する対応を行う。状態遷移マトリクスを利用した異常の原因の推定方法について、以下に例をあげて説明する。
　上記の手順により、システム異常と推定された場合は、完成した状態遷移マトリクスを利用して、以下に述べるような方法で、異常の原因の推定を行うことができる。上記の例で、「センサＡ」と「センサＢ」が２つとも異常値を示していると推定されたとする。このとき、相関関係情報２０２から「センサＡ」が因果関係の上流にあり「センサＢ」が因果関係の下流にあることがわかっていれば、「センサＡ」の変化により「センサＢ」が変化したことが推定されるので、原因は「センサＡ」まわりにあることが推定される。参考情報２０６も調べ、「センサＡ」まわりで部品の交換が長く行われておらず、交換時期が近づいていたなどの整備履歴があれば、この部品の劣化が原因である可能性が高いと推定されるため、異常検出の結果を受けて、最初にこの部品の点検を行うべきであるという判断ができる。直接相関や間接相関にあるセンサが複数あるときも同様にして、相関関係情報２０２、参考情報２０６などを利用して、状態遷移マトリクスから異常の原因を推定することができる。
　以上のように上記（１）～（４）の方法によって、閾値の設定が困難でセンサ単独での異常推定が難しい場合であっても、異常値を示しているセンサの推定、センサ故障かシステム異常かの推定、システム異常の場合には、異常の原因の推定が可能になる。
［処理フロー］
　以下、異常検出装置１の動作の内容について説明する。図１３は、異常検出装置１の処理の一連の流れを示すフローチャートである。例えば、取得部１０２は、異常検出対象のログ情報から検出結果情報３００を取得する（ステップＳ１００）。次に、解析部１０４、導出部１０６、非類似度算出部１０８を経て、取得部１０２によって取得された検出結果情報３００を解析した結果に基づいて、推定部１０９が異常値を示している特定センサを推定する（ステップＳ１０２）。ステップＳ１０２の処理の詳細については、後述する。出力部１１０は、推定部１０９によって特定センサであると推定されたセンサに係る情報を出力する（ステップＳ１０４）。
　以下、（１）多様な異常の特徴を抽出するためのモード情報取得処理を行って、モード情報取得情報を軌道形状を利用して可視化し、可視化した軌道形状を定量した軌道形状の数値化情報に基づいて非類似度を算出し、算出した非類似度に基づくクラスタリングの結果を利用して状態遷移マトリクスを生成する処理（２）状態遷移マトリクスを用いた詳細判定の順で、ステップＳ１０２の処理を実行する方法の一例について説明する。
［処理フロー：（１）多様な異常の特徴を抽出するためのモード情報取得処理を行って、モード情報取得情報を軌道形状を利用して可視化し、可視化した軌道形状を定量した軌道形状の数値化情報に基づいて非類似度を算出し、算出した非類似度に基づくクラスタリングの結果を利用して状態遷移マトリクスを生成する処理］
　図１４は、ステップＳ１０２に係る特定センサの推定処理の一連の流れを示すフローチャートである。以下、図１４を参照し、解析部１０４が各センサから得られた時系列情報に対して「第１種モード情報取得処理」により異常の特徴を抽出し、導出部１０６が「モード情報取得処理」が行われていない０次モードも含めた複数モードの中から適切と思われる組み合わせを選んで２次元位相面軌道を生成し、非類似度算出部１０８によって軌道形状を数値化し、これに基づいて算出された非類似度に基づくクラスタリングを行った結果を利用して、推定部１０９が状態遷移マトリクスを生成する場合を例にとって説明する。図１３の一例において、閾値が設定されておらず、単一センサの時系列情報からは、センサが異常値を示しているかどうかの推定が難しい場合について説明する。閾値が設定されていない場合は、縦軸と横軸とに示される要素が一致する対角要素の異常推定は、この対角要素より下の部分の状態遷移マトリクスの完成を待って行われる。なお上述したように、状態遷移マトリクスは、縦軸と横軸とに示される要素が一致する対角要素をはさんで上下対称であり、本実施例では「センサＡ－センサＢ」と「センサＢ－センサＡ」は区別しないので、図１２に示すように、対角要素を含む下半分だけを生成すればよい。また上記の実施例に示したように、図１２と同じように、「センサＡ」から「センサＦ」までの６個のセンサの計測値に対して、異常推定を実施する場合を想定して以下の説明を行う。例えば、「第１種モード情報取得処理」を行うため、解析部１０４は、取得部１０２が取得した検出結果情報３００を滑走窓法を利用してセンサ毎に時間ベクトルに変換する（ステップＳ３００）。解析部１０４は、この変換により多変量データに再構成されたデータセットが互いに独立であるとみなして、センサ毎に行列分解法を利用した多変量解析法を用いて、多様な異常の特徴を複数のモード情報として抽出する（ステップＳ３０２）。
　次に、導出部１０６は、解析部１０４によって取得されたセンサ毎の複数モード情報を利用して、位相面軌道を導出するために縦軸、及び横軸それぞれに使用するモードを選定縦軸、及び横軸センサの中から２つずつ組み合わせて、２次元位相平面上に２つのセンサ間の変化を表した位相面軌道を導出する（ステップＳ３０４）。導出部１０６が導出する位相面軌道は、例えば「センサＡ」のモード１結果ｘを第１軸（例えば、横軸）とし、「センサＢ」のモード１結果ｙを第２軸（例えば、縦軸）とした位相平面である。導出部１０６は、モード１結果ｘに含まれる各要素｛ｘ１、ｘ２、…、ｘｎ｝と、モード１結果ｙに含まれる各要素｛ｙ１、ｙ２、…、ｙｎ｝との中で検出タイミングが合致する要素によって示される座標（座標（ｘ１，ｙ１）、…、座標（ｘｎ，ｙｎ））を、検出時間順に直線によって結ぶことにより、軌道を導出する。
　非類似度算出部１０８は、導出された位相面軌道形状を複素自己回帰係数のような定量手法で数値化する（ステップＳ３０６）。非類似度算出部１０８は、数値化した位相面軌道形状と、軌道形状情報２０４に示される、数値化した複数の通常時の軌道形状との間の非類似度（例えば複素パワーケプストラム距離Ｄｃ）を計算して、通常情報との距離の大きさにより測定結果をクラスタリングし、図１１に示すデンドログラムのような可視化手法を用いてクラスタリング結果を可視化する（ステップＳ３０８）。推定部１０９は、この分類結果を用いて検出対象が、複数の通常事例がつくるグループとは、異なる事例として分類されたとき、通常事例に対する非類似度が大きいとして、当該センサ対のどちらかまたは両方が異常値を計測していると推定し、通常事例が創るグループと同じグループに分類されたとき、通常事例に対する非類似度が小さいとして当該センサ対のどちらも異常値を計測していないと推定する（ステップＳ３１０）。
　推定部１０９では、クラスタリングの結果を図１２に示すような状態遷移マトリクスに記入する。状態遷移マトリクスとは、上述したように、対象とするセンサ対を縦軸、及び横軸それぞれに対応付けて、異常推定結果を書き込んだものである。例えば「センサＡ」と「センサＢ」のセンサ対からなる軌道形状が通常時との非類似度が大きい（つまり異常値を示している）と推定された場合には、状態遷移マトリクスの該当箇所（横軸「センサＡ」、縦軸「センサＢ」の要素、「センサＡ－センサＢ」の欄）に「×」を記入する（ステップＳ３１２）。また異常値でないと推定された場合は、状態遷移マトリクスの該当箇所に「〇」を記入する（ステップＳ３１４）。
　推定部１０９は、ステップＳ３１２とステップＳ３１４において生成した状態遷移マトリクスが完成したか否かを判定する（ステップＳ３１６）。推定部１０９は、状態遷移マトリクスが完成するまでの間、ステップＳ３０４からＳ３１６の処理を繰り返す。以降の説明では、異常検知を実施する者が、異常判定を行う際に事前に、状態遷移マトリクス作成に使用するモードを複数種類選定しておき、推定部１０９が図１２に示すような「センサＡ」から「センサＦ」がつくる状態遷移マトリクスを生成する方法を説明する。この場合には、選定されたモード毎にひとつずつ状態遷移マトリクスが生成される。解析部１０４は、選定したすべてのモードで状態遷移マトリクスが完成したかを判定する（ステップＳ３１８）。推定部１０９は、必要とされるすべてのモードで状態遷移マトリクスが完成されるまでの間、ステップＳ３０４からステップＳ３１８の処理を繰り返す。状態遷移マトリクスの縦軸と横軸とに示されるセンサ名が一致しない要素のうち、一致する要素を対角要素とすると、この対角要素より下の部分の記入が終わった時点で、状態遷移マトリクスはいったん完成したものと判定し、次のステップＳ３２０でこの状態遷移マトリクスを用いて縦軸と横軸とに示されるセンサ名が一致する対角要素の異常推定などを行う。上述したように、閾値が設定されておらず、単一センサの時系列情報からは、異常値を示しているかどうかの判定が難しく、縦軸と横軸とに示されるセンサ名が一致した対角要素が空欄となっている場合について説明する。もし、検出したい異常が事前に決まっており、調べるべき最適モードが事前にわかっている場合は、必要なモードでのみ状態遷移マトリクスを用意すればよい。
　状態遷移マトリクスの各要素のうち、縦軸と横軸とに示されるセンサ名が一致する要素（例えば、上記で定義した「センサＡ－センサＡ」の要素）に関する推定部１０９の推定は、閾値が設定されている場合は、閾値を超えている場合に異常値（「×」）、閾値を超えていない場合に異常値でない（「〇」）と推定し、状態遷移マトリクスの該当箇所に結果を記入すればよい。一方、推定部１０９は、閾値が設定されておらず、単独のセンサ情報から、異常値か否かの判定が困難な場合は、状態遷移マトリクスのセンサ名が一致しない要素（例えば、「センサＡ－センサＢ」「センサＡ－センサＣ」などの要素）がすべて埋まってから、以下のセンサ故障かシステム異常かの総合的な判断とあわせて異常値を示しているセンサを特定する。
　推定部１０９は、縦軸と横軸とに示されるセンサ名が一致しない要素の記入が終わった状態推定マトリクスを用いて、相関関係情報２０２に格納されたセンサ間の相関関係に関する情報や、センサやその周辺部品の整備履歴などの情報が格納された参考情報２０６などに基づいて、異常値を示しているセンサの特定（縦軸と横軸とに示されるセンサ名が一致した要素の異常推定）や、センサ故障かシステム異常かの判定、システム故障と判定された場合は原因の推定を行う（ステップＳ３２０）。
［処理フロー：（２）状態遷移マトリクスを用いた詳細判定］
　以下、推定部１０９で状態遷移マトリクスを用いた詳細判定を行うステップＳ３２０の処理について説明する。図１５は、ステップＳ３２０に係る特定センサの推定処理の一連の流れを示すフローチャートである。図１５に示されるフローチャートにおいて、各種指標値や軌道の算出対象であるセンサが、「センサＡ」「センサＢ」「センサＣ」「センサＤ」「センサＥ」「センサＦ」の６個であり、相関関係情報２０２によれば、この中で相関関係にあるセンサは、「センサＡ」「センサＢ」のみであることがわかっている場合を例に説明する。この例では、複数のセンサ故障は、起きない場合について説明する。推定部１０９は、センサ名が一致しない要素の記入が終わった状態遷移マトリクス情報を取得し（ステップＳ４００）、「×」の付されたセンサ対があるかどうかを調べる（ステップＳ４０２）。推定部１０９は、「×」の付されたセンサ対がひとつもなければ、異常なしと推定し（ステップＳ４２８）、事前に選定された全モードの状態遷移マトリクスで検証を行ったかの判定へすすむ（ステップＳ４２０）。ステップＳ４２０以降の手順は、以下のステップを合わせて説明する。「×」の付されたセンサ対があれば、そのうちのひとつ選ぶ（例えば、縦軸「センサＡ」と横軸「センサＢ」の要素から「センサＡ」と「センサＢ」のセンサ対）（ステップＳ４０４）。選んだセンサ対の一方（例えば、「センサＢ」）を特定センサであるか否かを推定する推定対象のセンサ（以後第１センサとよぶ）とする（ステップＳ４０６）。このとき、同一の状態遷移マトリクスの同一のペアですでに第１センサに選ばれたセンサは第１センサには選ばない。
　推定部１０９は、第１センサを含まないセンサ対にも「×」があるかどうかを調べ（ステップＳ４０８）、第１センサを含まないセンサ対に「×」があれば、相関関係情報２０２に基づいて、第１センサと相関関係のあるセンサ（この例では「センサＡ」）を調べ、状態遷移マトリクス上で、念のためこのセンサが第１センサ以外のセンサとつくるセンサ対の欄がすべて「〇」であるか否かを調べる（ステップＳ４１０）。すべて「〇」である場合は、存在するべき相関関係（あるいは因果関係）が喪失しているので、推定部１０９は、第１センサの「センサ故障」と判定する（ステップＳ４１２）。このとき、推定部１０９は、判定結果を状態遷移マトリクスの該当欄（第１センサが「センサＢ」である場合には、状態遷移マトリクスの「センサＢ－センサＢ」の欄）に「×」を付し、センサ故障であることを付記する（ステップＳ４１４）。または上述のように、センサ周辺のシステム異常と区別して、センサ故障であることを示すために別な符号、例えば、「＃」などを用いてもよい。すべて「〇」ではなかった場合は（Ｓ４０８で第１センサを含むセンサ対にのみ「×」があるとしているのでこの場合はないはずであり）、Ｓ４０８のＮＯに該当するので、ステップＳ４２４に合流して詳しい解析を行う。状態遷移マトリクス上で、第１センサ以外のセンサとつくるセンサ対の欄がすべて「〇」で、第１センサが故障していると推定された場合は、Ｓ４１４で状態遷移マトリクスへの記録を終えたあと、ステップＳ４１６に進み、Ｓ４０４で選んだセンサ対の両センサを第１センサとして検証したかを確認し、両センサを第１センサとして検証が終わっていなければ、ステップＳ４０６に戻って検証を続ける。両センサで検証が終わっていれば、Ｓ４０４に戻り、別な「×」の付されたセンサ対を、選んで検証を続ける。これもすべて終わっていれば、ステップＳ４２０に進み、事前に選定した全モードで検証が終わっているかを判定し、終わっていなければＳ４００に戻る。終わっていれば、ステップＳ４２２に進み判定結果として完成した状態遷移マトリクスと判定結果を出力部１１０で出力する。ここでは、単一センサ故障を例に説明したが、複数センサが故障している場合も、推定部１０９は、同様の手順で推定を行う。
　この場合、出力部１１０は、特定センサであると推定されたセンサに係る情報をすべて出力してもよく、特定センサであると推定されたセンサのうち、一部のセンサに係る情報を出力してもよい。
　ステップＳ４０８で第１センサを含むセンサ対以外にも「×」が付されていれば、相関関係情報２０２に基づいて、その中に第１センサと相関関係のあるセンサが含まれているかを調べる（ステップＳ４２４）。第１センサと相関関係のあるセンサと第１センサ以外のセンサがつくるセンサ対の欄にも、「×」が付されている場合は、センサ間の相関関係の破綻が認められないので、センサ故障ではなく、第１センサまわりの「システム異常」であると推定し、ステップＳ４１４に移り、状態遷移マトリクスの対角要素の第１センサの欄（第１センサが「センサＡ」なら縦軸と横軸のセンサ名が一致する「センサＡ―センサＡ」の欄）に「×」を付すとともに、センサ故障ではなくシステム故障であることを記録して、ステップＳ４１６以降の処理を上述と同様に実行する。
　上記の処理により、「センサ故障」と推定された場合は、点検作業の際にセンサ故障が推定されるセンサに対してセンサ故障の原因を実際のセンサで確認し、必要に応じてセンサの部品の交換や、センサそのものの交換などの対処がなされる。「システム異常」と推定された場合は、状態遷移マトリクスを利用して、異常値を検出したセンサを特定するとともに、センサ間の相関関係情報２０２や異常値を示しているセンサまわりの装置やシステムの整備履歴などが納められた参考情報２０６に基づいて異常原因を総合的に推定することができる。例えば、システム異常と診断され、内部漏洩を測定する圧力センサと流量計が異常値を示しており、参考情報２０６からシール部品の交換期限が迫っていることがわかった場合には、シールの摩耗による内部漏洩が疑われるため、点検作業において、シール部分の点検から故障個所の特定を進めることで、効率よくシステムの整備を行うことができる。
［ステップＳ１０２の処理について］
　なお、上述では、特定センサを推定するステップＳ１０２において、解析部１０４における（１）「多様な異常の特徴を抽出するためのモード情報取得」（２）「軌道形状を利用したモード情報取得情報の可視化」（３－１）「モード情報取得情報により生成される軌道形状の定量化」（３－２）「軌道形状の数値化情報に基づく非類似度の算出とこれに基づくクラスタリング」（４）「非類似度算出部の計算結果を状態遷移マトリクスにより統合して実行する詳細な推定処理」を順番に実行する場合のフローチャートについて説明したが、本発明の実施例は、これに限定されない。ステップＳ１０２では、（１）～（４）の手順のうち、異常検出の目的に応じて、適切な手順を組み合わせてフローチャートを実行することができる。異常の性質によっては、高次モード変換を行わない０次モードの情報を利用して（あるいは高次モード情報と適当に組み合わせて）位相面軌道を生成するほうが適切な場合もある。また必ずしもすべてのモードについて状態遷移マトリクスを生成する必要もない。検知したい異常の特徴が分かっている場合には、適切なモードの組み合わせを選んで、異常検出を行うことで検知にかかる計算処理時間を短縮することができる。また必ずしもすべてのセンサ対について位相面軌道を生成する必要もない。ここで実施例として説明したのは記録された情報をオフラインで解析するのに適した手法であったが、例えば、オンラインでの素早い異常検出による緊急停止が目的であれば、対象とする故障モードを限定し、解析対象のセンサ対を絞ることで、異常検出にかかる時間を短縮することができる。０．０１秒を争う検知が必要でなければ効率的な整備を行うための異常検出やその原因の推定、あるいは異常の予兆の検知などが目的であれば、多少時間がかかっても、網羅的に解析を行う方法が望ましい。またモード情報取得処理の前に検出結果の平均と分散をそろえる標準化のようなモード情報取得処理以外の処理を加えてもよいし、モード情報取得処理後にこのような処理を行ってもよい。モード情報取得処理にも上述のように「第１種モード情報取得処理」、「第２種モード情報取得処理」、およびこれらを組み合わせた処理があり、位相面軌道の導出にもどのようなモード情報を選ぶかで任意の情報を組み合わせて軌道を導出できることや、位相面を用いない通常の時系列情報がつくる軌道からも同様の処理で非類似度が計算できることも述べている。
［相関関係の推定について］
　上記、２つのセンサ間の位相面軌道形状を利用した異常検出の処理フローは、２つのセンサ間の相関関係を推定する際にも用いることができる。例えば、図１１は、複数の通常事例間で「センサＡ」と「センサＢ」がつくる位相面軌道の非類似度の判定を行った結果を示したものだが、検出対象と複数の通常事例間の「センサＡ」と「センサＢ」のつくる位相面軌道形状の非類似度の比較ではなく、通常事例における全センサの組合せについて、位相面軌道形状を生成し、物理的な関係式などから相関関係が明らかなセンサ対をモデルペアとして、複数のモードでの非類似度を比較することにより、相関関係のあるグループとないグループにクラスタリングすることが可能である。このように上記異常検出装置を用いて解析した結果を、上記状態遷移マトリクス上に相関関係があると推定される場合に「〇」、相関関係がないと推定される場合に「×」を書き込むことにすれば、センサ間の相関関係マトリクスを上記異常検出装置を用いて生成することができる。複数モードで相関関係を調べることにより、ノイズが多く、判断が難しい場合も、相関関係の判定が可能になる。相関関係表の生成は、異常判定と並行して行うこともできるが、予め通常データを用いて実行しておき、相関関係情報２０２に記録しておいてもよい。
［モードの抽出について］
　また、上述では、「第１種モード情報取得処理」と名付けた、滑走窓毎の時系列データをベクトルの集まりに変換することにより得られたデータを、近傍の時間情報を変数としてもつ多変量データとして再構成することによって行列分解法を用いた多変量解析を行うことで複数モードにより多様な異常を抽出する例を説明したが、「第２種モード情報取得処理」も含めて、多様なモード情報を抽出する方法は、主成分分析に限られない。導出部１０６は、因子分析や、スパース主成分分析、ファジィ主成分分析、カーネル主成分分析、確率的主成分分析、ロバスト主成分分析など主成分分析と同様のモード情報取得が可能な行列分解手法を用いた多変量解析により複数モードの特徴抽出ができる方法を適用して、異常の特徴抽出を行ってもよい。また実施例の中で説明したように「第２種モード情報取得」と名付けた公知の時間領域での主成分分析や因子分析の手法に「第１種モード情報取得」の結果を利用することで、ラグ構造を導入できることから、センサの測定値がラグを伴って作用しあうような場合の異常検出にも力を発揮する。周波数の変化や、ラグを伴って発生する異常など検知する異常の性質にふさわしいモードを選定することで、見逃しや誤検知を減らし、正確な異常検出が可能になる。また「第１種モード情報取得処理」には、ラグ構造を導入して多変量データに再構成し、主成分分析をかける一連の処理には、ノイズフィルタのような効果が生まれるため、例えば窓幅などによって決まる多変量データに再構成する際のベクトル化の次元を検出結果に応じて調整することで、ノイズの多いデータからノイズを除去し、誤検知を減らすこともできる。
［特定センサの他の推定方法について］
　なお、上述では、解析部１０４、導出部１０６、非類似度算出部１０８、推定部１０９が（１）～（４）の手順によって特定センサを推定する場合について説明したが、本発明の実施例は、これに限られない。解析部１０４、導出部１０６、非類似度算出部１０８、推定部１０９は、例えば、検出結果情報３００のような情報を取得すると、通常時の検出結果と非類似度の高い異常値を計測しているセンサを、特定センサとして出力するように学習された学習モデルに基づいて、特定センサを推定してもよい。この場合、推定部１０９は、当該学習モデルに検出結果情報３００に示される検出結果を入力し、出力として得られたセンサを特定センサとして推定する。
　以上、本発明を実施するための形態について実施形態を用いて説明したが、本発明はこうした実施形態に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において種々の変形及び置換を加えることができる。
１…異常検出装置、１００…制御部、１０２…取得部、１０４…解析部、１０６…導出部、１０８…非類似度算出部、１０９…推定部、１１０…出力部、２００…記憶部、２０２…相関関係情報、２０４…軌道形状情報、２０６…参考情報、３００…検出結果情報、Ｄｃ…複素パワーケプストラム距離、Ｏｒｂ１…軌道、Ｏｒｂ２…軌道、ｘ、ｙ…モード１結果、ｚ１、ｚ２…複素自己回帰係数、ｚ１、ｚ２…複素パワーケプストラム距離

Claims (21)

  1. 　複数のセンサの各々が検出対象の状態を所定時間毎に検出した複数の検出結果を取得する取得部と、
    　前記取得部によって取得された複数の前記検出結果のそれぞれに対して、異常の性質に応じた特徴を抽出する複数のモード情報取得処理を行って得られた複数のモード情報に基づいて、前記複数のセンサのうち、異常の可能性がある特定センサを推定する推定部と、
    　を備える異常検出装置。
  2. 　前記取得部により所定時間毎に取得された前記検出結果の時系列情報を時間構造によって前記センサ毎の多変量データに再構成し、再構成した前記多変量データを多変量解析手法によって解析することにより、前記複数のモード情報を取得する解析部を更に備え、
    　前記推定部は、前記解析部によって取得された前記複数のモード情報に基づいて、前記特定センサを推定する、
    　請求項１に記載の異常検出装置。
  3. 　前記解析部は、滑走窓法によって前記多変量データに前記時間構造を導入する、
    　請求項２に記載の異常検出装置。
  4. 　前記解析部は、前記滑走窓法の窓幅を、所定値以上にすることにより、前記多変量データのノイズを除去する、
    　請求項３に記載の異常検出装置。
  5. 　前記解析部は、２つ以上の前記センサを選択したセンサグループ毎に、前記センサグループの同一時刻の検出結果を検出時刻毎に組み合わせることで、前記検出結果を多変量データに再構成し、再構成された前記多変量データを多変量解析手法によって解析することにより、前記複数のモード情報を取得する、
    　請求項２から請求項４のいずれか一項に記載の異常検出装置。
  6. 　前記解析部は、２つ以上の前記センサを選択したセンサグループ毎に、時間構造の導入によって前記センサ毎に再構成した前記多変量データから多変量解析により得られたモード情報を検出時刻毎に組み合わせることで、前記センサ毎のモード情報を多変量データに再構成し、再構成された前記多変量データを多変量解析手法によって解析することにより、前記複数のモード情報を取得する、
    　請求項５に記載の異常検出装置。
  7. 　前記解析部は、前記複数のセンサの中で、所定の数のセンサについて、所定の種類のセンサについて、又は互いに相関関係にあるセンサについて、前記モード情報を取得する、
    　請求項２から請求項６のうちいずれか一項に記載の異常検出装置。
  8. 　前記モード情報を２つ組合せ、一方のモード情報に含まれる要素をそれぞれ、２次元の位相平面の第１軸上の値とし、他方のモード情報に含まれる要素をそれぞれ前記位相平面の第２軸上の値とし、同時刻の要素を座標とした軌道を導出する導出部を更に備え、
    　前記推定部は、前記導出部によって導出された前記軌道に基づいて、前記特定センサを推定する、
    　請求項１から請求項７のうちいずれか一項に記載の異常検出装置。
  9. 　前記モード情報の計測時刻に係る情報を、座標平面の第１軸とし、前記モード情報の計測値に係る情報を、前記座標平面の第２軸とする２次元平面上の軌道を導出する導出部を更に備える、
    　請求項１から請求項８のうちいずれか一項に記載の異常検出装置。
  10. 　前記導出部は、
    　相関関係のある前記センサに関して前記軌道を導出する、
    　請求項８又は請求項９に記載の異常検出装置。
  11. 　前記推定部は、前記導出部によって導出された前記軌道の形状が、通常状態の軌道の形状と異なる場合、前記軌道に係る前記センサを、前記特定センサとして推定する、
    　請求項８から請求項１０のいずれか一項に記載の異常検出装置。
  12. 　前記軌道の形状に基づいて、非類似度を算出する非類似度算出部を更に備え、
    　前記推定部は、前記非類似度算出部によって算出された前記非類似度が通常時のデータ相互間の非類似度と比して大きい場合、当該非類似度に係る前記センサを、前記特定センサとして推定する、
    　請求項８から請求項１１のうちいずれか一項に記載の異常検出装置。
  13. 　前記非類似度算出部は、クラスタリング手法によって複数の通常データと検出対象データを前記非類似度に基づいて分類し、
    　前記推定部は、検出対象が通常データと異なるクラスタに分類された場合、前記センサを、前記特定センサとして推定する、
    　請求項１２に記載の異常検出装置。
  14. 　相関関係がある前記センサを示す相関関係情報を生成する相関関係情報生成部を更に備え、
    　前記相関関係情報生成部は、前記非類似度算出部によって算出された相関関係にあるセンサの組合せがつくる位相面軌道形状と相関関係の判定をするセンサの組合せがつくる位相面軌道形状を非類似度によりクラスタリングした結果に基づいて、相関関係の判定をするセンサの組合せが、前記推定部により相関関係のあるセンサグループに分類されていると推定された結果を、前記相関関係情報として生成する、
    　請求項１３に記載の異常検出装置。
  15. 　前記推定部は、前記センサの組合せ毎に、前記センサ、及びセンサの組合せが異常値を示しているか否かを推定した要素を有する状態遷移マトリクスを生成し、生成した前記状態遷移マトリクスに基づいて、前記特定センサの推定、特定センサに係る異常の原因、又は異常の種類を推定する、
    　請求項８から請求項１４のうちいずれか一項に記載の異常検出装置。
  16. 　前記推定部は、前記状態遷移マトリクスにおいて、前記センサ間の相関関係の破綻の有無に基づいて、前記特定センサに係る異常の原因、又は異常の種類を推定する、
    　請求項１５に記載の異常検出装置。
  17. 　前記推定部は、前記状態遷移マトリクスに異常値を示している要素があり、且つ当該要素に係るセンサと相関関係にあるセンサの要素が異常値を示していない場合、センサ間の相関関係が破綻していると推定し、異常値を示している要素に係るセンサが故障していると推定する、
    　請求項１５、又は請求項１６のうちいずれか一項に記載の異常検出装置。
  18. 　前記推定部は、前記状態遷移マトリクスに異常値を示している要素があり、且つ当該要素に係るセンサと相関関係にあるセンサの要素が異常値を示している場合、センサ間の相関関係が破綻していないと推定し、前記センサの検出対象のシステムが故障していると推定する、
    　請求項１５から請求項１７のうちいずれか一項に記載の異常検出装置。
  19. 　コンピュータに、
    　複数のセンサの各々が検出対象の状態を所定時間毎に検出した複数の検出結果を取得させ、
    　取得された複数の前記検出結果のそれぞれに対して、異常の性質に応じた特徴を抽出する複数のモード情報取得処理を行って得られた複数のモード情報に基づいて、前記複数のセンサのうち、異常の可能性がある特定センサを推定させる、
    　プログラム。
  20. 　複数のセンサの各々が検出対象の状態を所定時間毎に検出した複数の検出結果を取得する取得部と、
    　前記取得部により取得された前記検出結果の時系列情報を時間構造によって前記センサ毎の多変量データに再構成し、再構成した前記多変量データを主成分分析法によって解析することにより、前記センサ毎に特徴が抽出された主軸上のスコアの時系列情報であるモード情報を取得する解析部と、
    　前記解析部で取得された複数のモード情報から選択した２つのモード情報について、一方のモード情報に含まれる要素をそれぞれ、２次元の位相平面の第１軸上の値とし、他方のモード情報に含まれる要素をそれぞれ前記位相平面の第２軸上の値とし、同時刻の要素を座標とした軌道を導出する導出部と、
    　前記導出部によって導出された前記軌道の形状が、通常状態の軌道の形状と異なる場合、前記軌道に係る２つのセンサを、異常の可能性がある特定センサとして推定する推定部と、
    　を備える異常検出装置。
  21. 　前記導出部は、前記複数のセンサそれぞれについて、当該センサに係るモード情報および他のセンサに係るモード情報に基づく軌道を導出し、
    　前記推定部は、前記導出部によって導出されたそれぞれの軌道に基づいて、当該センサおよび他のセンサを特定センサとして推定し、
    　前記複数のセンサそれぞれについて、当該センサと他のセンサの組合せ毎に、前記センサの組合せが特定センサとして推定したか否かを示す要素を有する状態遷移マトリクスを生成し、
    生成した前記状態遷移マトリクスに基づいて、異常の可能性のある１のセンサを推定する、
    　請求項２０に記載の異常検出装置。
PCT/JP2019/042252 2018-10-30 2019-10-29 異常検出装置、異常検出方法、およびプログラム WO2020090770A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
EP19880340.5A EP3876056A4 (en) 2018-10-30 2019-10-29 FAULT DETECTION SYSTEM, FAULT DETECTION METHOD, AND PROGRAM
US17/288,611 US11669080B2 (en) 2018-10-30 2019-10-29 Abnormality detection device, abnormality detection method, and program
JP2020553912A JP7340265B2 (ja) 2018-10-30 2019-10-29 異常検出装置、異常検出方法、およびプログラム
CN201980071540.XA CN112955839A (zh) 2018-10-30 2019-10-29 异常检测装置、异常检测方法和程序

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2018-204515 2018-10-30
JP2018204515 2018-10-30

Publications (1)

Publication Number Publication Date
WO2020090770A1 true WO2020090770A1 (ja) 2020-05-07

Family

ID=70464517

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/042252 WO2020090770A1 (ja) 2018-10-30 2019-10-29 異常検出装置、異常検出方法、およびプログラム

Country Status (5)

Country Link
US (1) US11669080B2 (ja)
EP (1) EP3876056A4 (ja)
JP (1) JP7340265B2 (ja)
CN (1) CN112955839A (ja)
WO (1) WO2020090770A1 (ja)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023189022A1 (ja) * 2022-03-29 2023-10-05 三菱ケミカルエンジニアリング株式会社 診断装置、診断方法、及び診断プログラム
JP7491065B2 (ja) 2020-06-04 2024-05-28 株式会社豊田中央研究所 状態推定装置、及び状態推定方法、状態推定プログラム

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020166011A1 (ja) * 2019-02-14 2020-08-20 日本電気株式会社 時系列データ処理方法
CN112348644B (zh) * 2020-11-16 2024-04-02 上海品见智能科技有限公司 一种通过建立单调正相关过滤网的异常物流订单检测方法
CN114052676B (zh) * 2021-11-19 2024-05-07 南开大学 一种中医脉搏精简阵列传感器及其全阵列脉搏信息获取算法
US11743344B1 (en) * 2022-03-15 2023-08-29 International Business Machines Corporation Edge resource processing
CN115755752B (zh) * 2023-01-06 2023-04-14 山东鸿德电力科技有限公司 一种基于plc的自动化设备节能控制方法及系统
CN115979350A (zh) * 2023-03-20 2023-04-18 北京航天华腾科技有限公司 一种海洋监测设备数据采集系统
CN117250576B (zh) * 2023-11-16 2024-01-26 苏芯物联技术(南京)有限公司 一种基于多维传感数据的电流传感器实时异常检测方法
CN117609831A (zh) * 2023-11-21 2024-02-27 国网宁夏电力有限公司电力科学研究院 一种基于随机矩阵的变压器局放检测方法、介质及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010191556A (ja) * 2009-02-17 2010-09-02 Hitachi Ltd 異常検知方法及び異常検知システム
US20110270579A1 (en) * 2010-04-28 2011-11-03 Nellcor Puritan Bennett Ireland Systems and methods for signal monitoring using lissajous figures
WO2014034273A1 (ja) * 2012-08-29 2014-03-06 株式会社日立製作所 設備状態監視方法及び設備状態監視装置
JP2017156201A (ja) 2016-03-01 2017-09-07 旭化成エレクトロニクス株式会社 異常診断装置、異常診断方法、角度検出装置、およびプログラム
JP2018081523A (ja) * 2016-11-17 2018-05-24 株式会社日立製作所 設備診断装置及び設備診断方法
JP2018204515A (ja) 2017-06-02 2018-12-27 株式会社デンソー 送風装置

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3784248B2 (ja) * 2000-10-02 2006-06-07 株式会社ジェイテクト 回転角度検出装置、トルクセンサ及び舵取装置
JP4071449B2 (ja) 2001-03-27 2008-04-02 株式会社東芝 センサ異常検出方法及びセンサ異常検出装置
JP2004045041A (ja) * 2002-07-08 2004-02-12 Koyo Seiko Co Ltd 回転角検出装置及びトルク検出装置
JP2007026134A (ja) * 2005-07-19 2007-02-01 Matsushita Electric Works Ltd 異常判定装置
JP2008171074A (ja) * 2007-01-09 2008-07-24 National Institute Of Advanced Industrial & Technology 三次元形状モデル生成装置、三次元形状モデル生成方法、コンピュータプログラム、及び三次元形状モデル生成システム
JP5793299B2 (ja) * 2010-12-28 2015-10-14 株式会社東芝 プロセス監視診断装置
JP5538597B2 (ja) * 2013-06-19 2014-07-02 株式会社日立製作所 異常検知方法及び異常検知システム
US10281909B2 (en) * 2013-07-10 2019-05-07 Globiz Co., Ltd. Signal measuring/diagnosing system, and method for applying the same to individual devices
US10878385B2 (en) * 2015-06-19 2020-12-29 Uptake Technologies, Inc. Computer system and method for distributing execution of a predictive model
CN105487524B (zh) * 2015-12-29 2017-12-08 浙江中烟工业有限责任公司 具有多工况特性的超高速小盒包装机状态监测与诊断方法
US11327475B2 (en) * 2016-05-09 2022-05-10 Strong Force Iot Portfolio 2016, Llc Methods and systems for intelligent collection and analysis of vehicle data
JP7017861B2 (ja) * 2017-03-23 2022-02-09 株式会社日立製作所 異常検知システムおよび異常検知方法
CN108594790B (zh) * 2018-04-11 2019-12-10 浙江大学 一种基于结构化稀疏型主元分析的故障检测和分离方法
US11579588B2 (en) * 2018-07-30 2023-02-14 Sap Se Multivariate nonlinear autoregression for outlier detection

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010191556A (ja) * 2009-02-17 2010-09-02 Hitachi Ltd 異常検知方法及び異常検知システム
US20110270579A1 (en) * 2010-04-28 2011-11-03 Nellcor Puritan Bennett Ireland Systems and methods for signal monitoring using lissajous figures
WO2014034273A1 (ja) * 2012-08-29 2014-03-06 株式会社日立製作所 設備状態監視方法及び設備状態監視装置
JP2017156201A (ja) 2016-03-01 2017-09-07 旭化成エレクトロニクス株式会社 異常診断装置、異常診断方法、角度検出装置、およびプログラム
JP2018081523A (ja) * 2016-11-17 2018-05-24 株式会社日立製作所 設備診断装置及び設備診断方法
JP2018204515A (ja) 2017-06-02 2018-12-27 株式会社デンソー 送風装置

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7491065B2 (ja) 2020-06-04 2024-05-28 株式会社豊田中央研究所 状態推定装置、及び状態推定方法、状態推定プログラム
WO2023189022A1 (ja) * 2022-03-29 2023-10-05 三菱ケミカルエンジニアリング株式会社 診断装置、診断方法、及び診断プログラム

Also Published As

Publication number Publication date
JPWO2020090770A1 (ja) 2021-09-24
EP3876056A4 (en) 2022-08-03
CN112955839A (zh) 2021-06-11
EP3876056A1 (en) 2021-09-08
US20210397175A1 (en) 2021-12-23
US11669080B2 (en) 2023-06-06
JP7340265B2 (ja) 2023-09-07

Similar Documents

Publication Publication Date Title
WO2020090770A1 (ja) 異常検出装置、異常検出方法、およびプログラム
US10977568B2 (en) Information processing apparatus, diagnosis method, and program
US8630962B2 (en) Error detection method and its system for early detection of errors in a planar or facilities
CN107111309B (zh) 利用监督式学习方法的燃气涡轮机故障预测
US20120290879A1 (en) Method and device for monitoring the state of a facility
Lundgren et al. Data-driven fault diagnosis analysis and open-set classification of time-series data
WO2020090767A1 (ja) 異常診断装置、異常診断方法、及びプログラム
US20200143292A1 (en) Signature enhancement for deviation measurement-based classification of a detected anomaly in an industrial asset
WO2019043600A1 (en) ESTIMATOR OF REMAINING USEFUL LIFE ESTIMATOR
JP6714498B2 (ja) 設備診断装置及び設備診断方法
CN111964909A (zh) 滚动轴承运行状态检测方法、故障诊断方法及系统
JP2012018066A (ja) 異常検査装置
CN111678699B (zh) 一种面向滚动轴承早期故障监测与诊断方法及系统
CN115758200A (zh) 一种基于相似度度量的振动信号故障识别方法及系统
Dervilis et al. Robust methods for outlier detection and regression for SHM applications
US20080219336A1 (en) System and method for fault detection and localization in time series and spatial data
US11663815B2 (en) System and method for inspection of heat recovery steam generator
CN110057588B (zh) 基于奇异值与图论特征融合的轴承早期故障检测与诊断方法及系统
Popescu et al. Change detection in vibration analysis—A review of problems and solutions
JP7433648B2 (ja) 異常検出方法
Baggeröhr et al. Novel bearing fault detection using generative adversarial networks
Overbey et al. Damage assessment using generalized state-space correlation features
Oddan Multivariate statistical condition monitoring
Wang et al. Fault diagnosis of rolling bearings based on undirected weighted graph
JP2010191564A (ja) 特性解析方法および装置、特性分類方法および装置、上記特性解析方法または特性分類方法をコンピュータに実行させるためのプログラム、上記プログラムを記録したコンピュータ読み取り可能な記録媒体

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19880340

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2020553912

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

ENP Entry into the national phase

Ref document number: 2019880340

Country of ref document: EP

Effective date: 20210531