WO2024029387A1 - 時系列データ解析システム - Google Patents

時系列データ解析システム Download PDF

Info

Publication number
WO2024029387A1
WO2024029387A1 PCT/JP2023/027035 JP2023027035W WO2024029387A1 WO 2024029387 A1 WO2024029387 A1 WO 2024029387A1 JP 2023027035 W JP2023027035 W JP 2023027035W WO 2024029387 A1 WO2024029387 A1 WO 2024029387A1
Authority
WO
WIPO (PCT)
Prior art keywords
time
series data
time series
point
processing unit
Prior art date
Application number
PCT/JP2023/027035
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 Ｂｉｐｒｏｇｙ株式会社
Publication of WO2024029387A1 publication Critical patent/WO2024029387A1/ja

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D1/00—Measuring arrangements giving results other than momentary value of variable, of general application
    • G01D1/10—Measuring arrangements giving results other than momentary value of variable, of general application giving differentiated values
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D3/00—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups
    • G01D3/028—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups mitigating undesired influences, e.g. temperature, pressure
    • G01D3/032—Indicating or recording apparatus with provision for the special purposes referred to in the subgroups mitigating undesired influences, e.g. temperature, pressure affecting incoming signal, e.g. by averaging; gating undesired signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01K—MEASURING TEMPERATURE; MEASURING QUANTITY OF HEAT; THERMALLY-SENSITIVE ELEMENTS NOT OTHERWISE PROVIDED FOR
    • G01K3/00—Thermometers giving results other than momentary value of temperature
    • G01K3/08—Thermometers giving results other than momentary value of temperature giving differences of values; giving differentiated values
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations

