WO2014115615A1 - 異常診断方法およびその装置 - Google Patents

異常診断方法およびその装置 Download PDF

Info

Publication number
WO2014115615A1
WO2014115615A1 PCT/JP2014/050547 JP2014050547W WO2014115615A1 WO 2014115615 A1 WO2014115615 A1 WO 2014115615A1 JP 2014050547 W JP2014050547 W JP 2014050547W WO 2014115615 A1 WO2014115615 A1 WO 2014115615A1
Authority
WO
WIPO (PCT)
Prior art keywords
feature
abnormality
sensor
vector
calculated
Prior art date
Application number
PCT/JP2014/050547
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 US14/761,748 priority Critical patent/US9779495B2/en
Publication of WO2014115615A1 publication Critical patent/WO2014115615A1/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/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/0227—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions
    • G05B23/0235—Qualitative history assessment, whereby the type of data acted upon, e.g. waveforms, images or patterns, is not relevant, e.g. rule based assessment; if-then decisions based on a comparison with predetermined threshold or range, e.g. "classical methods", carried out during normal operation; threshold adaptation or choice; when or how to compare with the threshold
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F01—MACHINES OR ENGINES IN GENERAL; ENGINE PLANTS IN GENERAL; STEAM ENGINES
    • F01D—NON-POSITIVE DISPLACEMENT MACHINES OR ENGINES, e.g. STEAM TURBINES
    • F01D21/00—Shutting-down of machines or engines, e.g. in emergency; Regulating, controlling, or safety means not otherwise provided for
    • F01D21/003—Arrangements for testing or measuring
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F01—MACHINES OR ENGINES IN GENERAL; ENGINE PLANTS IN GENERAL; STEAM ENGINES
    • F01D—NON-POSITIVE DISPLACEMENT MACHINES OR ENGINES, e.g. STEAM TURBINES
    • F01D21/00—Shutting-down of machines or engines, e.g. in emergency; Regulating, controlling, or safety means not otherwise provided for
    • F01D21/14—Shutting-down of machines or engines, e.g. in emergency; Regulating, controlling, or safety means not otherwise provided for responsive to other specific conditions
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F01—MACHINES OR ENGINES IN GENERAL; ENGINE PLANTS IN GENERAL; STEAM ENGINES
    • F01D—NON-POSITIVE DISPLACEMENT MACHINES OR ENGINES, e.g. STEAM TURBINES
    • F01D25/00—Component parts, details, or accessories, not provided for in, or of interest apart from, other groups
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M15/00—Testing of engines
    • G01M15/14—Testing gas-turbine engines or jet-propulsion engines
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2413—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on distances to training or reference patterns
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2413—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on distances to training or reference patterns
    • G06F18/24133—Distances to prototypes
    • G06F18/24137—Distances to cluster centroïds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/20—Drawing from basic elements, e.g. lines or circles
    • G06T11/206—Drawing of charts or graphs
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0004—Industrial image inspection
    • G06T7/001—Industrial image inspection using an image reference approach
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/761—Proximity, similarity or dissimilarity measures
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2260/00—Function
    • F05D2260/80—Diagnostics
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2270/00—Control
    • F05D2270/70—Type of control algorithm
    • F05D2270/708—Type of control algorithm with comparison tables
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T50/00—Aeronautics or air transport
    • Y02T50/60—Efficient propulsion technologies, e.g. for aircraft

