Distributed optical fiber sensing signal processing method for safety monitoring of underground pipe network

Abstract

A distributed optical fiber sensing signal processing method for safety monitoring of underground pipe network, which belongs to infrastructure safety monitoring field, which is aimed to improve the intelligent ability of detection and identification of the existing distributed optical fiber sound/vibration sensing system under complex application conditions. The present invention utilizes the distributed optical fiber sound/vibration sensing system to pick up the sound or vibration signal of the whole line along the detection cable; and the customized short time feature and long time feature are respectively extracted from the relative quantity of the sound or the vibration signal at each spatial point in the whole monitoring range. The Bayesian identification and classification network at each spatial point is constructed and trained based on the prior knowledge of the collected signal features and their different background noises.

Claims

1. A distributed optical fiber sensing signal processing method for safety monitoring of an underground pipe network, comprising steps of: Step 1, picking up sound or vibration signal at each point along a line of a detection cable laid within a monitoring range; Step 2, extracting a customized short time feature according to a relative value of the sound or the vibration signal collected in a whole monitoring area at each spatial point; and Step 3, detecting an abnormal event occurs or not at each spatial point according to the short time feature, wherein then a long time feature is calculated and accumulated for the signal of the abnormal event, then values of these features are discretized after calculation and accumulation, then these discrete values of features are input into a Bayesian identification and classification network which is trained for each spatial point to realize identification and classification; wherein in the Step 3, detailed steps of construction of the Bayesian identification and classification network which is trained for the each spatial point are: Step 31, accumulating the sound or the vibration signal and then preprocessing it, wherein then a typical event data training set is built based on a long time signal which is obtained from pretreatment; Step 32, extracting the customized short time feature from the relative value of the signal for the typical event data training set, wherein the long time feature is calculated and accumulated, after the extraction, then the typical event data training set is transformed into a long time feature training set; and Step 33, constructing the Bayesian identification classification network topology and learning the parameter based on the probability distribution information of the long time feature in the long time data training set and combined with prior probability information of a background environment at each spatial point, thus obtaining the trained Bayesian identification classification network; wherein detailed steps of the Step 32 are: Step 321, extracting four customized short time features based on a relative quantity of the signals from the short time signal: a primary impulse intensity feature A.sub.1, a secondary impulse intensity feature A.sub.2, a primary and secondary impulse amplitude ratio feature A.sub.3, and a wavelet packet energy spectrum feature A.sub.4; Step 322, calculating and accumulating the short time features to get six customized long time features, which include a long time primary impulse intensity feature B.sub.1, a long time secondary impulse intensity feature B.sub.2, a long time amplitude ratio feature B.sub.3, a long time wavelet packet energy spectrum feature B.sub.4, a long time environmental background feature B.sub.5 and a long time strict impulse number feature B.sub.6; and Step 323, according to the long time features and corresponding event type tags thereof, transforming the typical event data training set into the long time feature training set.

2. A distributed optical fiber sensing signal processing method for safety monitoring of an underground pipe network, comprising steps of: Step 1, picking up sound or vibration signal at each point along a line of a detection cable laid within a monitoring range; Step 2, extracting a customized short time feature according to a relative value of the sound or the vibration signal collected in a whole monitoring area at each spatial point; and Step 3, detecting an abnormal event occurs or not at each spatial point according to the short time feature, wherein then a long time feature is calculated and accumulated for the signal of the abnormal event, then values of these features are discretized after calculation and accumulation, then these discrete values of features are input into a Bayesian identification and classification network which is trained for each spatial point to realize identification and classification; wherein detailed steps of the Step 3 are: Step 3.1, detecting abnormality of the sound or the vibration signal at each spatial point according to the short time feature, when threshold is given, and; Step 3.2, if the abnormal event is detected at certain spatial point, calculating and accumulating the long time feature at this point from a moment when the abnormal event occurs, then the values of these features are discretized after calculation and accumulation; and Step 3.3, inputting discretization results into the Bayesian identification and classification network which is trained for the corresponding spatial point to execute identification and classification; wherein in the Step 3, detailed steps of construction of the Bayesian identification and classification network which is trained for the each spatial point are: Step 31, accumulating the sound or the vibration signal and then preprocessing it, wherein then a typical event data training set is built based on a long time signal which is obtained from pretreatment; Step 32, extracting the customized short time feature from the relative value of the signal for the typical event data training set, wherein the long time feature is calculated and accumulated, after the extraction, then the typical event data training set is transformed into a long time feature training set; and Step 33, constructing the Bayesian identification classification network topology and learning the parameter based on the probability distribution information of the long time feature in the long time data training set and combined with prior probability information of a background environment at each spatial point, thus obtaining the trained Bayesian identification classification network; wherein detailed steps of the Step 32 are: Step 321, extracting four customized short time features based on a relative quantity of the signals from the short time signal: a primary impulse intensity feature A.sub.1, a secondary impulse intensity feature A.sub.2, a primary and secondary impulse amplitude ratio feature A.sub.3, and a wavelet packet energy spectrum feature A.sub.4; Step 322, calculating and accumulating the short time features to get six customized long time features, which include a long time primary impulse intensity feature B.sub.1, a long time secondary impulse intensity feature B.sub.2, a long time amplitude ratio feature B.sub.3, a long time wavelet packet energy spectrum feature B.sub.4, a long time environmental background feature B.sub.5 and a long time strict impulse number feature B.sub.6; and Step 323, according to the long time features and corresponding event type tags thereof, transforming the typical event data training set into the long time feature training set.