Definitions

  • the present invention relates to a time-series data analysis system and a time-series data analysis method for analyzing time-series data output from various measuring devices, various sensors, and the like.
  • an object of the present invention is to provide a data analysis system and a time-series data analysis method that make it possible to more effectively analyze data even when the data contains a lot of noise.
  • a time differentiation processing means for performing time differentiation on time series data and displaying a waveform diagram of the time series data after time differentiation as necessary; a range setting processing section that sets a time range on a time axis when the time differentiation processing section executes time differentiation processing on time series data to be processed; a peak position detection unit that detects a position where the time series data to be processed has a peak in the time range set by the range setting processing unit; (first invention).
  • the range setting processing unit allows the user to directly set at least one boundary of the time range on the time axis of time series data (second invention).
  • the range setting processing unit is capable of allowing the user to set a specific value on the vertical axis, and to set the corresponding value on the time axis as at least one boundary of the time range. 3 inventions).
  • the range setting processing unit searches for a peak value on the vertical axis of the time series data detected by the peak position detection unit, and calculates a value on the time axis corresponding to the peak value at the time. It can be one boundary of the range (fourth invention).
  • the moving average processing section further includes, as a noise removal means, a moving average processing unit that executes moving average processing on the time series data and displays a waveform diagram of the time series data after the moving average processing, if necessary; and an integral difference processing unit that executes integral difference processing on the time series data and displays a waveform diagram of the time series data after the integral difference processing as necessary. It is possible (fifth invention).
  • a time differentiation processing means for performing time differentiation on time series data and displaying a waveform diagram of the time series data after time differentiation as necessary;
  • the first pulse waveform has a pulse width that is the length of time that the value is continuously positive, a pulse height that is the length of time that the value is continuously positive, and a first pulse waveform that has a continuous positive value.
  • a pulse generation unit that displays a waveform diagram of both or one of the pulse waveforms
  • a subtraction processing unit that executes subtraction processing between a plurality of time series data or between a plurality of pulse waveforms, and displays the time series data or pulse waveform after the subtraction as necessary
  • a threshold processing unit that compares the value of the pulse waveform with a prespecified threshold and selects a pulse larger or smaller than the threshold
  • the moving average processing section further includes, as a noise removal means, a moving average processing unit that executes moving average processing on the time series data and displays a waveform diagram of the time series data after the moving average processing, if necessary; and an integral difference processing unit that executes integral difference processing on the time series data and displays a waveform diagram of the time series data after the integral difference processing as necessary. It is possible (10th invention).
  • the ninth invention described above may further include a user input processing section through which the user inputs user data (twelfth invention).
  • the method for finding the inflection point of time series data imported into a time series data analysis system is as follows.
  • the pulse generation unit included in the time series data analysis system defines the length of time during which the second time series data is continuously positive as a pulse width, and the length of time during which the second time series data is continuously positive as a pulse height.
  • a first pulse waveform in which the length of time during which the fourth time series data is continuously negative is the pulse width, and the length of time during which the fourth time series data is continuously negative is the pulse height.
  • the third time differentiation is performed after noise is removed by performing integral difference on the time series data obtained by performing the second time differentiation. This can be used as second time series data (fifteenth invention).
  • the time series data analysis method in the present invention includes: A time series data analysis method for analyzing time series data by obtaining an approximate curve of a waveform diagram of the time series data, the method comprising: provisionally arranging from a first dividing point to an nth dividing point (n is an integer of 2 or more) on the time series data between a prespecified starting point and an ending point; When the first dividing point is moved from the starting point to the second dividing point, the point where the first straight line connecting the starting point and the first dividing point best fits the time series data in that range is determined as the position of the first dividing point.
  • FIG. 1 is a diagram illustrating functional blocks of a time-series data analysis system according to an embodiment of the present invention.
  • This is time-series data of the output of a pressure sensor that measures changes in pressure from when resin is poured until it solidifies during injection molding.
  • 3 is a waveform diagram in which a part of the waveform diagram of FIG. 2 is expanded in the time direction.
  • FIG. 3 is a waveform diagram showing an example of time series data obtained by time-differentiating time series data.
  • FIG. 3 is a waveform diagram showing an example of time series data obtained by taking a moving average of time series data.
  • FIG. 3 is a diagram for explaining an integral difference with respect to time series data.
  • FIG. 3 is a diagram for explaining an integral difference with respect to time series data.
  • FIG. 3 is a diagram for explaining an integral difference with respect to time series data.
  • FIG. 3 is a diagram illustrating a problem when determining the time point when passing through zero by threshold processing.
  • FIG. 4 is a waveform diagram showing the operation of a pulse generation section that generates pulses based on a waveform diagram obtained by differentiation.
  • FIG. 7 is a waveform diagram showing an example of performing various processes on time series data to obtain an inflection point of the original time series data.
  • 10 is a waveform diagram illustrating that threshold processing is performed on the refraction points obtained in FIG. 9 to select pulses that are equal to or greater than a predetermined threshold;
  • FIG. FIG. 3 is a waveform diagram showing a method for correctly determining a peak point from a waveform diagram containing a lot of noise.
  • FIG. 6 is a waveform diagram showing an example of a method for finding one refraction point of an original waveform diagram.
  • FIG. 6 is a waveform diagram showing an example of a method for finding one refraction point of an original waveform diagram.
  • FIG. 2 is a diagram summarizing processes that can be executed by a system according to an embodiment of the present invention.
  • FIG. 1 is a diagram showing a user interface of a system according to an embodiment of the present invention.
  • FIG. 2 is a diagram illustrating a problem when difference processing is performed on the entire time series data.
  • FIG. 3 is a diagram illustrating that noise can be reduced by increasing the degree of moving average.
  • FIG. 3 is a diagram illustrating problems caused by increasing the number of moving average points.
  • FIG. 1 is a diagram showing a user interface of a system according to an embodiment of the present invention.
  • FIG. 2 is a diagram illustrating a problem when difference processing is performed on the entire time series data.
  • FIG. 3 is a diagram
  • FIG. 6 is a diagram illustrating a case where an operator manually sets a time range.
  • FIG. 3 is a diagram illustrating a method of determining a range by focusing on a value range on the vertical axis.
  • FIG. 6 is a diagram illustrating a method of setting a search range based on characteristics of feature points.
  • FIG. 3 is a diagram illustrating a method of setting a range based on the maximum point of speed data.
  • FIG. 3 is a diagram illustrating a process of setting a processing range based on the mutual relationship of a plurality of time series data.
  • FIG. 3 is a diagram illustrating a specific method of actually setting a processing range for time-series data.
  • FIG. 3 is a diagram illustrating a specific method of actually setting a processing range for time-series data.
  • FIG. 3 is a diagram illustrating a specific method of actually setting a processing range for time-series data.
  • FIG. 3 is a diagram illustrating a specific method of actually setting a processing range for time-series data.
  • FIG. 3 is a diagram for explaining a method for obtaining an approximate curve of original time series data.
  • FIG. 3 is a diagram for explaining a method for obtaining an approximate curve of original time series data.
  • FIG. 3 is a diagram for explaining a method for obtaining an approximate curve of original time series data.
  • FIG. 1 is a functional block diagram of the main parts of a time-series data analysis system according to an embodiment of the present invention.
  • Each processing unit shown in FIG. 1 includes a processor in a general computer including a personal computer, a temporary storage device such as a DRAM, a non-temporary storage device such as an HDD device, a keyboard and mouse as data input means, a display device, etc.
  • Various types of processing are performed by a combination of hardware and software such as a computer program read from a storage device and executed by a processor.
  • a data input processing unit 10 performs a process of importing time-series data from measurement devices, various sensors, etc. into the system.
  • time-series data refers to data that is captured at regular time intervals, or data that has undergone some processing on time-series data that has already been captured into the system.
  • the user input processing unit 12 performs a process of importing data input by a user from a keyboard, mouse, etc. into the system.
  • the image output processing unit 14 performs processing for outputting a user interface including characters, graphics, etc., actual time-series waveform data, etc. to a display or, if necessary, to an output device such as a printer.
  • the band-pass filter processing unit 16 performs band-pass filter processing on time-series data that has been imported into the system, or on time-series data that has already been imported into the system and has undergone some processing, and performs band-pass filter processing as necessary. Display the waveform diagram obtained. Note that the band-pass filter processing unit 16 not only directly performs band-pass filter processing on time-series data, but also performs fast Fourier transform as necessary, performs band-pass filter processing in frequency space, and then performs inverse processing. It is also possible to perform fast Fourier transform to return to time and space.
  • the moving average processing unit 18 averages the time series data taken into the system, or the time series data that has already been taken into the system and has undergone some processing, at a predetermined number of adjacent time positions. Execute processing. For this predetermined number, the default setting may be used as is, or a value set by the user via a user interface to be described later may be used.
  • the differential processing unit 20 performs a time differential operation on the time series data taken into the system or on the time series data that has already undergone some processing in the system, and generates the obtained waveform as necessary. Show diagram. For example, if the time interval (reciprocal of the sampling frequency) for capturing time series data is 1/1000th of a second, then the difference between two adjacent data values at this 1/1000th of a second interval is expressed as 1/1000th of a second. By dividing, we get the time derivative. However, it is also possible to execute the differentiation process by widening the time interval, such as every other or every second time interval, for each time series. Taking time differentiation at finite time intervals in this way is also called difference.
  • the integral difference processing unit 22 performs an integral operation on the time series data taken into the system or on the time series data that has already undergone some processing in the system, and further performs an integral operation on the obtained integral data. On the other hand, difference processing is executed, and the obtained waveform diagram is displayed as necessary.
  • the subtraction processing unit 24 executes subtraction processing between two certain waveform data, and displays the obtained waveform diagram as necessary.
  • the pulse generation unit 26 generates a first pulse waveform for the time series data, in which the length of time during which the value is continuously positive is the pulse width, and the length of time during which the value is continuously positive is the pulse height. and a second pulse waveform whose pulse width is the length of time during which the value is continuously negative and whose pulse height is the length of time during which the value is continuously negative. A waveform diagram of both or one of the pulse waveform and the second pulse waveform is displayed.
  • the threshold processing unit 28 executes a process of selecting a pulse that exceeds or is below the threshold from among a plurality of pulses by setting a threshold by the user, and selects the selected pulse as necessary. indicate.
  • the division point search processing unit 30 approximates a range from a certain start point to an end point of a certain time series data with a curve made of straight lines connecting several points (dividing points) included in the range, Processing is performed to determine the positions of these division points so that a sufficient approximation can be made to judge the characteristics of the time series data.
  • the processing range setting processing unit 32 does not process the entire time series data that has been imported into the system or has already undergone some processing within the system, but rather sets a time range based on certain rules. , and execute processing that targets only time-series data within the set time range. The specific contents of the processing performed by the processing range setting processing section 32 will be described in detail later.
  • the moving average processing section 18 and the integral difference processing section 22 shown in FIG. 1 perform noise removal processing, and are not necessarily required as part of the system configuration. Either one of them may be used as necessary. Systems can be configured to include one or both. Furthermore, it is not the case that the present invention cannot be realized unless all the processing units shown in FIG. 1 are included. For example, if only the differential processing section 20, range setting processing section 32, and peak position detection section 34 are included, it is possible to construct a system that performs the minimum necessary processing on time series data, and the differential processing section 20, Even if only the pulse generation section 26, the subtraction processing section 24, and the threshold processing section 28 are included, it is possible to construct a system that performs the minimum necessary processing on time series data. After that, the system can be constructed by adding other processing units shown in FIG. 1 according to actual processing needs.
  • Figure 2 shows how the pressure applied during injection molding (vertical axis) changes over time (horizontal axis) from when the resin is poured into the mold until it solidifies, and shows an actual pressure sensor. Shows output time series data. Incorporation of such data into the time series data analysis system is performed by the data input processing section 10 in FIG.
  • circles 10a, 10b, 10c, and 10d are shown at four locations superimposed on the signal waveform.
  • the leftmost circle 10a is the point at which the pressure starts to rapidly increase from almost zero, and shows the moment just before the resin spreads throughout the inside of the mold. Then, the pressure increases rapidly from point 10a and reaches the peak of circle 10b. This peak corresponds to the point when the entire mold is filled with resin and resin pouring is stopped.
  • the pressure rapidly decreases, and the pressure remains constant for a while from the point of circle 10c.
  • This constant pressure period corresponds to a period during which the temperature of the resin is cooled. Then, when the resin is sufficiently cooled and solidified, the pressure rapidly decreases again from the point of circle 10d.
  • the differential processing unit 20 performs time differentiation one more time on the curve in FIG. 4(a), which is the same as FIG. 3(c), to obtain the curve in FIG. 4(b). If we know the point at which the curve in Fig. 4(b) passes through zero (indicated by an upward arrow), this point will be the maximum point of the curve in Fig. 4(a), that is, the point at which the resin has spread throughout the inside of the mold. It can be determined that there is. However, the curve shown in Figure 3(a) is only a curve that assumes the change in pressure when the resin is spread throughout the inside of the mold, and the signal waveform of the actual pressure sensor does not follow such a smooth graph. Instead, it is output with considerable noise added to it.
  • FIG. 5 shows the problem when trying to find a local maximum point from time series data that contains a lot of noise.
  • FIG. 5(a) shows an example of a signal waveform corresponding to FIG. 4(a) for which the maximum point is to be determined, and it contains considerable noise when compared with the curve in FIG. 4(a). . Therefore, first, the moving average processing unit 18 of FIG. 1 takes a moving average.
  • the waveform in FIG. 5(b) is a 49-point moving average of the signal waveform in FIG. 5(a), and the waveform in FIG. 5(c) is a 99-point moving average.
  • the waveform in FIG. 5(d) is obtained by taking a 199-point moving average.
  • the multiple-point moving average corresponds to step difference, that is, performing differentiation, but in general, the more points the moving average has, the more the noise decreases.
  • the waveforms in (b), (c), and (d) in Figure 5 are obtained by taking moving averages at different points, and the waveforms in (b), (c), and (d) each pass through zero.
  • the time points t b , t c , and t d do not coincide, but are slightly shifted back and forth. This shows that the point at which the obtained waveform passes through zero changes depending on how the data is processed, that is, in this case, depending on how many points the moving average process is performed on. Therefore, with the method shown in FIG. 4 of taking the differential (difference) with respect to acceleration, it is not possible to accurately determine the point in time when the trend changes from the output signal containing noise.
  • a process called integral difference is executed. That is, for a waveform whose point in time passing through zero is desired to be determined, it is integrated at certain equal intervals, and the difference is calculated with respect to this integrated value. in particular, Perform the calculation.
  • the above integral difference equation will be explained with reference to FIGS. 6A and 6B.
  • the point of interest is i and the range of integration is w
  • the right side of the above equation The first term corresponds to the integrated value of waveform a from i to i+w, that is, the area under the curve of a from i to i+w (indicated by hatching).
  • the second term on the right side corresponds to the sum of the values of waveform a from iw to i, that is, the area under the curve from iw to i (shown in pear pattern).
  • the integral difference represents the difference between the area to the right of i and the area to the left of i.
  • w is 99
  • w is 199
  • 199-point integral difference it is called a 199-point integral difference.
  • the waveform in (b) shows a waveform diagram obtained by performing a 99-point integral difference on the waveform a in (a).
  • 6A corresponds to the case where the value of the integral difference is zero
  • FIG. 6B corresponds to the case where the difference between the first term and the second term on the right side is almost minimum.
  • the noise is reduced to the extent that it does not interfere with identifying the local maximum point, that is, the point in time when passing zero, without significantly losing the characteristics of the original waveform. You can see that it has been removed. Looking at the waveform in (a) of FIG.
  • the waveform shown in FIG. 6A (b) is actually time series data plotted against discrete times (horizontal axis) and connected by curves.
  • the upper waveform diagram corresponds to FIG. 6A (b)
  • the lower waveform diagram is obtained by time-differentiating the upper waveform diagram. Therefore, the point in time when the lower waveform diagram becomes zero is the point in time when the upper waveform diagram passes through zero.
  • the lower waveform diagram is also plotted against discrete times, similar to the upper waveform diagram. Since the time resolution is thus limited, it is rare for the time series data shown at the bottom of FIGS. 7(a) and 7(b) to contain complete zeros.
  • a threshold value of a constant width centered on zero on the vertical axis It is conceivable to set a band-shaped portion with a constant width centered on the time axis, and determine that data included in this range is zero.
  • data within the threshold range is shown by lines parallel to the vertical axis for ease of understanding. In this case, if the threshold width is narrow as in (a), there may be no data that falls within the range where the curve changes significantly, and if the threshold width is wide as in (b), there may be no corresponding data.
  • Another method is to search the time series data point by point and find the point at which the sign changes from plus to minus or from minus to plus. However, this method is inefficient because the processing is scalar processing.
  • FIG. 8(a) corresponds to the upper waveform diagram in FIGS. 7(a) and (b)
  • FIG. 8(b) corresponds to the lower waveform diagram in FIGS. 7(a) and (b).
  • FIG. 8(b) is a waveform obtained by time-differentiating the waveform of FIG. 8(a).
  • the pulse generating section 26 shown in FIG. 1 generates pulse waveforms as shown in FIGS. 8(c) and 8(d).
  • the pulse waveform of FIG. 8(c) is derived from the waveform of FIG.
  • FIG. 9 FIGS. 9(a) to 9(d) correspond to FIGS. 8(a) to 8(d)
  • the subtraction processing unit 24 of FIG. By subtracting and taking the difference in Fig. 9(c), the waveform diagram in Fig. 9(e) is obtained, and by further performing time differentiation on the waveform diagram in Fig. 9(e), Fig. 9(f) is obtained.
  • a waveform diagram as shown in which pulses remain only at the time of switching from positive to negative and the time of switching from negative to positive is obtained. Depending on whether this pulse is positive or negative, it can be determined whether the pulse has switched from positive to negative or from negative to positive.
  • the waveform diagram of FIG. 10(h) is obtained ((a) of FIG. 10 shows the pressure inside the mold detected. (b) to (g) correspond to (a) to (f) in FIG. 9).
  • this time point becomes the inflection point of the time series data in FIG. 10(a), i.e. It can be seen that the trend has changed.
  • FIG. 11 (a) is the output data of a pressure sensor that detects the actual pressure inside the mold, and the maximum pressure point is determined from this data that includes considerable noise.
  • the horizontal axis is time.
  • a fast Fourier transform (FFT) is performed on the output data in FIG. This is the waveform shown in FIG. 11(b). Note that in FIG. 11(b), "0-40 Hz” means that waveforms with frequencies higher than 40 Hz are cut and only waveforms with frequencies between 0 and 40 Hz are extracted.
  • the waveform shown in FIG. 11(c) is obtained by taking the moving average of nine points for the waveform in FIG. 11(b).
  • the maximum point of the waveform in FIG. 11(c) can be considered to be the point of maximum pressure.
  • the actual peak position in FIG. 11(a) is shifted to the left side from the peak position in FIG. 11(c). .
  • This is considered to be due to the influence of noise.
  • the influence of noise can be eliminated and the pressure can be maximized. You can find the point in time when The portion that performs the processing from (a) to (b) to (c) in FIG. 11 functionally corresponds to the peak position detection unit 34 shown in FIG. 1.
  • FIG. 12 shows a state in which the output waveform of the pressure sensor shown in FIG. 2 is expanded in the direction of the time axis (horizontal axis).
  • the waveform shown in FIG. 12(a) is subjected to band-pass filtering similar to that shown in FIG. 11(b) by the band-pass filter processing unit 16 shown in FIG. ) is the waveform shown.
  • the waveform shown in (c) is obtained by taking a nine-point moving average of the waveform in (b) by the moving average processing unit 18 in FIG. Furthermore, the waveform shown in (d) is obtained by time-differentiating the waveform in (c) by the differentiation processing unit 20, and the waveform (e) obtained by further performing time differentiation on this is obtained. ) is the waveform shown. In the waveform shown in (e), a peak can be seen slightly before the value of 600 on the horizontal axis, but if you compare it with the waveform in (a), this point is the start of a stable period. I understand.
  • FIG. 13 shows the output waveform of the pressure sensor shown in FIG. 2 extended in the direction of the time axis (horizontal axis). This corresponds to the point at which the pressure decreases as it cools and solidifies.
  • the waveform in FIG. 13(a) is subjected to band-pass filter processing similar to that in FIG. 11(b) by the band-pass filter processing unit 16 in FIG. 1, and the result is shown in FIG. 13(b).
  • the waveform is shown in
  • the waveform shown in FIG. 13(c) is obtained by taking a nine-point moving average of the waveform shown in FIG. 13(b) by the moving average processing unit 18 of FIG. Furthermore, the waveform shown in (d) is obtained by subjecting the waveform in (c) to time differentiation by the differential processing section 20.
  • the waveform shown in FIG. 13(d) it can be seen that the point in time when the value temporarily changes from a certain approximately constant positive value to a negative value corresponds to the portion shown by the circle in FIG. 13(a).
  • the waveform (e) obtained by further time-differentiating the waveform (d) that is, the waveform corresponding to acceleration, has a gradual downward movement, so this point can be successfully detected. It's difficult to do.
  • FIG. 14 is a diagram summarizing what can be done with the system of this embodiment.
  • band pass filter processing moving average processing
  • difference integration processing difference integration processing
  • differential processing differential processing
  • FIG. 15 shows an example of a user interface of the time series data analysis system shown in FIG. 1. Utilizing such a user interface improves convenience when the user actually performs processing on time-series data.
  • the above-described processes can be executed on arbitrary time-series data, and the results of the processes can be visually viewed immediately.
  • the "Data source" field at the top left of this screen lists the name of the time-series data to be processed, and the "Raw data” field immediately to the right of this field lists the data to be processed, such as the output data of pressure sensors and weight sensors.
  • the actual time series data input to the target computer will be displayed.
  • an "FFT spectrum” column is provided on the right side, and an FFT spectrum obtained by performing fast Fourier transform (FFT) on the waveform diagram displayed in the "Raw data” column is displayed here.
  • FFT fast Fourier transform
  • This FFT processing can be performed on all the input data displayed in the "Raw data” column, or the user can select whether or not to perform the FFT processing.
  • bandpass design is possible in this "FFT spectrum” column. Looking at the waveform in the "FFT spectrum” column of FIG. 15, it can be seen that the peaks are concentrated around 50 Hz on the horizontal axis. Therefore, 25 Hz is specified here as the upper limit of the bandpass (gray shading in the frequency portion higher than 25 Hz indicates this), and only waveforms below 25 Hz are extracted.
  • a column for selecting a sampling rate is provided below the "data source” column.
  • This sampling data corresponds to the time interval of adjacent individual data of the time series data, and here, "100Hz (1/100 second interval)", “1000Hz (1/1000 second interval)", It is possible to select one of "100 kHz (1/100,000 second interval)” (1000 Hz is selected in FIG. 15).
  • Bandpass Lower Limit (Hz)” and “Bandpass Upper Limit (Hz)” the upper and lower limits of bandpass processing can be graphically set on the screen.
  • the “Moving average range” column below is for setting the range for taking the moving average, that is, the number of consecutive time series data for taking the average, and can be set in the range from 1 to 19.
  • the “center” column below is a switch for selecting the center of the moving average.
  • the “refraction point detection” column below is set by checking or not checking the check column to determine whether or not to perform refraction point detection.
  • the "Peak Direction” column below is a column for selecting to detect either or both of an upward peak and a downward peak when detecting a peak.
  • the "Peak Width” field is a field where, when detecting a peak, it is set so that only peaks with a width equal to or less than the width set here are detected.
  • Peak threshold ratio to maximum peak
  • an inflection point where time-series data shows a characteristic change can be found as a peak in displacement, velocity, or acceleration.
  • the actual method of determining velocity from displacement and acceleration from velocity is based on time differentiation, that is, taking the step difference.
  • time differentiation that is, taking the step difference.
  • the process of calculating the difference is performed on the entire input time series data, the long-term trend in the time series data will be lost. Loss of long-term trends emphasizes short-term fluctuations, that is, noise, and hinders the extraction of feature points.
  • FIG. 16 is a diagram explaining this, in which (a) is time series data showing changes in a certain load, (b) is velocity data obtained by time differentiating the time series data in (a), ( c) shows acceleration data obtained by time-differentiating the velocity data in (b).
  • FIG. 16(c) when the process of calculating the time differential for the entire time series data in (a) is performed continuously, the noise is emphasized and the trend indicated by the circle in (a) is It becomes difficult to find changed feature points.
  • FIG. 17 shows that noise can be reduced by increasing the degree of moving average.
  • Figure 17 (a) shows a 5-point moving average of velocity and acceleration, (b) a 9-point moving average, and (c) a 19-point moving average. ing.
  • the upward arrow shown at the bottom of the figure indicates the point at which it was determined to be a minutiae, but in (a) it is not the minutiae point where the load begins to decrease, but it is shifted to the right of that point.
  • Time points are extracted as feature points.
  • the noise is reduced compared to (a), and the feature points can be detected, but there is a high possibility that they were detected by chance, and there may be a problem with robustness.
  • the noise has been considerably removed, and the noise has been further reduced and the feature has been detected.
  • the solid line is the actual curve showing the load fluctuation
  • the dotted line is the curve obtained by taking a 19-point moving average of the original load data
  • the dashed-dotted line is the curve obtained by moving 39 points relative to the original load data.
  • Each averaged curve is shown.
  • the curve will deviate significantly to the left (earlier side of time) from the actual rise. In this way, when the number of moving average points is increased, the waveform changes significantly from the original curve, and it becomes difficult to obtain accurate feature points simply by increasing the number of moving average points.
  • the top graph in Figure 19 shows how the pressure applied during injection molding (vertical axis) changes over time (horizontal axis) from when the resin is poured into the mold until it solidifies. It shows a graph of time series data of the output of the pressure sensor, and the second graph is a graph of time series data obtained by performing one differentiation on the top time series data to find the velocity.
  • the bottom graph shows time-series data obtained by differentiating the acceleration once more.
  • the value on the vertical axis rises rapidly when the value on the horizontal axis is around 250 (indicated by t 1 ), and then reaches a peak and declines just after 300 (indicated by t 2 ).
  • FIG. 20 shows a method of determining the range by focusing on the value range on the vertical axis instead of the time axis on the horizontal axis.
  • the graph shown in Figure 20 is the same as the graph shown in Figure 19, but in the speed graph shown second, a range in which the vertical axis value is less than -0.1 is set, and this setting range is This instructs the system to execute the processes shown in FIGS. 3 to 13 only. By doing this, it is possible to exclude from the feature point search process a portion where the acceleration changes significantly, where the horizontal axis is around 300. In this way, by specifying a threshold based on the value of the vertical axis and specifying a range where the original value is smaller or larger than the threshold, it is possible to reliably find the feature point that you really want to find.
  • FIG. 21 is the same graph as FIGS. 19 and 20, and here, an attempt is made to find the time point at which the acceleration is maximum before and after the time point corresponding to t 3 in FIG. 19. Therefore, the system is instructed to execute the processes shown in FIGS. 3 to 13, limited to a certain range (in this case, a range of 150 points) before and after time t3 .
  • a certain range in this case, a range of 150 points
  • the efficiency of the feature point extraction process can be improved, and accuracy can also be improved.
  • FIG. 22 shows a search range between the minimum acceleration point and the maximum speed point in the graphs showing the temporal changes in the sensors shown in FIGS. 19 to 21, and the point where the acceleration is maximum within this range is searched. This is an example of instructing the system to execute the processes shown in FIGS. 3 to 13.
  • the upper graph in FIG. 23 shows three time series data captured simultaneously. These are graphs showing how the output from pressure sensors provided at three locations on an injection mold changes over time before and after resin is applied to the mold. Three channels of data are generated simultaneously from three pressure sensors and input into the system. As can be seen from the graph shown at the top of FIG. 23, even within the same mold, the pressure changes are not the same depending on the location, but the overall tendency of the changes is similar. Setting the range based on the relationship between the data from these three sensors will be explained.
  • FIG. 24 shows a graph of the time series data shown in FIG. 19 at the bottom, and an example of a DSL for data input for the operator to specify a processing range for the time series data below at the top.
  • a gray band is displayed in the time series data in the range from 350 to 450 on the horizontal axis (time axis), and this range is designated as the processing range.
  • FIG. 25 shows a specific method for setting a processing range using DSL for the time series data shown in FIG. 20, focusing on the range on the vertical axis.
  • “ROLLING_MEAN 5" indicates that the moving average score is "5"
  • FIG. 26 shows how to set a processing range using DSL to search for the maximum acceleration point for the time series data shown in FIG. 22.
  • HORIZONTAL_LIMIT [IDXMAX(DST), IDXMAX(DST)+150]
  • DST range of the original waveform
  • TARGET IDXMAX(ACC)
  • FIG. 27 shows the waveform of output data from a sensor called a strain gauge used during press working.
  • the division point search processing unit 30 of the time series data analysis system shown in FIG. 1 executes the following processing.
  • FIG. 28 Although the waveform shown in FIG. 28 is not the same as the waveform shown in FIG. 27, the general tendency is the same, so the division point search process will be explained using FIG. 28.
  • three division points p 1 , p 2 , and p 3 are temporarily placed on the waveform data between two peaks, the starting point x 1 and the ending point x 2 , and the points are connected by a straight line.
  • the tentative positions of the division points p 1 , p 2 , and p 3 can be arranged evenly between x 1 and x 2 by the division point search processing unit 24, but other temporary locations can be arranged by the user. can also be selected and set.
  • FIG. 29 is a diagram for explaining how the positions of the three division points p 1 , p 2 , and p 3 are finally determined.
  • x 1 and the tentatively determined position of p 2 are fixed.
  • x 1 is located around 210
  • p 2 is located around 240.
  • p 1 is moved at predetermined intervals, for example, by 1 on the scale of the horizontal axis, and at each point, the straight line a connecting x 1 and p 1 best fits the waveform of the output data in the corresponding range. Determine the position of p 1 as a point.
  • the best-fitting point is determined, for example, as the point p 1 where the regression error is minimized by performing a regression calculation to approximate the straight line a and the output data waveform by the method of least squares.
  • the processing up to this point requires calculations at 30 points from 210 to 240 on the horizontal axis, but since there is only one division point p1 , the loop processing only needs to be done once, and the load on the computer is not that large. .
  • p 2 is moved by 1 between the determined p 1 and the tentatively determined p 3 , and a straight line b connecting p 1 and p 2 at each point.
  • the position of p 2 is determined as the point that best fits the waveform of the output data in the corresponding range.
  • p 3 is moved by 1 between the determined p 2 and x 2 , and the straight line c connecting p 2 and p 3 at each point is The position of p 3 is determined as the point that best fits the waveform of the output data.
  • the approximate curve obtained by straight lines a, b, c, and d connecting p 1 , p 2 , and p 3 determined as above is a curve that is approximated quite well compared to the waveform of the output data. . Thereby, it is possible to determine whether the press working was successful based on the approximate curve obtained by the straight lines a, b, c, and d. Moreover, since the processing described above only moves one point, the calculation load is not so high and it can be sufficiently executed by a normal personal computer.
  • a pressure sensor is provided inside the mold, and the pressure sensor detects the internal pressure during actual injection molding, and the waveform of the output data can be obtained.
  • the pressure inside the mold increases when liquid resin is injected, but as the resin injection ends and the resin inside cools and hardens, the pressure gradually decreases, and eventually the molded product leaves the mold. At the stage of extraction, the internal pressure decreases rapidly. Whether or not the injection molding process has been carried out properly can be confirmed by looking at the waveform of output data from such an output sensor. In such a case, the system of this embodiment can be utilized.
  • the number of dividing points was three points p 1 , p 2 , and p 3 , but if the characteristics of such pressure changes are known in advance, the number of dividing points can be increased or decreased based on that. The number can be adjusted to suit the application.
  • this embodiment can be applied to any data as long as it is time-series data. For example, since the output data of the pressure sensor shown in Fig. 2 is also time series data, a fitting curve is obtained using an appropriate number of dividing points in the second embodiment, and the steel plate cutting operation is performed appropriately based on the fitting curve. The system of this embodiment can be used to determine whether or not a request has been made.

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Automation & Control Theory (AREA)
  • Indication And Recording Devices For Special Purposes And Tariff Metering Devices (AREA)
  • Testing And Monitoring For Control Systems (AREA)
  • Complex Calculations (AREA)

