CN105607477A - Industrial control circuit oscillation detection method based on improved local mean value decomposition - Google Patents

Abstract

The invention discloses an industrial control loop oscillation detection method based on improved local mean value decomposition. The method comprises the following steps: in a control loop to be detected, pre-acquiring a group of process historical data; by use of a similar waveform matching method, performing end point extension on process output data; performing the improved local mean value decomposition on the process output data after the end point extension; calculating an autocorrelation function zero crossing point regularity index of each decomposition subsignal; and determining whether each regularity index exceeds a set threshold, and according to all determination results, obtaining a detection result. By using the method provided by the invention, time variation and multi-period oscillation behaviors of an industrial control loop can be quantitatively detected, such components as time-variation oscillation, multiple oscillation, intermittent oscillation, non-stationary signals and the like can be distinguished, at the same time, at the same time, the rule degree and period of each oscillation component can be obtained, and rich data support is provided for evaluation and fault source diagnosis of oscillation behaviors.

Description

A kind of Industry Control oscillation circuit detection method based on improving local mean value decomposition

Technical field

The present invention relates to the Performance Evaluation field in Industry Control, be specifically related to a kind of local based on improvingThe Industry Control oscillation circuit detection method that average is decomposed.

Background technology

Modern industry process device have scale greatly, comprehensively spend high, manipulation is complicated, variable is many, and longTime operates in the inferior feature of closed-loop control. The chemical process that industry is common, often comprises thousands ofUp to ten thousand control loops, and these control loops are because coupled relation interacts. Due to workThe characteristics such as industry control loop middle controller is crossed and adjusted, external disturbance and control valve nonlinear operation generalExist, the oscillation behavior of control loop occurs often, and this has greatly affected the operation of industrial flow equipmentEconomic benefit and stability.

Industrial flow equipment is carried out to oscillation test tentatively accurately and can reduce waste product output, reduceDisqualification rate, improves reliability, security in industrial flow equipment running process, reduces system simultaneouslyCause this. Many controllers can also keep good performance at initial operating stage, but As time goes on,Due to the impact of external interference factor or equipment self problem, the performance of controller can reduce even graduallyLost efficacy. Be embodied in control loop process all kinds of oscillation behaviors occur, wherein may comprise multiple shakingSwing, intermittent oscillation, the composition such as non-linear, thereby threaten the safe and stable operation of industrial process.

Meanwhile, because apparatus of load in real time environment and operating mode often change, industrial process is also shownReveal the one side of Non-stationary Data characteristic, be embodied in the local mean value transport phenomena of process data.For important control loop, find that in time the oscillating characteristic of its running contributes to engineering staff to enterRow fault diagnosis and investigation. Therefore,, in industrial control system Performance Evaluation process, design effectivelyOn-line monitoring means, in time, accurately detect all kinds of vibrations of non-stationary process data in control loopComposition, and distinguish different separately frequency ranges, for controller performance assessment and control loop eventBarrier is diagnosed all important in inhibitings.

Existing Industry Control oscillation circuit detection technique, the overwhelming majority is based on stationary process dataAnalytical method. In recent two decades, some oscillation test for non-stationary process data are also there areMethod. Can roughly be summarized as three kinds by its main thought: based on the analysis side of process data Time-domain StatisticsMethod; The analytical method of the auto-correlation function (ACF) based on process data; And based on process dataSignal decomposition method (comprise the conversion of empirical mode decomposition EMD and base decompose).

Detection method based on process data Time-domain Statistics or auto-correlation function domain analysis is in commercial ApplicationIn have three shortcomings: one, the method need to have certain priori to loop to be detected or process, certainA little parameters are also to determine according to experience; Two, the industrial process to non-stationary and many cycles of oscillation of existence,Cannot realize the full-automatic nothing intervention and detect, need to design wave filter targetedly and carry out data tranquilization placeReason separates with oscillationg component; Three, most detection algorithms cannot quantitatively calculate the regular degree of oscillationg component.

Signal decomposition method based on process data is compared above-mentioned detection method and is had progress at present, stillIts limitation is mainly reflected in: the subsignal number of plies redundancy that existing signal decomposition technical point solution obtains is various,Cause many subsignals to lack the support of actual physics meaning, do not there is good representativeness, and theseMethod is also poor to the degree of fitting of non-stationary signal trend, and computation complexity is also higher. In addition,Method based on decomposition technique can not process data discontinuous vibration, time the composition such as variable oscillation.