3. A distributed optical fiber sensing signal processing method for safety monitoring of an underground pipe network, comprising steps of: Step 1, picking up sound or vibration signal at each point along a line of a detection cable laid within a monitoring range; Step 2, extracting a customized short time feature according to a relative value of the sound or the vibration signal collected in a whole monitoring area at each spatial point; and Step 3, detecting an abnormal event occurs or not at each spatial point according to the short time feature, wherein then a long time feature is calculated and accumulated for the signal of the abnormal event, then values of these features are discretized after calculation and accumulation, then these discrete values of features are input into a Bayesian identification and classification network which is trained for each spatial point to realize identification and classification; wherein detailed steps of the Step 3 are: Step 3.1, detecting abnormality of the sound or the vibration signal at each spatial point according to the short time feature, when threshold is given, and; Step 3.2, if the abnormal event is detected at certain spatial point, calculating and accumulating the long time feature at this point from a moment when the abnormal event occurs, then the values of these features are discretized after calculation and accumulation; and Step 3.3, inputting discretization results into the Bayesian identification and classification network which is trained for the corresponding spatial point to execute identification and classification; wherein in the Step 3, detailed steps of construction of the Bayesian identification and classification network which is trained for the each spatial point are: Step 31, accumulating the sound or the vibration signal and then preprocessing it, wherein then a typical event data training set is built based on a long time signal which is obtained from pretreatment; Step 32, extracting the customized short time feature from the relative value of the signal for the typical event data training set, wherein the long time feature is calculated and accumulated, after the extraction, then the typical event data training set is transformed into a long time feature training set; and Step 33, constructing the Bayesian identification classification network topology and learning the parameter based on the probability distribution information of the long time feature in the long time data training set and combined with prior probability information of a background environment at each spatial point, thus obtaining the trained Bayesian identification classification network; wherein the step 33 comprises steps of: Step 331, discretizing the long time feature training set to construct a discrete long time feature training set; Step 332, according to mutual information of all long time feature attribute nodes and class nodes obtained from the discrete long time feature training set at each spatial point, and conditional mutual information between feature attribute nodes under a class node condition, constructing a Bayesian network topology; and Step 333, constructing the Bayesian network topology at each spatial point wherein parameters are further studied based on the constructed Bayesian network topology to obtain the trained Bayesian identification classification network wherein the step 332 comprises steps of: Step 3321, calculating the mutual information between the long time feature attribute node and the class node: wherein the mutual information between all the long time feature attribute node B and the class node C is calculated as follows: $\begin{array}{cc}I\left(B_j;C\right)=\underset{b{S}_b}{\overset{}{.Math.}}\underset{c\left\{1,2\right\}}{\overset{}{.Math.}}p\left(b_j^{},c\right)\log \frac{p\left(b_j^{},c\right)}{p\left(b_j^{}\right)p\left(c\right)}\left(j=1,2,\mathrm{.Math.},6\right);& \left(15\right)\end{array}$ in Eq. (15), S.sub.b denotes a set of all possible discrete values taken for the long time feature attribute node B.sub.j(j=1, 2, . . . , 6)p(b.sub.j,c) represents a joint probability when B.sub.j=B.sub.j, C=c; p(b.sub.j) represents a probability when B.sub.j=b.sub.j; values of p(b.sub.j,c) and p(b.sub.j) are obtained from the discrete long time feature training set, and a calculation method is: p(b.sub.j,c) equals a number of records in the discrete long time feature training set when B.sub.j=b.sub.j,C=c, divided by a total number of records in the discrete long time feature training sets; and p(b.sub.j) equals a number of records in the discrete long time feature training set when B.sub.j=b.sub.j, divided by the total number of records in the discrete long time feature training sets; Step 3322, calculating the conditional mutual information between each two of all the feature attribute nodes under the class node condition; wherein according to the discrete long time feature training set, calculation of the condition mutual information between each two feature attribute nodes B.sub.i,B.sub.j(ij, and i=1,2, . . . , 6, j=1,2, . . . , 6) under the condition of class node C is: $\begin{array}{cc}I\left(B_i;B_j|C\right)=\underset{b_i^{}=0}{\overset{n}{.Math.}}\underset{b_j^{}=0}{\overset{h}{.Math.}}\underset{c=1}{\overset{2}{.Math.}}p\left(b_i^{},b_j^{},c\right)\log \frac{p\left(b_i^{},b_j^{}|c\right)}{p\left(b_i^{}|c\right)p\left(b_j^{}|c\right)},\left(i=1,2,\mathrm{.Math.},6,j=1,2,\mathrm{.Math.},6;ij\right);& \left(16\right)\end{array}$ in Eq. (16), b.sub.i represents a discrete value of B.sub.i, and n represents a maximum of possible discrete values of B.sub.i; b.sub.i represents a discrete value of B.sub.j, and h represents a maximum of the discrete value of B.sub.j; c represents a value of a class node C, and p(b.sub.i,b.sub.j,c) represents a joint probability when B.sub.i=b.sub.i,B.sub.j=b.sub.j,C=c; p(b.sub.i,b.sub.j|c), p(b.sub.i|c) and p(b.sub.j|c) are corresponding conditional probabilities, whose calculation methods are: p(b.sub.i,b.sub.j,c) equals the number of records in the discrete long time feature training set when B.sub.i=b.sub.i,B.sub.j=b.sub.j,C=c, divided by the total number of records in discrete long time feature training sets; p(b.sub.i,b.sub.j|c) equals the number of records in the discrete long time feature training set when B.sub.i=b.sub.i,B.sub.j=b.sub.j,C=c, divided by the number of C=c records in the discrete long time feature training set; calculation of p(b.sub.i|c) and p(b.sub.j|c) is similar to that of p(b.sub.i,b.sub.j|c); Step 3323, initialization: defining a set S1=, and a set S2={B.sub.1, B.sub.2, . . . , B.sub.6}, wherein one of the feature attribute node B.sub.j(j=1,2, . . . , 6) in S2 is chosen to make I(B.sub.j;C) achieve a maximum value thereof; then an arc is added from the class node C to B.sub.j, that is, the node C is taken as a parent node of B.sub.j; and Step 3324, choosing one long time feature attribute node B.sub.j in S2 to make the I(B.sub.j;C) get to the maximum value; wherein then the arc is added from C to B.sub.j; let k=min(|S1|,2), top k long time feature attribute nodes in S1 to make I(B.sub.i;B.sub.j|C) get to a maximum; then k arcs are added from node B.sub.i to node B.sub.j in the network, and the long time feature attribute node B.sub.j is added to S1, and B.sub.j is deleted from S2; the above operation is cycled from choice of B.sub.j in S2 until the set S1 include all of the long time feature attribute nodes; when determining the Bayesian network topology, the B.sub.6 is deleted from the set S2 first, and add to the set S1; remaining nodes in the set S2 are {B.sub.1,B.sub.2,B.sub.3,B.sub.4,B.sub.5}; in a new set S2, B.sub.2 and the class node C have biggest mutual information, and k=min(|S1|,2)=1, then k nodes are selected in S1 as the parent nodes of B.sub.2 which makes I(B.sub.i;B.sub.2|C) the largest; it also means B.sub.6 is the parent node of B.sub.2, then B.sub.2 is deleted from the set S2, and add to the set S1; meanwhile, S2={B.sub.1,B.sub.3,B.sub.4,B.sub.5}, and S1={B.sub.2,B.sub.6}; then in the new set S2, B.sub.5 and the class node C have biggest mutual information, and k=2, then we choose k nodes in S1 as parent nodes of B.sub.5, which makes I(B.sub.i;B.sub.2|C) the largest; then B.sub.5 is deleted from the set S2, and added to the set S1; the above operation is cycled until the set S2 is empty.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a schematic diagram of a system configuration and method flow of the present invention.

(2) FIG. 2 is a distributed optical fiber sensing signal processing flow based on Bayesian identification classification network in the present invention.

(3) FIG. 3 is the accumulation and preprocessing of the distributed optical fiber sensing signal in the present invention;

(4) (a) an accumulation and preprocessing flow of the distributed optical fiber sensing signal;

(5) (b) a spatio-temporal response graph of the anomalous events after accumulation and preprocessing.

(6) FIG. 4 is a probability distribution of the feature values of the three typical events in the present invention;

(7) (a) normal events;

(8) (b) real abnormal events; and

(9) (c) environmental disturbance events.

(10) FIG. 5 is a probability distribution of the long time features B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5 of the long time signal at a certain spatial point in the present invention;

(11) (a) the probability distribution of B.sub.1;

(12) (b) the probability distribution of B.sub.2;