Abstract

ノイズが多く含まれるデータに対してもより有効にデータの解析を行うことを可能とするデータ解析システム及び時系列データ解析方法を提供する。時系列データ解析システムは、バンドパスフィルタ処理部１６、移動平均処理部１８、微分処理部２０、積分階差処理２２、減算処理部２４、パルス生成処理部２６、閾値処理部２８、分割点探索処理部３０、範囲設定処理部３２を含んでいる。これらの各処理部において、時系列データに対して各種の処理をさまざまに適用し実行することによって、時系列データの中の屈折点、すなわちトレンドが変化する点を高い精度で見いだすことができる。範囲設定処理部３２は、例えば、ユーザが時間軸に直接範囲を設定する処理を実行する。さらに、分割点探索処理部によって、コンピュータに対する負荷の少ない処理により時系列データによくフィットする近似曲線を求める。

Description

時系列データ解析システム
　本発明は、各種の測定装置や各種のセンサなどから出力される時系列データを解析するための、時系列データ解析システム及び時系列データ解析方法に関する。
　製造業の生産現場では、さまざまな測定装置やセンサによって各種のデータが取得されており、製品の品質の維持・向上のためには、これらのデータを解析することが必要不可欠とされている。このようなデータ解析のための従来の方法としては、これまでのところ、ＭＴ法などの波形全体の特徴を捉える方法が主流であった。
　しかしながら、実際のデータには多数のノイズが含まれており、このようなノイズを多く含むデータに対しては、従来のＭＴ法などでは、波形に含まれる「特定の形」を演算で取り出すことが難しかった。
　そこで、本発明は、ノイズが多く含まれるデータに対してもより有効にデータの解析を行うことを可能とするデータ解析システム及び時系列データ解析方法を提供することを目的とする。
　上記の課題を解決するために、本発明に係る時系列データ解析システムは、
　時系列データに対して時間微分を実行し、必要に応じて時間微分後の時系列データの波形図を表示する時間微分処理手段と、
　処理対象となる時系列データに対し、前記時間微分処理部による時間微分処理を実行するときに時間範囲を時間軸上で設定する範囲設定処理部と、
　前記範囲設定処理部によって設定された時間範囲において、処理対象となる時系列データがピークとなる位置を検出するピーク位置検出部と、
を備える（第１発明）。
　上記の第１発明において、前記範囲設定処理部は、時系列データの時間軸上に、前記時間範囲の少なくとも一方の境界をユーザが直接設定することができる（第２発明）。
　上記の第１発明において、前記範囲設定処理部は、ユーザが縦軸の特定の値を設定し、それに対応する時間軸上の値を前記時間範囲の少なくとも一方の境界とすることができる（第３発明）。
　上記の第１発明において、前記範囲設定処理部は、ピーク位置検出部によって検出された時系列データの縦軸のピーク値を検索し、当該ピーク値に対応する時間軸上の値を、前記時間範囲の一方の境界とすることができる（第４発明）。
　上記の第１発明において、さらに、ノイズ除去手段として、時系列データに対して移動平均処理を実行し、必要に応じて移動平均処理後の時系列データの波形図を表示する移動平均処理部、及び、時系列データに対して積分階差処理を実行し、必要に応じて積分階差処理後の時系列データの波形図を表示する積分階差処理部のいずれか一方又は両方を備えることができる（第５発明）。
　上記の第１発明において、時系列データに対してバンドパスフィルタ処理を実行するバンドパスフィルタ処理部をさらに有することができる（第６発明）。
　上記の第５発明において、ユーザがユーザデータを入力するユーザ入力処理部をさらに有することができる（第７発明）。
　上記の第５発明において、外部から時系列データをシステム内に取り込むデータ入力部をさらに有することができる（第８発明）。
　上記の課題を解決するために、本発明に係る時系列データ解析システムは、
　時系列データに対して時間微分を実行し、必要に応じて時間微分後の時系列データの波形図を表示する時間微分処理手段と、
　時系列データに対して、その値が連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高とする第１のパルス波形と、その値が連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高とする第２のパルス波形とを求め、必要に応じて第１のパルス波形及び第２のパルス波形の両方又はいずれか一方の波形図を表示するパルス生成部と、
　複数の時系列データの間又は複数のパルス波形の間で減算処理を実行し、必要に応じてその減算後の時系列データ又はパルス波形を表示する減算処理部と、
　パルス波形の値を予め指定された閾値と比較し、当該閾値より大きいパルス又は小さいパルスを選択する閾値処理部と、
　を備える（第９発明）。
　上記の第９発明において、さらに、ノイズ除去手段として、時系列データに対して移動平均処理を実行し、必要に応じて移動平均処理後の時系列データの波形図を表示する移動平均処理部、及び、時系列データに対して積分階差処理を実行し、必要に応じて積分階差処理後の時系列データの波形図を表示する積分階差処理部のいずれか一方又は両方を備えることができる（第１０発明）。
　上記の第９発明において、時系列データに対してバンドパスフィルタ処理を実行するバンドパスフィルタ処理部をさらに有することができる（第１１発明）。
　上記の第９発明において、ユーザがユーザデータを入力するユーザ入力処理部をさらに有することができる（第１２発明）。
　上記の第９発明において、外部から時系列データをシステム内に取り込むデータ入力部をさらに有することができる（第１３発明）。
　上記の課題を解決するために、時系列データ解析システムに取り込まれた時系列データの屈折点を求める方法は、
　前記時系列データ解析システムに含まれる微分処理部によって、取り込まれた第１の時系列データに対して３回の時間微分を行って第２の時系列データを得る第１ステップと、
　前記時系列データ解析システムに含まれるパルス生成部によって、前記第２の時系列データが連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高とする第１のパルス波形と、前記第４の時系列データが連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高とする第２のパルス波形とを求める第２ステップと、
　前記時系列データ解析システムに含まれる減算処理部によって、前記第１のパルス波形と前記第２のパルス波形との減算を行って第３のパルス波形を求める第３ステップと、
　前記微分処理部によって、前記第３のパルス波形に対して時間微分を行って第４のパルス波形を求める第４ステップと、
　前記時系列データ解析システムに含まれる閾値処理部によって、前記第４のパルス波形に対して閾値処理を行い、所定の閾値を超えるパルスのみを選択して、前記第１の時系列データの屈折点とする第５ステップと、
　を含む（第１４発明）。
　上記の第１４発明において、前記第１ステップで、２回目の時間微分を行って得られる時系列データに対して積分階差を行ってノイズを除去したあとに３回目の時間微分を実行して第２の時系列データとするとすることができる（第１５発明）。
　上記の課題を解決するために、本発明における時系列データ解析方法は、
　時系列データの波形図の近似曲線を求めることによって時系列データを解析する時系列データ解析方法であって、
　予め指定された始点と終点との間の前記時系列データ上に第１分割点から第ｎ分割点（ｎは２以上の整数）までを仮に配置するステップと、
　第１分割点を始点から第２分割点まで移動させたときに始点と第１分割点とを結ぶ第１直線がその範囲で時系列データと最も良くフィットする点を第１分割点の位置と決定し、
　第２分割点を第１分割点から第３分割点まで移動させたときに第１分割点と第２分割点とを結ぶ第２直線がその範囲で時系列データと最も良くフィットする点を第２分割点の位置と決定するステップと、
　以下同様の手順を続行し、第ｎ分割点を第ｎ－１分割点から終点まで移動させたときに第ｎ－１分割点と第２分割点とを結ぶ第ｎ直線がその範囲で時系列データと最も良くフィットする点を第ｎ分割点の位置と決定するステップと、
　前記始点と、第１分割点と、・・・第ｎ分割点と、終点とをこの順番で直線によって連結して得られる曲線を前記時系列データの近似曲線とするステップと、
　を含む（第１６発明）。
本発明の実施の一形態に係る時系列データ解析システムの機能ブロックズである。 射出成形時に樹脂が流し込まれてから固化するまでの圧力の変化を測定する圧力センサの出力の時系列データである。 図２の波形図の一部を時間方向に引き延ばした波形図である。 時系列データに時間微分を行って得られる時系列データの一例を示す波形図である。 時系列データに対して移動平均をとって得られる時系列データの一例を示す波形図である。 時系列データに対する積分階差について説明するための図である。 時系列データに対する積分階差について説明するための図である。 ゼロを通過する時点を閾値処理で判定する場合の問題点を示す図である。 微分によって得られた波形図に基づいてパルスを生成するパルス生成部の動作を示す波形図である。 時系列データに対して各種処理を実行して、元の時系列データの屈折点を求める例を示す波形図である。 図９で求めた屈折点に対して閾値処理を実行して所定の閾値以上のパルスを選択することを示した波形図である。 ノイズが多く含まれる波形図から正しくピークの時点を求める方法を示す波形図である。 元の波形図の一つの屈折点を求める方法の一例を示した波形図である。 元の波形図の一つの屈折点を求める方法の一例を示した波形図である。 本発明の実施の一形態のシステムで実行できる処理ことをまとめた図である。 本発明の実施の一形態のシステムのユーザインターフェイスを示す図である。 時系列データ全体に対して階差処理を行った場合の問題点を説明する図である。 移動平均の度合いを強めることによってノイズを低減できるとこを説明する図である。 移動平均の点数を増やすことによる問題点を説明する図である。 オペレータが手動で時間範囲を設定する場合について説明する図である。 縦軸の値域に着目して範囲を決める方法について説明する図である。 特徴点の特性に基づいて検索する範囲を設定する方法について説明する図である。 速度データの最大点に基づいて範囲を設定する方法について説明する図である。 複数の時系列データのお互いの関係性に基づいて処理範囲を設定する処理について説明する図である。 実際に時系列データに対して処理範囲を設定する具体的な方法について説明する図である。 実際に時系列データに対して処理範囲を設定する具体的な方法について説明する図である。 実際に時系列データに対して処理範囲を設定する具体的な方法について説明する図である。 元の時系列データの近似曲線を求める方法を説明するための図である。 元の時系列データの近似曲線を求める方法を説明するための図である。 元の時系列データの近似曲線を求める方法を説明するための図である。
　図１は、本発明の実施の一形態に係る時系列データ解析システムの主要部の機能ブロック図である。図１に示す各処理部は、パーソナルコンピュータを含む一般のコンピュータにおけるプロセッサ、ＤＲＡＭ等の一時的記憶装置、ＨＤＤ装置等の非一時的記憶装置、データ入力手段としてのキーボードやマウス、ディスプレイ装置等のハードウェアと、記憶装置から読み出されてプロセッサによって実行されるコンピュータプログラム等のソフトウェアとの組み合わせによって、各種処理が行われる。
　図１において、データ入力処理部１０は、測定装置や各種のセンサなどからの時系列データをシステム内に取り込む処理を行う。ここで、時系列データとは、一定の時間間隔で取り込まれたデータ、あるいは既にシステムに取り込まれた時系列データに対して何らかの処理が行われたデータを指す。ユーザ入力処理部１２は、ユーザがキーボードやマウス等から入力したデータをシステム内に取り込む処理を行う。画像出力処理部１４は、文字、図形等を含むユーザインターフェイスや、実際の時系列の波形データ等をディスプレイや、必要に応じてプリンタ等の出力装置に出力する処理を行う。
　バンドパスフィルタ処理部１６は、システムに取り込まれた時系列データに対して、あるいは既にシステムに取り込まれて何らかの処理がなされた時系列データに対して、バンドパスフィルタ処理を実行し、必要に応じて得られた波形図を表示する。なお、バンドパスフィルタ処理部１６は、時系列データに対して直接バンドパスフィルタ処理を実行するだけでなく、必要に応じて高速フーリエ変換を行い、周波数空間でバンドパスフィルタ処理を行ってから逆高速フーリエ変換を行って時間空間に戻す処理を実行することもできる。移動平均処理部１８は、システムに取り込まれた時系列データに対して、あるいは既にシステムに取り込まれて何らかの処理がなされた時系列データに対して、隣接する所定の数の時間位置で平均をとる処理を実行する。この所定の数は、デフォルトで設定されているものをそのまま使用することも、後述するユーザインターフェイスを介してユーザが設定した値を使用することもできる。
　微分処理部２０は、システムに取り込まれた時系列データに対して、あるいは既にシステム内で何らかの処理がなされた時系列データに対して、時間微分演算を実行し、必要に応じて得られた波形図を表示する。時系列データを取り込む時間間隔（サンプリング周波数の逆数）が例えば１０００分の１秒であったとすると、この１０００分の１秒の間隔で隣り合う２つのデータの値の差を１０００分の１秒で割ることによって時間微分が得られる。ただし、この時系列ごとに１つ置き又は２つ置き等時間間隔を広げて微分処理を実行することもできる。このように有限の時間間隔で時間微分をとることを、階差とも呼ぶ。積分階差処理部２２は、システムに取り込まれた時系列データに対して、あるいは既にシステム内で何らかの処理がなされた時系列データに対して、積分演算を実行し、さらに得られた積分データに対して、階差処理を実行して、必要に応じて得られた波形図を表示する。
　減算処理部２４は、ある二つの波形データの間で減算処理を実行し、必要に応じて得られた波形図を表示する。パルス生成部２６は、時系列データに対して、その値が連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高とする第１のパルス波形と、その値が連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高とする第２のパルス波形とを求め、必要に応じて第１のパルス波形及び第２のパルス波形の両方又はいずれか一方の波形図を表示する。閾値処理部２８は、ユーザが閾値を設定することによって、複数のパスルの中からその閾値を超える、又はその閾値を下回るパルスを選択するという処理を実行し、必要に応じて選択されたパルスを表示する。分割点探索処理部３０は、ある時系列データのある始点から終点までの範囲を、その範囲に含まれるいくつかの点（分割点）を直線でつないだ直線からなる曲線で近似するときに、その時系列データの特性を判断するのに十分な近似ができるように、これらの分割点の位置を決めるという処理を実行する。
　処理範囲設定処理部３２は、システムに取り込まれた時系列データ、あるいは既にシステム内で何らかの処理がなされた時系列データの全体を処理対象とするのではなく、ある一定のルールに基づいて時間範囲を設定し、その設定された時間範囲の時系列データだけを処理対象とする処理を実行する。処理範囲設定処理部３２によって行う具体的な処理の内容については後に詳述する。
　図１に示す移動平均処理部１８及び積分階差処理部２２は、ノイズ除去処理を行うものであり、システムの構成として必ず必要なものというわけではなく、必要に応じてこれらのうちのいずれか一方又は両方を備えるようシステムを構成することができる。また、図１に示したすべての処理部を備えていなければ本発明が実現できないというものではない。例えば、微分処理部２０、範囲設定処理部３２、ピーク位置検出部３４だけを含んでいれば、時系列データに対する最低限必要な処理を行うシステムを構築することができるし、微分処理部２０、パルス生成部２６、減算処理部２４、閾値処理部２８だけを含んでいても、時系列データに対する最低限必要な処理を行うシステムを構築することができる。あとは、実際の処理の必要性に応じて図１に示すその他の処理部を追加してシステムを構築することができる。
　図２は、射出成形において、樹脂が型に流し込まれてから固化するまでの間に加わる圧力（縦軸）が時間（横軸）とともにどのように変化するかを示した、実際の圧力センサの出力の時系列データを示している。このようなデータの時系列データ解析システムへの取り込みは、図１のデータ入力処理部１０によって行われる。
　図２には、信号波形に重ねて４カ所に円１０ａ、１０ｂ、１０ｃ、１０ｄが示されている。図２に示す信号波形において、一番左の円１０ａは、圧力がほぼゼロの状態から急激に増大し始める時点であり、樹脂が型の内部全体に行き渡る直前の瞬間を示している。そして、１０ａの時点から圧力が急激に増大して円１０ｂのピークに到達する。このピークは、樹脂が型全体に充填され、樹脂の流し込みが停止された時点に対応する。樹脂の流し込みが停止されると圧力は急速に低下するとともに、円１０ｃの時点からしばらく圧力一定の状態が続く。この圧力一定の期間は、樹脂の温度が冷やされている期間に対応する。そして、樹脂が十分に冷やされて固化すると、円１０ｄの時点から再び急速に圧力が低下する。
　図２に示した型の内部の圧力の時系列的な変化の仕方は、型の大きさや種類、樹脂の組成等によって変わってくる。したがって、例えば図２の円１０ａ、１０ｂ、１０ｃ、１０ｄに対応する特徴的な変化が時系列の中で現れる時点（屈折点）を捉えることで、成形品の不良の発生や、射出成形作業の不具合等を知ることができる。すなわち、射出成形作業そのものや成形品の状態を、圧力センサによってモニタすることができる。
　ところで、射出成形作業や鋼板の切り出し作業を圧力センサや加重センサでモニタする場合には、大量のデータが得られる。このような大量のデータから上述のような特徴的な屈折点を正確に抽出するには、迅速な処理が必要となるので、処理の高速化が求められ、それだけ計算機に対する負荷が増大する。そこで、本実施形態では、以下のようなデータ処理を行う。
　上述の射出成形の例で、樹脂が型の内部全体に行き渡る瞬間（図２の円１０ａ）の前後の変化を、時間を大きく引き延ばすと、図３（ａ）のようになると考えられる。図３（ａ）に円１１ａで示す、樹脂が型の内部全体に行き渡る時点から圧力が上昇し始める。すなわち、この時点が屈折点となる。屈折点を知りたい場合に、まず図１に示した微分処理部２０により、図３（ａ）に示したグラフに対して時間微分をとると、図３（ｂ）のような時系列データの曲線が得られる。図３（ａ）の変化を変位（displacement）と考えれば、図３（ｂ）の波形は速度（velocity）に対応する。そして図３（ｂ）の曲線に対し、微分処理部２０によりさらに微分をとると、図３（ｃ）に示す曲線が得られる。これは加速度（acceleration）に対応する。この図３（ｃ）の曲線において、円１１ｂのような極大点が見つかれば、その時点が、樹脂が型の内部全体に行き渡った瞬間の屈折点であると判断することができる。
　そこで、図３（ｃ）と同じ図４（ａ）の曲線に対し、微分処理部２０によりさらにもう１回時間微分を行って、図４（ｂ）の曲線を得る。この図４（ｂ）の曲線がゼロを通過する時点（上向き矢印で示す）が分かれば、この時点が図４（ａ）の曲線の極大点、すなわち樹脂が型の内部全体に行き渡った時点であると判断することができる。しかしながら、図３（ａ）に示した曲線は、樹脂が型の内部全体に行き渡った場合の圧力の変化を想定した曲線に過ぎず、実際の圧力センサの信号波形はこのような滑らかなグラフにはならず、これにかなりのノイズが載った状態で出力される。
　図５は、ノイズが多く含まれる時系列データから極大点を求めようとした場合の問題点を示している。図５（ａ）は、図４（ａ）に対応する、極大点を求めようとする信号波形の一例を示しており、図４（ａ）の曲線と比較するとかなりのノイズが含まれている。そこで、まず、図１の移動平均処理部１８により移動平均をとる。この図５（ａ）の信号波形に対して、４９点移動平均をとったものが図５（ｂ）の波形であり、９９点移動平均をとったものが図５（ｃ）の波形であり、１９９点移動平均をとったものが図５（ｄ）の波形である。ここで、複数点移動平均は、階差すなわち微分を行ったことに相当するが、一般には、移動平均の点数が多くなればそれだけノイズは減少する。
　図５の（ｂ）、（ｃ）、（ｄ）の波形は、異なる点数で移動平均をとったものであるが、（ｂ）、（ｃ）、（ｄ）の波形がゼロを通過するそれぞれの時点ｔb、ｔc、ｔdは一致してはおらず、前後に多少ずれている。このことは、データの処理の仕方によって、すなわちこの場合だと、何点で移動平均処理を行うかによって、得られた波形がゼロを通過する時点が変わってしまうことを示している。したがって、上記図４に示した、加速度に対して微分（階差）をとるという方法では、ノイズが含まれる出力信号から、トレンドが変わる時点を正確に求めることはできない。
　そこで、上記の問題を改善するために、本実施形態では、積分階差という処理を実行する。すなわち、ゼロを通過する時点を求めたい波形に対して、ある等間隔の時間で積分を行い、この積分値に対して差分をとる。具体的には、
Figure JPOXMLDOC01-appb-I000001
という計算を実行する。
　上記の積分階差の式について、図６Ａ及び図６Ｂを参照しながら説明する。まず、図６Ａ及び図６Ｂの（ａ）の波形（図５（ａ）の波形と同じである）において、いま注目している時点をｉとし、積分範囲をｗとすると、上記の式の右辺第１項は、波形ａの値をｉからｉ＋ｗまで積算したもの、すなわちｉからｉ＋ｗまでのａの曲線の下側の面積（ハッチングで示す）に対応する。同様に右辺第２項は、波形ａの値をｉ－ｗからｉまで積算したもの、すなわちｉ－ｗからｉまでの曲線の下側の面積（梨子地で示す）に対応する。したがって、積分階差は、ｉの右側の面積とｉの左側の面積との差を表している。ここで、ｗが９９であれば９９点積分階差、ｗが１９９であれば１９９点積分階差という。
　図６Ａ及び図６Ｂにおいて、（ｂ）の波形は、（ａ）の波形ａに対して、９９点積分階差を実行して得られた波形図を示している。そして、図６Ａでは、積分階差の値がゼロとなる場合に対応し、図６Ｂでは、右辺第１項と第２の差がほぼ極小となっている場合に対応している。図６Ａ及び図６Ｂにおける（ｂ）の波形を見ると分かるように、元の波形の特徴を大きく失うことなく、極大点、すなわちゼロを通過する時点を特定するのに支障がない範囲でノイズが除去されていることが分かる。図６Ａの（ａ）の波形を見ると、ｉの時点よりもやや左側にこの波形が実際にピークとなる時点がある。人がこのピークを見れば、このピークがノイズに起因するものだと判断することができるが、このような判断をコンピュータに実行させようとすることは非常に難しい。これに対して、上記の積分階差処理を実行して得られる図６Ａの（ｂ）の波形をみると、この波形がゼロを通過する時点が、図６Ａの（ａ）の波形曲線の全体的な極大点と一致していることが分かる。
　次に、図６Ａの（ｂ）に示すような波形がゼロを通過する時点をコンピュータによって判定する具体的な方法について説明する。図６Ａの（ｂ）に示す波形も、実際には飛び飛びの時間（横軸）に対してプロットし、これを曲線で結んだ時系列のデータである。図７（ａ）、（ｂ）は、上の波形図が、図６Ａの（ｂ）に対応するものであり、下の波形図は上の波形図を時間微分したものである。したがって、下の波形図がゼロになる時点が、上の波形図がゼロを通過する時点となる。しかし下の波形図も上の波形図と同様に飛び飛びの時間に対してプロットされるものである。このように時間解像度が有限であるため、図７（ａ）、（ｂ）の下に示す時系列データに完全なゼロが含まれるとは稀である。
　データがゼロとなる点を探索する場合には、図７（ａ）、（ｂ）の下の波形図に示すように、縦軸のゼロを中心とした一定幅のしきい値（横軸の時間軸を中心とする一定幅の帯状の部分）を設定し、この範囲に含まれるデータをゼロと判定することが考えられる。図７（ａ）、（ｂ）の下の図で、閾値の範囲にあるデータについて、分かりやすいように縦軸に平行な線で示している。この場合、閾値の幅が（ａ）のように狭いと、曲線が大きく変化する範囲に該当するデータが存在しない場合があったり、また、閾値の幅が（ｂ）のように広いと、該当するデータが多数存在して本来見つけたいデータを見つけることができないといった事態が生じる。このように、アプリケーションごとに適切な閾値を設定することは容易ではなく、効率も悪いという問題がある。別の方法として、時系列データを１点ずつ探索し、符号がプラスからマイナスへ、あるいはマイナスからプラスへ変わる時点を見つける方法も考えられる。しかしながら、この方法では処理がスカラー処理となるため、効率が悪い。
　そこで本実施形態では、次のようにして波形図がゼロを通過する時点を見いだす。図８（ａ）は、図７（ａ）、（ｂ）の上の波形図に対応し、図８（ｂ）は、図７（ａ）、（ｂ）の下の波形図に対応するもので、図８（ｂ）は図８（ａ）の波形を時間微分した波形である。この図８（ｂ）の波形について、図１に示したパルス生成部２６によって、図８（ｃ）、（ｄ）のようなパルス波形を導く。図８（ｃ）のパルス波形は、図８（ｂ）の波形において、連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高として導いたものである。例えば、図８（ｂ）の中央の大きく変化する部分でＴ1の時間幅に対応する図８（ｃ）の部分にそのＴ1の値に対応する高さのパルスが示されている。一方、図８（ｄ）のパルス波形は、図８（ｂ）の波形において、連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高として導いたものである。例えば、図８（ｂ）の中央の大きく変化する部分でＴ2の時間幅に対応する図８（ｄ）の部分にそのＴ2の値に対応する高さのパルスが示されている。図８の（ｃ）と（ｄ）は、縦軸の値は別として、時間軸で見ると互い違いになっている。
　次に、図９に示すように（図９（ａ）～（ｄ）は、図８（ａ）～（ｄ）に対応する）、図１の減算処理部２４によって、図９（ｄ）と図９（ｃ）を減算して差分をとると、図９（ｅ）の波形図が得られ、さらに図９（ｅ）の波形図に対して時間微分を行うと、図９（ｆ）に示すような、正から負へ切り換わった時点と負から正に切り換わった時点のみにパルスが残る波形図が得られる。そしてこのパルスが正か負かによって、正から負に切り換わったのか、負から正に切り換わったのかが分かる。図９（ｆ）の波形図の各ピークは、極大点の時間的な位置に対応し、各ピークの高さは各階差が同じ符号で連続する長さを反映している。なお、ここまでの処理はすべてベクトル処理であり、データを一つひとつ検索する必要がないため、計算処理の負荷が大幅に軽減される。
　さらに、図９（ｆ）の波形図に対して、前述の積分階差を実行すると、図１０（ｈ）の波形図が得られる（図１０の（ａ）は、型の内部の圧力を検出する圧力センサの出力波形であり、（ｂ）～（ｇ）は、図９の（ａ）～（ｆ）に対応する）。ここで、図１の閾値処理部２６により、ある一定の閾値を設定しておいて、これを超えるパルスのみを選択すれば、この時点が図１０（ａ）の時系列データの屈折点、すなわちトレンドが変化した時点であることが分かる。
　次に、ノイズの含まれるデータから実際のピークとなる時点を探索する例について説明する。図１１において、（ａ）は実際の型の内部の圧力を検出する圧力センサの出力データであり、このかなりのノイズが含まれるデータから圧力の最大点を求める。図１１において、横軸は時間である。この図１１（ａ）の出力データに対して、高速フーリエ変換（ＦＦＴ）を実行し、得られたスペクトルに対してバンドパスフィルタで不要な成分を除去し、さら逆ＦＦＴを実行して得られたのが、図１１（ｂ）に示す波形である。なお、図１１（ｂ）において、「０－４０Ｈｚ」は、４０Ｈｚより大きい周波数の波形をカットして０～４０Ｈｚの周波数の波形だけを取り出すことを意味する。
　そして、図１１（ｂ）の波形に対して、９点の移動平均をとったものが図１１（ｃ）に示した波形である。この図１１（ｃ）の波形の最大点が、圧力の最大時点と考えることができる。ところで、図１１（ａ）の波形と図１１（ｃ）の波形を比較すると、図１１（ａ）の実際のピークの位置は、図１１（ｃ）のピークの位置よりも左側にずれている。これは、ノイズの影響によるものと考えられる。このようにノイズによってピークの位置がずれるような場合であっても、図１１の（ａ）→（ｂ）→（ｃ）の処理を行うことによって、ノイズの影響を排除して圧力が最大となる時点を求めることができる。この図１１の（ａ）→（ｂ）→（ｃ）の処理を行う部分は、機能的には図１に示すピーク位置検出部３４に対応する。
　次に、図２の圧力センサの出力波形の円１０ｃに当たる時点、すなわち樹脂が充填され圧力がピークになった直後に、圧力が急激に低下したあとにほぼ一定の安定な期間が始まる時点を見いだす処理について説明する。図１２において、（ａ）は図２に示した圧力センサの出力波形を、横軸の時間軸の方向に伸ばした状態を示している。図１２の（ａ）の波形に対して、図１のバンドパスフィルタ処理部１６によって、図１１の（ｂ）と同様のバンドパスフィルタ処理を施して得られたのが、図１２の（ｂ）に示す波形である。
　そして、（ｂ）の波形に対して、図１の移動平均処理部１８によって９点の移動平均をとって得られたのが（ｃ）に示す波形である。そしてさらに、（ｃ）の波形に対して微分処理部２０によって時間微分を行って得られたのが（ｄ）に示す波形であり、これにさらに時間微分を行って得られたのが（ｅ）に示す波形である。この（ｅ）に示す波形には、横軸の値が６００のやや手前にピークが見られるが、（ａ）の波形と対照させると、この時点が、安定な期間の開始する時点であることが分かる。
　続いて、図２の圧力センサの出力波形の円１０ｄに当たる時点、すなわち樹脂が十分に冷やされて固化して圧力が低下する時点を見いだす処理について説明する。図１３において、（ａ）は図２に示した圧力センサの出力波形を、横軸の時間軸の方向に伸ばした状態を示しており、（ａ）の円で示した部分が、樹脂が十分に冷やされて固化して圧力が低下する時点に該当する。図１３（ａ）の波形に対して、図１のバンドパスフィルタ処理部１６によって、図１１の（ｂ）と同様のバンドパスフィルタ処理を施して得られたのが、図１３の（ｂ）に示す波形である。
　そして、図１３の（ｂ）の波形に対して、図１の移動平均処理部１８によって９点の移動平均をとって得られたのが（ｃ）に示す波形である。そしてさらに、（ｃ）の波形に対して微分処理部２０によって時間微分を行って得られたのが（ｄ）に示す波形である。ここでは、この（ｄ）に示す波形から、あるほぼ一定の正の値から一時的に負になる時点が、図１３（ａ）の円で示した部分に対応することが分かる。なお、この例では、（ｄ）の波形に対してさらに時間微分を行って得られる（ｅ）の波形、すなわち加速度に対応する波形では、下降の動きが緩やかであるため、この時点をうまく検出することは難しい。
　これまで説明してきたように、図１に示した各処理部により、時系列データに対して各種の処理をさまざまに適用し実行することによって、時系列データの中の屈折点、すなわちトレンドが変化する点を高い精度で見いだすことができる。図１４は、本実施形態のシステムでできることをまとめた図である。データ処理としてできる基本的な処理は、バンドパスフィルタ処理、移動平均処理、階差積分処理、微分処理の４つで、これらを組み合わせることによって各種のデータ処理が可能となる。そして、処理がなされたデータに対して、最大値、最大値の位置（時間）、最小値の位置（時間）、閾値処理、積分のうちの１つ又はいくつかを求めることによって屈折点、すなわちトレンドの変化点を特定することができる。
　図１５は、図１に示した時系列データ解析システムのユーザインターフェイスの一例を示している。このようなユーザインターフェイスを利用することにより、ユーザが実際に時系列データに対する処理を実行する際の利便性が向上する。このユーザインターフェイスを介して、任意の時系列データに対して上述の各処理を実行し、その処理結果を直ちに視覚的に見ることができる。この画面の左上の「データソース」欄には、処理対象となる時系列データの名称が記載され、このすぐ右側の「Ｒａｗ ｄａｔａ」欄に、圧力センサや加重センサの出力データなどの、処理の対象となるコンピュータに入力された実際の時系列データが表示される。
　さらに、その右側には「ＦＦＴ ｓｐｅｃｔｒｕｍ」欄が設けられ、ここには「Ｒａｗ ｄａｔａ」欄に表示されている波形図に対して高速フーリエ変換（ＦＦＴ）が行われたＦＦＴスペクトルが表示される。このＦＦＴ処理は、「Ｒａｗ ｄａｔａ」欄に表示される入力データすべてに対して実行することもできるし、ユーザがＦＦＴ処理を実行するかどうかを選択するようにもできる。なお、この「ＦＦＴ ｓｐｅｃｔｒｕｍ」欄では、バンドパスの設計が可能である。図１５の「ＦＦＴ ｓｐｅｃｔｒｕｍ」欄の波形を見ると、横軸の５０Ｈｚ前後にピークが集中していることが分かる。そこで、ここではバンドパスの上限として２５Ｈｚが指定されており（２５Ｈｚより高い周波数部分のグレーの網かけがそのことを示している）、２５Ｈｚ以下の波形だけが抽出されている。
　また、図１５において、「データソース」欄の下には、サンプリングレートを選択する欄が設けられている。このサンプリングデータは、時系列データの隣り合う個々のデータの時間間隔に対応するものであり、ここでは「１００Ｈｚ（１００分の１秒間隔）」、「１０００Ｈｚ（１０００分の１秒間隔）」、「１００ｋＨｚ（１０万分の１秒間隔）」のうちのいずれかを選択することができるようになっている（図１５では１０００Ｈｚが選択されている）。その下の「検出モード」欄では、時系列データに対して変位（displacement）、速度（velocity）、加速度（acceleration）のいずれかを選択して表示させることができる。さらにその下の「バンドパス下限（Ｈｚ）」及び「バンドパス上限（Ｈｚ）」では、バンドパス処理の上限値及び下限値を、画面上でグラフィック的に設定できるようになっている。
　その下の「移動平均範囲」欄は、移動平均をとる範囲、すなわち平均をとる連続する時系列データの数を設定するものであり、１から１９までの範囲で設定できる。その下の「ｃｅｎｔｅｒ」欄は、移動平均の中央を選択するスイッチである。その下の「屈折点検出」欄は、屈折点検出を実行するかどうかをチェック欄にチェックするかしないかで設定する。その下の「ピーク方向」欄は、ピークを検出する場合に、上方向のピーク及び下方向のピークのいずれか又は両方を検出することを選択する欄である。「ピーク幅」欄は、ピークを検出する場合に、ここで設定した幅以下の幅のピークだけ検出するように設定する欄である。「ピーク閾値（最大ピークに対する割合）」は、ここで設定した閾値を超える値を有するもののみをピークとして検出するように設定する欄である。そして、図１５の右下で最も広い領域には、これまでの各欄で設定した内容に基づいて入力された時系列データを処理して得られる波形図が示されている。
　以上の説明から理解されるように、時系列データが特徴的な変化を示す屈折点は、変位、速度、加速度のいずれかのピークとして見いだすことができる。そして、変位から速度を求める方法、及び、速度から加速度を求める実際の方法は、実装上は時間微分、すなわち階差をとることによる。しかしながら、入力された時系列データの全体に対して階差をとる処理を行うと、その時系列データにおける長期のトレンドが失われることになる。長期のトレンドが失われることは、短期の変動、すなわちノイズを強調することになって、特徴点の抽出の妨げとなる。図１６はこのことを説明する図であり、（ａ）は、ある荷重の変化を示す時系列データ、（ｂ）は、（ａ）の時系列データを時間微分して得られる速度データ、（ｃ）は（ｂ）の速度データを時間微分して得られた加速度データを示している。図１６の（ｃ）に示すように、（ａ）の時系列データ全体に対して時間微分を求める処理を連続して行うと、ノイズが強調されて、（ａ）に円で示したトレンドが変化した特徴点を見つけることが困難となる。
　また、図１７は、移動平均の度合いを強めることによってノイズを低減できるとこを示している。図１７の（ａ）は、速度及び加速度に対して５点移動平均をとったもの、（ｂ）は９点移動平均をとったもの、（ｃ）は１９点移動平均をとったものを示している。これらの図で、図の下に示した上向きの矢印は、特徴点と判定された時点を示しているが、（ａ）では、荷重が低下し始める特徴点ではなくそれよりも右側にずれた時点を特徴点として抽出している。（ｂ）では、（ａ）に比べてノイズは低減されており、特徴点を検出できているが、たまたま検出できた可能性が高く、堅牢性に問題がありうる。（ｃ）では、ノイズはかなり除去されており、ノイズはさらに低減されて特徴件が検出されている。
　しかしながら、移動平均の点数を増やすと、図１８に示すような問題が生じうる。図１８において、実線が荷重の変動を示す実際の曲線を、点線は元の荷重のデータに対して１９点移動平均を取った曲線を、一点鎖線は元の荷重のデータに対して３９点移動平均を取った曲線をそれぞれ示している。図１８において、例えば荷重のかかり始めを示す曲線の立ち上がりを知りたいときに、３９点移動平均を取ってしまうと、実際の立ち上がりの時点から左側（時間の早い側）に大きくずれてします。このように、移動平均を取る点数を増やしてゆくと、元の曲線から大きく波形が変化してしまい、単純に移動平均の点数を増やしてゆくだけでは正確な特徴点が難しくなる。
　図１９の一番上のグラフは、射出成形において、樹脂が型に流し込まれてから固化するまでの間に加わる圧力（縦軸）が時間（横軸）とともにどのように変化するかを示した圧力センサの出力の時系列データのグラフを示しており、２番目のグラフは一番上の時系列データに対し、一回の微分を行って速度を求めた時系列データのグラフであり、一番下のグラフは、さらにもう一回微分して加速度を求めた時系列データを示している。一番上のグラフでは、横軸の値が２５０の辺り（ｔ1で示す）で縦軸の値が急激に上昇し、３００を少し過ぎた辺り（ｔ2で示す）でピークに達して低下し始め、４００を少し超えた辺り（ｔ3で示す）で低下速度が急速に減少してからしばらくほぼ一定の値が続いている。このように、横軸が２００から５００までの間に大きく３つのトレンドの変化が認められ、それぞれに意味が異なっている。
　圧力センサからの出力が、図１９の一番上のグラフのような挙動を示すことを分かっていれば、例えばｔ3の特徴点を抽出しようとするオペレータは、横軸が４００を少し超えた辺りの前後の適当な範囲に限定して、図３から図１３に示したような処理を実行するようにすれば、特徴点の抽出処理を効率化することができるとともに、精度も向上する。そこで、本実施形態では、図１９のような時系列データに対して特徴点の抽出処理を実行するときに、図１９の時系列データを表示させた状態で、横軸が４００を少し超えた辺りの前後にグレーの時間幅を手動で直接設定することによって、その設定範囲についてのみ図３から図１３に示したような処理を実行するようシステムに指示することができる。　
　図２０は、横軸の時間軸ではなく、縦軸の値域に着目して範囲を決める方法を示している。図２０に示すグラフは、図１９に示したグラフと同じであるが、２番目に示した速度のグラフにおいて、縦軸の値がマイナス０．１以下になる範囲を設定し、この設定範囲についてのみ図３から図１３に示したような処理を実行するようシステムに指示する。このようにすることによって、横軸が３００前後になる加速度が大きく変化する部分を、特徴点の探索処理から排除することができる。このように、縦軸の値に基づいて閾値を指定し、それよりも元の値が小さい範囲又は大きい範囲を指定することで、本当に見つけたい特徴点を確実に見つけることができる。
　図２１は、図１９、図２０と同じグラフであり、ここでは、図１９のｔ3に対応する時点の前後で加速度が最大となる時点を見つけようとしている。そこで、時刻ｔ3の前後のある範囲（ここでは１５０ポイントの範囲）に限定し図３から図１３に示したような処理を実行するようシステムに指示する。このようにして特徴点の特性に基づいて検索する範囲を設定することによっても、特徴点の抽出処理を効率化することができ、精度も向上する。
　図２２は、図１９～図２１に示したセンサセンサの時間的変化を示すグラフのうち、加速度最小点と速度最大点の間を検索範囲とし、この範囲で加速度が最大となる点を検索するたに、図３から図１３に示したような処理を実行するようシステムに指示する例である。
　次に、複数の時系列データのお互いの関係性に基づいて処理範囲を設定する処理について説明する。図２３の上のグラフには、同時に取り込まれた３つの時系列データが示されている。これらは、射出成形の型の３カ所に設けられた圧力センサからの出力が、樹脂が型に重点される前後で圧力が時間的にどのように変化するかを示したグラフである。３つの圧力センサから同時に３チャネルのデータが発生し、システムに入力される。図２３の上に示すグラフから分かるように、同じ一つの型の内部でも、場所によって圧力の変化が同一ではないが、全体としての変化の傾向は類似していることが分かる。この３つセンサからのデータの関係性をみて、範囲を設定することについて説明する。
　図２３の上に示した３チャネルの時系列データから、圧力上昇が開始される点を見つける場合を考える。図２３の下に示すグラフは、上に示す３つのデータについての分散（３つの時系列データの各時点における平均値とそれぞれのデータとの差分をとり、さらに各差分を２乗して合計したもの）を求めたものである。図２３の下のグラフを見ると、樹脂が充填される前はほとんど変化がない。圧力センサが設けられた位置は互いに異なるため、圧力が上昇し始める時刻が異なるため、３つの圧力が上昇を始める直前に分散がなり最大となり、そしていったん樹脂が充填されると、圧力センサの出力は再び類似した挙動に戻る。そこで、図２３の下に示す例では、分散が最初にピークとなる時点の前後の－５０～＋１０ポイントの範囲を検索対象として設定すれば、３つの圧力センサの立ち上がり点を検索範囲とすることができる。
　次に、図２４～２６を参照して、ＤＳＬを用いて、実際に時系列データに対して処理範囲を設定する具体的な方法について説明する。図２４は、下部に図１９に示した時系列データのグラフを示し、上部に、下の時系列データに対してオペレータが処理範囲を指定するためのデータ入力のためのＤＳＬの一例を示している。この例では、オペレータが、図の下部の時系列データを見ながら「HORIZONTAL_LIMIT = 」の右側の括弧内に「３８０」及び「４５０」と入力した状態を示している。この操作により、時系列データには、横軸（時間軸）が３５０から４５０までの範囲にグレーの帯が表示され、この範囲が処理範囲として指定される。このグレーの帯の幅を見て適切でないと判断したときは、「HORIZONTAL_LIMIT = 」の右側の括弧内の数字を変更して適切な範囲を再度設定することができる。なお、「ROLLING_MEAN = 5」は、移動平均の点数が「５」であることを示しており、「VERTICAL_LIMIT = [None, None]」は、垂直方向の範囲指定をしないことを示しており、「TARGET = IDXMIN(ACC)」は、加速度（ＡＣＣ）が最小になる位置を検索することを指示したことを示している
　図２５は、図２０に示した時系列データに対して、縦軸の値域に着目してＤＳＬを用いて処理範囲を設定する具体的な方法を示している。図２５では、「HORIZONTAL_LIMIT = [None, None]」となっており、横軸について範囲設定は行わないが、代わりに「VERTIAL_LIMIT = [None, -0.1，VCT]」と設定されて、２番目の速度の縦軸の値が－０．１以下の範囲について検索することを指示したこと示している。なお、「ROLLING_MEAN = ５」は、移動平均の点数が「５」であることを示しており、「TARGET = IDXMAN(ACC)」は、加速度が最大にある位置を検索することを指示したことを示している。
　図２６は、の下部には、図２２に示した時系列データに対して、加速度最大点を検索するための、ＤＳＬを用いた処理範囲の設定の仕方を示している。図２６において、「HORIZONTAL_LIMIT = [IDXMAX(DST), IDXMAX(DST)+150]」とあるのは、元波形（ＤＳＴ）の値域が最大となる点を検索範囲の左端とし、元波形の縦軸の値が＋１５０である点を検索範囲の右端として設定したことを示している。また、「TARGET = IDXMAX(ACC)」とあるのは、加速度が最大となる位置を検索することを指示したことを示している。
　以上のように、オペレータが時系列データの特性を知っている場合には、その特性に合わせて容易に範囲を設定して、時系列データの全体ではなく設定された一部の範囲だけをシステム検索させることができる。このようにすることによって、検索処理効率及び検索精度を向上させることができ、合わせて図１８に関連して説明した問題点を回避することができる。
　次に、時系列データに対し少ない数の直線を分割点で連結して時系列データにより良くフィットする近似曲線を求める実施形態につい手説明する。図２７は、プレス加工の際に使用する歪みゲージというセンサからの出力データの波形を示している。図２７の出力データの波形に矢印ｘ1、ｘ2で示した２つのピークの間の波形の様子を見ることによって、プレス加工がうまくいったか、あるいは何か不都合が生じたか等を知ることができることが経験的に知られている。そこで、プレス加工がうまくいったかどうかを自動的に判定するために、２つのピークであるｘ1とｘ2の間をなるべく少ない数で分割し、その分割点を直線でつなぐことによって最も良く出力データの波形にフィットする近似曲線を求め、この近似曲線に基づいて上記のような判断を行う。そのために、本実施形態では、図１に示した時系列データ解析システムの分割点探索処理部３０により、以下のような処理を実行する。
　図２８に示した波形は、図２７に示した波形と同じではないが、大体の傾向は一致しているので、図２８を用いて、分割点探索処理について説明する。図２８において、始点ｘ1と終点ｘ2という２つのピークの間にｐ1、ｐ2、ｐ3という３つの分割点を波形データの上に仮に置き、その間を直線でつなぐ。分割点ｐ1、ｐ2、ｐ3の仮の位置は、分割点探索処理部２４により、例えばｘ1とｘ2の間に均等に配置することもできるが、それ以外の仮の配置をユーザが選択して設定することもできる。
　ｐ1、ｐ2、ｐ3の最初の位置が決まれば、ｘ1とｐ1とをつなぐ直線ａ、ｐ1とｐ2とをつなぐ直線ｂ、ｐ2とｐ3とをつなぐ直線ｃ、ｐ3とｘ2とをつなぐ直線ｄという４つの直線ａ、ｂ、ｃ、ｄが得られる。分割点ｐ1、ｐ2、ｐ3はあくまでも仮に置かれた点なので、この時点の４つの直線ａ、ｂ、ｃ、ｄからなる曲線が、図２７の出力データの波形を最も良く表している可能性は低い。しかしながら、分割点ｐ1、ｐ2、ｐ3の数がこのように少数であっても、これらをｘ1、ｘ2の２点間のすべての点で任意に移動させて最小二乗法により直線で近似する回帰演算は、アルゴリズム自体は短いが、複数の分割点を組み合わせることで、計算量は大幅に増大する。
　そこで、３つの分割点ｐ1、ｐ2、ｐ3を移動させて、図２７の出力データの波形を最も良く表す４つの曲線ａ、ｂ、ｃ、ｄを求めることを考える。図２９は、３つの分割点ｐ1、ｐ2、ｐ3の位置を最終的にどのように決めるかを説明するための図である。まず、図２９（ａ）では、ｘ1と、仮に決められたｐ2の位置を固定する。図２９（ａ）の横軸の目盛りでは、ｘ1は２１０の辺り、ｐ2は２４０の辺りに位置している。この間でｐ1を決められた間隔、例えば横軸の目盛りで１ずつ移動させ、各点において、ｘ1とｐ1とを結ぶ直線ａが、対応する範囲の出力データの波形と最も良くフィットする点としてｐ1の位置を決定する。
　ここで、最もよくフィットする点とは、例えば直線ａと出力データ波形とを最小二乗法により近似する回帰演算を行って、回帰の誤差が最小となるｐ1の点として決定する。ここまでの処理では、横軸の２１０から２４０までの３０点で計算が必要となるが、分割点ｐ1は１つだけなので、ループ処理は１回で済み、コンピュータの負荷はそれほど大きくはない。
　次に、図２９（ｂ）に示すように、決定したｐ1から、仮に決められたｐ3までの間でｐ2を１ずつ移動させ、各点においてｐ1とｐ2とを結ぶ直線ｂが、対応する範囲の出力データの波形と最も良くフィットする点としてｐ2の位置を決定する。さらに、図２９（ｃ）に示すように、決定したｐ2からｘ2までの間でｐ3を１ずつ移動させ、各点においてｐ2とｐ3とを結ぶ直線ｃが、対応する範囲の出力データの波形と最も良くフィットする点としてｐ3の位置を決定する。以上のようにして決定したｐ1、ｐ2、ｐ3を結ぶ直線ａ、ｂ、ｃ、ｄによって得られる近似曲線は、出力データの波形と比較してかなりよく近似された曲線となっている。これにより、直線ａ、ｂ、ｃ、ｄによって得られる近似曲線に基づいて、プレス加工がうまくいったかどうかを判定することができる。しかも、上で説明した処理は、移動させる点が１つだけなので、計算負荷はそれほど高くなく、通常のパーソナルコンピュータでも十分に実行可能である。
　なお、図２８及び図２９に示した処理を一回行っただけでは十分にフィットする曲線が得られない場合には、このような処理を複数回行って、回帰誤差が最も小さくなる回の曲線を採用することもできる。そのような場合であっても、一回の処理にかかるコンピュータに対する負荷及び要する時間が大きくないので、全体としてコンピュータに対する負荷及び要する時間を削減することができる。
　例えば、射出成形では金型内部に圧力センサを設けて、実際に射出成形を行う際に内部の圧力を圧力センサで検知し、その出力データの波形を得ることができる。金型内部の圧力は、液体状の樹脂を注入している段階では高まるが、樹脂の注入が終わり、内部の樹脂が冷えて固まるに従って圧力は徐々に低下し、最終的に金型から成形品を取り出す段階では内部の圧力は急減する。射出成形の工程が適正に行われたかどうかを、このような出力センサの出力データの波形を見ることによって確認することができる。このような場合に、本実施形態のシステムを活用することができる。
　上述の例では、分割点の数はｐ1、ｐ2、ｐ3の３点だったが、このような圧力変化の特性が予め分かっていれば、それに基づいて、分割点の数を増減させてその用途に合う適切な数とすることができる。また、本実施形態を適用できるデータは、時系列のデータであれば、どのようなものにも適用できる。例えば、図２に示した圧力センサの出力データも時系列データであることら、第２実施形態における適切な数の分割点を使ってフィットする曲線を求め、それに基づいて鋼板の切断作業が適切に行われたかどうかを判定する、という用途に本実施形態のシステムを用いることができる。
１０　データ入力部
１２　ユーザ入力処理部
１４　画像出力処理部
１６　バンドパスフィルタ処理部
１８　移動平均処理部
２０　微分処理部
２２　積分階差処理
２４　減算処理部
２６　パルス生成部
２８　閾値処理部
３０　分割点探索処理部
３２　範囲設定処理部