Definitions

  • the present invention relates to an abnormality diagnosis method and apparatus for detecting an abnormality at an early stage based on multi-dimensional time-series data output from a plant or equipment and identifying a sensor related to the detected abnormality.
  • Electric power companies use waste heat from gas turbines to supply hot water for district heating and supply high-pressure steam and low-pressure steam to factories.
  • Petrochemical companies operate gas turbines and other power sources. As described above, in various plants and facilities using a gas turbine or the like, abnormality detection for detecting a malfunction of the facility or its sign is extremely important for minimizing damage to society.
  • Patent Document 1 the similarity between observation data and past learning data is calculated by an original method and estimated by linear combination of data weighted according to the similarity. A method is disclosed in which a value is calculated and an abnormality is detected based on the degree of deviation between the estimated value and the observed data.
  • Patent Document 2 discloses an abnormality detection method for detecting the presence or absence of an abnormality based on an abnormality measure calculated by comparison with a model created from past normal data. It is disclosed to create a model by a local subspace method.
  • Patent Document 3 detects an abnormality based on a T 2 statistic and a Q statistic calculated from a plurality of measurement variables, A process state monitoring device that enumerates abnormal factor candidates based on the contribution amount of each measurement variable to these statistics after detection is disclosed.
  • the sensor exceeding the monitoring standard is a sensor related to the abnormality.
  • an abnormality as a combination of two or more types of sensor signals is detected in the first place. There is a problem that you can not.
  • Patent Document 1 and Patent Document 2 when observation data that is not in learning is observed by providing normal data as learning data, these can be detected as abnormal. Since a plurality of sensor signals are calculated together, an abnormality as a combination of sensor signals can also be detected. However, there is no description on a method for correctly diagnosing a sensor related to an abnormality.
  • an abnormality as a combination of two or more types of sensor signals can be detected by using T 2 statistics and Q statistics of a plurality of measurement items as an abnormality index. It is possible to diagnose a sensor related to an abnormality based on the contribution amount of each measurement item to.
  • the abnormality detection as an indicator of T 2 statistic and Q statistic when applied to equipment and plant switches of switches and operating conditions of stopping the operation is frequently and irregularly, becomes many false alarms, It is difficult to apply as it is.
  • the local subspace method described in Patent Document 2 can suppress false alarms even in the above-described case by collecting learning data exhaustively.
  • the local subspace method uses a projection distance onto an affine subspace created using learning data (x1, x2, x3) close to observation data (q) as an index of abnormality. That is, the magnitude of each element of the vector from b to q is the contribution amount with reference to the foot (b) of the perpendicular from the observation data (q) to the affine subspace.
  • the reference (b) follows this, and the contribution amount of other sensors not related to the abnormality may increase.
  • FIGS. 3A and 3B A simplified example is shown in FIGS. 3A and 3B.
  • FIG. 3A normal data is distributed as shown in the abc three-dimensional sensor signal, and when the observation data changes from q ⁇ q ′ ⁇ q, the reference data changes from b ⁇ b ′ ⁇ b. To do. Since it is as shown in FIG. 3B when represented in a time series graph, the contribution amount of a and b is large, but since c is considered to be related to an abnormality in practice, it is erroneously diagnosed. Become. This applies not only to the local subspace method but also to all case-based methods whose criteria change according to observation data, such as a cluster-based method such as ART2 and a method called similarity-based modeling described in Patent Document 1. .
  • an object of the present invention is to provide an abnormality diagnosis method and apparatus for solving the above-described problems, detecting an abnormality with high sensitivity based on a multidimensional time series sensor signal, and correctly diagnosing a sensor related to the detected abnormality. There is.
  • the present invention provides a method and apparatus for diagnosing abnormality of equipment or apparatus based on a plurality of sensor signals output from a plurality of sensors attached to the equipment or apparatus.
  • a multidimensional feature vector is obtained from the feature values of a plurality of sensor signals in the learning period, and a multidimensional feature vector is extracted at each time from the feature values of the plurality of sensor signals at the time of diagnosis.
  • a reference feature vector for each time is calculated based on a set of multidimensional feature vectors and a multidimensional feature vector for each time extracted at the time of diagnosis, and the extracted multidimensional feature vector for each time and the calculated reference feature for each time Calculate the anomaly measure based on the difference from the vector, detect the anomaly by comparing the calculated anomaly measure with a predetermined threshold, and the time when the anomaly was detected
  • a sensor for abnormality associated detected based on the 2-dimensional distribution density of feature values definitive plurality of sensor signals and to identify from the plurality of sensors.
  • a process of creating and storing learning data based on a plurality of sensor signals output from a plurality of sensors attached to the facility or apparatus, and the facility or apparatus In the method and apparatus for diagnosing abnormality of equipment, including the step of diagnosing abnormality of a plurality of sensor signals newly output from a plurality of sensors attached to the sensor, the step of creating and storing learning data includes a plurality of steps Extracting a feature vector from the sensor signal, storing the extracted feature vector as learning data, and a predetermined number of learning data stored according to the feature vector for each of the feature vectors stored as learning data Selecting a feature vector, creating a reference vector for learning using a predetermined number of selected feature vectors, and creating learning A step of calculating an abnormality measure for each of feature vectors stored as learning data based on a reference vector of the method, a step of calculating a threshold value based on the calculated abnormality measure, and a two-dimensional feature distribution of all combinations from the
  • a step of creating and storing learning data based on a plurality of sensor signals output from a plurality of sensors attached to the facility or apparatus, and the facility includes A step of performing mode division for each operating state based on an event signal output from the facility or apparatus, a step of extracting a feature vector from a plurality of sensor signals, a step of storing the extracted feature vector as learning data, For each feature vector stored as learning data, a process for selecting a predetermined number of feature vectors from the learning data according to the feature vector And a step of creating a learning reference vector using the selected predetermined number of feature vectors, a step of calculating an abnormality measure of the feature vector accumulated as learning data based on the created learning reference vector, and
  • the steps to be performed are the step of extracting feature vectors from a plurality of sensor signals output from a plurality of sensors attached to the equipment or the apparatus to obtain observation vectors, and the learning data accumulated according to the extracted observation vectors.
  • the process of creating the network the process of calculating the abnormality measure of the observation vector based on the created reference vector for abnormality diagnosis, and the observation based on the calculated abnormality measure and the threshold value calculated for each divided mode and divided mode
  • a step of identifying a sensor related to the device A step of determining whether a vector is abnormal or normal, and an abnormality based on a two-dimensional feature distribution density calculated for each mode and mode at the time when a sensor signal corresponding to the observation vector determined to be abnormal is detected.
  • the abnormality detection of the data at each time is performed based on the distance from the learning data and the reference data calculated from the data at that time, it is possible to detect the abnormality with high sensitivity.
  • the sensor having a large contribution amount is not used as it is as the related sensor, but the sensor related to the abnormality is specified based on the two-dimensional distribution density of the sensor signal. can do.
  • wrist which shows the example of a sensor signal. It is a flowchart which shows the flow of a process from a sensor signal input to threshold value setting.
  • the present invention relates to a method and apparatus for diagnosing an abnormality of a facility or apparatus based on a plurality of sensor signals output from a plurality of sensors attached to the facility or apparatus, and a multidimensional feature vector for each time from the sensor signal.
  • a reference vector at each time based on a set of feature vectors for a learning period specified in advance and a feature vector at each time, and an anomaly measure based on the difference between the feature vector at each time and the reference feature vector
  • the abnormality is detected by comparing the abnormality measure with a predetermined threshold value, and the sensor related to the abnormality can be identified based on the two-dimensional distribution density of the feature value at the time when the abnormality is detected.
  • FIG. 1 shows a configuration example of a system for realizing the abnormality diagnosis method of the present invention.
  • This system includes a sensor signal storage unit 103 that stores the sensor signal 102 output from the equipment 101, a feature vector extraction unit 104 that extracts a feature vector based on the sensor signal 102, and a feature vector of a learning period specified in advance.
  • a reference vector calculation unit 105 that calculates a reference feature vector at each time based on the set and a feature vector at each time; an abnormal measure calculation unit that calculates an abnormal measure based on a difference between the feature vector at each time and the reference feature vector;
  • a threshold value calculation unit 107 that calculates a threshold value based on an abnormality measure for a specified learning period, an abnormality detection unit 108 that detects an abnormality by comparing the abnormality measure with the calculated threshold value, and a sensor signal for the learning period
  • a distribution density calculation unit 109 for calculating a two-dimensional distribution density of the sensor, and identifies a sensor related to the abnormality at the time when the abnormality is detected. Configured with an associated sensor identification unit 110.
  • the operation of this system has two phases: “learning” that generates and stores learning data using accumulated data, and “abnormality diagnosis” that detects abnormalities based on input signals and identifies related sensors. is there. Basically, the former is offline processing and the latter is online processing. However, the latter can be processed offline. In the following explanation, they are distinguished by the terms of learning and abnormality diagnosis.
  • the equipment 101 subject to state monitoring is equipment or a plant such as a gas turbine or a steam turbine.
  • the facility 101 outputs a sensor signal 102 indicating the state.
  • the sensor signal 102 is stored in the sensor signal storage unit 103.
  • An example in which the sensor signals 102 are listed and displayed in a tabular form is shown in FIG.
  • the sensor signal 102 is a multi-dimensional time series signal acquired at regular intervals, and a table in which the sensor signal 102 is listed is a date / time column 401 and a plurality of sensor values provided in the equipment 101 as shown in FIG. It consists of a data column 402.
  • the number of types of sensors may be several hundreds to thousands, for example, the temperature of cylinders, oil, cooling water, etc., the pressure of oil or cooling water, the rotational speed of the shaft, the room temperature, the operating time, and the like.
  • a control signal for controlling something to a certain value.
  • feature vectors are extracted using data for a specified period of data stored in the sensor signal storage unit 103.
  • the data extracted here is called learning data.
  • An abnormal measure of learning data is calculated by cross-validation, and an abnormality determination threshold value is calculated based on the abnormal measure of all learning data.
  • the distribution density calculation unit 109 calculates a two-dimensional distribution density of feature values.
  • the flow of processing in the feature vector extraction unit 104, the reference vector calculation unit 105, the abnormality measure calculation unit 106, and the threshold value calculation unit 107 will be described with reference to FIG.
  • the sensor signal 102 of the period specified as the learning period is input from the sensor signal storage unit 103 (S501), and after normalization for each sensor signal (S502), the feature vector is extracted. Is performed (S503).
  • the reference vector calculation unit 105 assigns a group number for intersection verification to the extracted feature vector (S504).
  • a group is determined by specifying a period, for example, one day is one group, or is divided into equal numbers of groups specified in advance.
  • the first feature vector is noticed from the extracted feature vectors (S505), and a predetermined number of feature vectors are selected from the feature vectors in a group different from the notice vector in the order closest to the notice vector (S506). ), A reference vector is created using the selected feature vector (S507).
  • the anomaly measure calculation unit 106 calculates an anomaly measure based on the distance from the feature vector of interest to the reference vector (S508). If the abnormal measure calculation for all vectors has been completed (S509), the threshold value calculation unit 107 sets a threshold value based on the abnormal measure for all vectors (S511). If all the abnormal measure calculations have not been completed in step S509, the next feature vector is noticed (S510), and the processes in steps S506 to S509 are repeated.
  • each sensor signal is normalized. For example, using the average and standard deviation of a specified period, conversion is performed so that the average becomes 0 and the variance becomes 1.
  • the average and standard deviation of each sensor signal is stored so that the same conversion can be performed at the time of abnormality diagnosis.
  • conversion is performed so that the maximum is 1 and the minimum is 0 using the maximum and minimum values of the specified period of each sensor signal.
  • a preset upper limit value and lower limit value may be used instead of the maximum value and the minimum value.
  • the maximum value and minimum value or the upper limit value and the lower limit value of each sensor signal are stored so that the same conversion can be performed at the time of abnormality diagnosis.
  • the canonicalization of sensor signals is for simultaneously handling sensor signals of different units and scales.
  • step S503 feature vectors are extracted at each time. It is conceivable to arrange sensor signals in canonical order as they are. However, there is a window of ⁇ 1, ⁇ 2, ... for a certain time, and features of [Window width (3, 5, ...) x number of sensors] It is also possible to extract features representing temporal changes in data by using vectors. Alternatively, discrete wavelet transform (DWT: “Discrete” Wavelet “Transform”) may be applied to be decomposed into frequency components. In step S503, feature selection is performed. As a minimum processing, it is necessary to exclude sensor signals with very small variance and monotonously increasing sensor signals.
  • DWT discrete wavelet transform
  • LSC Local subspace
  • PDM projection distance method
  • the local subspace method is a method of creating a k ⁇ 1 dimensional affine subspace using the k ⁇ neighbor vector of the vector of interest q.
  • the k-neighbor vector is the feature vector selected in step S506.
  • the number to be specified is the value of k.
  • the point b on the affine subspace closest to the vector of interest q may be obtained.
  • the projection distance method is a method of creating a subspace having a unique origin for a selected feature vector, that is, an affine subspace (space of maximum dispersion).
  • the number specified in step S506 may be any number, but if it is too large, both vector selection and subspace calculation take time, so it is preferable to set it to several tens to several hundreds.
  • a method for calculating the affine subspace will be described. First, the average ⁇ of the selected feature vectors and the covariance matrix ⁇ are obtained, and then the eigenvalue problem of ⁇ is solved, and a matrix U in which eigenvectors corresponding to r eigenvalues designated in advance from the larger value are arranged is affine. Let it be an orthonormal basis of the subspace. r is a number smaller than the dimension of the feature vector and smaller than the number of selected data. Alternatively, r may not be a fixed number, and may be a value when the contribution rate accumulated from the larger eigenvalue exceeds a predesignated ratio. The anomaly measure is the projection distance of the vector of interest onto the affine subspace.
  • a local average distance method in which the distance from the vector of interest q to the average vector of k-neighbor vectors is an abnormal measure, a Gaussian process, or the like may be used.
  • the method for setting the threshold value in step S511 will be described.
  • the abnormal measures of all the feature vectors in the learning period are sorted in ascending order, and a value reaching a ratio close to 1 specified in advance is obtained.
  • the threshold value is calculated by processing such as adding an offset or multiplying by a constant. If the offset is 0 and the magnification is 1, this value becomes the threshold value.
  • the calculated threshold value is not shown, it is recorded in association with the learning data.
  • the threshold value is calculated based on the abnormality measure calculated by cross-validation of the learning data, it is possible to suppress false alarms.
  • the processing range of the distribution density calculation is calculated by expanding the range outward from the minimum value and the maximum value (S605). For example, MIN is changed to MIN ⁇ S ⁇ M and MAX is changed to MAX + S ⁇ M.
  • M is a predetermined integer of 1 or more.
  • the bin number (BNO) is calculated from the feature value (F) by the following equation (S606).
  • BNO INT ((F-MIN) ⁇ N / (MAX-MIN))
  • INT (X) represents the integer part of X.
  • step S607 If unprocessed features remain in the above processing (S607), pay attention to the next feature (S608), and perform the processing from step S603 to S606. If all the features have been processed in step S607, the process proceeds to step S609.
  • step S609 focus on the first combination of two features.
  • the two features may be the same.
  • a two-dimensional array for calculating distribution density is secured, and 0 is set to all elements (S610).
  • the size of the array is N + 2M.
  • 1 is added to the elements of the array corresponding to the bin numbers of the two features (S611).
  • a two-dimensional frequency distribution hoverogram
  • This frequency distribution is converted into an image and stored (S612). The conversion method will be described later. If the above process has been performed for all combinations of two features (S613), the process ends (S615). Otherwise, the next combination is noted (S614), and the processes from step S610 to S612 are performed. Do.
  • the size of the two-dimensional array and the minimum and maximum values of each feature calculated in step S605 are recorded.
  • the maximum value of array elements that is, the maximum frequency is obtained.
  • the image size is the same as the array size, and the pixel value of the corresponding coordinate from the value of each element is, for example, [255 ⁇ element value of array / maximum frequency].
  • 255 is the maximum value when the pixel value is represented by 8 bits, and if this value is used, it can be stored in the bitmap format as it is.
  • the pixel value is [255 ⁇ LOG (element value of array + 1) / LOG (maximum frequency + 1)].
  • LOG (X) represents the logarithm of X.
  • the image obtained by the above processing is referred to as a distribution density image because a portion with a high density in the two-dimensional feature space is represented by a high pixel value.
  • the method of creating an image is not limited to the above method. For example, instead of a simple frequency distribution, a Gaussian distribution or another weighted filter may be assigned to one piece of data and superimposed. Alternatively, a maximum value filter of a predetermined size may be applied to the image obtained by the above method, or an average filter or other weighted filter may be applied. In addition, it is not always necessary to save in the image format, and the two-dimensional array may be saved in the text format.
  • an abnormality measure of data for a specified period or newly observed data among the data accumulated in the sensor signal accumulation unit 103 is calculated, and it is determined whether it is normal or abnormal. Identify the relevant sensor for time.
  • FIG. 7 is a diagram for explaining the flow of processing at the time of abnormality diagnosis in the feature vector extraction unit 104, the reference vector calculation unit 105, the abnormality measure calculation unit 106, and the abnormality detection unit 108.
  • a sensor signal is input from the sensor signal storage unit 103 or the facility 101 (S701), and after normalization for each sensor signal (S702), the feature vector is extracted (S703).
  • the same parameters as those used for the canonicalization of the learning data are used in the process shown in step S502 of FIG.
  • the feature vector is extracted by the same method as that in step S503. Therefore, when the feature selection is executed in step S503, the same feature is selected.
  • the feature vector extracted here is referred to as an observation vector in order to distinguish it from learning data.
  • the reference vector creation unit 105 selects a predetermined number of feature vectors from the feature vectors of the learning data in the order closest to the observation vector (S704), and creates a reference vector using these feature vectors. (S705).
  • the anomaly measure calculation unit 106 calculates an anomaly measure based on the distance of the observation vector to the reference vector (S706).
  • the abnormality detection unit 108 compares the threshold value calculated at the time of learning with the abnormality measure, and determines that it is abnormal if the abnormality measure is larger than the threshold value, and normal otherwise (S707).
  • FIG. 8 is a diagram for explaining an embodiment of the flow of processing in the related sensor specifying unit 110. This process is performed for the data determined to be abnormal in step S707. First, the difference between the observed value and the reference value of each feature is calculated from the observed vector and the reference vector (S801), and the feature having the maximum value is searched for (S802). This is referred to as feature A. Next, from the observation vector, the bin number of each feature is calculated using the size of the two-dimensional array used in the distribution density calculation process shown in FIG. 6 and the minimum and maximum values of each feature calculated in step S605 ( S803).
  • This is a feature B, and corresponding points of the observation vector are plotted on the distribution density images of the features A and B (S806).
  • a color image is converted and an unused color is assigned to the coordinates corresponding to the bin numbers of the feature A and the feature B. If the pixel value in S804 is 0 or smaller than a predetermined value, feature A and feature B are set as related sensors (S807).
  • FIG. 9 is a diagram for explaining another embodiment of the process flow in the related sensor identification unit 110. Similar to the processing of FIG. 8, in step S707, the data determined to be abnormal are executed. First, the bin number of each feature is calculated from the observation vector using the size of the two-dimensional array used in the distribution density calculation process shown in FIG. 6 and the minimum and maximum values of each feature calculated in step S605. (S901).
  • a feature that minimizes s i is searched (S904), and this is set as a feature C.
  • the pixel value of the coordinate corresponding to each feature and the bin number of the feature C is read (S905), and the feature that minimizes this is searched (S906).
  • the corresponding points of the observation vector are plotted on the distribution density image of the feature C and the feature D in the same manner as in step S806 (S907). If the pixel value in S905 is 0 or smaller than a predetermined value, feature C and feature D are set as related sensors (S908).
  • the bin number of each feature is calculated in step S901, and the pixel value of the coordinate corresponding to the observation vector is read for all distribution density images in step S902. Look for these as feature E and feature F.
  • the corresponding points of the observation vector are plotted on the distribution density images of the features E and F in the same manner as in step S806, and the pixel value read from the distribution density images of the features E and F is 0 or predetermined. If it is smaller than the above value, the feature E and the feature F are used as related sensors.
  • any one of the related sensor identification methods described above may be performed, or may be performed sequentially until the related sensor can be identified, or all of the methods may be performed to display a distribution density image in which observation vectors are plotted. Then, the user may decide.
  • a method of analyzing in detail when the related sensor cannot be identified by the above method will be described with reference to FIG.
  • a feature vector for a learning period is input (S1001).
  • the bin number of each feature is determined using the size of the two-dimensional array used in the distribution density calculation process shown in FIG. 6 and the minimum and maximum values of each feature calculated in step S605.
  • Calculate (S1002) Similarly, the bin number of each feature is calculated for the observation vector (S1003).
  • the pixel value of the coordinate corresponding to each feature and the bin number of the feature G is read (S1009), and the feature having the smallest value is searched for as the feature H (S1010). .
  • the corresponding points of the observation vector are plotted on the partial distribution density images of the feature G and the feature H (S1011). If the pixel value in S1009 is 0 or smaller than a predetermined value, the feature A, feature G, and feature H are set as related sensors (S1012). In these processes, the above-described feature C may be used instead of the feature A.
  • the above method is a method of specifying a related sensor at each time when an abnormality is detected.
  • the method is applied when abnormality is continuously detected, it becomes a redundant process. Therefore, it is conceivable that information is stored in a buffer and processing is performed collectively while abnormalities are continuously detected.
  • the target for sensor identification may be narrowed down by the duration of abnormality or the cumulative abnormality degree.
  • a section in which abnormality is continuously detected is specified based on the result of abnormality determination in step S707 in the abnormality detection unit 108 (S1101). Meanwhile, the observation vector and the reference vector are stored in the buffer. Note that “continuous” does not necessarily have a strict meaning, and may include interruption of a predetermined length.
  • the sum of the absolute values of the difference between the observation value and the reference value is calculated for each feature (S1102), and the feature having the maximum value is searched (S1103). This is referred to as feature A.
  • the size of the two-dimensional array used in the distribution density calculation process shown in FIG. 6 and the minimum and maximum values of each feature calculated in step S605 are used.
  • the bin number of the feature is calculated (S1104).
  • a distribution density image whose one feature is the feature A is read, and pixel values of coordinates corresponding to the bin numbers of the features and the feature A are read for all the observation vectors included in the section (S1105).
  • Non-zero pixels are counted for each feature for all observation vectors included in the section (S1106).
  • the sum of pixel values is calculated for each feature for all observation vectors included in the section (S1107).
  • the value calculated for feature i is sum (i).
  • a feature having the smallest cnt (i) is searched for. If there are a plurality of features, a feature having the smallest sum (i) is searched for (S1108). This is a feature B, and corresponding points of observation vectors included in the abnormal section are plotted on the distribution density images of the features A and B (S1109). If the number of non-zero pixels in S1106 is smaller than the number of observation vectors included in the section, that is, if the number of 0 pixels is 1 or more, feature A and feature B are set as related sensors (S1110).
  • a feature with high correlation of learning data may be searched for as feature B.
  • the level of correlation is determined based on, for example, the correlation value between the feature A and the other feature.
  • a correlation density image having a small number of non-zero pixels is considered to have a high correlation.
  • FIG. 11 The process flow of FIG. 11 described above is a process in which the process shown in FIG. 8 is extended to a plurality of observation vectors, but the process shown in FIG. 9 or 10 is extended to a plurality of observation vectors. Or they may be used in combination.
  • the horizontal axis and the vertical axis represent the relative values of the respective features, and dark portions are portions where the learning data density is high.
  • An asterisk represents observation data included in the abnormal section. Actually, it may be displayed with a color.
  • FIG. 12A shows the distribution density of features a and b. From the distribution of the learning data, it can be seen that there is a steady state where both the features a and b are low and the features a and b are both high, and the feature b rises behind the feature a during the state transition. The relationship of time lag is stable. It is detected as an abnormality that the rise of the feature b occurs earlier than the timing of the rise of the feature a. Alternatively, the shift of the descending timing may be detected, but it can be confirmed from the time series graph of the sensor signal.
  • FIG. 12B represents the distribution density of the features c and d. There is a state where both the features c and d are low, a feature c is high and a feature d is low, a feature c and d are both high, and an intermediate state between them, but the feature d is high and the feature c is not low. Such a state is detected as abnormal.
  • FIG. 12C shows the value of the feature e on both the horizontal axis and the vertical axis. Although it is not substantially two-dimensional, it is possible to simplify display and sensor specific processing by representing the images in the same format.
  • the learning data is distributed on a straight line. It is detected as abnormal that the value of the feature e is higher than any value of the learning data.
  • the two-dimensional feature distribution density of all combinations is calculated in advance using learning data, and the sensor related to the abnormality is specified based on the distribution density of the corresponding points of the observation vector when the abnormality is detected. Therefore, it is possible to specify the correct sensor even in the situation shown in FIGS. 3A and 3B.
  • the tendency of the distribution of learning data can be known. For example, as shown in FIG. 12D, when the low-density part of the learning data exists in a wide range, even if the abnormality indicated by the asterisk 1201 is detected, it must be doubted whether it is really abnormal. In that case, the learning period is increased so as to increase the density, or the combination is not processed together.
  • the feature g is preferably excluded from the processing target.
  • FIG. 1 An example of a GUI for setting the learning period and processing parameters is shown in FIG. In the following description, this setting is simply referred to as recipe setting.
  • the past sensor signal 102 is assumed to be stored in the database in association with the equipment ID and time.
  • the recipe setting screen 1301 the target device, learning period, use sensor, reference calculation parameter, and threshold setting parameter are input.
  • the facility ID input window 1302 the ID of the target facility is input.
  • a list of device IDs of data stored in the database is displayed by pressing the equipment list display button 1303, it is selected and input from the list.
  • the learning period input window 1304 the start date and end date of the period for which learning data is to be extracted are input.
  • the sensor selection window 1305 a sensor to be used is input. By clicking the list display button 1306, a sensor list 1307 is displayed. It is possible to select multiple items from the list.
  • the reference calculation parameter input window 1308 parameters specified in normal model creation are input. The figure shows an example in which a local subspace is adopted as a normal model, and the number of neighboring vectors used for model creation and regularization parameters are input.
  • the regularization parameter is a small number added to the diagonal component in order to prevent the inverse matrix of the correlation matrix C from being obtained in Equation 2.
  • a radio button is used to select how to determine a cross-validation group in the processing shown in FIG. Choose whether one day is divided into one group or a specified number of groups. In the latter case, the number of groups is entered.
  • a value to be applied to the cumulative histogram in the threshold value setting process in step S511 is input.
  • the recipe name input window 1310 a unique name associated with the input information is input. When all the information is input, the test target period is input to the test period input window 1311 and the test of the recipe is performed by pressing the test button 1312.
  • the serial number of the test executed with the same recipe name is assigned.
  • the feature vector extraction unit 104 extracts feature vectors from the sensor signal 102 in the designated learning period.
  • the average and standard deviation are obtained using all data in the designated learning period.
  • the average and standard deviation values are also stored in association with the recipe name and test number for each sensor.
  • the reference vector calculation unit 105 generates a reference vector at each time by cross validation, and the abnormal measure calculation unit 106 calculates an abnormal measure corresponding to step S508. Further, in step S511, the threshold value calculation unit 107 calculates an abnormality determination threshold value and stores it in association with the recipe name and test number. In the distribution density calculation unit 109, a distribution density image of a two-dimensional feature is created and stored by the process shown in FIG. In addition, the processing range calculated in step S605 is stored for all features. In addition, device ID information, sensor usage information, learning period, parameters used for feature vector extraction, and reference creation parameters are stored in association with recipe names and test numbers.
  • the abnormality detection process shown in FIG. 7 is performed using the sensor signal 102 of the designated test period as an input.
  • the related sensor specifying process is performed according to the procedure shown in any of FIGS.
  • a serial number is added to each detected abnormality.
  • a serial number, a section start time, a section end time, and a distribution density image name in which corresponding points of a plurality of related sensor names and observation vectors are plotted are recorded.
  • FIGS. 14A, 14B, and 14C Examples of GUI for showing the test result to the user are shown in FIGS. 14A, 14B, and 14C. By selecting a tab displayed at the top of each screen, the entire result display screen 1401, the enlarged result display screen 1402, and the abnormality diagnosis result display screen 1403 can be switched.
  • FIG. 14A shows the entire result display screen 1401.
  • the entire result display screen 1401 displays a time series graph of the abnormality measure, threshold value, determination result, and sensor signal for the entire specified period.
  • the period display window 1404 displays the designated learning period and test period.
  • the anomaly measure display window 1405 displays the anomaly measure, threshold value, and determination result for the designated learning period and test period.
  • the sensor signal display window 1406 displays the output value of the specified sensor for the specified period.
  • the designation of the sensor is performed by inputting to the sensor name selection window 1407. However, the head sensor is selected before the user designates it.
  • a cursor 1408 represents the starting point for enlarged display and can be moved by a mouse operation.
  • the display day designation window 1409 although not used on this screen, the number of days from the start point to the end point of the enlarged display on the result enlarged display screen 1402 is displayed. You can also enter on this screen.
  • the date display window 1410 displays the date at the cursor position. When the end button 1411 is pressed, the entire result display screen 1401, the enlarged result display screen 1402, and the abnormality diagnosis result display screen 1403 are erased and the processing ends.
  • FIG. 14B shows a result enlarged display screen 1402.
  • the result enlarged display screen 1402 displays a time series graph of the abnormality measure, threshold value, determination result, and sensor signal of the specified number of days starting from the date indicated by the cursor 1408 in the entire result display screen 1401. .
  • the period display window 1404 displays the same information as the entire result display screen 1401.
  • the abnormality measure display window 1405 and the sensor signal display window 1406 the same information as that on the entire result display screen 1401 is displayed in an enlarged manner.
  • the display day designation window 1409 the number of days from the start point to the end point of the enlarged display is designated.
  • the starting date of the enlarged display is displayed. It is also possible to change the starting point of display with the scroll bar 1412, and this change is reflected in the position of the cursor 1408 and the display of the date display window 1410.
  • the entire length of the scroll bar display area 1413 corresponds to the entire period displayed on the entire result display screen 1401.
  • the length of the scroll bar 1412 corresponds to the number of days specified in the display day specification window 1409, and the left end of the scroll bar 1412 corresponds to the starting point of the enlarged display.
  • the abnormal section number 1414 is displayed in a balloon at the corresponding position in the abnormality measure graph. The process ends when the end button 1411 is pressed.
  • FIG. 14C shows an example of the abnormality diagnosis result display screen 1403.
  • a two-dimensional feature distribution density of the sensor related to the abnormality and a time series graph are displayed together.
  • the abnormality with the youngest number among the abnormalities displayed on the result enlarged display screen 1402 is set as a display target.
  • the date display window 1415 an abnormal date to be displayed is displayed.
  • the abnormal section number display window 1416 displays an abnormal section number to be displayed.
  • an abnormal section number can be designated by numerical input.
  • the time display window 1417 displays the section start time and section end time of the abnormality to be displayed.
  • the distribution density display window 1418 for example, the image stored in step S1109 of FIG. 11 is displayed.
  • This image represents the two-dimensional frequency distribution of the two sensors identified in step S1110.
  • the horizontal axis is the normalized value of “SensorX”, and the vertical axis is the normalized value of “SensorY”.
  • 1419 represents the distribution of the learning period, and the darker the display, the higher the frequency.
  • 1420 is data corresponding to the observation vector of the abnormal section to be displayed.
  • the check box can be unchecked so that it is not displayed, and the distribution density of data during the learning period of the feature value can be confirmed.
  • the processing ranges calculated in step S605 of the two sensors are displayed. These values are the “SensorX” value corresponding to the left edge of the distribution density image, the “SensorX” value corresponding to the right edge, the “SensorY” value corresponding to the lower edge, and the “SensorY” value corresponding to the upper edge. is there.
  • the anomaly measure display window 1423 displays an anomaly measure, a threshold value, and a determination result for a period including an anomaly section to be displayed.
  • the display period is determined in advance so as to be sufficiently enlarged so that a change in signal can be observed.
  • the output value of the same period as the abnormality measure display window 1423 of one of the two sensors specified in step S1110 is displayed in a graph.
  • the second related sensor signal display window 1425 the other sensor signal output value is displayed in a graph.
  • the reference values can be displayed in different colors in the first related sensor signal display window 1424 and the second related sensor signal display window 1425 (not shown). it can.
  • Buttons 1427 and 1428 are buttons for switching the abnormality diagnosis result to be displayed.
  • button 1427 “Next” is clicked, the abnormality of the next lower number in the abnormality diagnosis result display screen 1403 than the abnormality data currently displayed is displayed. Is displayed.
  • button 1428 “Previous” is clicked, abnormality data having a number one lower than that of the abnormality data currently displayed is displayed on the abnormality diagnosis result display screen 1403.
  • the distribution density display window 1418 may display the images of the number of tables stored in step S1109 side by side.
  • time series graphs of the two sensors are displayed separately, but they may be displayed in different colors in one window.
  • the process ends when the end button 1411 is pressed.
  • the registration button 1315 When the registration button 1315 is pressed, the information stored in association with the recipe name and the test number being displayed in the test number display window 1313 is registered in association with the recipe name, and the process ends. If a cancel button 1316 is pressed, the process ends without saving anything.
  • test result list display screen 1501 shown in FIG. 15 is displayed.
  • the test result list 1502 displays recipe information such as the learning period, test period, selected sensor number, reference creation parameter, and threshold setting parameter for all tests, and test result information such as the threshold and the number of abnormal sections. To do. There is a selection check button at the left end of the list, and only one of them can be selected.
  • the detail display button 1503 is pressed, information stored in association with the recipe name and test number is loaded, and the entire result display screen 1401 is displayed.
  • the result enlarged display screen 1402 or the abnormality diagnosis result display screen 1403 can be displayed by switching the tabs.
  • the display returns to the test result list display screen 1501 by pressing the end button 1411.
  • the registration button 1504 By pressing the registration button 1504, the information stored in association with the selected test number is registered in association with the recipe name, and the display of the test result list display screen 1501 and the display of the recipe setting screen 1301 are terminated.
  • the return button 1505 is pressed, the display returns to the recipe setting screen 1301 without registering the recipe.
  • Registered recipes are managed with an active or inactive label, and for newly observed data, information on active recipes with matching device IDs is used with reference to FIGS. Processing from any of the described feature vector extraction to abnormality detection and related sensor identification is performed, and the result is stored in association with the recipe name.
  • FIG. 16 shows an example of a GUI for showing the result of the above abnormality detection diagnosis processing to the user.
  • FIG. 16 is an example of a GUI for designating a display target.
  • the display target specification screen 1601 specifies the display target equipment, recipe, and period.
  • a device ID is selected by a device ID selection window 1602.
  • a recipe to be displayed is selected from a list of recipes targeted for the device ID by means of a recipe name selection window 1603.
  • the data recording period display unit 1604 displays the start date and end date of the period that is processed using the input recipe and remains recorded.
  • the result display period designation window 1605 the start date and end date of the period for which the result is to be displayed are input.
  • the display sensor designation window 1606 the name of the sensor to be displayed is input.
  • the display period designated in the result display period designation window 1605 is displayed as shown in FIG.
  • the window corresponding to the abnormal measure display window 1405 displays the abnormal measure, the threshold value, and the determination result of the designated display period.
  • the output value of the sensor designated by the display sensor designation window 1606 for the designated period is displayed.
  • the display target sensor can be changed by inputting to a window corresponding to the sensor name selection window 1408.
  • the above embodiment processes learning data setting offline, abnormality diagnosis processing in real time, and result display offline, but the result display can also be performed in real time.
  • the length of the display period, the recipe to be displayed, and the information to be displayed may be determined in advance, and the latest information may be displayed at regular intervals.
  • the sensor related to the detected abnormality can be correctly diagnosed. Became possible.
  • the two-dimensional feature distribution density obtained as a result of the inspection is displayed as an image, it becomes easy to distinguish the sensor related to the abnormality. In addition, the trend of the distribution of learning data can be confirmed on the screen.
  • FIG. 18 shows a second embodiment of the system for realizing the abnormality diagnosis method of the present invention.
  • the configuration in this embodiment is obtained by adding a mode dividing unit 1802 to the configuration shown in FIG. 1 described in the first embodiment.
  • the same components as those in the first embodiment are denoted by the same reference numerals.
  • a mode division unit 1802 receives an event signal 1801 from the facility, and divides the event signal 1801 into modes representing the operation state of the facility based on the event signal 1801.
  • the result of the mode division is input to the threshold value calculation unit 107 ′, and the processing flow described in FIG. 5 is executed for each divided mode, and the threshold setting in step S511 is performed for each mode.
  • the processing flow described in FIG. 7 is also executed for each divided mode, and the abnormality determination in step S707 in the abnormality detection unit 108 ′ is performed using the threshold value of the corresponding mode. Further, the distribution density calculation unit 109 ′ also inputs the result of mode division, creates a distribution density image for each mode, and performs sensor specification processing in the related sensor specification unit 110 ′ for each mode.
  • FIG. 19A An example of the event signal is shown in FIG. 19A.
  • a list 1920 shown in FIG. 19A is a signal representing the operation / failure / warning of the equipment that is output irregularly, and includes a time string and a character string or code representing the operation / failure / warning.
  • this event signal is input to the mode dividing unit 1802 (S1901), and a start sequence and a stop sequence are cut out by searching for a predetermined character string or code (S1902). Based on the result, “steady OFF” mode 1911 from the end time of the stop sequence to the start time of the start sequence, “start” mode 1912 in the start sequence, from the end time of the start sequence to the start time of the stop sequence The operation is divided into four operating states of “steady ON” mode 1913 and “stop” mode 1914 in the stop sequence (S1903). An example is shown in FIG. 19C.
  • a start event and an end event of the sequence are designated in advance, and the event signal is cut out while scanning from the beginning to the end in the following manner.
  • search for a start event When it is found, the sequence starts.
  • search for an end event If found, end the sequence.
  • the end event is a specified end event, a failure, a warning, or a specified start event.
  • the event signal it is possible to accurately divide various operating states, and by setting a threshold value for each mode, the transition period between the “start” mode 1912 and the “stop” mode 1914 Even when it is necessary to lower the sensitivity due to lack of learning data, the “steady OFF” mode 1911 and the “steady ON” mode 1913 enable highly sensitive abnormality detection.
  • the present embodiment it is possible to set a threshold value according to various operating states of the facility, enable highly sensitive abnormality detection, and correctly diagnose a sensor related to the detected abnormality. Became possible. Further, by calculating the distribution density of the learning data for each mode, data in different states is not included, so diagnosis is facilitated and understanding of the distribution density image is facilitated.