(13) (c) the probability distribution of B.sub.3;

(14) (d) the probability distribution of B.sub.4;

(15) (e) the probability distribution of B.sub.5.

(16) FIG. 6 is a probability distribution of the long time feature B.sub.6 of the long time signal at a certain spatial point in the present invention.

(17) FIG. 7 is a Bayesian identification classification network constructed at a certain spatial point in the present invention.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

(18) In order to make the purpose, technical scheme and advantages of this invention more clearly understood, the following description will be combined with drawings and implementation examples to further elaborate the present invention. It is to be understood that the specific embodiments described herein are merely illustration of the present invention and not intended to limit the present invention.

(19) According to a first preferred embodiment of the present invention, the structure diagram and the sensing principle of the distributed optical fiber sensing signal processing method for safety monitoring of the underground pipe network are shown in FIG. 1. The hardware of the distributed optical fiber sound/vibration sensing system based on linear demodulation phase-sensitive OTDR (O-OTDR) technology mainly consists of three parts, an optical detection fiber, an optical signal demodulator, and a signal processing host. The optical detection fiber usually uses ordinary single-mode communication optical fiber. It is generally buried along the underground pipelines, transmission cables, or the urban roads, and it can also directly take a spare core of the communication cable already laid along the pipelines or the urban roads. The optical signal demodulator is the core of the system. Its internal components include optical devices and electrical devices. A continuous coherent optical signal produced by the ultra-narrow line-width lasers is modulated into optical pulse signal by acousto-optic or electro-optical modulator (AOM/EOM). The optical pulse signal is amplified by the erbium-doped fiber amplifier (EDFA), and the amplified light pulse signal is then injected into the detection fiber through the isolator, the port 1, port 2 of the circulator; The optical pulse signal generates Rayleigh scattering along the transmission in the fiber, and the backward Rayleigh scattering signal returns along the fiber, which is received by the port 2, port 3 of the circulator. The noise signal is filtered by the filter and coupled by the first coupler and then injected into a non-balanced Mach-Zehnder or Michelson interferometer, etc., as is determined by a specific demodulation method. Three demodulated phase signals with a difference of 120 degrees are output via a 3*3 second coupler. The phase changes introduced by the external disturbances demonstrate the force, the sound waves and vibration operating on the fiber. The demodulated optical signals are transformed into electrical signals by the photoelectric detector. The signal is acquired by the analog-to-digital converter synchronously triggered and controlled by the waveform generation card. Finally, the digital signals are transmitted to the signal processing host through the network interface at real time. The signal processing host is generally the ordinary computer host (PC) or FPGA/DSP embedded motherboard for fiber sensing signal analysis and processing. The event information inducing by the sound waves or vibration can be obtained by a specific signal processing algorithm, and intelligently identified and classified, and its location can be determined from the optical time domain reflection principle.

(20) Distributed optical fiber sensing signal processing method for the underground pipe network safety monitoring is detailed as follows. The optical detection fiber itself is equivalent to a series of sensor nodes evenly distributed along the fiber, and each sensor node corresponds to a data collection node in each spatial point. These distributed spatial sensor nodes or the acquisition nodes cooperates to pick up the sound or vibration signals on the whole line and extract distinguishable features of different sound and vibration signals. The features mainly denote the particularly defined short time features and long time feature extraction based on the relative amounts of signals in the present invention. Combined with a priori knowledge of the event occurring probability in different background environment of each spatial point, a distributed Bayesian identification classification network is constructed. Finally, the distributed detection, identification and classification of the abnormal events are carried out by the constructed Bayesian identification classification network, based on the short time and long time features extracted from the sensing signal at each spatial point. Specifically, the short time feature extraction method for the acquired sound or vibration signal in the whole monitoring range is the same as the method in constructing the Bayesian identification and classification network.

(21) According to a second preferred embodiment of the present invention, the distributed optical fiber sensor signal processing flow based on the distributed Bayesian identification classification network is shown in FIG. 2, which is divided into two parts:

(22) In the first part, the Bayesian identification and classification network at each spatial point is constructed based on the training database. Firstly, the sound/vibration sensing signals are accumulated and preprocessed. The typical event data training set is constructed from the preprocessed long time signal. The short time feature is extracted based on the relative quantity of the signal for the typical event data training set, and the long time feature is calculated by the extracted short time feature. After discretization, based on the probability distribution information of the long time feature of the discrete long time feature training set, combined with the prior probability information of the background environment at each point, the Bayesian identification network topology is constructed and the parameters are learned. Then the trained Bayesian identification classification network is obtained.

(23) In the second part, online detection, identification and classification are carried out based on the on-site event signal: Given a judgment threshold, the short time feature based on the relative quantity of the signal is used to detect the abnormal sound or vibration signal at each point. If any abnormal event is detected at a spatial point, the long time features can be cumulatively calculated from the short time features when the abnormal events occur at the point. The cumulative calculation result will be discretely processed. After the discretization, the result will be input into the constructed distributed Bayesian identification classification network for online identification and classification of the field event signal.

(24) According to a third preferred embodiment of the present invention, the sound/vibration sensing signal accumulation and preprocessing process at each spatial point is shown in FIG. 3 (a). The optical signal demodulator takes a trigger pulse as a signal acquisition period. The signal collected in this time unit is normally an exponentially decaying OTDR trajectory along the fiber. Assume the signal collected in the entire monitoring range at the k th trigger pulse period as:
X.sub.k={x.sub.ki(i=1,2, . . . ,N)}=[x.sub.k1,x.sub.k2, . . . ,x.sub.kN](1);

(25) In Eq. (1), i represents the spatial index, N is the total number of data collected in the entire monitoring range, which can be seen as N spatial sensor nodes. With the light pulse triggered periodically, the original OTDR trajectory collected in each light pulse trigger period along the time axis is accumulated. After continuous accumulation of T collected original OTDR trajectories, a temporal-spatial signal response matrix with N spatial dimensions and T temporal dimensions is constructed,
XX={x.sub.ki(k=1,2, . . . ,T;i=1,2, . . . ,N)}(2);

(26) T is the T th trigger pulse period. The graphic display of the temporal-spatial response matrix XX is shown in FIG. 3 (b), that is, the temporal-spatial response graph of the anomalous event obtained after the accumulation and preprocessing of the distributed optical fiber sensor signal. In the embodiment of the present invention, the short time feature extraction and the long time feature accumulative calculation based on the relative amount of the signal are performed for each record in the typical event data training set. The accumulation time of the short time and the long time signals are 1 s and 30 s respectively. The short time signal accumulated in the i th spatial point or sampling node is defined as:
w.sub.i={x.sub.ki,k=1,2, . . . ,M}(3);

(27) In the embodiment of the present invention, M is the number of data points accumulated for 1 s on the time axis. As the cycle trigger frequency is 508 Hz, the number of data points accumulated M is 508 within 1 s. m short time signals can be accumulated on the time axis and constitute a long time signal at each spatial point,
W.sub.i={x.sub.ki(k=1,2, . . . ,L)}
L=mM m1(4);