In the practical application of process oscillation test algorithm, can effectively detect whether tool of Industry Control loopThere is oscillation behavior, and a rule degree index of qualitative assessment oscillation behavior, and become while being generally applicable to existThe process data of vibration, intermittent oscillation, non-stationary and non-linear component, for Accurate Diagnosis industry mistakeThe existence of journey vibration has very important Practical significance, is also conducive to determining of industrial process control performanceAmount assessment.

Summary of the invention

The invention provides a kind of Industry Control oscillation circuit based on improving local mean value decomposition detectsMethod, the Industry Control circuit process of the behavior such as variable oscillation, multicycle vibration can be applicable to exist time,Detection method is generally applicable to non-stationary or process data stably, only need obtain conventional operation data,Without process mechanism knowledge.

Based on an Industry Control oscillation circuit detection method of improving local mean value decomposition, comprising:

Step 1, the output of process signal of one group of control loop to be detected of collection;

Step 2, utilizes similar waveform matching process, and the output of process signal is carried out to end points continuation;

Step 3, carries out improved local mean value decomposition to the output of process signal after continuation, obtains pointSeparate subsignal;

Step 4, calculates the regular index of auto-correlation function zero cross point of each decomposition subsignal(Regularityofauto-covariancefunctionzero-crossingintervals)；

Step 5, judges whether the regular index of each auto-correlation function zero cross point exceedes threshold value, if superCross threshold value, in the decomposition subsignal that control loop is corresponding, have concussion.

If there are multiple decomposition subsignals to have oscillation behavior in the output of process signal gathering, judgementThere is multicycle oscillation behavior in this control loop.

Accuracy in detection and the reliability of the behavior such as change, multicycle vibration when the present invention can improve,The aspect of increasing economic efficiency has important practical value.