Abstract

　多次元時系列センサ信号に基づく異常検知において、対策、調査など次のアクションを決めるために、どのセンサが異常に関連しているのかを特定する必要があるが、基準値と観測値の差が大きいものを関連センサとする従来の方法では、誤る場合があったために、本発明では、センサ信号から時刻毎に多次元特徴ベクトルを抽出し、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準特徴ベクトルを算出し、各時刻の特徴ベクトルと基準特徴ベクトルの差に基づき異常測度を算出し、異常測度と所定のしきい値との比較により異常を検出し、異常が検出された時刻について、特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定するようにした。

Description

異常診断方法およびその装置
　本発明は、プラントや設備などの出力する多次元時系列データをもとに異常を早期に検知し、検知された異常に関連するセンサを特定する異常診断方法およびその装置に関する。
　電力会社では、ガスタービンの廃熱などを利用して地域暖房用温水を供給したり、工場向けに高圧蒸気や低圧蒸気を供給したりしている。石油化学会社では、ガスタービンなどを電源設備として運転している。このようにガスタービンなどを用いた各種プラントや設備において、設備の不具合あるいはその兆候を検知する異常検知は、社会へのダメージを最小限に抑えるためにも極めて重要である。
　ガスタービンや蒸気タービンのみならず、水力発電所での水車、原子力発電所の原子炉、風力発電所の風車、航空機や重機のエンジン、鉄道車両や軌道、エスカレータ、エレベータ、機器・部品レベルでも、搭載電池の劣化・寿命など、上記のような予防保全を必要とする設備は枚挙に暇がない。
　このため、対象設備やプラントに複数のセンサを取り付け、センサ毎の監視基準に従って正常か異常かを判断することが行われている。米国特許第６，９５２，６６２号（特許文献１）には、観測データと過去の学習データとの類似度を独自の方法で計算し、類似度に応じて重み付けされたデータの線形結合により推定値を算出して、推定値と観測データのはずれ度合をもとに異常を検出する方法が開示されている。また、特開２０１１－７０６３５号公報（特許文献２）には、過去の正常データから作成されたモデルとの比較によって算出される異常測度に基づいて異常の有無を検知する異常検知方法において、正常モデルを局所部分空間法によって作成することが開示されている。
　しかし、異常を検知するのみでは対策や調査などの次のアクションが決められないため、どのセンサが異常に関連しているのかを診断したいというニーズがある。
　このようなニーズに対応する技術として、特開２０１２－１３８０４４号公報（特許文献３）には、複数の測定変数から算出されるＴ２統計量およびＱ統計量に基づいて異常を検出し、異常検出後にこれらの統計量に対する各測定変数の寄与量に基づいて異常要因候補を列挙するプロセス状態監視装置が開示されている。
米国特許第６，９５２，６６２号明細書 特開２０１１－７０６３５号公報 特開２０１２－１３８０４４号公報
　センサ毎の監視基準に従って正常か異常かを判断する方法では、監視基準を超えたセンサがすなわち異常に関連するセンサであるが、この方法では２種以上のセンサ信号の組合せとしての異常はそもそも検知できないという問題がある。
　特許文献１や特許文献２に記載の方法によれば、学習データとして正常時のデータを与えることにより、学習にない観測データが観察されると、これらを異常として検出することができる。複数のセンサ信号をまとめて演算するため、センサ信号の組合せとしての異常も検出できる。しかし、異常に関連するセンサを正しく診断する方法については記載がない。
　特許文献３に記載の方法によれば、複数の測定項目のＴ２統計量やＱ統計量を異常の指標とすることにより、２種以上のセンサ信号の組合せとしての異常も検知でき、統計量に対する各測定項目の寄与量に基づいて異常に関連するセンサを診断することが可能である。しかし、Ｔ２統計量やＱ統計量を指標とした異常検出は、運転と停止の切り替りや運転条件の切り替りが頻繁かつ不規則にある設備やプラントに適用すると、誤報が多くなってしまい、そのままでは適用困難である。特許文献２に記載の局所部分空間法は、学習データを網羅的に集めることにより、上記のようなケースでも誤報を抑制することが可能である。
　しかし、特許文献１に記載の異常検知手法に、特許文献２に記載の寄与量に基づいて異常に関連するセンサを診断する方法を適用すると、誤る場合がある。この理由について簡単に説明する。局所部分空間法は、図２に示す通り、観測データ(q)に近い学習データ(x1，x2，x3)を用いて作られるアフィン部分空間への投影距離を異常の指標としている。つまり、観測データ(q)からアフィン部分空間への垂線の足(b)を基準として、bからqのベクトルの各要素の大きさが寄与量となる。あるセンサ信号に異常がある場合、これに基準(b)が追随して異常とは関係ない他のセンサの寄与量が大きくなることが起こりうる。
　単純化した例を図３Ａ及び図３Ｂに示す。図３Ａにおいて、abcの３次元のセンサ信号において正常なデータが図示するように分布しており、観測データがq→q’→qと変化した場合、基準データはb→b’→bと変化する。時系列のグラフで表すと図３Ｂに示す通りとなるため、aとbの寄与量が大きくなるが、実際にはcが異常に関連していると考えられるため、誤って診断されることになる。このことは、局所部分空間法のみでなく、ＡＲＴ２のようなクラスタに基づく方法や特許文献１に記載の類似度ベースモデリングと呼ばれる手法など、観測データによって基準が変化する事例ベースの手法すべてにあてはまる。
　そこで、本発明の目的は、上記課題を解決し、多次元時系列センサ信号に基づいて高感度に異常を検出し、検出した異常に関連するセンサを正しく診断する異常診断方法および装置を提供することにある。
　上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて設備または装置の異常を診断する方法及びその装置において、予め指定された学習期間における複数のセンサ信号の特徴値から多次元の特徴ベクトルを求め、診断時に複数のセンサ信号の特徴値から時刻毎に多次元特徴ベクトルを抽出し、求めた予め指定された学習期間の多次元の特徴ベクトルの集合と診断時に抽出した時刻毎の多次元の特徴ベクトルに基づいて時刻毎の基準特徴ベクトルを算出し、抽出した時刻毎の多次元特徴ベクトルと算出した時刻毎の基準特徴ベクトルとの差に基づき異常測度を算出し、算出した異常測度と所定のしきい値とを比較することにより異常を検出し、異常が検出された時刻における複数のセンサ信号の特徴値の２次元の分布密度に基づいて検出した異常に関連するセンサを複数のセンサの中から特定するようにした。
　又、上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法及びその装置において、学習データを作成して蓄積する工程は、複数のセンサ信号から特徴ベクトルを抽出する工程と、抽出した特徴ベクトルを学習データとして蓄積する工程と、学習データとして蓄積した特徴ベクトルの各々について特徴ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、作成した学習用の基準ベクトルに基づいて学習データとして蓄積した特徴ベクトルの各々について異常測度を算出する工程と、算出した異常測度に基づきしきい値を算出する工程と、特徴ベクトルから全組合せの２次元の特徴分布密度を算出する工程とを有し、センサ信号の異常診断を行う工程は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、抽出した観測ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、作成した異常診断用の基準ベクトルに基づいて観測ベクトルの異常測度を算出する工程と、算出した異常測度と学習データを作成して蓄積する工程において算出したしきい値に基づき観測ベクトルが異常か正常かを判定する工程と、判定する工程において異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における観測ベクトルと２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有して構成した。
　さらにまた、上記課題を解決するために、本発明では、設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法及びその装置において、学習データを作成して蓄積する工程は、設備または装置から出力されるイベント信号に基づいて稼動状態別のモード分割を行う工程と、複数のセンサ信号から特徴ベクトルを抽出する工程と、抽出した特徴ベクトルを学習データとして蓄積する工程と、学習データとして蓄積した特徴ベクトルのおのおのについて特徴ベクトルに応じて学習データの中から所定数の特徴ベクトルを選択する工程と、選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、作成した学習用の基準ベクトルに基づき学習データとして蓄積した特徴ベクトルの異常測度を算出する工程と、算出した異常測度に基づきモード毎にしきい値を算出する工程と、学習データとして蓄積した特徴ベクトルから全組合せの２次元の特徴分布密度をモード毎に算出する工程を有し、センサ信号の異常診断を行う工程は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、抽出した観測ベクトルに応じて蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、作成した異常診断用の基準ベクトルに基づいて観測ベクトルの異常測度を算出する工程と、算出した異常測度と分割したモードと分割したモード別に算出したしきい値に基づき観測ベクトルが異常か正常かを判定する工程と、判定により異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における観測ベクトルとモードとモード別に算出した２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有して構成した。
　本発明によれば、各時刻のデータの異常検出を、学習データとその時刻のデータから算出される基準データからの距離に基づいて行うため、高感度な異常検出が可能である。異常が検出された時刻の関連センサ特定において、寄与量が大きいセンサをそのまま関連センサとするのではなく、センサ信号の２次元の分布密度に基づいて異常に関連するセンサを特定するため、正しく診断することができる。
　以上の手法を適用したシステムにより、ガスタービンや蒸気タービンなどの設備のみならず、水力発電所での水車、原子力発電所の原子炉、風力発電所の風車、航空機や重機のエンジン、鉄道車両や軌道、エスカレータ、エレベータ、そして機器・部品レベルでは、搭載電池の劣化・寿命など、種々の設備・部品において、対象の異常を早期に検出するのみならず、関連センサに応じた調査あるいは対策を迅速に行うことが可能となる。
本発明の実施例１に関する設備状態監視システムの概略の構成を示すブロック図である。 局所部分空間法を説明する図である。 寄与量に基づくセンサ信号特定における誤りの例を説明する図で、正常データの分布を示す３次元のグラフである。 寄与量に基づくセンサ信号特定における誤りの例を説明する図で、３つのセンサ信号の基準値に対する観測値の変化を示すセンサ信号のグラフである。 センサ信号の例を示す信号リストの表である。 センサ信号入力からしきい値設定までの処理の流れを示すフロー図である。 ２次元の特徴の分布密度を算出する処理の流れを示すフロー図である。 センサ信号入力から異常判定までの処理の流れを示すフロー図である。 異常に関連するセンサを特定する処理の流れを示すフロー図の例である。 異常に関連するセンサを特定する処理の流れを示すフロー図の別の例である。 異常に関連する３個のセンサを特定する処理の流れを示すフロー図である。 異常が連続する区間毎に関連するセンサを特定する処理の流れを示すフロー図である。 分布密度画像に異常区間の対応点をプロットした例を表す図である。 分布密度画像に異常区間の対応点をプロットした例を表す図である。 分布密度画像に異常区間の対応点をプロットした例を表す図である。 分布密度画像に異常区間の対応点をプロットした例を表す図である。 レシピ設定のためのＧＵＩの１例を表す表示画面の正面図である。 結果表示画面の例で複数の時系列データを表示した画面の正面図である。 結果表示画面の例で複数の時系列データを拡大表示した画面の正面図である。 結果表示画面の例で異常関連センサの分布密度と複数の時系列データを表示した画面の正面図である。 レシピ設定のためのテスト結果の一覧を表示した画面の正面図である。 結果表示対象指定のためのＧＵＩの１例を表す表示画面の正面図である。 結果表示画面に含まれる期間表示ウィンドウの表示例である。 本発明の実施例２に関する設備状態監視システムの概略の構成を示すブロック図である。 本発明の実施例２に関するイベント信号の例を示す信号リストの表である。 本発明の実施例２に関するイベント信号に基づくモード分割処理の流れを示すフロー図である。 本発明の実施例２に関する設備の可動状態を分割して４種のモードの何れかに分類した状態を示すイベント信号の模式図である。
　本発明は、設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて設備または装置の異常を診断する方法及びその装置において、センサ信号から時刻毎に多次元の特徴ベクトルを抽出し、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準ベクトルを算出し、各時刻の特徴ベクトルと基準特徴ベクトルとの差に基づいて異常測度を算出し、異常測度と所定のしきい値との比較により異常を検出し、異常が検出された時刻における特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定できるようにしたものである。
　以下に、本発明の実施例を図を用いて説明する。
　以下、図面を用いて本発明の内容を詳細に説明する。
　図１に、本発明の異常診断方法を実現するシステムの一構成例を示す。
　本システムは、設備１０１から出力されるセンサ信号１０２を蓄積するセンサ信号蓄積部１０３、センサ信号１０２をもとに特徴ベクトルを抽出する特徴ベクトル抽出部１０４、予め指定された学習期間の特徴ベクトルの集合と各時刻の特徴ベクトルに基づいて各時刻の基準特徴ベクトルを算出する基準ベクトル算出部１０５、各時刻の特徴ベクトルと基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出部１０６、予め指定された学習期間の異常測度に基づきしきい値を算出するしきい値算出部１０７、異常測度と算出されたしきい値との比較により異常を検出する異常検出部１０８、学習期間のセンサ信号の２次元の分布密度を算出する分布密度算出部１０９、異常が検出された時刻について、異常に関連するセンサを特定する関連センサ特定部１１０を備えて構成される。
　本システムの動作には、蓄積されたデータを用いて学習データの生成、保存を行う「学習」と入力信号に基づき異常を検出し、関連センサの特定を行う「異常診断」の二つのフェーズがある。基本的に前者はオフラインの処理、後者はオンラインの処理である。ただし、後者をオフラインの処理とすることも可能である。以下の説明では、それらを学習時、異常診断時という言葉で区別する。
　状態監視の対象とする設備１０１は、ガスタービンや蒸気タービンなどの設備やプラントである。設備１０１は、その状態を表すセンサ信号１０２を出力する。センサ信号１０２はセンサ信号蓄積部１０３に蓄積されている。センサ信号１０２をリスト化して表形式に表した例を図４に示す。センサ信号１０２は一定間隔毎に取得される多次元時系列信号であり、それをリスト化した表は、図４に示すように、日時の欄４０１と設備１０１に設けられた複数のセンサ値のデータの欄４０２からなる。センサの種類は、数百から数千といった数になる場合もあり、例えば、シリンダ、オイル、冷却水などの温度、オイルや冷却水の圧力、軸の回転速度、室温、運転時間などである。出力や状態を表すのみならず、何かをある値に制御するための制御信号の場合もある。
　学習時の処理の流れを図５を用いて説明する。学習時はセンサ信号蓄積部１０３に蓄積されたデータのうち指定された期間のデータを用いて、特徴ベクトルを抽出する。ここで抽出されたデータを学習データと呼ぶ。交差検証によって学習データの異常測度を算出し、全学習データの異常測度に基づき異常判定のしきい値を算出する。また、分布密度算出部１０９において特徴値の２次元の分布密度を算出する。
　まず、特徴ベクトル抽出部１０４、基準ベクトル算出部１０５、異常測度算出部１０６、しきい値算出部１０７における処理の流れを図５を用いて説明する。はじめに、特徴ベクトル抽出部１０４において、センサ信号蓄積部１０３から学習期間として指定された期間のセンサ信号１０２を入力し（Ｓ５０１）、センサ信号毎に正準化した後（Ｓ５０２）、特徴ベクトルの抽出を行う（Ｓ５０３）。次に、基準ベクトル算出部１０５において、抽出された特徴ベクトルに交差検証用のグループ番号を付与する（Ｓ５０４）。グループは、１日分を１グループとするなど、期間を指定して決めるか、予め指定されたグループ数に等分に分けて決める。
　抽出された特徴ベクトルから、１個目の特徴ベクトルに注目し（Ｓ５０５）、注目ベクトルと異なるグループの特徴ベクトルの中から、注目ベクトルに近い順に予め指定された数の特徴ベクトルを選択し（Ｓ５０６）、選択した特徴ベクトルを用いて基準ベクトルを作成する（Ｓ５０７）。異常測度算出部１０６において、注目した特徴ベクトルの基準ベクトルまでの距離に基づいて異常測度を算出する（Ｓ５０８）。全ベクトルの異常測度算出が終了していれば（Ｓ５０９）、しきい値算出部１０７において、全ベクトルの異常測度に基づいてしきい値を設定する（Ｓ５１１）。ステップＳ５０９において全ての異常測度算出が終了していなければ、次の特徴ベクトルに注目し（Ｓ５１０）、ステップＳ５０６からＳ５０９の処理を繰り返す。
　次に、各ステップについて詳細に説明する。  
  ステップＳ５０２においては、各センサ信号の正準化を行う。例えば、指定された期間の平均と標準偏差を用いて、平均を０、分散を１となるように変換する。異常診断時に同じ変換ができるよう、各センサ信号の平均と標準偏差を記憶しておく。あるいは、各センサ信号の指定された期間の最大値と最小値を用いて最大が１、最小が０となるように変換する。あるいは、最大値と最小値の代わりに予め設定した上限値と下限値を用いてもよい。異常診断時に同じ変換ができるよう、各センサ信号の最大値と最小値または上限値と下限値を記憶しておく。センサ信号の正準化は、単位およびスケールの異なるセンサ信号を同時に扱うためのものである。