Claims (16)

  1. 　時系列データに対して時間微分を実行し、必要に応じて時間微分後の時系列データの波形図を表示する時間微分処理手段と、
    　処理対象となる時系列データに対し、前記時間微分処理部による時間微分処理を実行するときに時間範囲を時間軸上で設定する範囲設定処理部と、
    　前記範囲設定処理部によって設定された時間範囲において、処理対象となる時系列データがピークとなる位置を検出するピーク位置検出部と、
    を備える、時系列データ解析システム。
  2. 　前記範囲設定処理部は、時系列データの時間軸上に、前記時間範囲の少なくとも一方の境界をユーザが直接設定するものである、請求項１に記載の時系列データ解析システム。
  3. 　前記範囲設定処理部は、ユーザが縦軸の特定の値を設定し、それに対応する時間軸上の値を前記時間範囲の少なくとも一方の境界とするものである、請求項１に記載の時系列データ解析システム。
  4. 　前記範囲設定処理部は、ピーク位置検出部によって検出された時系列データの縦軸のピーク値を検索し、当該ピーク値に対応する時間軸上の値を、前記時間範囲の一方の境界とするものである、請求項１に記載の時系列データ解析システム。
  5. 　さらに、ノイズ除去手段として、時系列データに対して移動平均処理を実行し、必要に応じて移動平均処理後の時系列データの波形図を表示する移動平均処理部、及び、時系列データに対して積分階差処理を実行し、必要に応じて積分階差処理後の時系列データの波形図を表示する積分階差処理部のいずれか一方又は両方を備える、請求項１に記載の時系列データ解析システム。
  6. 　時系列データに対してバンドパスフィルタ処理を実行するバンドパスフィルタ処理部をさらに有する請求項１に記載の時系列データ解析システム。
  7. 　ユーザがユーザデータを入力するユーザ入力処理部をさらに有する請求項５に記載の時系列データ解析システム。
  8. 　外部から時系列データをシステム内に取り込むデータ入力部をさらに有する請求項５に記載の時系列データ解析システム。
  9. 　時系列データに対して時間微分を実行し、必要に応じて時間微分後の時系列データの波形図を表示する時間微分処理手段と、
    　時系列データに対して、その値が連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高とする第１のパルス波形と、その値が連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高とする第２のパルス波形とを求め、必要に応じて第１のパルス波形及び第２のパルス波形の両方又はいずれか一方の波形図を表示するパルス生成部と、
    　複数の時系列データの間又は複数のパルス波形の間で減算処理を実行し、必要に応じてその減算後の時系列データ又はパルス波形を表示する減算処理部と、
    　パルス波形の値を予め指定された閾値と比較し、当該閾値より大きいパルス又は小さいパルスを選択する閾値処理部と、
    　を備える時系列データ解析システム。
  10. 　さらに、ノイズ除去手段として、時系列データに対して移動平均処理を実行し、必要に応じて移動平均処理後の時系列データの波形図を表示する移動平均処理部、及び、時系列データに対して積分階差処理を実行し、必要に応じて積分階差処理後の時系列データの波形図を表示する積分階差処理部のいずれか一方又は両方を備える、請求項９に記載の時系列データ解析システム。
  11. 　時系列データに対してバンドパスフィルタ処理を実行するバンドパスフィルタ処理部をさらに有する請求項９に記載の時系列データ解析システム。
  12. 　ユーザがユーザデータを入力するユーザ入力処理部をさらに有する請求項９に記載の時系列データ解析システム。
  13. 　外部から時系列データをシステム内に取り込むデータ入力部をさらに有する請求項９に記載の時系列データ解析システム。
  14. 　時系列データ解析システムに取り込まれた時系列データの屈折点を求める方法であって、
    　前記時系列データ解析システムに含まれる微分処理部によって、取り込まれた第１の時系列データに対して３回の時間微分を行って第２の時系列データを得る第１ステップと、
    　前記時系列データ解析システムに含まれるパルス生成部によって、前記第２の時系列データが連続して正である時間の長さをパルス幅、連続して正である時間の長さをパルス高とする第１のパルス波形と、前記第４の時系列データが連続して負である時間の長さをパルス幅、連続して負である時間の長さをパルス高とする第２のパルス波形とを求める第２ステップと、
    　前記時系列データ解析システムに含まれる減算処理部によって、前記第１のパルス波形と前記第２のパルス波形との減算を行って第３のパルス波形を求める第３ステップと、
    　前記微分処理部によって、前記第３のパルス波形に対して時間微分を行って第４のパルス波形を求める第４ステップと、
    　前記時系列データ解析システムに含まれる閾値処理部によって、前記第４のパルス波形に対して閾値処理を行い、所定の閾値を超えるパルスのみを選択して、前記第１の時系列データの屈折点とする第５ステップと、
    　を含む、時系列データ解析方法。
  15. 　前記第１ステップで、２回目の時間微分を行って得られる時系列データに対して積分階差を行ってノイズを除去したあとに３回目の時間微分を実行して第２の時系列データとする、請求項１４に記載の時系列データ解析方法。
  16. 　時系列データの波形図の近似曲線を求めることによって時系列データを解析する時系列データ解析方法であって、
    　予め指定された始点と終点との間の前記時系列データ上に第１分割点から第ｎ分割点（ｎは２以上の整数）までを仮に配置するステップと、
    　第１分割点を始点から第２分割点まで移動させたときに始点と第１分割点とを結ぶ第１直線がその範囲で時系列データと最も良くフィットする点を第１分割点の位置と決定し、
    　第２分割点を第１分割点から第３分割点まで移動させたときに第１分割点と第２分割点とを結ぶ第２直線がその範囲で時系列データと最も良くフィットする点を第２分割点の位置と決定するステップと、
    　以下同様の手順を続行し、第ｎ分割点を第ｎ－１分割点から終点まで移動させたときに第ｎ－１分割点と第２分割点とを結ぶ第ｎ直線がその範囲で時系列データと最も良くフィットする点を第ｎ分割点の位置と決定するステップと、
    　前記始点と、第１分割点と、・・・第ｎ分割点と、終点とをこの順番で直線によって連結して得られる曲線を前記時系列データの近似曲線とするステップと、
    　を含む時系列データ解析方法。