(28) wherein L is the data length of the long time signal, m is a positive integer. In the embodiment of the present invention, the cumulative length of the long time signal is 30 seconds, so m is 30.

(29) According to a forth preferred embodiment of the present invention, the typical event data training set is constructed based on the accumulated long time signal of the typical event at each spatial point. The operation process is detailed as follows. The long time signals with the length of L at each spatial point are accumulated and labeled as the event tag according to the types of events that actually occurs. During the underground pipe network safety monitoring process, the typical events at each spatial point include: the stationary environmental noise without abnormal events, artificial excavation, the traffic and factory interference which are easy to misjudge. In the present embodiment, the stationary environmental noise is the normal event and the event tag is set as 0; all the events such as artificial excavation, traffic and factory interference are referred as abnormal events, where the artificial excavation is a real anomaly event and the event tag is set as 1, while traffic and factory disturbances are environmental interferences and the event label is set as 2. Long time signals are added to the database according to the event tags to complete the three types of typical event data training set building.

(30) According to a fifth preferred embodiment of the present invention, each record of the typical event data training set is subjected to the short time feature extraction and the long time feature accumulative calculation based on the relative quantity of the signal in the present invention, and the typical event data training set is transformed into the corresponding long time feature training set, which is detailed as follows:

(31) 1. The Short Time Feature Extraction Method of Short Time Signal

(32) In order to avoid the influence on the detection and recognition caused by the fiber-sensing signal amplitude unevenness at different locations in a wide-range monitoring, the relative quantity features are extracted in the present invention. For the short time signal of is, four customized short time features based on the relative quantity of the signals are extracted: primary impulse intensity feature, A.sub.1, secondary impulse intensity feature, A.sub.2, primary and secondary impulse amplitude ratio feature, A.sub.3 and wavelet packet energy spectrum feature A.sub.4. Among them, the primary impulse intensity feature A.sub.1 extracts the maximum impulse features of the signal; the secondary impulse intensity feature A.sub.2 extracts the impulse features of the residual signal after removing the maximum impulse part of the signal; primary and secondary impulse amplitude ratio feature A.sub.3 extracts impulse amplitude ratio of the maximum impulse part and the residual signal after removing the maximum impulse part of the signal; the wavelet packet energy spectrum feature A.sub.4 extracts the energy spectrum feature of the important frequency band of the impulse signal.

(33) (1) The Extraction Method of the Primary Impulse Intensity Feature A.sub.1

(34) Generally, normal events are usually stationary random noise signals, and typical abnormal event signals are usually accompanied by great vibration shocks. The first short time feature extracted in the present invention is the primary impulse intensity feature A.sub.1. The short time signal w.sub.i is taken as the processing object for the short time feature extraction. First, histogram statistics is carried out for the short time signal amplitude values. The amplitude range of the short time signal w is divided into 10 consecutive intervals. max(w.sub.i) represents the maximum value of w.sub.i min(w.sub.i) represents the minimum value of w.sub.i, and each interval is recorded as s.sub.l(l=1, 2, . . . , 10), that is:
s.sub.l=[(l1)averagelength+min(w.sub.i),laveragelength+min(w.sub.i)],(l=1,2, . . . ,10)(5);

(35) In Eq. (5), averagelength is the width of each interval, that is:

(36) $\begin{array}{cc}\mathrm{averlength}=\frac{\max \left(w_i\right)-\min \left(w_i\right)}{10};& \left(6\right)\end{array}$

(37) Calculate the corresponding sampling point number distributed in each interval for the short time period signal w.sub.i respectively and generate a frequency distribution histogram, namely compute the percentage of each interval, percent:

(38) $\begin{array}{cc}\{\begin{array}{c}\mathrm{percent}=\left\{{\mathrm{percent}}_l,l=1,2,\mathrm{.Math.},10\right\}\\ {\mathrm{percent}}_l=\frac{{}_l}{M}\end{array};& \left(7\right)\end{array}$

(39) In Eq. (7), n.sub.i is the sampling point number of the short time signal, whose amplitudes range in s.sub.l, M is the data length of the short time signal. Assuming that the maximum of percent.sub.l(l[1,10]) is got at l=q, where q is the amplitude interval corresponding to the highest point of the histogram, the primary impulse intensity feature A.sub.1 is computed as:

(40) $\begin{array}{cc}{A}_1=\{\begin{array}{c}{p}_q+P_{q+1}\left(q\left[2,9\right],p_{q-1}{p}_{q+1}\right)\\ p_q+P_{q-1}\left(q\left[2,9\right],p_{q-1}>p_{q+1}\right)\\ p_q+P_{q+1}\left(q=1\right)\\ p_q+P_{q-1}\left(q=10\right)\end{array};& \left(8\right)\end{array}$

(41) For the stationary ambient noise i.e. the short time signal of the normal event, the horizontal amplitude distribution interval is small, and the sampling point number, i.e. the height of the histogram difference is not obvious in each interval. At this time, the calculated value of A.sub.1 is very small. While the horizontal amplitude distribution interval will be abruptly increased for the short time signal of the abnormal event, simultaneously, the absolute amplitude value of most sampling points is small, and most of the sampling points are centralized in the two smaller amplitude distribution intervals in the middle. At this time, the calculated value of A.sub.1 is larger. So it can be used as the first short time distinguishable feature in the extraction operation, which is denoted as the primary impulse intensity feature. As shown in FIG. 4, there is a significant difference in the probability distributions of the value of the primary impulse intensity feature A.sub.1 for the three typical events. The value of A.sub.1 for the normal event signal in FIG. 4 (a) is small and the distribution center is near 0.4, while the value of A.sub.1 for the real anomaly events in FIG. 4 (b) and that for the environmental interference events in FIG. 4 (c) are both large and the distribution centers are near 0.9.

(42) (2) The Extraction Method of Secondary Impulse Intensity Feature A.sub.2

(43) As shown in FIG. 4, two types of abnormal events, such as real anomaly events and environmental disturbances, usually have a large primary impulse intensity, and thus we continue to extract a second short time feature, the secondary impulse intensity feature A.sub.2. Based on the data analysis, the short time signal of real abnormal events, i.e. artificial excavation, generally only have a large impulse intensity. After the impulse is removed, the remaining short time signal is close to the stationary environmental noise; while the short time signal of traffic interference and other environmental interference events is normally accompanied by several impulses with different intensity due to uneven road or stone shocks. After the primary impulse is removed, the remaining signal still contains more than one impulse. The extraction method of the secondary impulse intensity feature A.sub.2 is detailed as follows: assuming that the absolute maximum value of the short time signal w.sub.i={x.sub.ki, k=1, 2, . . . , M} is obtained at k=t, where t is the time when the absolute value of the short time signal is max. Remove M/4 sampling points near k=t from w.sub.i to obtain the remaining 3M/4 sampling points, which constitute the time sequence w.sub.i, where M is the short period of signal length:

(44) $\begin{array}{cc}{w}_i^{}=\{\begin{array}{c}{x}_{\mathrm{ki}}\left(\begin{array}{c}k=1,2,\mathrm{.Math.},t-\frac{M}{8},\\ t+\frac{M}{8},t+\frac{M}{8}+1,\mathrm{.Math.},M\end{array}\right)\left(\frac{M}{8}<t<\frac{7M}{8}\right)\\ x_{\mathrm{ki}}\left(k=\frac{M}{4}+1,\frac{M}{4}+2,\mathrm{.Math.},M\right)\left(t\frac{M}{8}\right)\\ x_k\left(k=1,2,\mathrm{.Math.},\frac{3M}{4}\right)\left(t\frac{7M}{8}\right)\end{array};& \left(9\right)\end{array}$

(45) Finally, the secondary impulse intensity feature A.sub.2 is extracted for w.sub.i, by the same method as that of extracting the primary impulse feature A.sub.1, as shown in Eq. (5), (6), (7) and (8). In general, the value of the secondary impulse intensity feature A.sub.2 is small for the short time signal of the artificial excavation event, while that for the short time signal of the traffic and factory disturbance event is larger.

(46) (3) The Extraction Method of the Primary and Secondary Impulse Amplitude Ratio Feature A.sub.3

(47) Due to the fact that the short time signal of real abnormal event after the primary impulse is removed is similar to the stationary noise, its amplitude is very small compared with the maximum impulse amplitude. While the short time signal of the traffic and factory interference events after the primary impulse is removed have several impulses with different intensity, whose amplitudes are always large compared with the primary impulse amplitude. Therefore, the calculation method of the primary and secondary impulse amplitude ratio feature A.sub.3 as a third short time feature, is:

(48) $\begin{array}{cc}{A}_3=\frac{\max \left(.Math.w_i^{}.Math.\right)}{\max \left(.Math.w_i^{}.Math.\right)};& \left(10\right)\end{array}$

(49) In Eq. (10), max(|w.sub.i|) represents the maximum absolute value of the amplitude of the time series signal w.sub.i; max(|w.sub.i|) represents the maximum absolute value of the amplitude of the time series signal w.sub.i.

(50) (4) The Extraction Method of Wavelet Packet Energy Spectrum Feature A.sub.4

(51) In general, the low- and middle-frequency information is relatively rich in the instantaneous impulses for real abnormal events; while low-frequency information is rich in that of the traffic and factories or other environmental interference events. Based on this point, a three-layer db6 wavelet packet decomposition is carried out on the short time signal w.sub.i. We sort the decomposed sub-bands from low to high according to the frequency component, which is marked as s.sub.0s.sub.7. Finally, the wavelet packet coefficient energy of each frequency band is calculated to form the wavelet packet coefficient energy feature vector E.sub.wp=[E.sub.0, E.sub.1, . . . , E.sub.7]. The wavelet packet energy spectrum feature A.sub.4 is obtained as:

(52) $\begin{array}{cc}{A}_4=\frac{E_2+E_3}{E_0+E_1};& \left(11\right)\end{array}$

(53) In Eq. (11), E.sub.2, E.sub.3 represents the energy of the medium and low frequency sub-bands, E.sub.0, E.sub.1 represents the energy of the low frequency sub-bands. The low- and medium-frequency information of the short time signal is rich for real-time event, and the value of A.sub.4 is large, while the low-frequency information of the short time signal is rich for the traffic and factory interferences or other environmental events, and its value of A.sub.4 is generally small.

(54) 2. The Cumulative Calculation Method of Long Time Feature of Long Time Signal

(55) Short time signal is a local observation of the event and the event classification based on the short time features always cannot obtain a high identification accuracy. So after extraction of the short time features, we continue to accumulate a long time signal and extract its features to observe the event for a longer time, to get an accurate identification based on the rich event information. In the present invention, there are six customized long time features for the accumulated long time signals, including the long time primary impulse intensity feature B.sub.1, the long time secondary impulse intensity feature B.sub.2, the long time amplitude ratio feature B.sub.3, the long time wavelet packet energy spectrum feature B.sub.4, long time environmental background feature B.sub.5 and long time strict impulse number feature B.sub.6, and all of them are obtained from the cumulative calculation result based on the extracted customized short time features, where the long time primary impulse intensity feature B.sub.1 represents the average of the primary impulse intensities of all the short time signals of the abnormal events during the accumulation time period; the long time secondary impulse intensity feature B.sub.2 represents the average of the secondary impulse intensities of all the short time signals of the abnormal events during the accumulation time period; the long time amplitude ratio feature B.sub.3 represents the average of the amplitude ratio feature of all the short time signals of the abnormal events during the accumulation time period; the long time wavelet packet energy spectrum feature B.sub.4 represents the average of the wavelet packet energy spectrum feature of all the short time signals of the abnormal events during the accumulation time period; the long time environmental background feature B.sub.5 represents the background fluctuation of the long time signal; the long time strict impulse number feature B.sub.6 represents the detected number of strict impulses within long time signal. The long time features of the long time signal don't need to be re-extracted. They can be obtained by accumulating or averagely calculating of m short time features corresponding to the long time signal, and the detailed steps are:

(56) (a) initialization

(57) Set the variable t=1; denote the set T.sub.1=, T.sub.2=, T.sub.3=, T.sub.4=, T.sub.5=, T.sub.6=; and consider the long time signal as m short time signal w.sub.i,t (t=1, 2, . . . , m), where w.sub.i,t represents the t th short time signal.
B.sub.1=average(T1),B.sub.2=average(T2),B.sub.3=average(T3),B.sub.4=average(T4)(12);

(58) (b) Repeat the Following Operations Until t=m

(59) The values of the short time features A.sub.1, A.sub.2, A.sub.3, A.sub.4 of the short time signal are denoted as a1, a2, a3, a4 respectively. Observe the four types of short time feature of m short time signals within the window according to the long time signal and put the feature values in different sets T1, T2, T3, T4, T5, T6 respectively, which is aimed to prepare for the extraction of six long time features. For example, if a1>thr, the short time signal is judged to be an abnormal event, a1, a2, a3, a4 will be added to the sets T1, T2, T3, T4 respectively; if a1<thr, the short time signal is judged to be a normal event, a1 will be added to the set T5; if a1>thr and a2<thr1, a3<thr2, which meets the strict impulse conditions, a1 will be put into the set T6. thr, thr1, thr2 are the empirical threshold to judge the impulse intensity strong or weak respectively. When the m short time feature values of the short time signal are all placed in the set T1, T2, T3, T4, T5, T6, the six long time feature values of the long time signal are calculated. Among them, averaging the short time feature values in the set T1, T2, T3, T4 can obtain the value of B.sub.1, B.sub.2, B.sub.3, B.sub.4 in turn.

(60) The set T5 contains a1 of all normal events during long time period. Generally, the artificial excavation and other real abnormal events appear intermittently, and the short time signal nearby is relatively stationary, thus the values of a1 are always small; while traffic vehicles, factory interference and other environmental interference events are more complex, and the short time signal nearby also contains more small fluctuations, thus the values of a1 are always larger. In the end, all the a1 values in the set T5 are averaged, and the background fluctuation of the long time signal can be obtained. Therefore, the long-term environmental background feature B.sub.5 is calculated as:
B.sub.5=average(T5)(13).

(61) Generally, the impulse of the environmental interference event signal and that of the real abnormal event signal is significantly different. In order to minimize the false alarms induced by the environmental interference, strictly selection and exclusion for the impulse signal generated by environmental interference events is needed. Then we take the strict impulse number feature B.sub.6 as one of the criterion for the abnormal event determination. The impulse simply judged by the single empirical threshold thr actually contain two categories, pseudo-impulse and strict impulse; Then another two empirical threshold thr1 and thr2 can be added to select the strict impulses, which can be used to exclude the pseudo-impulses and detect the real abnormal events. Specifically, when the value of the short time primary impulse intensity feature A.sub.1 satisfies a1>thr, if and only if the value of the short time secondary impulse intensity feature A.sub.2 satisfies a2<thr1, and when the short time amplitude ratio feature A.sub.3 satisfies a3<thr2, it is then judged that the short time signal w.sub.i,t is a strict impulse signal; wherein the selection criteria of the empirical threshold thr1 and thr2 are determined by the training signals in the training event database. In the present embodiment, when the real abnormal event signal caused by artificial excavation is processed by the thresholds thr, thr1 and thr2, it is judged to be a strict impulse signal at a probability of 70% or more. After the threshold processing with thr, thr1, thr2, only 10% of environmental disturbances such as traffic and factory interference can be judged as a strict impulse signal. Finally, the number of long time impulse feature B.sub.6 is calculated as:
B.sub.6=length(T6)(14);

(62) Eq. (14) represents the number of data in the set T6, that is, the impulse number meeting short time strict impulse conditions in the long time period.

(63) 3. Construction of the Long Time Feature Training Set Corresponding to the Typical Event Data Training Set

(64) According to the feature extraction method above, the short time feature extraction and the long time feature accumulation calculation are carried out for each data record in the typical event data training set. The long time features attached with the event type tag of the record will eventually constitute a long time feature vector corresponding to data record, which is denoted as the Train_Feature=[class, b1, b2, b3, b4, b5, b6]; where in class represents the event type tag for the record, b1,b2,b3,b4,b5,b6 are respectively the values of the long time feature B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5, B.sub.6. Finally, the whole typical event data training set is transformed into the corresponding long time feature training set TrainingFeature.

(65) According to a sixth preferred embodiment of the present invention, the long time feature training set of the typical event data training set is discretized, and the specific steps are as follows:

(66) 1. Discretization of Long Time Feature Attribute Value

(67) Since the first five long time feature attribute values of the feature set {B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5, B.sub.6} are continuous and the value of each attribute node in the Bayesian recognition classification network should be discrete, it is necessary to discretize the long time feature attribute value first. The method comprises of the following steps: taking a certain spatial point as an example, the probability distribution of the long time feature is calculated from the continuous long time feature training set of the point, and then the interval division and discretization are carried out according to the probability distribution.

(68) (1) Discretization of the Attribute Value of Long Time Feature B.sub.1

(69) The probability distribution of the long time feature B.sub.1 of a long time signal at a certain spatial point is calculated based on the continuous long time feature training set. As shown in FIG. 5(a), the probability density of the environmental interference event is large in the interval [0.80, 0.91]; and the probability density of the real anomaly event is larger in the interval [0.92,1.0]; while the probability density of the two kinds of abnormal events is not different in the interval [0.91, 0.92]. Therefore, the discrete values of the long time feature B.sub.1 in intervals (0,0.8), [0.80, 0.91], (0.91, 0.92], and (0.92,1.0] are set to 0, 1, 2, 3 respectively.

(70) (2) Discretization of the Attribute Value of Long Time Feature B.sub.2

(71) The probability distribution of the long time feature B.sub.2 of a long time signal at a certain spatial point is calculated based on the continuous long time feature training set. As shown in FIG. 5 (b), the probability density of the real anomaly event is larger than that of the environmental interference event in the interval [0.63,0.68); while the probability density of the two kinds of abnormal events is not different in the interval [0,0.63); Therefore, the discrete values of the long time feature B.sub.2 in intervals [0,0.63), [0.63,0.68), and [0.68,1.0] are set to 0, 1, 2.

(72) (3) Discretization of the Attribute Value of Long Time Feature B.sub.3

(73) The probability distribution of the long time feature B.sub.3 of a long time signal at a certain spatial point is calculated based on the continuous long time feature training set. As shown in FIG. 5 (c), the probability density of the real anomaly event is large in the interval [0,0.32); the probability density of the environmental interference event is large in the interval [0.39,1.0]; and the probability density of the two kinds of abnormal events is not different in the interval [0.32,0.39). Therefore, the discrete values of the long time feature B.sub.3 in intervals [0,0.32), [0.32,0.39), and [0.39,1.0] are set to 0, 1, 2.

(74) (4) Discretization of the Attribute Value of Long Time Feature B.sub.4

(75) The probability distribution of the long time feature B.sub.4 of a long time signal at a certain spatial point is calculated based on the continuous long time feature training set. As shown in FIG. 5 (d), the probability density of the environmental interference event is large in the interval [0,0.6), the probability density of the real anomaly event is large and that of the environmental interference event is small in the interval [0.8,3.0], the probability density of the two kinds of abnormal events is not very different in the interval [0.6,0.8), Therefore, the discrete values of the long time feature of B.sub.4 in intervals [0,0.6), [0.6,0.8), [0.8,3.0], and (3.0,.sub.) are set to 0, 1, 2, 3.

(76) (5) Discretization of the Attribute Value of Long Time Feature B.sub.5

(77) The probability distribution of the long time feature B.sub.5 of a long time signal at a certain spatial point is calculated based on the continuous long time feature training set. As shown in FIG. 5 (e), the probability density of the real anomaly event is large in the interval [0,0.59), while the probability density of the environmental interference event is large in the interval [0.67,1.0], and the probability density of the two kinds of abnormal events is not different in the interval [0.59,0.67). Thus in this inventive embodiment, the discrete values of the long time feature B.sub.5 in the interval [0,0.59), [0.59,0.67), [0.67,1.0] are set to 0, 1, 2, respectively.

(78) Long time feature B.sub.6 is a count number of the strict impulses of m short time signal. As the value of itself is discrete, there is no need for further discretization. The probability distribution of the long time feature B.sub.6 as is shown in FIG. 6, can be used directly to construct the Bayesian identification classification network.

(79) 2. The Construction of Discrete Long Time Feature Training Set

(80) After the feature attribute value discretization as above, the continuous long time feature training set TrainingFeature corresponding to each data record, can be transformed into the discrete long time feature training set TrainingFeature, as stated in the fifth embodiment. The discrete long time feature training set is denoted as, Train_Feature=[class, b1, b2, b3, b4, b5, b6], where b1, b2, b3, b4, b5, b6 are the corresponding discrete feature values of the long time feature B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5, B.sub.6 after the discretization.

(81) According to a seventh preferred embodiment of the present invention, the distributed Bayesian identification classification network at each spatial point (corresponding to each spatial node) is constructed, which includes Bayesian network topology and network parameter learning at that spatial point.

(82) Taking the Bayesian identification and classification network construction process at a certain point as an example, according to the mutual information of all long time feature attribute nodes and class nodes obtained from the discrete long time feature training set at this point, and the conditional mutual information between the feature attribute nodes under the class node condition, the Bayesian network topology is constructed; and the parameters are further studied based on the constructed Bayesian network topology. On the basis of the network topology, the conditional probability of each feature attribute node is calculated, and the detailed steps are as follows:

(83) 1. The Construction of the Bayesian Network Topology Based on the Mutual Information and the Conditional Mutual Information

(84) (1) calculation of the mutual information between the long time feature attribute node and the class node

(85) The mutual information between all the long time feature attribute node B.sub.j and the class node C is calculated as follows:

(86) $\begin{array}{cc}I\left(B_j;C\right)=\underset{b{S}_b}{.Math.}\underset{c\left\{1,2\right\}}{.Math.}p\left(b_j^{},c\right)\log \frac{p\left(b_j^{},c\right)}{p\left(b_j^{}\right)p\left(c\right)}\left(j=1,2,\mathrm{.Math.},6\right);& \left(15\right)\end{array}$

(87) In Eq. (15), S.sub.b denotes a set of all the possible discrete values taken for the long time feature attribute node B.sub.j(j=1, 2, . . . , 6), p(b.sub.j, c) represents the joint probability when B.sub.j=b.sub.j, C=c; p(b.sub.b) represents the probability when B.sub.j=b.sub.j; the values of p(b.sub.j, c) and p(b.sub.j) can be obtained from the discrete long time feature training set, and the calculation method is: p(b.sub.j, c) equals the number of records in the discrete long time feature training set when B.sub.j=b.sub.j, C=c, divided by the total number of records in discrete long time feature training sets; and p(b.sub.j) equals the number of records in the discrete long time feature training set when B.sub.j=b.sub.j, divided by the total number of records in discrete long time feature training sets. In the present embodiment, the mutual information I(B.sub.j; C) between all the discrete long time characteristic attribute nodes B.sub.j and the class nodes C is calculated according to the method described as above.

(88) (2) The Calculation of Conditional Mutual Information Between Each Two of all the Feature Attribute Nodes Under the Class Node Condition

(89) According to the discrete long time feature training set, the calculation of the condition mutual information between each two feature attribute nodes B.sub.i, B.sub.j(ij, and i=1, 2, . . . , 6, j=1, 2, . . . , 6) under the condition of class node C is:

(90) $\begin{array}{cc}I\left(B_i;B_j|C\right)=\stackrel{n}{\underset{b_i^{}=0}{.Math.}}\underset{b_j^{}=0}{\overset{h}{.Math.}}\underset{c=1}{\overset{2}{.Math.}}p\left(b_1^{},b_j^{},c\right)\log \frac{p\left(b_i^{},b_j^{}|c\right)}{p\left(b_i^{}|c\right)p\left(b_j^{}|c\right)},\left(i=1,2,\mathrm{.Math.},6,j=1,2,\mathrm{.Math.},6;i1\right);& \left(16\right)\end{array}$

(91) In Eq. (16), b.sub.i represents the discrete value of B.sub.i, and n represents the maximum of the possible discrete values of B.sub.i; b.sub.j represents the discrete value of B.sub.j, and h represents the maximum of the discrete value of B.sub.j; C represents the value of class node C, and p(b.sub.i, b.sub.j, c) represents the joint probability when B.sub.i=b.sub.i, B.sub.j=b.sub.j, C=c; p(b.sub.i, b.sub.j|c), p(b.sub.i|c) and p(b.sub.j|c) are the corresponding conditional probabilities, whose calculation methods are: p(b.sub.i, b.sub.j, c) equals the number of records in the discrete long time feature training set when B.sub.i=b.sub.i, B.sub.j=b.sub.j, C=c, divided by the total number of records in discrete long time feature training sets; p(b.sub.i, b.sub.j|c) equals the number of records in the discrete long-term feature training set when B.sub.i=b.sub.i, B.sub.j=b.sub.j, C=c, divided by the number of C=c records in the discrete long time feature training set; calculation of p(b.sub.i|c) and p(b.sub.j|c) is similar to that of p(b.sub.i, b.sub.j|c). In the present embodiment, the conditional mutual information I(B.sub.i; B.sub.j|C) between each two feature attribute nodes under the class nodes condition is calculated according to the method described as above.

(92) (3) The Construction of Bayesian Network Topology

(93) According to the mutual information between the long time feature attribute node and the class node, and the conditional mutual information between each two feature attribute nodes under the class node condition, the Bayesian network topology is determined as follows:

(94) (a) Initialization

(95) Let the set S1=, and the set S2={B.sub.1, B.sub.2, . . . , B.sub.6}, one of the feature attribute node B.sub.j(j=1, 2, . . . , 6) in S2 is chosen to make I(B.sub.j; C) achieve its maximum value of. Then an arc is added from the class node C to B.sub.j, that is, the node C is taken as a parent node of B.sub.j;

(96) (b) The Bayesian Network Topology is Determined by the Repetition of the Following Operation

(97) One long time feature attribute node B.sub.j in S2 is chosen to make the I(B.sub.j; C) get to its maximum value. Then an arc is added from C to B.sub.j; let k=min(|S1|, 2), the top k long time feature attribute nodes in S1 to make I(B.sub.i; B.sub.j|C) get to its maximum. Then k arcs are added from node B.sub.i to node B.sub.j in the network, and the long time feature attribute node B.sub.j is added to S1, and delete B.sub.j from S2. The above operation is cycled from the choice of B.sub.j in S2 until the set S1 include all of the long time feature attribute nodes.

(98) For example, when we determine the Bayesian network topology in this embodiment, since the mutual information between B.sub.6 and the class node C is the largest, we delete the B.sub.6 from set S2 first, and add it to set S1; the remaining nodes in set S2 are {B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5}; in the new set S2, B.sub.2 and the class node C have the biggest mutual information, and k=min(|S1|,2)=1, then we select k nodes in S1 as parent nodes of B.sub.2, which can make I(B.sub.1; B.sub.2|C) the largest. It also means we make B.sub.6 as the parent node of B.sub.2, then we delete B.sub.2 from the set S2, and add it to the set S1; at this time, S2={B.sub.1, B.sub.3, B.sub.4, B.sub.5}, and S1={B.sub.2,B.sub.6}. Then in the new set S2, B.sub.5 and the class node C have the biggest mutual information, and k=2, then we choose k nodes in S1 as parent nodes of B.sub.5, which can make I(B.sub.1; B.sub.2|C) the largest. Then B.sub.5 is deleted from the set S2, and added to the set S1; the above operation is cycled until the set S2 is empty. In the present embodiment, the Bayesian network topology is constructed as shown in FIG. 7.

(99) 2. The Learning of Bayesian Network Parameter

(100) The network parameters is learned based on the above constructed Bayesian network topology. Specifically, the conditional probability p(b.sub.i|Pa.sub.i,c) (i=1, 2, . . . , 6) is calculated for each feature attribute node of the Bayesian network topology, where b.sub.i represents the value of the long time feature attribute node B.sub.i, and Pa.sub.i obtained by the training process above is the parent node of the attribute node B.sub.i. The calculation method is detailed as follows: p(b.sub.i|Pa.sub.i, c) equals the number of records in the discrete long time feature training set when B.sub.i=b.sub.i, Pa.sub.i=pa.sub.i, C=c, divided by the number of records in the discrete long time feature training set when Pa.sub.i=pa.sub.i, C=c, where Pa.sub.i is the value of B.sub.i's parent node Pa.sub.i.

(101) Here we just take the construction and training of the Bayesian identification classification network at one spatial point as a specific example. The Bayesian identification and classification network at all the other spatial points along the fiber can be distributedly constructed and trained as above.

(102) According to an eighth preferred embodiment of the present invention, the detection and recognition process of an on-line anomalous event is carried out based on the distributed Bayesian identification and classification network constructed and trained for all the spatial point. Firstly, the abnormality of signal is judged from the short time feature. When an anomaly event is detected at a spatial point, its long time feature accumulation calculation and discretization process is carried out. The discretization results is then input into the constructed distributed Bayesian identification classification network, then the abnormal event can be identified and classified at this point, that is to say it can be determined whether this abnormality is a true abnormal event caused by the artificial excavation, or the environmental interference events caused by traffic, factory and other typical interference. The detailed detection process is as follows: the short time signal is accumulated for one second(s) at each spatial point along the time axis, the short time primary impulse feature A.sub.1 is then extracted from a unit of short time signal for abnormal events detection. If the value of the short time feature A.sub.1 is smaller than the empirical threshold thr, then the short time signal is judged as a kind of normal event. Then we can move to the next spatial point to perform similar short time feature extraction and abnormal event detection; while when A.sub.1 exceeds the empirical threshold thr, it is judged that an abnormal event occurs at that spatial point.

(103) According to a ninth preferred embodiment of the present invention, based on the Bayesian identification and classification network trained at each spatial point, the abnormal events detected in the field application can be identified and classified on line as follows: when an anomaly event is detected at a spatial point, the signal is accumulated for a period of 30 s from that moment; when the accumulation time is out, the corresponding long time feature B.sub.1, B.sub.2, B.sub.3, B.sub.4, B.sub.5, B.sub.6 is accumulated and calculated which is similar to the training phase; then these features are discretized, and input into the distributed Bayesian identification classification network prepared at that point, and the abnormal event is identified and classified, and the abnormal event type such as the real abnormal event, the environment interference event and etc can be output from the network. The detailed steps are as follows:

(104) 0$\begin{array}{cc}{\mathrm{PS}}_1=\frac{p\left(c=1\right)p\left(b_1^{},\mathrm{.Math.},b_6^{}|c=1\right)}{p\left(b_1^{},\mathrm{.Math.},b_6^{}\right)}p\left(c=1\right)\underset{i=1}{\overset{6}{.Math.}}p\left(b_i^{}|{\mathrm{pa}}_i,c=1\right);& \left(17\right)\\ {\mathrm{PS}}_2=\frac{p\left(c=2\right)p\left(b_1^{},\mathrm{.Math.},b_6^{}|c=2\right)}{p\left(b_1^{},\mathrm{.Math.},b_6^{}\right)}p\left(c=2\right)\underset{i=1}{\overset{6}{.Math.}}p\left(b_i^{}|{\mathrm{pa}}_i,c=2\right);& \left(18\right)\end{array}$

(105) In Eq. (17), p(c=1) denotes the probability when the event category is a real abnormal event, which is a priori probability value; and p(b.sub.1, . . . , b.sub.6|c=1) is a conditional probability, indicating the probability that the long time feature takes the current value under the condition when the event category is a real abnormal event; and p(b.sub.1, . . . , b.sub.6) indicates the probability when the long time feature takes the current value; and p(b.sub.i|pa.sub.i, c=1) represents the probability when the discrete value of the long time feature attribute node B.sub.i is b.sub.i, under the condition that the discrete value of the parent node of long time feature attribute node B.sub.i is pa.sub.i and the category is a real abnormal event. The physical meaning of each symbol in Eq. (18) is similar to that of Eq. (17). Eq. (17) is computed for the posterior probability of the real anomaly event, and Eq. (18) is for the posterior probability of the environmental disturbance event. Since the value of p(b.sub.1, . . . , b.sub.6) of PS.sub.1 and PS.sub.2 in actual calculation process are the same, computation of p(b.sub.1, . . . , b.sub.6) is equivalent to that of

(106) $\underset{i=1}{\overset{6}{.Math.}}p\left(b_i^{}|{\mathrm{pa}}_i,c=1\right),$
where p(b.sub.i|pa.sub.i, c=1) is the parameter obtained in the training process. Then the computation of PS.sub.1 and PS.sub.2 is equivalent to PS.sub.1 and PS.sub.2 respectively:

(107) $\begin{array}{cc}{\mathrm{PS}}_1^{}=p\left(c=1\right)\underset{i=1}{\overset{6}{.Math.}}p\left(b_i^{}|{\mathrm{pa}}_i,c=1\right);& \left(19\right)\\ {\mathrm{PS}}_2^{}=p\left(c=2\right)\underset{i=1}{\overset{6}{.Math.}}p\left(b_i^{}|{\mathrm{pa}}_i,c=2\right);& \left(20\right)\end{array}$

(108) The event category is thus determined according to the calculated posterior probability values as shown in Eq. (19) and (20). For example, if a posterior probability at a spatial point is calculated as PS.sub.1>PS.sub.2, the abnormal event is judged as a true abnormal event with probability PS.sub.1, and its classification type is class=1; otherwise, it is judged as an environmental interference event with probability PS.sub.2 and its classification type is class=2.

(109) The present invention discloses a kind of distributed optical fiber sensing signal processing method for safety monitoring of underground pipe network, which is applicable not only to the distributed optical fiber acoustic/vibration sensing system based on the phase sensitive optical time domain reflection technology, but also to MZ, Micelson, Sagnac interferometer and other types of distributed optical fiber sound/vibration sensing systems. And it is also applicable to the quasi-distributed fiber optic acoustic/vibration sensing system based on grating and Febry-perot array; the short time and the long time feature extraction and discretization in this method can be adapted to the practical applications. The distributed identification network construction method according to the sensing signals of the real abnormal events and the environment interference at each spatial point can be directly applied to the safety monitoring of underground pipe network such as oil and gas pipelines, water supply pipelines, heat pipes and so on, and also applies to communication Fiber optic cables, power cables, perimeter security, structural health monitoring and other distributed optical fiber sensing applications. What is described above is provided only as a preferred embodiment of the present invention and not intended to limit the present invention, thus any modifications, equivalent substitutions and improvements within the spirit and principles of the present invention are intended to be protected by the present invention within the range.