　ステップＳ５０３においては、時刻毎に特徴ベクトル抽出を行う。センサ信号を正準化したものをそのまま並べることが考えられるが、ある時刻に対して±1，±2，…のウィンドウを設け，[ウィンドウ幅（3，5，…）×センサ数]の特徴ベクトルにより、データの時間変化を表す特徴を抽出することもできる。また、離散ウェーブレット変換(DWT: Discrete Wavelet Transform)を施して、周波数成分に分解してもよい。さらに、ステップＳ５０３において、特徴選択を行う。最低限の処理として、分散が非常に小さいセンサ信号および単調増加するセンサ信号を除く必要がある。
　また、相関解析による無効信号を削除することも考えられる。これは、多次元時系列信号に対して相関解析を行い、相関値が１に近い複数の信号があるなど、極めて類似性が高い場合に、これらは冗長だとして、この複数の信号から重複する信号を削除し、重複しないものを残す方法である。このほか、ユーザが指定するようにしてもよい。また、長期変動が大きい特徴を除くことも考えられる。長期変動が大きい特徴を用いることは正常状態の状態数を多くすることにつながり、学習データの不足を引き起こすためである。例えば、１周期期間毎の平均と分散を算出し、それらのばらつきによって長期変動の大きさを推定できる。
　基準ベクトル作成手法としては、局所部分空間法(LSC: Local Sub-space Classifier)や投影距離法(PDM: Projection Distance Method)が考えられる。
　局所部分空間法は、注目ベクトルqのk－近傍ベクトルを用いてｋ－１次元のアフィン部分空間を作成する方法である。図２にｋ＝３の場合の例を示す。k-近傍ベクトルはすなわちステップＳ５０６において選択された特徴ベクトルである。指定する数は、ｋの値である。
　図２に示すように、異常測度は図に示す投影距離で表されるため、注目ベクトルqに最も近いアフィン部分空間上の点bを求めればよい。評価データqとそのk－近傍ベクトルxi( i＝1，…，k )からbを算出するには、qをk個並べた行列Qとxiを並べた行列Xから
Figure JPOXMLDOC01-appb-M000001
 により相関行列Cを求め、