PCT/JP2023/027035 2022-08-01 2023-07-24 時系列データ解析システム WO2024029387A1 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2022122546A JP2024019830A (ja) 2022-08-01 2022-08-01 時系列データ解析システム
JP2022-122546 2022-08-01

Publications (1)

Publication Number Publication Date
WO2024029387A1 true WO2024029387A1 (ja) 2024-02-08

Family

ID=89848944

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2023/027035 WO2024029387A1 (ja) 2022-08-01 2023-07-24 時系列データ解析システム

Country Status (2)

Country Link
JP (1) JP2024019830A (ja)
WO (1) WO2024029387A1 (ja)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005216144A (ja) * 2004-01-30 2005-08-11 Omron Corp センサシステム、処理装置および拡張ユニット
WO2005111641A1 (ja) * 2004-05-13 2005-11-24 Mitsubishi Denki Kabushiki Kaisha 状態把握装置およびこの状態把握装置を使用した電力開閉機器の開閉制御装置
JP2008310723A (ja) * 2007-06-18 2008-12-25 Ono Sokki Co Ltd 時系列データ処理システム及びコンピュータプログラム
JP2018069026A (ja) * 2016-10-20 2018-05-10 パナソニックＩｐマネジメント株式会社 脈波計測装置、及び、脈波計測方法
JP2018159618A (ja) * 2017-03-23 2018-10-11 中山水熱工業株式会社 波形分析補助装置、及び波形分析補助システム
JP2022108079A (ja) * 2021-01-12 2022-07-25 日本光電工業株式会社 生体情報処理装置、生体情報処理方法、プログラム及び記憶媒体

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005216144A (ja) * 2004-01-30 2005-08-11 Omron Corp センサシステム、処理装置および拡張ユニット
WO2005111641A1 (ja) * 2004-05-13 2005-11-24 Mitsubishi Denki Kabushiki Kaisha 状態把握装置およびこの状態把握装置を使用した電力開閉機器の開閉制御装置
JP2008310723A (ja) * 2007-06-18 2008-12-25 Ono Sokki Co Ltd 時系列データ処理システム及びコンピュータプログラム
JP2018069026A (ja) * 2016-10-20 2018-05-10 パナソニックＩｐマネジメント株式会社 脈波計測装置、及び、脈波計測方法
JP2018159618A (ja) * 2017-03-23 2018-10-11 中山水熱工業株式会社 波形分析補助装置、及び波形分析補助システム
JP2022108079A (ja) * 2021-01-12 2022-07-25 日本光電工業株式会社 生体情報処理装置、生体情報処理方法、プログラム及び記憶媒体