The present invention directly adopts the measurable variable of chemical process as the output of process signal, this output of processSignal obtains by field real-time acquisition, and along with passage of time, constantly gather and renewal process defeatedGo out signal to monitoring system. First collected Process History data are carried out to end points continuation processing,Then utilize improved local mean value to decompose and obtain decomposing subsignal set { x^k, calculate each and decompose sonSignal x^kThe regular index η of corresponding auto-correlation function zero cross point^k，η^kComputation complexity minimum,Also can carry out large batch of multi-group data simultaneously.

As preferably, in step 2, get one section of waveform, searching process at the left end of the output of process signalIn output signal, with the left side waveform of the highest waveform of this section of Waveform Matching degree, utilize left side waveform to mistakeThe left end of journey output signal carries out continuation; Get one section of waveform at the right-hand member of the output of process signal, foundIn journey output signal, with the right waveform of the highest waveform of this section of Waveform Matching degree, utilize the right waveform pairThe right-hand member of the output of process signal carries out continuation.

As preferably, in step 3, while carrying out improved local mean value decomposition, adopt adaptive windowIt is average that mouth size Selection strategy carries out sliding window, adopts the judgement form of two norms to adjust as pure frequencyThe criterion of signal processed.

Adaptive window size selection strategy refers to, in each iterative process, window size all canTo be adjusted according to the feature of active procedure output signal, this strategy allows longer sequential to criticizeProcess, preferably, adaptive window size computing formula is as follows:

$W=\frac{T_{max}}{{\mathrm{\σ}}_{T_i}+C}+\frac{S_{\mathrm{\ω}}}{S_{\mathrm{\ω}}+{10}^{n}C}{\mathrm{\μ}}_{T_i};$

Wherein, W represents the window size of moving average, T_iRepresent in the output of process signal all continuousThe time interval of extreme point, T_max＝max{T_i}，Represent respectively T_iAverage,Represent T_iStandardPoor, S_ωFor the data length of active procedure output signal, n is S_ωThe order of magnitude, C is a constant.C default value is 3.

In iterative process, the criterion of pure frequency modulated signal refers to, for the modulation after iteration n timeSignal S_kn(t), judge its condition as pure frequency adjusted signal, preferably, in step 3, pure frequencyThe criterion of rate modulation signal is s_kn(t) local envelope function meets with lower inequality:

$norm\{a_{k(n+1)}-1\}\≤\sqrt{length\left\{x\left(t\right)\right\}}*\mathrm{\δ};$

Wherein, s_kn(t) be the modulation signal after iteration n time; a_k(n+1)For s_kn(t) local envelope function,Norm{} represents to get two norms, length{x (t) } be the data length of active procedure output signal x (t), δ isConstant, span is δ=0.001～0.01.

As preferably, in step 4, calculate the auto-correlation function zero cross point rule of each decomposition subsignalThe step of property index is as follows:

Step 4-1, calculates each decomposition subsignal x^kAuto-correlation function, k decomposes the numbering of subsignal;

Step 4-2, the interval in statistics auto-correlation function between all continuous zero cross pointsRespectivelyCounting periodAverageAnd standard deviation

Step 4-3, utilizes following formula to calculate the regular index η of auto-correlation function zero cross point^k：

${\mathrm{\η}}^k=\frac{1}{3}\×\frac{{\mathrm{\μ}}_{Z_i}^k}{{\mathrm{\σ}}_{Z_i}^k}.$

As preferably, the threshold value in step 5 is 1.

The beneficial effect that the present invention compared with prior art has:

1, without outside additional signal excitation, also can not introduce additional disturbance to control system, canRealize detection and the diagnosis of non-intrusion type.

2, computation complexity is low, convenient operation, and algorithm is write simple and easy, is beneficial in existing DCS workStand or control system host computer on implement.

3, the signal decomposition method adopting has been realized the automatic separation of non-stationary component in process data,Than existing other decomposition techniques, decomposition efficiency is higher, and computation complexity is lower.

4, the moving average window size selection strategy proposing has adaptive characteristic, changes simultaneouslyEntering stopping criterion for iteration is later more suitable in process control loop oscillation test.

5, can to Industry Control loop time become, multicycle oscillation behavior carries out quantizating index detection,For assessment and the source of trouble diagnosis of loop performance to be detected provide abundant Data support.

6, adopt the method for data driven type completely, without process priori, without design in advanceWave filter, does not need to carry out manual intervention yet.

Brief description of the drawings

Fig. 1 is the schematic flow sheet of chemical process in embodiment;

Fig. 2 is the output of process signal of one group of furnace outlet temperature control loop gathering in embodiment;

Fig. 3 is the data of the output of process signal after the continuation of similar waveform matching process end points in embodiment;

Fig. 4 is the schematic diagram of the output of process signal after local mean value is decomposed in embodiment;

Fig. 5 is method flow diagram of the present invention.

Detailed description of the invention

Comment with the performance of main heating furnace in domestic certain large petrochemical plant delayed coking production process belowEstimate for example, the time variable oscillation behavior detection method of the chemical process that has control valve viscosity property is done in detailThin description.

As shown in Figure 1, petrochemical process heating furnace is important step and the main energy consumption list in production procedureOne of unit, the steady control of heater outlet temperature has important meaning for improving product quality and reducing energy consumptionJustice.

Heating furnace is supplied heat-obtaining by device in Gas, and gas amount changes and fluctuates according to upstream oiliness, needsControl air intake and make device in Gas fully burn to obtain maximum heat, should ensure certain sky simultaneouslyGas surplus, but too much Cryogenic air can be taken away furnace heat, causes waste of fuel, the economic effect of lossBenefit, therefore, using furnace outlet temperature as controlled variable, fuel device in Gas aperture becomes as operationAmount is carried out circuit controls, and process exists random perturbation simultaneously.

Device in Gas degree adjustment valve (control valve) belongs to the executing agency of this control loop, moves one sectionAfter time, occur certain nonlinear characteristic, because controller is crossed the reason such as adjust, control loop is easyThere is persistent oscillation behavior. And external disturbance also can be introduced this loop by coupling circuit, leadsCause loop and produce other hunting of frequencys.

As shown in Figure 5, a kind of Industry Control oscillation circuit detection side based on improving local mean value decompositionMethod, comprising:

Step 1, the output of process signal of one group of control loop to be detected of collection.

The method of gatherer process output signal is within default each sampling interval, to record to be detectedControl loop in process data, and the process data collecting in each sampling interval is added onThe process data end previously gathering.

Sampling interval refers to the sampling interval of performance evaluation system. Process data is along with passage of time is continuousUpgrade, every time span through a sampling interval, all has new process data to add to previously and adoptsThe end of the process data of collection. In general and industrial control system of sampling interval of performance evaluation systemControl cycle is identical, also can be chosen as the integral multiple of control cycle, specifically according to performance monitoring and workThe requirement of real-time at industry scene and memory data output restriction are determined.

The output of process signal that the present embodiment gathers is in device in Gas control valve viscous situation, between having againFurnace outlet temperature data when the external disturbance of having a rest is introduced. Furnace outlet temperature after centralizationAs shown in Figure 2, in Fig. 2, abscissa is sampled point ordinal number to degrees of data, and unit is Sample (1 SampleThe sampling interval of corresponding data), ordinate is heating furnace under the nominal situation after centralizationOutlet temperature, unit is DEG C.

Step 2, utilizes similar waveform matching process, and the output of process signal is carried out to end points continuation.

Original procedure output signal is carried out to left and right end points continuation processing, and object is to eliminate input signal pairThe end effect impact of decomposition method, its detailed description of the invention is:

Step 2-1, with the left end point x (t of the output of process signal₀) be starting point, get one of original signal to the rightPortion waveshape w (t₀:t₀+l)＝[x(t₀),x(t₀+1),…x(t₀+ l)], wherein l meets w (t₀:t₀+ l) onlyThe maximum that comprises a zero cross point;

Step 2-2, finds w (t₀:t₀+ mid point l)Wherein[] representsGet toward the nearest integer of positive infinity;

Step 2-3, to right search in next sequence x (t) withBe worth identical point, be designated asx(t₁), with x (t₁) be one section of intermediate point extraction and w (t₀:t₀+ l) the identical waveform of length, is designated asw(t₁-[l/2]:t₁+[l/2]-1)；

Step 2-4, calculates w (t₀:t₀+ l) with w (t₁-[l/2]:t₁+ [l/2]-1) Waveform Matching degree m₁。Wherein Waveform Matching degree can according to prior art " Zhu Xiaojun, Lv Shiqin, Wang Yanfei, etc. improvedLMD algorithm and the application in EEG signal characteristic abstraction [J] thereof. Institutes Of Technology Of Taiyuan's journal, 2012,43 (3): 339-343. " calculate.

Step 2-5, repeating step 2-3 and step 2-4, until the output of process signal search finishes, like thisCan obtain a series of matching degree m=[m₁,m₂,…m_n]；

Step 2-6, finds minimum of a value m in m_b, corresponding wave band w (t_b-[l/2]:t_b+[l/2]-1)Be (the t with w₀:t₀+ l) a section of mating most of waveform, now the left end of signal x (t) is availablew(t_b-l-[l/2]:t_b-[l/2]) carry out continuation;

Step 2-7, adopts the right endpoint using the same method to signal x (t) and carries out continuation.

The data of the present embodiment after to the continuation of the output of process signal as shown in Figure 3.

Step 3, carries out improved local mean value decomposition to the output of process signal after continuation, obtains pointSeparate subsignal.

The present invention decomposes and improves local mean value of the prior art, has retained former methodical instituteThere are mathematics and calculated characteristics, just at moving average window size selection strategy and pure frequency modulated signalCriterion (end condition) on modify, improved decomposition method is for same process data,Than former method, capacity of decomposition is stronger, and still less, iteration time is faster for the subsignal quantity of acquisition,Be more suitable in analyzing original signal oscillation behavior. Original local mean value is decomposed can be according to prior art“SmithJS.ThelocalmeandecompositionanditsapplicationtoEEGperceptiondata[J].JournaloftheRoyalSocietyInterface,2005,2(5):443-454. " carry out.

During improved local mean value is decomposed, the selection strategy of moving average window size refers to, at every turnIn iterative process, window size all can be adjusted according to the feature of current demand signal, simultaneously this strategyAllow longer sequential to carry out batch processing. Its concrete formula is as follows:

$W=\frac{T_{\max }}{{\mathrm{\σ}}_{T_i}+C}+\frac{S_{\mathrm{\ω}}}{S_{\mathrm{\ω}}+{10}^{n}C}{\mathrm{\μ}}_{T_i};$

Wherein, W represents the window size of moving average, T_iRepresent in sequence between all continuous threshold pointsThe time interval, T_max＝max{T_i}，WithRepresent respectively T_iAverage and standard deviation, S_ωFor currentSequence data length, n is S_ωThe order of magnitude, C is a constant, default value is 3.

S_ω/(S_ω+10ⁿC) be the contraction factor of formula, it has ensured that the recruitment of W can not exceedOtherwise window size will be excessive. From formula, also can find out, the selection of window size is also examined with currentThe data length S considering_ωRelevant, S_ωLess, W is less. One largerT is often described_iChange brightAobvious, now must select less S_ωTo ensure less T_iBe not left in the basket, on the contrary, one largerExplanation can be selected larger S_ωTo ensure normally carrying out of moving average. C is a constant,Default value is 3, but in the time analyzing fast-changing signal, better effect can be obtained in 3≤C≤5.

During improved local mean value is decomposed, in iterative process, the criterion of pure frequency modulated signal isRefer to, for the modulation signal s after iteration n time_kn(t) condition that, determines that it is pure frequency modulated signal iss_kn(t) local envelope function meets with lower inequality:

$norm\{a_{k(n+1)}-1\}\≤\sqrt{length\left\{x\left(t\right)\right\}}*\mathrm{\δ};$

Wherein, a_k(n+1)For s_kn(t) local envelope function, norm{} represents to get two norms, length{x (t) }For the data length of current sequence x (t), δ is a constant between 0 and 1, its common valueScope is δ=0.001～0.01.

The process data x that the present embodiment obtains continuation carries out improved local mean value decomposition, obtains pointSeparate subsignal arrangement set { x¹,x²,x³, as shown in Figure 4.

Step 4, calculates the regular index of auto-correlation function zero cross point of each decomposition subsignal, concreteStep is as follows:

Step 4-1, calculates each decomposition subsignal x^kAuto-correlation function, k decomposes the numbering of subsignal;

Step 4-2, the interval in statistics auto-correlation function between all continuous zero cross pointsRespectivelyCounting periodAverageAnd standard deviation

Step 4-3, utilizes following formula to calculate the regular index η of auto-correlation function zero cross point^k：

${\mathrm{\η}}^k=\frac{1}{3}\×\frac{{\mathrm{\μ}}_{Z_i}^k}{{\mathrm{\σ}}_{Z_i}^k}.$

The regular index η of correlation function zero cross point^kCan according to prior art " ThornhillNF,HuangB,ZhangH.Detectionofmultipleoscillationsincontrolloops[J].JournalofProcessControl, 2003,13 (1): 91-100. " calculate. Correlation function zero cross point ruleThe Analysis of Deep Implications of rule property index is that, for standard oscillator signal, its all waveforms that comprise are all answeredThere is identical wavelength span, therefore η^k→ ∝, and in actual motion, due to environment and measurement mistakeThe unfavorable factor impacts such as difference, the method regulation η^k> 1 can judge in original signal and greatly may shake in additionSwing behavior signal. In addition, the regular index of correlation function zero cross point is also subject to zero cross point in signalNumber impact, if zero cross point number is very few in original signal, the estimation of average, standard deviation all mayThere is relatively large deviation, therefore the present invention's regulation, in former the output of process signal, all-wave number must be more than or equal to 11。

In the present embodiment, x³Comprise zero cross point number lower than 11, calculating correlation function zero crossingWhen the regular index of some, this is cast out, calculate remaining x¹And x²Corresponding auto-correlation function zeroRegular index { the η in crosspoint¹,η²}。

Two regular indexs of auto-correlation function zero cross point of decomposing subsignal are respectively{η¹＝0.23,η²=4.70}. Obviously decompose subsignal x²Obtained correlation function zero cross point ruleProperty index, considerably beyond given threshold value Ω=1, can be decomposed subsignal x certainly²For oscillating component. And it is rightDecompose subsignal x in ground floor¹, its value of statistical indicant, for being far smaller than given threshold value, should be thought notThere is the subsignal of vibration, as shown in Figure 4.

Step 5, if the regular index η of one of them correlation function zero cross point^kExceed threshold value Ω,Judge the decomposition subsignal x that this control loop is corresponding^kThere is vibration, if in the process data gatheringThere are multiple decomposition subsignals to have oscillation behavior, judge that this control loop exists multicycle vibration rowFor. Described in the present embodiment, threshold value Ω is 1, works as η^k> 1, illustrates x^kIn there is oscillation behavior.

Utilize the inventive method, carrying out on the basis of single vibration, multicycle oscillation test, can alsoEnough time variable oscillation behaviors to Industry Control loop quantitatively detect, the rule of variable oscillation component when acquisitionDegree and the cycle, for evaluation and the source of trouble diagnosis of oscillation behavior provide abundant Data support.

Claims (7)

1. the Industry Control oscillation circuit detection method based on improving local mean value decomposition, its spyLevy and be, comprising:

Step 1, the output of process signal of one group of control loop to be detected of collection;

Step 2, utilizes similar waveform matching process, and the output of process signal is carried out to end points continuation;

Step 3, carries out improved local mean value decomposition to the output of process signal after continuation, obtains pointSeparate subsignal;

Step 4, calculates the regular index of auto-correlation function zero cross point of each decomposition subsignal;

Step 5, judges whether the regular index of each auto-correlation function zero cross point exceedes threshold value, if superCross threshold value, in the decomposition subsignal that control loop is corresponding, have vibration.

2. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 1Survey method, is characterized in that, in step 2, gets one section of waveform at the left end of the output of process signal, seeksThe left side waveform of looking in the output of process signal the waveform the highest with this section of Waveform Matching degree, utilizes left side rippleShape is carried out continuation to the left end of the output of process signal; Get one section of waveform at the right-hand member of the output of process signal,The right waveform of finding in the output of process signal the waveform the highest with this section of Waveform Matching degree, utilizes the rightWaveform carries out continuation to the right-hand member of the output of process signal.

3. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 2Survey method, is characterized in that, in step 3, while carrying out improved local mean value decomposition, adopts adaptiveIt is average that the window size selection strategy of answering carries out sliding window, adopts the judgement form of two norms as pureThe criterion of frequency modulated signal.

4. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 3Survey method, is characterized in that, adaptive window size computing formula is as follows:

$W=\frac{T_{\max }}{{\mathrm{\σ}}_{T_i}+C}+\frac{S_{\mathrm{\ω}}}{S_{\mathrm{\ω}}+{10}^n}{\mathrm{\μ}}_{T_i};$

Wherein, W represents the window size of moving average, T_iRepresent in the output of process signal all continuousThe time interval of extreme point, T_max＝max{T_i}，Represent respectively T_iAverage,Represent T_iStandardPoor, S_ωFor the data length of active procedure output signal, n is S_ωThe order of magnitude, C is a constant.

5. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 4Survey method, is characterized in that, in step 3, the criterion of pure frequency modulated signal is s_kn(t) officePortion's envelope function meets with lower inequality:

$norm\{a_{k(n+1)}-1\}\≤\sqrt{length\left\{x\left(t\right)\right\}}*\mathrm{\δ};$

Wherein, s_kn(t) be the modulation signal after iteration n time; a_k(n+1)For s_kn(t) local envelope function,Norm{} represents to get two norms, length{x (t) } be the data length of active procedure output signal x (t), δ isConstant, span is δ=0.001～0.01.

6. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 5Survey method, is characterized in that, in step 4, calculates the auto-correlation function zero crossing of each decomposition subsignalThe step of the regular index of point is as follows:

Step 4-1, calculates each decomposition subsignal x^kAuto-correlation function, k decomposes the numbering of subsignal;

Step 4-2, the interval in statistics auto-correlation function between all continuous zero cross pointsRespectivelyCounting periodAverageAnd standard deviation

Step 4-3, utilizes following formula to calculate the regular index η of auto-correlation function zero cross point^k：

${\mathrm{\η}}^k=\frac{1}{3}\×\frac{{\mathrm{\μ}}_{Z_i}^k}{{\mathrm{\σ}}_{Z_i}^k}.$

7. the Industry Control oscillation circuit inspection based on improving local mean value decomposition as claimed in claim 6Survey method, is characterized in that, the threshold value in step 5 is 1.