Figure JPOXMLDOC01-appb-M000002
 によりbを計算する。ここまでが、ステップＳ５０７における基準ベクトル作成にあたる。
異常測度ｄはｑとｂの間の距離であるから次式で表される。　
Figure JPOXMLDOC01-appb-M000003
　なお、図２ではｋ＝３の場合を説明したが、特徴ベクトルの次元数より十分小さければいくつでもよい。ｋ＝１の場合は、最近傍法と等価の処理になる。
　投影距離法は、選択された特徴ベクトルに対し独自の原点をもつ部分空間すなわちアフィン部分空間（分散最大の空間）を作成する方法である。ステップＳ５０６において指定する数はいくつでも良いが、大きすぎるとベクトルの選択、部分空間の算出ともに時間がかかってしまうので数十から数百とするとよい。
　アフィン部分空間の算出方法について説明する。まず、選択された特徴ベクトルの平均μと共分散行列Σ を求め、次にΣの固有値問題を解いて値の大きい方から予め指定したr個の固有値に対応する固有ベクトルを並べた行列Ｕをアフィン部分空間の正規直交基底とする。ｒは特徴ベクトルの次元より小さくかつ選択データ数より小さい数とする。あるいはｒを固定した数とせず、固有値の大きい方から累積した寄与率が予め指定した割合を超えたときの値としてもよい。異常測度は、注目ベクトルのアフィン部分空間への投影距離とする。
　この他、注目ベクトルqのk-近傍ベクトルの平均ベクトルまでの距離を異常測度とする局所平均距離法や、ガウシアンプロセスなどを用いてもよい。
　ステップＳ５１１におけるしきい値を設定する方法について説明する。学習期間の全特徴ベクトルの異常測度を昇順にソートし、予め指定した１に近い比率に到達する値を求める。この値を基準としてオフセットを加える、定数倍するなどの処理によりしきい値を算出する。オフセット０、倍率を１とすれば、この値がしきい値となる。算出したしきい値は、図示していないが、学習データと対応付けて記録しておく。
　上記に示したように、注目ベクトルの近傍ベクトルを用いて基準ベクトルを作成することにより、状態が複雑に変化する設備に対しても適切な基準により精度の高い異常測度を算出することが可能である。学習データの交差検証により算出された異常測度に基づいてしきい値を算出するため、誤報を抑制することができる。
　次に、分布密度算出部１０９における処理の流れを図６を用いて説明する。　
　始めに、学習期間の特徴ベクトルを入力する（Ｓ６０１）。最初の特徴に注目し（Ｓ６０２）、学習期間のデータの最大値（ＭＡＸ）と最小値（ＭＩＮ）を求める（Ｓ６０３）。次に、最小値から最大値を指定された数Ｎで分割する際の刻み幅Ｓを算出する（Ｓ６０４）。Ｓ＝（ＭＡＸ－ＭＩＮ）／Ｎ　で計算できる。次に、最小値と最大値から外側に範囲を広げて分布密度算出の処理範囲を算出する（Ｓ６０５）。広げる範囲は例えばＭＩＮをＭＩＮ－Ｓ×Ｍ、ＭＡＸをＭＡＸ＋Ｓ×Ｍに変更する。ここでＭは予め決められた１以上の整数とする。
　次に学習期間の全データについて、特徴値（Ｆ）からビン番号（ＢＮＯ）を次式で算出する（Ｓ６０６）。　
　ＢＮＯ＝ＩＮＴ（（Ｆ－ＭＩＮ）・Ｎ／（ＭＡＸ－ＭＩＮ））
ただしＩＮＴ（Ｘ）はＸの整数部を表す。
　以上の処理について未処理の特徴が残っている場合は（Ｓ６０７）、次の特徴に注目し（Ｓ６０８）、ステップＳ６０３からＳ６０６の処理を行う。ステップＳ６０７において全特徴について処理していれば、ステップＳ６０９に進む。
　ステップＳ６０９においては、特徴２個の最初の組合せに注目する。２個の特徴は同じものであってもよい。次に分布密度算出用の２次元配列を確保し、すべての要素に０をセットする（Ｓ６１０）。配列のサイズはＮ＋２Ｍである。学習期間の全データについて、２個の特徴のビン番号に対応する配列の要素に１を加算する（Ｓ６１１）。この処理により、特徴２個による２次元の頻度分布（ヒストグラム）が算出される。この頻度分布を画像に変換して保存する（Ｓ６１２）。変換方法については後述する。上記処理が特徴２個の全ての組合せについて実施済みであれば（Ｓ６１３）、処理を終了し（Ｓ６１５）、そうでなければ次の組合せに注目し（Ｓ６１４）、ステップＳ６１０からＳ６１２までの処理を行う。図示はしていないが、２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を記録しておく。
　ステップＳ６１２における、画像変換方法の例を説明する。始めに配列要素の最大値すなわち最大頻度を求める。画像サイズは配列サイズと同じとし、各要素の値から対応する座標の画素値を例えば、[２５５×配列の要素値／最大頻度]とする。２５５は画素値を８ビットで表す場合の最大値であり、この値を用いれば、そのままビットマップ形式で保存できる。あるいは、画素値を[２５５×ＬＯＧ（配列の要素値＋１）／ＬＯＧ（最大頻度＋１）]とする。ただしＬＯＧ（Ｘ）はＸの対数を表す。このような変換式を用いれば、最大頻度が大きい場合も非ゼロの頻度に非ゼロの画素値を対応させることが可能になる。
　上記処理により得られた画像は、２次元の特徴空間上で密度が高いところが高い画素値で表されているため、分布密度画像と呼ぶこととする。画像の作り方は、上記方法に限定されない。例えば単純な頻度分布ではなく、１個のデータにガウス分布や他の重みつきフィルタを割り当て、それを重畳するようにしてもよい。あるいは、上記方法で得られた画像に所定サイズの最大値フィルタをかけたり、平均フィルタ、その他の重みつきフィルタをかけたりしてもよい。また必ずしも画像形式で保存する必要はなく２次元配列をテキスト形式で保存してもよい。
　異常診断時の処理の流れを図７ないし図９を用いて説明する。異常診断時はセンサ信号蓄積部１０３に蓄積されたデータのうち指定された期間のデータあるいは新たに観測されたデータの異常測度を算出し、正常か異常かの判定を行い、異常と判定された時刻について関連するセンサを特定する。
　図７は特徴ベクトル抽出部１０４、基準ベクトル算出部１０５、異常測度算出部１０６、異常検出部１０８における異常診断時の処理の流れを説明する図である。始めに、特徴ベクトル抽出部１０４において、センサ信号蓄積部１０３または設備１０１からセンサ信号を入力し（Ｓ７０１）、センサ信号毎に正準化した後（Ｓ７０２）、特徴ベクトルの抽出を行う（Ｓ７０３）。センサ信号の正準化には、図５のステップＳ５０２に示す処理において、学習データの正準化に用いたものと同じパラメータを用いる。特徴ベクトルの抽出は、ステップＳ５０３の処理と同じ方法で行う。したがって、ステップＳ５０３において特徴選択を実行した場合は同じ特徴を選択する。ここで抽出された特徴ベクトルを、学習データと区別するために観測ベクトルと呼ぶこととする。
　次に、基準ベクトル作成部１０５において、学習データの特徴ベクトルの中から、観測ベクトルに近い順に予め指定された数の特徴ベクトルを選択し（Ｓ７０４）、それらの特徴ベクトルを用いて基準ベクトルを作成する（Ｓ７０５）。異常測度算出部１０６において、観測ベクトルの基準ベクトルまでの距離に基づいて異常測度を算出する（Ｓ７０６）。異常検出部１０８において、学習時に算出したしきい値と異常測度を比較して、異常測度がしきい値より大きければ異常、そうでなければ正常と判定する（Ｓ７０７）。
　図８は、関連センサ特定部１１０における処理の流れの一実施例を説明する図である。この処理は、ステップＳ７０７において、異常と判定されたデータについて実施する。始めに、観測ベクトルと基準ベクトルから、各特徴の観測値と基準値の差を算出し（Ｓ８０１）、これが最大となる特徴を探す（Ｓ８０２）。これを特徴Ａとしておく。次に観測ベクトルから、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ８０３）。一方の特徴が特徴Ａである分布密度画像を読み込み、各特徴と特徴Ａのビン番号に対応する座標の画素値を読み（Ｓ８０４）、これが最小となる特徴を探す（Ｓ８０５）。これを特徴Ｂとし、特徴Ａと特徴Ｂの分布密度画像に観測ベクトルの対応点をプロットする（Ｓ８０６）。例えば、カラー画像に変換して、特徴Ａと特徴Ｂのビン番号に対応する座標に使用していない色をつける。Ｓ８０４における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ａと特徴Ｂを関連センサとする（Ｓ８０７）。
　図９は、関連センサ特定部１１０における処理の流れの別の一実施例を説明する図である。図８の処理と同様、ステップＳ７０７において、異常と判定されたデータについて実施する。始めに、観測ベクトルから、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ９０１）。
　次に全ての分布密度画像を読み込み、それぞれの分布密度画像について２個のビン番号に対応する座標の画素値を読み（Ｓ９０２）、特徴毎にその値の合計siを算出する（Ｓ９０３）。つまり、特徴iと特徴jの分布密度画像から読み取った値をI(i，j)とすると、siは次式で計算される。　
Figure JPOXMLDOC01-appb-M000004
　このsiが最小となる特徴を探し（Ｓ９０４）、これを特徴Ｃとしておく。次に一方の特徴が特徴Ｃである分布密度画像について、各特徴と特徴Ｃのビン番号に対応する座標の画素値を読み（Ｓ９０５）、これが最小となる特徴を探す（Ｓ９０６）。これを特徴Ｄとし、ステップＳ８０６と同様の方法で特徴Ｃと特徴Ｄの分布密度画像に観測ベクトルの対応点をプロットする（Ｓ９０７）。Ｓ９０５における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ｃと特徴Ｄを関連センサとする（Ｓ９０８）。
　また別の実施例では、ステップＳ９０１で各特徴のビン番号を算出し、ステップＳ９０２で全ての分布密度画像について観測ベクトルに対応する座標の画素値を読んだ後、これが最小となる特徴の組合せを探し、これらを特徴Ｅおよび特徴Ｆとする。ステップＳ８０６と同様の方法で特徴Ｅと特徴Ｆの分布密度画像に観測ベクトルの対応点をプロットし、特徴Ｅと特徴Ｆの分布密度画像から読んだ画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ｅと特徴Ｆを関連センサとする。
　以上、説明した関連センサ特定の方法はいずれかひとつを実施してもよいし、関連センサが特定できるまでシーケンシャルに実施してもよいし、全て実施して観測ベクトルをプロットした分布密度画像を表示して、ユーザに決めさせてもよい。
　さらに、上記方法で関連センサが特定できない場合に詳細に解析する方法について、図１０を用いて説明する。始めに、学習期間の特徴ベクトルを入力する（Ｓ１００１）。次に学習期間の全データについて、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ１００２）。同様に、観測ベクトルについて各特徴のビン番号を算出する（Ｓ１００３）。
　次に、学習期間の特徴ベクトルから前述の特徴Ａのビン番号が観測ベクトルと同じデータのみ抽出し（Ｓ１００４）、これらを用いてステップＳ６１０からＳ６１２に従い、特徴Ａを除く組合せの部分的分布密度画像を作成する（Ｓ１００５）。全ての部分的分布密度画像について２個のビン番号に対応する座標の画素値を読み（Ｓ１００６）、特徴毎にその値の合計を算出する（Ｓ１００７）。これが最小となる特徴を探し（Ｓ１００８）、これを特徴Ｇとしておく。
　一方の特徴が特徴Ｇである部分的分布密度画像について、各特徴と特徴Ｇのビン番号に対応する座標の画素値を読み（Ｓ１００９）、これが最小となる特徴を探し特徴Ｈとする（Ｓ１０１０）。特徴Ｇと特徴Ｈの部分的分布密度画像に観測ベクトルの対応点をプロットする（Ｓ１０１１）。Ｓ１００９における画素値が０であれば、あるいは予め決められた値より小さければ、特徴Ａと特徴Ｇと特徴Ｈを関連センサとする（Ｓ１０１２）。なお、これらの処理では特徴Ａのかわりに前述の特徴Ｃを用いてもよい。
　上記方法は、異常検出した時刻毎に関連センサを特定する方法であるが、連続して異常検出された場合に適用すると冗長な処理になってしまう。そこで、異常が連続して検出されている間はバッファに情報をためておき、まとめて処理を行うことが考えられる。この場合、異常の継続時間や累積異常度でセンサ特定を行う対象を絞り込んでもよい。
　連続して検出された一連の異常区間について関連するセンサを特定する方法を、図１１を用いて説明する。まず、異常検出部１０８においてステップＳ７０７で異常判定された結果をもとに、連続して異常検出された区間を特定する（Ｓ１１０１）。その間、観測ベクトルと基準ベクトルはバッファに蓄積しておく。なお「連続」というのは必ずしも厳密な意味ではなく、予め決められた長さの中断を含んでいてもよい。次に、上記区間に含まれる全ての観測ベクトルについて、特徴毎に観測値と基準値の差の絶対値の和を算出し（Ｓ１１０２）、これが最大となる特徴を探す（Ｓ１１０３）。これを特徴Ａとしておく。
　次に、上記区間に含まれる全ての観測ベクトルについて、図６で示した分布密度算出処理で使用した２次元配列のサイズおよびステップＳ６０５で算出した各特徴の最小値と最大値を用いて、各特徴のビン番号を算出する（Ｓ１１０４）。一方の特徴が特徴Ａである分布密度画像を読み込み、上記区間に含まれる全ての観測ベクトルについて、各特徴と特徴Ａのビン番号に対応する座標の画素値を読む（Ｓ１１０５）。上記区間に含まれる全ての観測ベクトルについて特徴毎に非０の画素をカウントする（Ｓ１１０６）。特徴iについてカウントした値をcnt(i)とする。また、上記区間に含まれる全ての観測ベクトルについて特徴毎に画素値の合計を算出する（Ｓ１１０７）。特徴iについて算出した値をsum(i)とする。cnt(i)が最小となる特徴を探し、複数ある場合はそのうちsum(i)が最小となる特徴を探す（Ｓ１１０８）。これを特徴Ｂとし、特徴Ａと特徴Ｂの分布密度画像に異常区間に含まれる観測ベクトルの対応点をプロットする（Ｓ１１０９）。Ｓ１１０６における非０の画素数が区間に含まれる観測ベクトル数より小さければ、つまり０の画素数が１以上あれば、特徴Ａと特徴Ｂを関連センサとする（Ｓ１１１０）。
　あるいは、cnt(i)が最小となる特徴が複数ある場合は、学習データの相関の高い特徴を探して特徴Ｂとしてもよい。相関の高さは例えば特徴Ａともう一方の特徴の相関値に基づいて決める。または、分布密度画像の非０の画素の数が少ないものを相関が高いと考える。
　なお、上記に説明した図１１の処理フローは、図８に示す処理を複数の観測ベクトルに拡張した処理であるが、図９または図１０に示す処理を複数の観測ベクトルに拡張して用いてもよいし、それらを組み合わせて用いてもよい。
　ステップＳ１１０９にて作成される分布密度画像の例を図１２Ａ～Ｃに示す。横軸と縦軸はそれぞれの特徴の相対的な値を表し、色の濃い部分が学習データの密度が高い部分である。星印は異常区間に含まれる観測データを表す。実際には色をつけて表示すればよい。
　図１２Ａは、特徴ａ、ｂの分布密度を表す。学習データの分布からは、特徴ａ、ｂとも低い状態と特徴ａ、ｂとも高い状態の定常状態があり状態遷移時は特徴ａに遅れて特徴ｂが上昇していることがわかる。時間ずれの関係は安定している。特徴ａ上昇のタイミングに対して特徴ｂの上昇が早く起こっていることが異常として検出されている。あるいは下降タイミングのずれを検出しているかもしれないが、どちらであるかはセンサ信号の時系列グラフから確認できる。
　図１２Ｂは、特徴ｃ、ｄの分布密度を表す。特徴ｃ、ｄとも低い状態、特徴ｃが高く特徴ｄが低い状態、特徴ｃ、ｄとも高い状態とそれらの中間の状態はあるが特徴ｄが高く特徴ｃが低い状態にはならない。そのような状態が異常として検出されている。
　図１２Ｃは、横軸も縦軸も特徴eの値を表したものである。実質的には２次元ではないが、同じ形式の画像として表すことにより、表示やセンサ特定の処理を簡便にすることができる。学習データは直線上に分布する。特徴eの値が学習データのどの値よりも高いことが異常として検出されている。
　上記の通り、予め学習データを用いて全ての組合せの２次元特徴分布密度を算出しておき、異常が検出されたときの観測ベクトルの対応点の分布密度に基づいて異常に関連するセンサを特定するため、図３Ａ及び図３Ｂに示すような状況である場合でも正しいセンサを特定することが可能である。また、２次元特徴分布密度を画像で表し、これをユーザに提示することにより、それらのセンサが異常に関連することが理解されやすくなる。また、学習データの分布の傾向を知ることができる。例えば、図１２Ｄに示すように、学習データの密度の低い部分が広い範囲に存在する場合に、星印１２０１で示す異常を検出したとしてもそれが本当に異常かどうかは疑わなくてはならない。その場合、学習期間を増やして密度が高くなるようにするか、この組合せでは一緒に処理しないことにする、具体的には特徴ｇを処理対象からはずすようにするとよい。
　以上の方法を実現するシステムのＧＵＩの実施例を説明する。  
  学習期間および処理パラメータ設定のためのＧＵＩの例を、図１３に示す。以下の説明ではこの設定のことを単にレシピ設定と呼ぶことにする。また、過去のセンサ信号１０２は設備ＩＤおよび時刻と対応付けられてデータベースに保存されているものとする。レシピ設定画面１３０１では、対象装置、学習期間、使用センサ、基準算出パラメータ、しきい値設定パラメータを入力する。設備ＩＤ入力ウィンドウ１３０２には、対象とする設備のＩＤを入力する。設備リスト表示ボタン１３０３押下により図示はしていないがデータベースに保存されているデータの装置ＩＤのリストが表示されるので、リストから選択入力する。