Also Published As

Publication number Publication date
JP2024019830A (ja) 2024-02-14

Similar Documents

Publication Publication Date Title
US5586041A (en) Method and system for real-time statistical process monitoring
US5311759A (en) Method and system for real-time statistical process monitoring
JP6420885B1 (ja) 電磁振動成分の除去方法、回転機械診断方法、及び回転機械診断装置
US20040059760A1 (en) Waveform detector and state monitoring system using it
CN111898443B (zh) 一种fdm型3d打印机送丝机构流量监测方法
CN117150217B (zh) 一种基于大数据的数据压缩处理方法
US6721445B1 (en) Method for detecting anomalies in a signal
WO2024029387A1 (ja) 時系列データ解析システム
JP2018149183A5 (ja)
KR102026069B1 (ko) 반도체 설비의 센서 데이터 분할 시스템 및 그 방법
CN111176226A (zh) 一种基于运行工况的设备特征参数报警阈值自动分析方法
CN112474815B (zh) 一种控制轧制过程的方法及装置
JP2018153269A5 (ja)
US6309571B2 (en) Method and apparatus for the control of injection molding
CN111538755B (zh) 一种基于归一化互相关与单位根检验的设备运行状态异常检测方法
JP6794283B2 (ja) 生産工程分析装置及びこれを用いる生産管理システム
JP5751606B2 (ja) 機械設備における異常診断システム
WO2020124934A1 (zh) 一种伺服电机负载惯量的测定方法
Manhertz et al. Evaluation of short-time fourier-transformation spectrograms derived from the vibration measurement of internal-combustion engines
CN113531838B (zh) 一种智能化调整空间温度的制冷控制方法及系统
JP7255242B2 (ja) 情報処理装置、情報処理方法、プログラム、及び異常診断装置
US4334282A (en) Testing apparatus
JPH04276539A (ja) 回転機械の損傷・劣化の寿命予測方法
CN109798970A (zh) 异常检测装置、异常检测方法、异常检测系统及存储介质
CN108830845B (zh) 一种慢压射过程稳定性定量分析方法

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: 23849942

Country of ref document: EP

Kind code of ref document: A1