　学習期間入力ウィンドウ１３０４には、学習データを抽出したい期間の開始日と終了日を入力する。センサ選択ウィンドウ１３０５には、使用するセンサを入力する。リスト表示ボタン１３０６のクリックによりセンサリスト１３０７が表示されるので、リストから選択入力する。リストから複数選択することも可能である。基準算出パラメータ入力ウィンドウ１３０８には、正常モデル作成において指定するパラメータを入力する。図は正常モデルとして局所部分空間を採用した場合の例であり、モデル作成に使う近傍ベクトル数と正則化パラメータを入力する。正則化パラメータは、数２において相関行列Cの逆行列が求められないことを防ぐため、対角成分に加算する小さい数である。
　しきい値設定パラメータ入力ウィンドウ１３０９には、図５に示す処理における交差検証のグループをどのように決めるかをラジオボタンで選択する。１日を１グループとするか、指定された数のグループに分割するかを選び、後者の場合はグループ数を入力する。また、ステップＳ５１１のしきい値設定処理において累積ヒストグラムに適用する比率を値入力する。レシピ名入力ウィンドウ１３１０には、入力された情報に対応付けるユニークな名前を入力する。全ての情報を入力したらテスト期間入力ウィンドウ１３１１にテスト対象期間を入力し、テストボタン１３１２の押下により、レシピのテストを行う。
　この操作により、同じレシピ名で実行したテストの通し番号が採番される。次に、特徴ベクトル抽出部１０４において、指定した学習期間のセンサ信号１０２から特徴ベクトルの抽出を行う。図５で説明したステップＳ５０２の正準化においては、指定した学習期間の全データを用いて平均と標準偏差を求める。この平均と標準偏差の値もセンサ毎にレシピ名およびテスト番号に対応付けて保存しておく。
　図５で説明したステップＳ５０７に対応して、交差検証により基準ベクトル算出部１０５において各時刻の基準ベクトルを作成し、ステップＳ５０８に対応して、異常測度算出部１０６において異常測度を算出する。更に、ステップＳ５１１に対応して、しきい値算出部１０７において異常判定しきい値を算出し、レシピ名およびテスト番号と対応付けて保存する。分布密度算出部１０９において、図６に示した処理により２次元の特徴の分布密度画像を作成し、保存する。併せて、ステップＳ６０５で算出した処理範囲を全ての特徴について保存しておく。その他、装置ＩＤ情報、使用センサ情報、学習期間、特徴ベクトル抽出に用いるパラメータ、基準作成パラメータをレシピ名およびテスト番号と対応付けて保存する。
　次に、指定したテスト期間のセンサ信号１０２を入力として、図７に示した異常検知の処理を行う。検出された異常について、図８乃至図１１の何れかに示した手順で関連センサ特定の処理を行う。検出された異常は区間毎に通し番号が付加される。区間毎に、通し番号、区間開始時刻、区間終了時刻、複数の関連センサ名と観測ベクトルの対応点がプロットされた分布密度画像名が記録される。
　テストの結果をユーザに示すためのＧＵＩの例を図１４Ａおよび図１４Ｂおよび図１４Ｃに示す。各画面の上部に表示されたタブを選択することにより、結果全体表示画面１４０１と結果拡大表示画面１４０２と異常診断結果表示画面１４０３を切り換えることができる。
　図１４Ａには、結果全体表示画面１４０１を示す。結果全体表示画面１４０１には、指定された全期間の異常測度、しきい値、判定結果とセンサ信号の時系列グラフを表示する。期間表示ウィンドウ１４０４には、指定された学習期間およびテスト期間が表示される。異常測度表示ウィンドウ１４０５には、指定された学習期間およびテスト期間の異常測度としきい値と判定結果が表示される。センサ信号表示ウィンドウ１４０６には、指定された期間の指定されたセンサの出力値が表示される。センサの指定は、センサ名選択ウィンドウ１４０７への入力によって行う。ただし、ユーザが指定する前は、先頭のセンサが選択されている。カーソル１４０８は、拡大表示の時の起点を表し、マウス操作により移動できる。
　表示日数指定ウィンドウ１４０９には、この画面では使用しないが、結果拡大表示画面１４０２での、拡大表示の起点から終点までの日数が表示される。この画面で入力することもできる。日付表示ウィンドウ１４１０には、カーソル位置の日付が、表示される。終了ボタン１４１１押下により結果全体表示画面１４０１、結果拡大表示画面１４０２、異常診断結果表示画面１４０３とも消去し終了する。
　図１４Ｂには、結果拡大表示画面１４０２を示す。結果拡大表示画面１４０２には、結果全体表示画面１４０１において、カーソル１４０８で示された日付を起点として、指定された日数の異常測度、しきい値、判定結果とセンサ信号の時系列グラフを表示する。期間表示ウィンドウ１４０４には、結果全体表示画面１４０１と同じ情報が表示される。異常測度表示ウィンドウ１４０５およびセンサ信号表示ウィンドウ１４０６には、結果全体表示画面１４０１と同様の情報が、拡大表示される。表示日数指定ウィンドウ１４０９で、拡大表示の起点から終点までの日数を指定する。
　日付表示ウィンドウ１４１０には、拡大表示の起点の日付が表示されている。スクロールバー１４１２で表示の起点を変更することも可能であり、この変更はカーソル１４０８の位置と日付表示ウィンドウ１４１０の表示に反映される。スクロールバー表示領域１４１３の全体の長さは結果全体表示画面１４０１に表示されている全期間に相当する。また、スクロールバー１４１２の長さは表示日数指定ウィンドウ１４０９で指定された日数に相当し、スクロールバー１４１２の左端部が拡大表示の起点に対応する。異常が検出された時刻には、異常測度グラフの対応する位置に異常区間番号１４１４がバルーン表示される。終了ボタン１４１１押下により終了する。
　図１４Ｃには、異常診断結果表示画面１４０３の例を示す。異常に関連するセンサの２次元の特徴分布密度と時系列グラフを併せて表示する。異常診断結果表示画面１４０３に切り替ったときには、結果拡大表示画面１４０２に表示されている異常のうち最も番号の若い異常を表示対象とする。日付表示ウィンドウ１４１５には、表示対象の異常の日付が表示される。異常区間番号表示ウィンドウ１４１６には、表示対象の異常の区間番号が表示される。ウィンドウ１４１６では数値入力により、異常区間番号を指定することができる。時刻表示ウィンドウ１４１７には、表示対象の異常の区間開始時刻と区間終了時刻が表示される。
　分布密度表示ウィンドウ１４１８には、例えば図１１のステップＳ１１０９において保存された画像を表示する。この画像は、ステップＳ１１１０において特定された２個のセンサの２次元頻度分布を表したものであり、横軸は”SensorX”の規格化された値、縦軸は”SensorY” の規格化された値を表す。１４１９は、学習期間の分布を表しており、濃く表示されているほど高頻度であることを意味する。１４２０は、表示対象の異常区間の観測ベクトルに対応するデータである。異常表示指示ウィンドウ１４２２で、チェックボックスのチェックをはずすことにより、表示しないようにすることもでき、その特徴値の学習期間のデータの分布密度を確認できる。処理範囲表示ウィンドウ１４２１には２個のセンサのステップＳ６０５において算出された処理範囲が表示される。これらの値は、それぞれ分布密度画像の左端に対応する”SensorX”の値、右端に対応する”SensorX”の値、下端に対応する”SensorY”の値、上端に対応する”SensorY”の値である。
　異常測度表示ウィンドウ１４２３には、表示対象の異常区間を含む期間の異常測度としきい値と判定結果が表示される。表示期間は信号の変化が観測可能なよう十分に拡大される長さを予め決めておく。第一関連センサ信号表示ウィンドウ１４２４には、ステップＳ１１１０において特定された２個のセンサのうち一方の、異常測度表示ウィンドウ１４２３と同じ期間の出力値がグラフ表示される。第二関連センサ信号表示ウィンドウ１４２５には、もう一方のセンサ信号出力値がグラフ表示される。基準表示指示ウィンドウ１４２６のチェックボックスにチェックを入れることにより、図示していないが第一関連センサ信号表示ウィンドウ１４２４および第二関連センサ信号表示ウィンドウ１４２５に基準値を異なる色で重ねて表示することができる。
　ボタン１４２７及び１４２８は、表示する異常診断結果を切り替えるためのボタンで、ボタン１４２７「次」をクリックすると、現在表示されている異常のデータよりも異常診断結果表示画面１４０３において次に若い番号の異常のデータが表示される。また、ぼたん１４２８「前」をクリックすると、現在表示されている異常のデータよりも異常診断結果表示画面１４０３において１つ若い番号の異常のデータが表示される。
　又、図示はしていないが、分布密度表示ウィンドウ１４１８には、ステップＳ１１０９において保存された表数の画像を並べて表示するようにしてもよい。
　本図では２個のセンサの時系列グラフを別々に表示したが、１個のウィンドウに異なる色で重ねて表示してもよい。終了ボタン１４１１押下により終了する。
　図１４Ａ～Ｃに示すいずれかの画面で、終了ボタン１４１１押下により異常検出結果および診断結果の確認が終了したら、図１３に示すレシピ設定画面１３０１の表示に戻る。テスト番号表示ウィンドウ１３１３には、上記のテストで採番された番号が表示されている。確認した内容に問題があれば、学習期間や選択センサ、パラメータなどを変更し、テストボタン１３１２の押下により、再度テストを行う。あるいは、一度行ったテストの結果を再度確認することもできる。テスト番号表示ウィンドウ１３１３からテスト番号を選択入力し、表示ボタン１３１４を押下する。この操作により、レシピ名とテスト番号に対応付けて保存された情報をロードし結果全体表示画面１４０１を表示する。タブの切り替えにより結果拡大表示画面１４０２または異常診断結果表示画面１４０３を表示させることもできる。確認が済んだら終了ボタン１４１１押下により、レシピ設定画面１３０１の表示に戻る。
　登録ボタン１３１５の押下により、上記レシピ名とテスト番号表示ウィンドウ１３１３に表示中のテスト番号に対応付けて保存されている情報をレシピ名と対応付けて登録し、終了する。キャンセルボタン１３１６が押下された場合は、何も保存しないで終了する。
　また、テスト結果一覧ボタン１３１７が押下された場合は、図１５に示す、テスト結果一覧表示画面１５０１を表示する。テスト結果リスト１５０２には、全てのテストの学習期間、テスト期間、選択センサ番号、基準作成パラメータ、しきい値設定パラメータなどのレシピ情報と、しきい値、異常区間数などのテスト結果情報を表示する。リストの左端に選択チェックボタンがあり、いずれか一つのみ選択することができる。詳細表示ボタン１５０３押下により、レシピ名とテスト番号に対応付けて保存された情報をロードし、結果全体表示画面１４０１を表示する。タブの切り替えにより結果拡大表示画面１４０２または異常診断結果表示画面１４０３を表示させることもできる。確認が済んだら終了ボタン１４１１押下により、テスト結果一覧表示画面１５０１の表示に戻る。登録ボタン１５０４の押下により、選択中のテスト番号に対応付けて保存されている情報をレシピ名と対応付けて登録し、テスト結果一覧表示画面１５０１の表示およびレシピ設定画面１３０１の表示を終了する。戻るボタン１５０５が押下された場合は、レシピの登録は行わずにレシピ設定画面１３０１の表示に戻る。
　登録されたレシピは、活性か不活性かのラベルをつけて管理され、新しく観測されたデータに対しては、装置ＩＤが一致する活性なレシピの情報を用いて図７ないし図１１を用いて説明した何れかの特徴ベクトル抽出から異常検出、関連センサ特定までの処理を行い、結果をレシピ名と対応付けて保存しておく。
　以上の異常検診断処理の結果をユーザに示すためのＧＵＩの例を、図１６に示す。　
　図１６は、表示対象を指定するＧＵＩの例である。表示対象指定画面１６０１から表示対象の設備、レシピおよび期間を指定する。初めに、装置ＩＤ選択ウィンドウ１６０２により装置ＩＤを選択する。次に、レシピ名選択ウィンドウ１６０３により、装置ＩＤを対象としたレシピのリストから表示対象のレシピを選択する。データ記録期間表示部１６０４には、入力されたレシピを用いて処理され、記録が残されている期間の開始日と終了日が表示される。結果表示期間指定ウィンドウ１６０５には、結果を表示したい期間の開始日と終了日を入力する。表示センサ指定ウィンドウ１６０６には、表示したいセンサの名を入力する。表示ボタン１６０７押下により図１４Ａに示す結果全体表示画面１４０１を表示する。終了ボタン１６０８押下により終了する。
　結果表示にかかわるＧＵＩの画面および操作は、図１４Ａ～Ｃに示すテスト結果表示にかかわるＧＵＩとほぼ同じであるため、異なる部分のみ説明する。結果全体表示画面１４０１、結果拡大表示画面１４０２における期間表示ウィンドウ１４０４に対応する表示ウィンドウとしては、図１７に示すように結果表示期間指定ウィンドウ１６０５で指定された表示期間が表示される。結果全体表示画面１４０１に対応する画面において、異常測度表示ウィンドウ１４０５に対応するウィンドウには、指定された表示期間の異常測度としきい値と判定結果が表示される。センサ信号表示ウィンドウ１４０６に対応するウィンドウには、指定された期間の、表示センサ指定ウィンドウ１６０６により指定されたセンサの出力値が表示される。表示対象センサは、センサ名選択ウィンドウ１４０８に対応するウィンドウへの入力によって変更することも可能である。
　上記実施例は学習データ設定をオフライン、異常診断処理をリアルタイム、結果表示をオフラインでそれぞれ処理するものであるが、結果表示もリアルタイムに行うことが可能である。その場合、表示期間の長さ、表示対象とするレシピ、表示対象とする情報を予め定めておき、一定時間毎に最新の情報を表示するよう構成すればよい。
　逆に、任意の期間を設定し、レシピを選択して、オフラインで異常診断処理を行う機能を付加したものも本発明の範囲に含まれる。
　本実施例によれば、状態が複雑に偏けする設備に対しても精度の高い異常測度を算出することができる。又、誤報の発生を抑制することができる。
　更に、本実施例によれば、正常なセンサの信号から得られるデータが異常なセンサの信号のデータの影響を受けるような場合であっても、検出した異常に関連するセンサを正しく診断することが可能になった。
　また、検査の結果得られる２次元特徴分布密度を画像で表示するので、異常に関連するセンサの見分けが容易になる。さらに、学習データの分布の傾向を画面上で確認することができるようになった。
　以上、設備から出力されるセンサ信号に基づき異常診断する方法の実施例を説明したが、別の実施例として、さらに、設備から出力されるイベント信号も利用して異常診断する方法を説明する。
　図１８に、本発明の異常診断方法を実現するシステムの第二の実施例を示す。本実施例における構成は、実施例１で説明した図１に示した構成に、モード分割部１８０２を加えたものとなっている。実施例１と同じ構成については同じ番号を付してある。モード分割部１８０２は、設備からイベント信号１８０１を入力し、これに基づき設備の稼動状態を表すモードに分割する。モード分割の結果は、しきい値算出部１０７´に入力され、図５で説明した処理フローが分割されたモードごとに実行されてステップＳ５１１のしきい値設定がモード別に実施される。
　また、図７で説明した処理フローも分割されたモードごとに実行されて、異常検出部１０８´におけるステップＳ７０７の異常判定を、対応するモードのしきい値を用いて行う。また、分布密度算出部１０９´でもモード分割の結果を入力し、分布密度画像作成をモード毎に行い、関連センサ特定部１１０´におけるセンサ特定の処理をモード毎に行う。
　上記以外は、全て前述の実施例１で説明した方法と同様であるため、実施例１と異なるイベント信号に基づくモード分割方法の実施例を図１９Ａ～Ｃを用いて説明する。イベント信号の例を図１９Ａに示す。図１９Ａに示したリスト１９２０は、不定期に出力される設備の操作・故障・警告を表す信号であり、時刻と操作・故障・警告を表す文字列またはコードからなる。
　図１９Ｂに示すように、このイベント信号をモード分割部１８０２に入力し（Ｓ１９０１）、所定の文字列またはコードの検索により起動シーケンスと停止シーケンスの切り出しを行う（Ｓ１９０２）。その結果をもとに、停止シーケンスの終了時刻から起動シーケンスの開始時刻までの「定常OFF」モード１９１１、起動シーケンス中の「起動」モード１９１２、起動シーケンスの終了時刻から停止シーケンスの開始時刻までの「定常ON」モード１９１３、停止シーケンス中の「停止」モード１９１４の４つの稼動状態に分割する（Ｓ１９０３）。図１９Ｃに例を示す。
　シーケンス切り出しのためには、予めシーケンスの開始イベントおよび終了イベントを指定しておき、イベント信号の先頭から最後まで以下の要領でスキャンしながら切り出していく。
(1)シーケンスの途中でない場合は、開始イベントを探索する。見つかったらシーケンスの開始とする。
(2)シーケンスの途中の場合は、終了イベントを探索する。見つかったらシーケンスの終了とする。ここで終了イベントとは、指定の終了イベントのほか、故障、警告、指定の開始イベントとする。
　以上のように、イベント信号を利用することにより、多様な稼動状態を正確に分けることができ、モード別にしきい値を設定することにより、「起動」モード１９１２および「停止」モード１９１４の過渡期において学習データ不足により感度を落とす必要がある場合でも、「定常OFF」モード１９１１および「定常ON」モード１９１３では高感度な異常検知が可能になる。
　即ち、本実施例によれば、設備の多様な稼働状態に応じてしきい値を設定することができ、高感度な異常検出を可能にし、更に検出した異常に関連するセンサを正しく診断することが可能になった。また、学習データの分布密度をモード毎に算出することにより、異なる状態のデータが含まれなくなるため、診断が容易になり、分布密度画像の理解も容易になる。
　１０１・・・設備　１０２・・・センサ信号　１０３・・・センサ信号蓄積部　１０４・・・特徴ベクトル抽出部　１０５・・・基準ベクトル作成部　１０６・・・異常測度算出部　１０７・・・しきい値算出部　１０８・・・異常検出部　１０９・・・分布密度算出部　１１０・・・関連センサ特定部　１３０１・・・レシピ設定画面　　１４０１・・・結果全体表示画面　１４０２・・・結果拡大表示画面　１４０３・・・異常診断結果表示画面　　１５０１・・・テスト結果一覧表示画面　　１６０１・・・表示対象指定画面　　１８０１・・・イベント信号　１８０２・・・モード分割部。

Claims (13)

  1. 　設備または装置に装着された複数のセンサから出力される複数のセンサ信号に基づいて前記設備または装置の異常を診断する方法であって、予め指定された学習期間における前記複数のセンサ信号の特徴値から多次元の特徴ベクトルを求め、前記診断時に前記複数のセンサ信号の特徴値から時刻毎に多次元特徴ベクトルを抽出し、前記求めた予め指定された学習期間の前記多次元の特徴ベクトルの集合と前記診断時に抽出した時刻毎の前記多次元の特徴ベクトルに基づいて前記時刻毎の基準特徴ベクトルを算出し、前記抽出した時刻毎の多次元特徴ベクトルと前記算出した時刻毎の基準特徴ベクトルとの差に基づき異常測度を算出し、前記算出した異常測度と所定のしきい値とを比較することにより異常を検出し、前記異常が検出された時刻における前記複数のセンサ信号の特徴値の２次元の分布密度に基づいて前記検出した異常に関連するセンサを前記複数のセンサの中から特定することを特徴とする異常診断方法。
  2. 　上記異常診断方法において、前記所定のしきい値を予め指定された学習期間の異常測度に基づき算出することを特徴とする請求項１に記載の異常診断方法。
  3. 　設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法であって、
    　前記学習データを作成して蓄積する工程は、前記複数のセンサ信号から特徴ベクトルを抽出する工程と、前記抽出した特徴ベクトルを学習データとして蓄積する工程と、前記学習データとして蓄積した特徴ベクトルの各々について前記特徴ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、前記作成した学習用の基準ベクトルに基づいて前記学習データとして蓄積した特徴ベクトルの各々について異常測度を算出する工程と、前記算出した前記異常測度に基づきしきい値を算出する工程と、前記特徴ベクトルから全組合せの２次元の特徴分布密度を算出する工程とを有し、
    　前記センサ信号の異常診断を行う工程は、前記設備または装置に装着された前記複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、前記抽出した観測ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、前記作成した異常診断用の基準ベクトルに基づいて前記観測ベクトルの異常測度を算出する工程と、前記算出した異常測度と前記学習データを作成して蓄積する工程において算出したしきい値に基づき前記観測ベクトルが異常か正常かを判定する工程と、前記判定する工程において異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における前記観測ベクトルと前記２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有する
    ことを特徴とする異常診断方法。
  4. 　前記学習用の基準ベクトルを作成する工程と、前記異常診断用の基準ベクトルを作成する工程とは、何れも局所部分空間法を利用すること特徴とする請求項３に記載の異常診断方法。
  5. 　前記学習データを作成して蓄積する工程において算出される前記２次元の特徴分布密度は、前記特徴ベクトルの集合の２次元の特徴の頻度分布に基づき算出され、濃淡画像に変換して保存されることを特徴とする請求項３に記載の異常診断方法。
  6. 　前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、前記観測ベクトルと前記基準ベクトルに基づき観測値と基準値の差を算出する工程と、前記算出された差が最大となる特徴を構成するセンサを第一の関連センサとする工程と、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
  7. 　前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、全ての前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度を算出する工程と、前記特徴毎に前記算出した密度を合計した値が最小となる特徴を構成するセンサを第一の関連センサとする工程と、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
  8. 　前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、前記学習データのうち前記第一の関連センサのセンサ信号が異常と判定された観測ベクトルと近いデータを用いて前記第一の関連センサを除く全組合せの２次元の部分的特徴分布密度を算出する工程と、全ての前記２次元の部分的特徴分布密度に基づき、前記観測ベクトルに対応する値の密度を算出する工程と、特徴毎に前記密度を合計した値が最小となる特徴を構成するセンサを第二の関連センサとする工程と、前記第二の関連センサにかかわる前記２次元の部分的特徴分布密度に基づき、前記観測ベクトルに対応する値の密度が最小となる前記第二の関連センサのペアとなる特徴を構成するセンサを第三の関連センサとする工程とをさらに含む請求項６に記載の異常診断方法。
  9. 　前記センサ信号の異常診断を行う工程における前記異常に関連するセンサを特定する工程は、連続して異常検出された区間を特定する工程と、前記区間に含まれる全ての観測ベクトルについて、前記観測ベクトルと前記基準ベクトルに基づき観測値と基準値の差の絶対値の和を算出する工程と、前記算出された和が最大となる特徴を構成するセンサを第一の関連センサとする工程と、前記区間に含まれる全ての観測ベクトルについて、前記第一の関連センサにかかわる前記２次元の特徴分布密度に基づき、前記観測ベクトルに対応する値の密度の和および密度が非０となる回数を算出する工程と、前記算出された回数が最小かつそのうち前記算出された和が最小となる前記第一の関連センサのペアとなる特徴を構成するセンサを第二の関連センサとする工程からなる請求項３に記載の異常診断方法。
  10. 　前記異常測度および前記しきい値および前記判定結果の時系列グラフと、前記特定された異常に関連するセンサから出力されたセンサ信号の時系列グラフと、前記特定された２次元の特徴の学習データの分布密度を表す画像に、前記異常と判定された観測ベクトルに対応する点を重ねてプロットした画像を表示する工程をさらに含むことを特徴とする請求項３に記載の異常診断方法。
  11. 　設備または装置に装着された複数のセンサから出力される複数のセンサ信号をもとに学習データを作成して蓄積する工程と、前記設備または装置に装着された複数のセンサから新たに出力された複数のセンサ信号の異常診断を行う工程とを含む設備の異常を診断する方法であって、
    　前記学習データを作成して蓄積する工程は、設備または装置から出力されるイベント信号に基づいて稼動状態別のモード分割を行う工程と、前記複数のセンサ信号から特徴ベクトルを抽出する工程と、前記抽出した特徴ベクトルを学習データとして蓄積する工程と、前記学習データとして蓄積した特徴ベクトルのおのおのについて前記特徴ベクトルに応じて前記学習データの中から所定数の特徴ベクトルを選択する工程と、前記選択された所定数の特徴ベクトルを用いて学習用の基準ベクトルを作成する工程と、前記作成した学習用の基準ベクトルに基づき前記学習データとして蓄積した特徴ベクトルの異常測度を算出する工程と、前記算出した異常測度に基づき前記モード毎にしきい値を算出する工程と、前記学習データとして蓄積した特徴ベクトルから全組合せの２次元の特徴分布密度を前記モード毎に算出する工程を有し、
    　前記センサ信号の異常診断を行う工程は、前記設備または装置に装着された前記複数のセンサから出力される複数のセンサ信号から特徴ベクトルを抽出して観測ベクトルとする工程と、前記抽出した観測ベクトルに応じて前記蓄積した学習データの中から所定数の特徴ベクトルを選択する工程と、前記学習データの中から選択された所定数の特徴ベクトルを用いて異常診断用の基準ベクトルを作成する工程と、前記作成した異常診断用の基準ベクトルに基づいて前記観測ベクトルの異常測度を算出する工程と、前記算出した異常測度と前記分割したモードと前記分割したモード別に算出した前記しきい値に基づき前記観測ベクトルが異常か正常かを判定する工程と、前記判定により異常と判定された観測ベクトルに対応するセンサ信号が検出された時刻における前記観測ベクトルと前記モードと前記モード別に算出した前記２次元の特徴分布密度に基づき異常に関連するセンサを特定する工程とを有する
    ことを特徴とする異常診断方法。
  12. 　設備または装置に装着されたセンサから出力されるセンサ信号に基づいて前記設備または装置の異常を診断する装置であって、前記センサ信号を蓄積する生データ蓄積手段と、前記センサ信号の特徴値から特徴ベクトルを抽出する特徴ベクトル抽出手段と、予め指定された学習期間において前記特徴値から抽出した前記特徴ベクトルの集合と各時刻において前記特徴値から抽出した前記特徴ベクトルに基づいて前記各時刻における基準特徴ベクトルを算出する基準ベクトル算出手段と、前記各時刻の前記特徴ベクトルと前記基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出手段と、前記予め指定された学習期間における前記異常測度に基づいてしきい値を算出するしきい値算出手段と、前記異常測度と前記しきい値との比較により異常を検出する異常検出手段と、前記予め指定された学習期間における前記特徴値の２次元の分布密度を算出する分布密度算出手段と、前記異常検出手段により検出された異常に対応するセンサ信号が検出された時刻における前記特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定する関連センサ特定手段とを備えたことを特徴とする異常診断装置。
  13. 　設備または装置に装着されたセンサから出力されるセンサ信号およびイベント信号に基づいて前記設備または装置の異常を診断する装置であって、前記センサ信号を蓄積する生データ蓄積手段と、前記イベント信号に基づいて稼動状態別のモード分割を行うモード分割手段と、前記センサ信号の特徴値から特徴ベクトルを抽出する特徴ベクトル抽出手段と、予め指定された学習期間において前記特徴値から抽出した前記特徴ベクトルの集合と各時刻において前記特徴値から抽出した前記特徴ベクトルに基づいて前記各時刻における基準特徴ベクトルを算出する基準ベクトル算出手段と、前記各時刻の前記特徴ベクトルと前記基準特徴ベクトルの差に基づき異常測度を算出する異常測度算出手段と、前記予め指定された学習期間における前記異常測度に基づいて前記モード分割手段で分割したモード毎にしきい値を算出するしきい値算出手段と、前記異常測度と前記算出したモード毎のしきい値との比較により異常を検出する異常検出手段と、前記予め指定された学習期間における前記特徴値の２次元の分布密度を前記モード毎に算出する分布密度算出手段と、前記異常検出手段により検出された異常に対応するセンサ信号が検出された時刻における前記モードに対応する前記特徴値の２次元の分布密度に基づいて異常に関連するセンサを特定する関連センサ特定手段とを備えたことを特徴とする異常診断装置。
PCT/JP2014/050547 2013-01-22 2014-01-15 異常診断方法およびその装置 WO2014115615A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US14/761,748 US9779495B2 (en) 2013-01-22 2014-01-15 Anomaly diagnosis method and apparatus

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2013-009316 2013-01-22
JP2013009316A JP6076751B2 (ja) 2013-01-22 2013-01-22 異常診断方法およびその装置

Publications (1)

Publication Number Publication Date
WO2014115615A1 true WO2014115615A1 (ja) 2014-07-31

Family

ID=51227407

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2014/050547 WO2014115615A1 (ja) 2013-01-22 2014-01-15 異常診断方法およびその装置

Country Status (3)

Country Link
US (1) US9779495B2 (ja)
JP (1) JP6076751B2 (ja)
WO (1) WO2014115615A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017150286A1 (ja) * 2016-02-29 2017-09-08 日本電気株式会社 システム分析装置、システム分析方法、及び、コンピュータ読み取り可能な記録媒体
CN109507990A (zh) * 2018-12-25 2019-03-22 中南大学 一种故障溯源方法及系统
CN111061235A (zh) * 2019-12-20 2020-04-24 中核控制系统工程有限公司 一种具有故障预警功能的dcs设备诊断方法
CN111539482A (zh) * 2020-04-28 2020-08-14 三峡大学 基于rbf核函数的空间多维风电功率数据降维及重构方法
CN113027658A (zh) * 2021-03-24 2021-06-25 华中科技大学 一种水轮机转轮实时状态评估方法及其应用

Families Citing this family (57)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3015757B1 (fr) * 2013-12-23 2019-05-31 Electricite De France Procede d'estimation quantitative du colmatage des plaques d'un generateur de vapeur
JP6223936B2 (ja) * 2014-09-12 2017-11-01 株式会社日立ハイテクノロジーズ 異常傾向検出方法およびシステム
FR3032273B1 (fr) * 2015-01-30 2019-06-21 Safran Aircraft Engines Procede, systeme et programme d'ordinateur pour phase d'apprentissage d'une analyse acoustique ou vibratoire d'une machine
JP5849167B1 (ja) * 2015-04-09 2016-01-27 株式会社日立パワーソリューションズ 異常検知方法およびその装置
KR101647826B1 (ko) * 2015-05-27 2016-08-11 한국에너지기술연구원 전력시스템의 안정도를 판별하는 장치 및 그 전력시스템의 안정도를 관리하는 방법
GB2560643B (en) * 2015-09-30 2021-07-21 Schlumberger Technology Bv Downhole tool analysis using anomaly detection of measurement data
KR102384742B1 (ko) * 2015-10-27 2022-04-07 삼성에스디에스 주식회사 센서 데이터를 이용한 이상 감지 장치 및 방법
US10073421B2 (en) * 2015-11-17 2018-09-11 Rockwell Automation Technologies, Inc. Predictive monitoring and diagnostics systems and methods
GB2554038B (en) * 2016-05-04 2019-05-22 Interactive Coventry Ltd A method for monitoring the operational state of a system
EP3258333A1 (en) * 2016-06-17 2017-12-20 Siemens Aktiengesellschaft Method and system for monitoring sensor data of rotating equipment
US10565513B2 (en) * 2016-09-19 2020-02-18 Applied Materials, Inc. Time-series fault detection, fault classification, and transition analysis using a K-nearest-neighbor and logistic regression approach
US10677943B2 (en) * 2016-12-16 2020-06-09 Smiths Detection, Llc System and method for monitoring a computed tomography imaging system
KR101971553B1 (ko) * 2017-03-21 2019-04-23 (주)심플랫폼 사물인터넷 기반 기기 관리 시스템 및 방법
JP6749488B2 (ja) * 2017-05-26 2020-09-02 三菱電機ビルテクノサービス株式会社 異常重要度算出システム、異常重要度算出装置、及び異常重要度算出プログラム
JP6772963B2 (ja) * 2017-06-05 2020-10-21 トヨタ自動車株式会社 異常診断装置及び異常診断方法
US11669771B2 (en) 2017-07-13 2023-06-06 Nec Corporation Learning system, analysis system, learning method, and storage medium
JP7010641B2 (ja) 2017-09-27 2022-01-26 パナソニック インテレクチュアル プロパティ コーポレーション オブ アメリカ 異常診断方法および異常診断装置
JP7017363B2 (ja) * 2017-10-06 2022-02-08 株式会社日立パワーソリューションズ 異常検知装置および異常検知方法
JP6693938B2 (ja) * 2017-11-17 2020-05-13 ファナック株式会社 外観検査装置
JP6856167B2 (ja) * 2018-02-16 2021-04-07 日本電気株式会社 異常音判定基準作成装置および異常音検知装置
EP3759456A4 (en) * 2018-02-27 2021-12-01 Falkonry, Inc. SYSTEM AND METHOD OF EXPLANATION OF PREDICTIONS IN COMPLEX SYSTEMS
KR101967339B1 (ko) * 2018-03-06 2019-04-09 단국대학교 산학협력단 심층학습 기반의 ａｄａｓ 센서 고장 진단 및 백업을 위한 장치 및 방법
JP7013993B2 (ja) * 2018-03-26 2022-02-01 トヨタ自動車株式会社 診断装置及び診断方法
JP7043320B2 (ja) * 2018-03-30 2022-03-29 株式会社小松製作所 状態分析装置および状態分析方法
US10553046B2 (en) * 2018-04-05 2020-02-04 GM Global Technology Operations LLC Vehicle prognostics and remedial response
CN109002833B (zh) * 2018-06-12 2019-08-27 国家卫生健康委科学技术研究所 一种微液滴数据分析方法及系统
WO2019244203A1 (ja) 2018-06-18 2019-12-26 三菱電機株式会社 診断装置、診断方法及びプログラム
CN109005513B (zh) * 2018-06-26 2021-03-19 北京酷云互动科技有限公司 手机终端关联方法和手机终端关联系统
JP7032261B2 (ja) * 2018-07-24 2022-03-08 日立グローバルライフソリューションズ株式会社 異常検知システムおよび異常検知方法
US11122065B2 (en) * 2018-08-14 2021-09-14 Vmware, Inc. Adaptive anomaly detection for computer systems
CN112567306A (zh) * 2018-08-31 2021-03-26 东芝三菱电机产业系统株式会社 制造过程监视装置
KR101971278B1 (ko) * 2018-12-13 2019-04-26 주식회사 알고리고 인공신경망을 이용한 비정상 데이터 구분 장치
CN113228006A (zh) * 2018-12-17 2021-08-06 华为技术有限公司 检测连续事件中的异常的装置和方法及其计算机程序产品
JP7463055B2 (ja) * 2018-12-27 2024-04-08 株式会社Ihi 異常診断装置、異常診断方法、異常診断プログラム、及び記録媒体
US10922906B2 (en) 2019-03-28 2021-02-16 GM Global Technology Operations LLC Monitoring and diagnosing vehicle system problems using machine learning classifiers
WO2020202567A1 (ja) * 2019-04-05 2020-10-08 株式会社Ｉｈｉ原動機 振動音響解析方法及び装置と機器異常部位推定方法及び装置
KR102119891B1 (ko) * 2019-04-11 2020-06-05 주식회사 알고리고 미세 변화 데이터 및 공간 데이터를 기초로 한 인공신경망을 이용한 비정상 데이터 구분 장치
KR102198224B1 (ko) * 2019-04-11 2021-01-05 주식회사 알고리고 공간 데이터를 기초로 한 인공신경망을 이용한 비정상 데이터 구분 장치
JP7218867B2 (ja) * 2019-05-24 2023-02-07 Ｋｙｂ株式会社 異常検知装置及び異常検知方法
KR20200141812A (ko) 2019-06-11 2020-12-21 삼성전자주식회사 뉴럴 네트워크를 이용하여 이상 신호를 감지하는 방법 및 장치
JP7344015B2 (ja) * 2019-06-13 2023-09-13 株式会社日立ハイテクソリューションズ 異常検知装置及び異常検知方法
DE102019209797A1 (de) * 2019-07-03 2021-01-07 Thyssenkrupp Ag Verfahren und Einrichtung zur Ermittlung des fahrbetriebsbedingten Zustandes von Karosseriekomponenten eines Fahrzeugs sowie eines entsprechenden Fahrerverhaltens
EP4035008A4 (en) * 2019-09-27 2023-10-11 Tata Consultancy Services Limited METHOD AND SYSTEM FOR DIAGNOSING ANOMALIES IN A MANUFACTURING PLANT
CN110866908B (zh) * 2019-11-12 2021-03-26 腾讯科技（深圳）有限公司 图像处理方法、装置、服务器及存储介质
CN110831050B (zh) * 2019-11-21 2022-09-30 武汉宝久智控科技有限公司 一种传感器节点控制方法及系统
DE102019219858A1 (de) * 2019-12-17 2021-06-17 Zf Friedrichshafen Ag Verfahren zur Bestimmung eines Betriebszustands von Systemkomponenten eines Fahrmischers
CN115135358A (zh) * 2020-02-27 2022-09-30 美国西门子医学诊断股份有限公司 使用机器学习的自动传感器追踪验证
US11954309B2 (en) * 2020-05-04 2024-04-09 Adobe Inc. Systems for predicting a terminal event
CN111912638B (zh) * 2020-06-13 2021-12-21 宁波大学 一种在线故障根源识别的精馏塔故障诊断方法
DE102020123721A1 (de) 2020-09-11 2022-03-17 Volkswagen Aktiengesellschaft Zustandsbestimmungsverfahren für eine Fahrzeugkomponente, Schaltung, Kraftfahrzeug
CN112257754B (zh) * 2020-09-24 2023-07-28 北京航天测控技术有限公司 航天器运行状态的分析方法和装置
CN113442691A (zh) * 2021-06-18 2021-09-28 科大讯飞股份有限公司 智能车膜的控制方法、装置和存储介质及电子设备
WO2023002977A1 (ja) * 2021-07-19 2023-01-26 株式会社東芝 情報処理装置、情報処理方法およびプログラム
JP2023042945A (ja) * 2021-09-15 2023-03-28 株式会社東芝 監視装置、方法およびプログラム
CN114228637B (zh) * 2021-12-02 2024-02-20 科大讯飞股份有限公司 一种车辆断电保护方法、装置、存储介质及设备
CN116311374B (zh) * 2023-03-27 2023-10-20 淮阴工学院 一种化工厂工人异常行为识别与预警方法及系统
CN116957543B (zh) * 2023-09-19 2023-12-22 成都秦川物联网科技股份有限公司 基于大数据的智慧燃气设备管理方法和物联网系统

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11510898A (ja) * 1995-08-11 1999-09-21 フィッシャー−ローズマウント システムズ，インコーポレイテッド プロセスに於ける欠陥センサの検出と特定を行うための方法及び装置
JP2004213273A (ja) * 2002-12-27 2004-07-29 Tokyo Gas Co Ltd 故障検出装置及び故障検出方法
JP2012009064A (ja) * 2011-09-05 2012-01-12 Toshiba Corp 学習型プロセス異常診断装置、およびオペレータ判断推測結果収集装置
JP2013143009A (ja) * 2012-01-11 2013-07-22 Hitachi Ltd 設備状態監視方法およびその装置

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6952662B2 (en) 2000-03-30 2005-10-04 Smartsignal Corporation Signal differentiation system using improved non-linear operator
JP5431235B2 (ja) 2009-08-28 2014-03-05 株式会社日立製作所 設備状態監視方法およびその装置
JP5813317B2 (ja) 2010-12-28 2015-11-17 株式会社東芝 プロセス状態監視装置

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11510898A (ja) * 1995-08-11 1999-09-21 フィッシャー−ローズマウント システムズ，インコーポレイテッド プロセスに於ける欠陥センサの検出と特定を行うための方法及び装置
JP2004213273A (ja) * 2002-12-27 2004-07-29 Tokyo Gas Co Ltd 故障検出装置及び故障検出方法
JP2012009064A (ja) * 2011-09-05 2012-01-12 Toshiba Corp 学習型プロセス異常診断装置、およびオペレータ判断推測結果収集装置
JP2013143009A (ja) * 2012-01-11 2013-07-22 Hitachi Ltd 設備状態監視方法およびその装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017150286A1 (ja) * 2016-02-29 2017-09-08 日本電気株式会社 システム分析装置、システム分析方法、及び、コンピュータ読み取り可能な記録媒体
JPWO2017150286A1 (ja) * 2016-02-29 2018-10-25 日本電気株式会社 システム分析装置、システム分析方法、及び、プログラム
CN109507990A (zh) * 2018-12-25 2019-03-22 中南大学 一种故障溯源方法及系统
CN109507990B (zh) * 2018-12-25 2021-06-15 中南大学 一种故障溯源方法及系统
CN111061235A (zh) * 2019-12-20 2020-04-24 中核控制系统工程有限公司 一种具有故障预警功能的dcs设备诊断方法
CN111539482A (zh) * 2020-04-28 2020-08-14 三峡大学 基于rbf核函数的空间多维风电功率数据降维及重构方法
CN113027658A (zh) * 2021-03-24 2021-06-25 华中科技大学 一种水轮机转轮实时状态评估方法及其应用

Also Published As

Publication number Publication date
JP2014142697A (ja) 2014-08-07
JP6076751B2 (ja) 2017-02-08
US20150363925A1 (en) 2015-12-17
US9779495B2 (en) 2017-10-03

Similar Documents

Publication Publication Date Title
JP6076751B2 (ja) 異常診断方法およびその装置
JP6216242B2 (ja) 異常検知方法およびその装置
EP2905665B1 (en) Information processing apparatus, diagnosis method, and program
JP5301717B1 (ja) 設備状態監視方法およびその装置
JP5945350B2 (ja) 設備状態監視方法およびその装置
JP5342708B1 (ja) 異常検知方法及びその装置
JP5431235B2 (ja) 設備状態監視方法およびその装置
JP7017363B2 (ja) 異常検知装置および異常検知方法
JP5331774B2 (ja) 設備状態監視方法およびその装置並びに設備状態監視用プログラム
JP5538597B2 (ja) 異常検知方法及び異常検知システム
JP5501903B2 (ja) 異常検知方法及びそのシステム
WO2013011745A1 (ja) 設備状態監視方法およびその装置
JP6223936B2 (ja) 異常傾向検出方法およびシステム
JP2013143009A (ja) 設備状態監視方法およびその装置
JP5849167B1 (ja) 異常検知方法およびその装置
JP6831729B2 (ja) 異常検知装置
JP2014170269A (ja) 時系列データの異常監視装置、異常監視方法及びプログラム
JP2023106472A (ja) 異常検知装置及び異常検知方法
JP6213655B2 (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: 14743264

Country of ref document: EP

Kind code of ref document: A1

WWE Wipo information: entry into national phase

Ref document number: 14761748

Country of ref document: US

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 14743264

Country of ref document: EP

Kind code of ref document: A1