Distributed acoustic sensing road cavity detection method based on active excitation source

Abstract

The invention discloses a distributed acoustic sensing road cavity detection method based on an active excitation source. The method is based on distributed acoustic sensing (DAS), and comprises the following steps of: firstly, acquiring a DAS signal sampling data sequence to be processed; preprocessing the data of each DAS channel; calculating a short-time zero-crossing rate to remove channels with abnormal data records; then, signal space-time positioning is carried out through sliding window detection, signal envelopes in a signal area are extracted, and space-time distribution information of signals is obtained; and finally, extracting spatial correlation features and energy distribution features of the signals, and realizing automatic detection of the road cavity based on feature judgment. According to the method, all-weather sensing and high-resolution precision characteristics of the DAS array are fully utilized, accurate estimation of road cavity detection is achieved with small calculation amount, large-range sensing can be achieved, and the method is suitable for real-time engineering application occasions.

Description

Technical Field

[0001] The present invention belongs to the technical field of distributed acoustic sensing, and in particular relates to a distributed acoustic sensing road cavity detection method based on an active excitation source. Background Art

[0002] With the accelerated pace of urbanization, the monitoring and management of urban road voids face significant challenges. The intertwined impacts of various municipal projects and the complexity of underground pipeline networks in current urban development are gradually revealing structural defects in road subsurface structures. These hidden, sudden, and highly destructive potential hazards have become a significant source of risk to urban public safety. Existing monitoring methods are primarily implemented in two scenarios: the first is during the emergency response phase following a collapse incident, when the risk has already translated into actual losses, resulting in not only economic losses but also public concern; the second is during periodic road maintenance inspections. While this approach can provide early warning, it faces practical challenges such as high technical costs and intensive resource investment. This coexistence of passive response and costly preventative measures highlights the urgency and systemic challenges of urban underground space safety management. Due to the strong electromagnetic interference resistance, excellent electrical insulation, and light transmission characteristics of optical fiber, distributed sensing systems using optical fiber as both a sensitive element and a signal transmission medium have recently attracted increasing attention in applications such as urban road monitoring and intrusion detection. Therefore, applications such as urban road subsurface safety monitoring based on DAS are of great significance.

[0003] The existing methods for detecting cavities in urban roads mainly include: (1) ground penetrating radar method, which uses radar echoes for imaging to reflect the underground space structure and can perceive the underground space structure with high precision. However, this method cannot monitor the underground space structure in real time and has a limited detection range. At the same time, the deployment of ground penetrating radar is time-consuming and cannot meet the requirements of monitoring the underground space structure of a large number of existing urban roads; (2) passive noise imaging method based on DAS, which uses DAS to receive background noise signals such as vehicle vibration and subway running signals, and after long-term accumulation and superposition processing, it can reflect the structure of the underground space to a certain extent. However, this method requires a long time to record the background noise signal and is only of reference significance for underground structure monitoring under real-time conditions; (3) DAS transient surface wave imaging method based on active excitation source, which is widely used because the active excitation source signal characteristics are more obvious and can be controlled manually. However, the performance of surface wave imaging method deteriorates sharply under low signal-to-noise ratio conditions and cannot truly reflect the underground space structure information. At the same time, the optical fibers under urban roads are irregularly arranged, and the signals received by each DAS sensor cannot be effectively combined, and the imaging results often have large deviations. Summary of the Invention

[0004] The present invention aims to provide a distributed acoustic sensing road cavity detection method based on an active excitation source to solve the technical problems of poor real-time monitoring of underground space structures and limited detection range.

[0005] In order to solve the above technical problems, the specific technical solutions of the present invention are as follows:

[0006] A method for detecting road potholes using distributed acoustic sensing based on an active excitation source comprises the following steps:

[0007] Step 1: Initialize parameters and obtain the DAS signal sampling data sequence x(l,n) to be processed, l = 1, 2, ..., L, n = 1, 2, ..., N, where l is the discrete time index, n is the DAS channel number, L is the number of sampling points, and N is the number of sensor channels;

[0008] Step 2: pre-process the DAS signal sampling data sequence x(l,n) to be processed;

[0009] Step 3: Calculate the short-time zero-crossing rate Z(n) of all channels and remove the DAS channels with abnormal short-time zero-crossing rate;

[0010] Step 4: Extract the signal energy in different areas through the sliding window and locate all active source signal areas D m , m=1,2,...,M, m is the active source signal area index, M is the total number of active source signal areas;

[0011] Step 5: For the current active source signal area D m , extract the signal envelope, estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, determine whether it is a valid active source signal area, if so, estimate the signal source channel n s , otherwise go directly to step 7;

[0012] Step 6: For the effective active source signal area output in step 5, calculate the spatial correlation feature r(n), energy distribution kurtosis feature K, and energy symmetry feature S of the active source signal, and determine whether there is a hole in each channel based on the joint features;

[0013] Step 7: Determine whether the hole detection of all areas is completed. If it is completed, output all channels with detected holes. Otherwise, set m=m+1 and return to step 5.

[0014] Furthermore, in step 1, the following method is used to initialize parameters to obtain the DAS signal sampling data sequence x(l,n) to be processed, l=1,2,...,L,n=1,2,...,N, which specifically includes the following steps:

[0015] Step 1.1: Initialize the various perception parameters for hole detection, including the initialization of the following parameters:

[0016] Sampling frequency f s Initialization: 2<f s <8kHz;

[0017] Bandpass filter upper limit frequency f H Initialized to 50<f H <80Hz;

[0018] Bandpass filter lower limit frequency f L Initialized to 1<f L <15Hz;

[0019] Signal amplitude minimum ε x Initialized to 10 -6 <ε x <10 -4 real number;

[0020] Short-term zero-crossing rate abnormal lower limit threshold z L Initialized to 0.05<z L Real numbers < 0.5;

[0021] Signal area detection sliding window channel window length W c Initialized to 40≤W c Integer ≤ 60;

[0022] Signal area detection sliding window time window length W t Initialized to: 2f s ≤W t ≤5f s integer;

[0023] Signal area detection sliding window channel step S c Initialized to 8≤S c Integer ≤ 12;

[0024] Signal area detection sliding window time step S t Initialized to 0.25f s ≤S t ≤f s integer;

[0025] Signal area energy threshold E S Initialized to 0.01<E S Real numbers < 5;

[0026] Spatial correlation feature hard threshold th hard Initialize 0<th hard Real numbers < 0.5;

[0027] Energy kurtosis feature determination threshold th K Initialized to 0.5<th K Real numbers less than 0.9;

[0028] Energy symmetry feature judgment threshold th S Initialized to 0.1<th S A real number less than 0.45.

[0029] Step 1.2, obtain the DAS signal sampling data sequence x(l,n) to be processed, l = 1, 2, ..., L, n = 1, 2, ..., N: the real-time collected data of L sampling points simultaneously received from N sensors in an array is used as the data sequence x(l,n) to be processed, or the DAS signal data of L sampling points collected by N sensors is extracted from the memory as the data sequence to be processed.

[0030] Furthermore, in step 2, the DAS signal sampling data sequence x(l,n) is preprocessed using the following method, which specifically includes the following steps:

[0031] Step 2.1, design the frequency range to be [f L ,f H ] is used as a bandpass filter to filter the DAS signal sampling data sequence x(l,n) and obtain the filtered data x b (l,n);

[0032] Step 2.2: For the filtered data x b (l,n) is normalized to obtain the preprocessed data x p (l,n):

[0033] x p (l,n)=x b (l,n) / max(x max (n),ε x )

[0034] Among them, x max (n) represents all time frames x of channel n b (l,n),l=1,2,...,the maximum absolute value of L,ε x is the minimum value of the signal amplitude initialized in step 1, and max(·) means taking the larger value of the two.

[0035] 4. The method for detecting road potholes using distributed acoustic sensing based on an active excitation source according to claim 1, wherein in step 3, the short-time zero-crossing rate Z(n) of all channels is calculated using the following method, and DAS channels with abnormal data records are removed, specifically comprising the following steps:

[0036] Step 3.1, calculate the preprocessed data x for each channel n = 1, 2, ..., N p The short-time zero-crossing rate Z(n) of (l,n):

[0037]

[0038] Where sign(·) is the sign function, defined as:

[0039]

[0040] Among them, |·| is the modulus operation;

[0041] Step 3.2: Calculate the upper limit of the short-time zero-crossing rate anomaly Z H and the abnormal lower limit Z L :

[0042]

[0043] Among them, z H and z L are the upper and lower thresholds of the short-time zero-crossing rate anomaly initialized in step 1 respectively;

[0044] Step 3.3: For each channel in n=1, 2, ..., N, determine:

[0045]

[0046] If P(n)=1, then keep the channel data x p (l,n), otherwise the channel data is discarded and the data x after removing the abnormal channel is obtained. d (l,n).

[0047] Furthermore, in step 4, the following method is used to extract the signal energy in each area of the DAS data and locate the active source signal area:

[0048] Step 4.1: Signal region detection: Sliding window traverses the entire data segment and obtains the average energy E(p,q) in each window;

[0049]

[0050] in, is the time domain window index, is the channel window index, is the floor function, W t The length of the sliding window for detecting the signal region initialized in step 1, W c The signal area initialized in step 1 detects the sliding window channel window length; S tDetect the sliding window time step for the signal region initialized in step 1; S c Detect the sliding window channel step for the signal region initialized in step 1;

[0051] Step 4.2: Traverse the energy of each window and select the energy E(p,q)>E S The M windows are used as the active source signal detection area D m ,m=1,2,...,M,E S is the signal area energy threshold initialized in step 1;

[0052] Step 4.3: Initialize the active source signal detection area index m=1.

[0053] Furthermore, in step 5, the following method is used to extract the signal envelope in the current area, and estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, the signal source channel n s , specifically including the following steps:

[0054] Step 5.1: Use smoothing filtering to extract the current signal detection area D m The preprocessed data x d (l,n) each channel signal envelope;

[0055] Step 5.2: Perform peak detection on the extracted signal envelopes of each channel, estimate the start time t1, end time t2, start channel C1 and end channel C2 of the signal, and calculate the start time frame of the signal. End time frame And store the signal time range in Ω t , that is, Ω t =[L1,L2], the signal channel range is stored in Ω c , that is, Ω c =[C1,C2];

[0056] Step 5.3: Determine whether L1≠L2 and C1≠C2 are true. If so, proceed to step 5.4; otherwise, proceed to step 7.

[0057] Step 5.4: Estimate the signal source channel n based on the signal spatiotemporal detection results s :

[0058]

[0059] Furthermore, in step 6, based on the spatiotemporal position of the signal extracted in step 5, the following method is used to extract the spatial correlation feature r(n), the energy distribution kurtosis feature K, and the energy symmetry feature S of the active source signal, specifically comprising the following steps:

[0060] Step 6.1: Calculate the spatial correlation characteristic r(n) of the active source signal:

[0061]

[0062] in, Indicates that channel n is in Ω t The time domain average amplitude within the range is calculated as:

[0063]

[0064] Step 6.2: Extract Ω c The energy distribution kurtosis characteristic K in the range is:

[0065]

[0066] in, is the kurtosis of the energy distribution, calculated as:

[0067]

[0068] Among them, E ave (n) represents the time-domain average energy of the active source signal of channel n, represents the average energy of the active source signal, which is calculated as:

[0069]

[0070] Step 6.3: Extract signal source channel n s The energy distribution symmetry characteristic S nearby:

[0071]

[0072] where ξ is the mean square error about the symmetry of the energy distribution, calculated as:

[0073]

[0074] Step 6.4: For the signal channel range Ω c For each channel in , determine:

[0075]

[0076] Among them, soft is the soft threshold of spatial correlation feature, calculated as Ω c The average value of r(n) of all spatial correlation feature valley channels in ,th hard is the hard threshold of the spatial correlation feature initialized in step 1, th K The energy kurtosis feature determination threshold initialized in step 1, th Sis the energy symmetry feature determination threshold initialized in step 1, and min(·) represents the smaller value of the two. If H(n) = 1, it is determined that there is a hole in channel n and it is recorded, otherwise it is not recorded.

[0077] Furthermore, in step 7, the following method is used to determine whether the data processing is complete, specifically including the following steps: Step 7.1, determine whether the active source signal area index m=M is established. If so, the DAS data processing is completed and all channels with detected holes are output. Otherwise, the active source signal area index m=m+1 is set and the process returns to step 5.

[0078] The present invention provides a distributed fiber optic acoustic sensing road cavity detection method based on an active excitation source, which has the following advantages:

[0079] 1. This invention utilizes the characteristic that when the DAS channel coupling is abnormal or there is sustained abnormal shaking, the main frequency component of the signal deviates. Through short-term zero-crossing rate analysis, abnormal DAS channels are effectively identified and eliminated, as shown in step 3, thereby improving the robustness of road cavity detection.

[0080] 2. The present invention performs preliminary positioning of the signal area based on sliding window energy detection and performs spatiotemporal positioning of the active source signal only within the selected target area, as shown in steps 4 and 5. This avoids the extraction of signal parameters across the entire spatiotemporal range, significantly reduces the computational complexity of the method, and improves the application efficiency of the method in large-scale urban road cavity detection.

[0081] 3. The present invention makes full use of the propagation characteristics of active source signals in underground space. As shown in the corresponding processing steps in step 6, the extracted characteristic parameters have good discrimination between road cavity structures and non-cavity structures, and can still well reflect the structural characteristics of underground space at low signal-to-noise ratios. By combining spatial correlation characteristics and energy distribution characteristics, the robustness and accuracy of road cavity detection are further improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 Schematic diagram of the process of the present invention;

[0083] Figure 2 This is a local channel time domain waveform diagram of the DAS signal sampling data sequence of Example 1;

[0084] Figure 3 This is a time domain waveform diagram of a local channel of the preprocessed signal of Example 1;

[0085] Figure 4 This is a full-channel short-time zero-crossing rate detection diagram of Example 1;

[0086] Figure 5This is the local channel signal energy distribution diagram of Example 1;

[0087] Figure 6 This is the time domain waveform diagram of the active source signal of Example 1;

[0088] Figure 7 This is the spatial correlation characteristic diagram of the active source signal in Example 1;

[0089] Figure 8 This is a local channel time domain waveform diagram of the DAS signal sampling data sequence of Example 2;

[0090] Figure 9 This is a time domain waveform diagram of a local channel of the preprocessed signal of Example 2;

[0091] Figure 10 This is a full-channel short-time zero-crossing rate detection diagram of Example 2;

[0092] Figure 11 This is the local channel signal energy distribution diagram of Example 2;

[0093] Figure 12 This is the time domain waveform diagram of the active source signal of Example 2;

[0094] Figure 13 This is the spatial correlation characteristic diagram of the active source signal in Example 2. DETAILED DESCRIPTION

[0095] In order to better understand the purpose, structure and function of the present invention, the following is a further detailed description of a distributed acoustic sensing road cavity detection method based on an active excitation source of the present invention in conjunction with the accompanying drawings.

[0096] like Figure 1 As shown, the present invention proposes a distributed acoustic sensing road cavity detection method based on an active excitation source, comprising the following steps:

[0097] Step 1: Initialize parameters and obtain the DAS signal sampling data sequence x(l,n) to be processed, l = 1, 2, ..., L, n = 1, 2, ..., N, where l is the discrete time index, n is the DAS channel number, L is the number of sampling points, and N is the number of sensor channels;

[0098] Step 2: pre-process the DAS signal sampling data sequence x(l,n) to be processed;

[0099] Step 3: Calculate the short-time zero-crossing rate Z(n) of all channels and remove the DAS channels with abnormal short-time zero-crossing rate;

[0100] Step 4: Extract the signal energy in different areas through the sliding window and locate all active source signal areas D m, m=1,2,...,M, m is the active source signal area index, M is the total number of active source signal areas;

[0101] Step 5: For the current active source signal area D m , extract the signal envelope, estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, determine whether it is a valid active source signal area, if so, estimate the signal source channel n s , otherwise go directly to step 7;

[0102] Step 6: For the effective active source signal area output in step 5, calculate the spatial correlation feature r(n), energy distribution kurtosis feature K, and energy symmetry feature S of the active source signal, and determine whether there is a hole in each channel based on the joint features;

[0103] Step 7: Determine whether the hole detection of all areas is completed. If it is completed, output all channels with detected holes. Otherwise, set m=m+1 and return to step 5.

[0104] Furthermore, in step 1, the following method is used to initialize parameters to obtain the DAS signal sampling data sequence x(l,n) to be processed, l=1,2,...,L,n=1,2,...,N, which specifically includes the following steps:

[0105] Step 1.1: Initialize the various perception parameters for hole detection, including the initialization of the following parameters:

[0106] Sampling frequency f s Initialization: 2<f s <8kHz;

[0107] Bandpass filter upper limit frequency f H Initialized to 50<f H <80Hz;

[0108] Bandpass filter lower limit frequency f L Initialized to 1<f L <15Hz;

[0109] Signal amplitude minimum ε x Initialized to 10 -6 <ε x <10 -4 real number;

[0110] Short-term zero-crossing rate abnormal lower limit threshold z L Initialized to 0.05<z L Real numbers < 0.5;

[0111] Signal area detection sliding window channel window length W cInitialized to 40≤W c Integer ≤ 60;

[0112] Signal area detection sliding window time window length W t Initialized to: 2f s ≤W t ≤5f s integer;

[0113] Signal area detection sliding window channel step S c Initialized to 8≤S c Integer ≤ 12;

[0114] Signal area detection sliding window time step S t Initialized to 0.25f s ≤S t ≤f s integer;

[0115] Signal area energy threshold E S Initialized to 0.01<E S Real numbers < 5;

[0116] Spatial correlation feature hard threshold th hard Initialize 0<th hard Real numbers < 0.5;

[0117] Energy kurtosis feature determination threshold th K Initialized to 0.5<th K Real numbers less than 0.9;

[0118] Energy symmetry feature judgment threshold th S Initialized to 0.1<th S A real number less than 0.45.

[0119] Step 1.2, obtaining a DAS signal sampling data sequence x(l,n) to be processed, where l=1,2,...,L,n=1,2,...,N: real-time data collected at L sampling points simultaneously received from N sensors in an array is used as the data sequence x(l,n) to be processed, or DAS signal data of L sampling points collected by N sensors is extracted from a memory as the data sequence to be processed;

[0120] In step 1, in order to take into account both the computational complexity and the estimation accuracy of the present invention, the sampling frequency f s The preferred value is 4kHz, and the upper limit frequency of the bandpass filter f H The preferred value is 60Hz, and the lower limit frequency f of the bandpass filter is L The preferred value is 5Hz, and the minimum signal amplitude ε x The preferred value is 10 -5, short-term zero-crossing rate abnormal upper limit threshold z H The preferred value of is 2, and the short-term zero-crossing rate abnormal lower limit threshold z L The preferred value is 0.5, and the signal area detection sliding window channel size W c The preferred value is 50, and the signal area detection sliding window time size W t The preferred value is 12, and the signal area detection sliding window channel step window length S c The preferred value is 20, and the signal area detection sliding window time step is S t The preferred value of is 3, and the signal area energy threshold E S The preferred value is 0.5, and the hard threshold of spatial correlation feature th hard The preferred value is 0.3, and the energy kurtosis feature determination threshold th K The preferred value of is 0.6, and the energy symmetry feature determination threshold th S The preferred value of is 0.3.

[0121] Furthermore, in step 2, the DAS signal sampling data sequence x(l,n) is preprocessed using the following method, which specifically includes the following steps:

[0122] Step 2.1, design the frequency range to be [f L ,f H ] is used as a bandpass filter to filter the DAS signal sampling data sequence x(l,n) and obtain the filtered data x b (l,n);

[0123] Step 2.2: For the filtered data x b (l,n) is normalized to obtain the preprocessed data x p (l,n):

[0124] x p (l,n)=x b (l,n) / max(x max (n),ε x )

[0125] Among them, x max (n) represents all time frames x of channel n b (l,n),l=1,2,...,the maximum absolute value of L,ε x is the minimum value of the signal amplitude initialized in step 1, and max(·) means taking the larger value of the two.

[0126] 4. The method for detecting road potholes using distributed acoustic sensing based on an active excitation source according to claim 1, wherein in step 3, the short-time zero-crossing rate Z(n) of all channels is calculated using the following method, and DAS channels with abnormal data records are removed, specifically comprising the following steps:

[0127] Step 3.1, calculate the preprocessed data x for each channel n = 1, 2, ..., N p The short-time zero-crossing rate Z(n) of (l,n):

[0128]

[0129] Where sign(·) is the sign function, defined as:

[0130]

[0131] Among them, |·| is the modulus operation;

[0132] Step 3.2: Calculate the upper limit of the short-time zero-crossing rate anomaly Z H and the abnormal lower limit Z L :

[0133]

[0134] Among them, z H and z L are the upper and lower thresholds of the short-time zero-crossing rate anomaly initialized in step 1 respectively;

[0135] Step 3.3: For each channel in n=1, 2, ..., N, determine:

[0136]

[0137] If P(n)=1, then keep the channel data x p (l,n), otherwise the channel data is discarded and the data x after removing the abnormal channel is obtained. d (l,n).

[0138] Furthermore, in step 4, the following method is used to extract the signal energy in each area of the DAS data and locate the active source signal area:

[0139] Step 4.1: Signal region detection: Sliding window traverses the entire data segment and obtains the average energy E(p,q) in each window;

[0140]

[0141] in, is the time domain window index, is the channel window index, is the floor function, W t The length of the sliding window for detecting the signal region initialized in step 1, W c The signal area initialized in step 1 detects the sliding window channel window length; S t Detect the sliding window time step for the signal region initialized in step 1; S c Detect the sliding window channel step for the signal region initialized in step 1;

[0142] Step 4.2: Traverse the energy of each window and select the energy E(p,q)>E S The M windows are used as the active source signal detection area D m ,m=1,2,...,M,E S is the signal area energy threshold initialized in step 1;

[0143] Step 4.3: Initialize the active source signal detection area index m=1.

[0144] Furthermore, in step 5, the following method is used to extract the signal envelope in the current area, and estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, the signal source channel n s , specifically including the following steps:

[0145] Step 5.1: Use smoothing filtering to extract the current signal detection area D m The preprocessed data x d (l,n) each channel signal envelope;

[0146] Step 5.2: Perform peak detection on the extracted signal envelopes of each channel, estimate the start time t1, end time t2, start channel C1 and end channel C2 of the signal, and calculate the start time frame of the signal. End time frame And store the signal time range in Ω t , that is, Ω t =[L1,L2], the signal channel range is stored in Ω c , that is, Ω c =[C1,C2];

[0147] Step 5.3: Determine whether L1≠L2 and C1≠C2 are true. If so, proceed to step 5.4; otherwise, proceed to step 7.

[0148] Step 5.4: Estimate the signal source channel n based on the signal spatiotemporal detection results s :

[0149]

[0150] Furthermore, in step 6, based on the spatiotemporal position of the signal extracted in step 5, the following method is used to extract the spatial correlation feature r(n), the energy distribution kurtosis feature K, and the energy symmetry feature S of the active source signal, specifically comprising the following steps:

[0151] Step 6.1: Calculate the spatial correlation characteristic r(n) of the active source signal:

[0152]

[0153] in, Indicates that channel n is in Ω t The time domain average amplitude within the range is calculated as:

[0154]

[0155] Step 6.2: Extract Ω c The energy distribution kurtosis characteristic K in the range is:

[0156]

[0157] in, is the kurtosis of the energy distribution, calculated as:

[0158]

[0159] Among them, E ave (n) represents the time-domain average energy of the active source signal of channel n, represents the average energy of the active source signal, which is calculated as:

[0160]

[0161] Step 6.3: Extract signal source channel n s The energy distribution symmetry characteristic S nearby:

[0162]

[0163] where ξ is the mean square error about the symmetry of the energy distribution, calculated as:

[0164]

[0165] Step 6.4: For the signal channel range Ω c For each channel in , determine:

[0166]

[0167] Among them, soft is the soft threshold of spatial correlation feature, calculated as Ω cThe average value of r(n) of all spatial correlation feature valley channels in ,th hard is the hard threshold of the spatial correlation feature initialized in step 1, th K The energy kurtosis feature determination threshold initialized in step 1, th S is the energy symmetry feature determination threshold initialized in step 1, and min(·) represents the smaller value of the two. If H(n) = 1, it is determined that there is a hole in channel n and it is recorded, otherwise it is not recorded.

[0168] Furthermore, in step 7, the following method is used to determine whether the data processing is complete, specifically including the following steps: Step 7.1, determine whether the active source signal area index m=M is established. If so, the DAS data processing is completed and all channels with detected holes are output. Otherwise, the active source signal area index m=m+1 is set and the process returns to step 5.

[0169] Example 1:

[0170] Field verification of DAS data collection and cavity detection was carried out on Youhu Road, Jiangning District, Nanjing.

[0171] The following DAS data segment consists of 1200 DAS channels and lasts for 18 seconds for hole detection:

[0172] According to step 1, set the sampling frequency f s =4kHz, bandpass filter upper limit frequency f H =60Hz, bandpass filter lower limit frequency f L =5Hz, minimum signal amplitude ε x =10 -5 , short-term zero-crossing rate abnormal upper limit threshold z H =2, short-term zero-crossing rate abnormal lower limit threshold z L =0.5, signal area detection sliding window channel size W c =50, signal area detection sliding window time size W t =2f s , signal area detection sliding window channel step window length S c =20, signal area detection sliding window time step S t =0.5f s , signal area energy threshold E S =0.1, spatial correlation feature hard threshold th hard =0.4, energy kurtosis feature determination threshold th K =0.6, energy symmetry feature judgment threshold th S =0.3, the local channel time domain waveform of the DAS signal sampling data sequence read is as follows Figure 2 As shown;

[0173] According to step 2, the design frequency range is [f L ,f H ] is used as a bandpass filter to filter the DAS signal sampling data sequence x(l,n) and to filter the filtered data x b (l,n) is normalized to obtain the preprocessed data x p (l,n) Figure 3 As shown;

[0174] According to step 3, calculate the short-time zero-crossing rate Z(n) of all channels. The calculation results are as follows: Figure 4 As shown, the short-term zero-crossing rate is greater than the abnormal upper limit Z H The sum is less than the abnormal lower limit Z L The channel data is updated to get the updated data x d (l,n);

[0175] According to step 4, the sliding window extracts the signal energy in each area of the DAS data. The extracted local signal energy results are as follows: Figure 5 As shown, an active source signal area D1 is located as Figure 6 As shown, its channel range is [770,820] and time range is [1s,3s];

[0176] According to step 5, the signal envelope of the signal detection area is extracted, and the signal start time t1 = 1.7s, the end time t2 = 2.0s, the signal start channel C1 = 720, the end channel C2 = 770, and the signal source channel n are estimated. s =745;

[0177] According to step 6, the spatial correlation feature r(n), energy distribution kurtosis feature K and energy symmetry feature S of the active source signal are extracted, and it is determined whether there is a hole in the channel. The extracted spatial correlation feature r(n) of the active source signal is as follows: Figure 7 As shown, the spatial correlation feature soft threshold th is calculated soft =0.53, extract Ω c The energy distribution kurtosis characteristic K in the range is K = 0.23, and the signal source channel n s The energy distribution symmetry characteristic near S = 0.08, for the signal channel range Ω c Make a hollow judgment for each channel in;

[0178] According to step 7, it is determined whether the data processing is completed. If the processing is completed, all channels with detected holes are output. Otherwise, the active source signal area index m=m+1 is set and the process returns to step 5.

[0179] There is one road hole in the experimental section of this data, located in channel 799. This method successfully detects the channel where the hole is located.

[0180] Example 2:

[0181] DAS data collection and field verification of void-free features were carried out on Daishan Road in Yuhuatai District, Nanjing.

[0182] The following DAS data segment consists of 400 DAS channels and lasts for 18 seconds for hole detection:

[0183] According to step 1, set the sampling frequency f s =4kHz, bandpass filter upper limit frequency f H =60Hz, bandpass filter lower limit frequency f L =5Hz, minimum signal amplitude ε x =10 -5 , short-term zero-crossing rate abnormal upper limit threshold z H =2, short-term zero-crossing rate abnormal lower limit threshold z L =0.5, signal area detection sliding window channel size W c =30, signal area detection sliding window time size W t =2f s , signal area detection sliding window channel step window length S c =10, signal area detection sliding window time step S t =0.5f s , signal area energy threshold E S =2.5, hard threshold of spatial correlation feature th hard =0.5, energy kurtosis feature determination threshold th K =0.6, energy symmetry feature judgment threshold th S =0.3, the local channel time domain waveform of the DAS signal sampling data sequence read is as follows Figure 8 As shown;

[0184] According to step 2, the design frequency range is [f L ,f H ] is used as a bandpass filter to filter the DAS signal sampling data sequence x(l,n) and to filter the filtered data x b (l,n) is normalized to obtain the preprocessed data x p (l,n) Figure 9 As shown;

[0185] According to step 3, calculate the short-time zero-crossing rate Z(n) of all channels. The calculation results are as follows: Figure 4 As shown, the short-term zero-crossing rate is greater than the abnormal upper limit ZH The sum is less than the abnormal lower limit Z L The channel data is updated to get the updated data x d (l,n);

[0186] According to step 4, the sliding window extracts the signal energy in each area of the DAS data. The extracted local signal energy is as follows: Figure 11 As shown, an active source signal area D1 is located as Figure 12 As shown, its channel range is [240,270] and time range is [8s,10s];

[0187] According to step 5, the signal envelope of the signal detection area is extracted, and the signal start time t1 = 8.5s, the end time t2 = 8.7s, the signal start channel C1 = 240, the end channel C2 = 270, and the signal source channel n are estimated. s =255;

[0188] According to step 6, the spatial correlation feature r(n), energy distribution kurtosis feature K and energy symmetry feature S of the active source signal are extracted, and it is determined whether there is a hole in the channel. The extracted spatial correlation feature r(n) of the active source signal is as follows: Figure 7 As shown, the spatial correlation feature soft threshold th is calculated soft =0.64, extract Ω c The energy distribution kurtosis characteristic K in the range is K = 0.89, and the signal source channel n s The energy distribution symmetry characteristic near S = 0.77, for the signal channel range Ω c Make a hollow judgment for each channel in;

[0189] According to step 7, it is determined whether the data processing is completed. If the processing is completed, all channels with detected holes are output. Otherwise, the active source signal area index m=m+1 is set and the process returns to step 5.

[0190] This experiment is a closed road test verification, and a road section without holes is selected for feature verification. This method does not detect the existence of holes in this area.

[0191] It will be understood that the present invention is described by way of some embodiments, and it will be appreciated by those skilled in the art that various changes or equivalent substitutions may be made to these features and embodiments without departing from the spirit and scope of the present invention. In addition, under the teachings of the present invention, these features and embodiments may be modified to adapt to specific circumstances and materials without departing from the spirit and scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are intended to be protected by the present invention.

Claims

1. A distributed acoustic sensing road cavity detection method based on active excitation source, characterized in that: The following steps are involved: Step 1: Initialize parameters and obtain the DAS signal sampling data sequence x(l,n) to be processed, l = 1, 2, ..., L, n = 1, 2, ..., N, where l is the discrete time index, n is the DAS channel number, L is the number of sampling points, and N is the number of sensor channels; Step 2: pre-process the DAS signal sampling data sequence x(l,n) to be processed; Step 3: Calculate the short-time zero-crossing rate Z(n) of all channels and remove the DAS channels with abnormal short-time zero-crossing rate; Step 4: Extract the signal energy in different areas through the sliding window and locate all active source signal areas D m , m=1,2,...,M, m is the active source signal area index, M is the total number of active source signal areas; Step 5: For the current active source signal area D m , extract the signal envelope, estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, determine whether it is a valid active source signal area, if so, estimate the signal source channel n s , otherwise go directly to step 7; Step 6: For the effective active source signal area output in step 5, calculate the spatial correlation feature r(n), energy distribution kurtosis feature K, and energy symmetry feature S of the active source signal, and determine whether there is a hole in each channel based on the joint features; Step 7: Determine whether the hole detection of all areas is completed. If it is completed, output all channels with detected holes. Otherwise, set m=m+1 and return to step 5.

2. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 1, the following method is used to initialize parameters and obtain the DAS signal sampling data sequence x(l,n) to be processed, l=1,2,...,L,n=1,2,...,N, which specifically includes the following steps: Step 1.1: Initialize the various perception parameters for hole detection, including the initialization of the following parameters: Sampling frequency f s Initialization: 2<f s <8kHz; Bandpass filter upper limit frequency f H Initialized to 50<f H <80Hz; Bandpass filter lower limit frequency f L Initialized to 1<f L <15Hz; Signal amplitude minimum ε x Initialized to 10 -6 <ε x <10 -4 real number; Short-term zero-crossing rate abnormal upper limit threshold z H Initialized to 1.5<z H Real numbers < 5; Short-term zero-crossing rate abnormal lower limit threshold z L Initialized to 0.05<z L Real numbers < 0.5; Signal area detection sliding window channel window length W c Initialized to 40≤W c Integer ≤ 60; Signal area detection sliding window time window length W t Initialized to: 2f s ≤W t ≤5f s integer; Signal area detection sliding window channel step S c Initialized to 8≤S c Integer ≤ 12; Signal area detection sliding window time step S t Initialized to 0.25f s ≤S t ≤f s integer; Signal area energy threshold E S Initialized to 0.01<E S Real numbers < 5; Spatial correlation feature hard threshold th hard Initialize 0<th hard Real numbers < 0.5; Energy kurtosis feature determination threshold th K Initialized to 0.5<th K Real numbers less than 0.9; Energy symmetry feature judgment threshold th S Initialized to 0.1<th S Real numbers <0.45; Step 1.2, obtaining the DAS signal sampling data sequence x(l,n) to be processed, l = 1, 2, ..., L, n = 1, 2, ..., N: real-time data of L sampling points simultaneously received from N sensors in an array is used as the data sequence x(l,n) to be processed, or DAS signal data of L sampling points collected by N sensors is extracted from a memory as the data sequence to be processed.

3. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 2, the DAS signal sampling data sequence x(l,n) is preprocessed using the following method, which specifically includes the following steps: Step 2.1, design the frequency range to be [f L ,f H ] is used as a bandpass filter to filter the DAS signal sampling data sequence x(l,n) and obtain the filtered data x b (l,n); Step 2.2: For the filtered data x b (l,n) is normalized to obtain the preprocessed data x p (l,n): x p (l,n)=x b (l,n) / max(x max (n),ε x ) Among them, x max (n) represents all time frames x of channel n b (l,n),l=1,2,...,the maximum absolute value of L,ε x is the minimum value of the signal amplitude initialized in step 1, and max(·) means taking the larger value of the two.

4. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 3, the short-time zero-crossing rate Z(n) of all channels is calculated using the following method, and the DAS channels with abnormal data records are removed. Specifically, the following steps are included: Step 3.1, calculate the preprocessed data x for each channel n = 1, 2, ..., N p The short-time zero-crossing rate Z(n) of (l,n): Where sign(·) is the sign function, defined as: Among them, |·| is the modulus operation; Step 3.2: Calculate the upper limit of the short-time zero-crossing rate anomaly Z H and the abnormal lower limit Z L : Among them, z H and z L are the upper and lower thresholds of the short-time zero-crossing rate anomaly initialized in step 1 respectively; Step 3.3: For each channel in n=1, 2, ..., N, determine: If P(n)=1, then keep the channel data x p (l,n), otherwise the channel data is discarded and the data x after removing the abnormal channel is obtained. d (l,n).

5. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 4, the following method is used to extract the signal energy in each area of the DAS data and locate the active source signal area: Step 4.1: Signal region detection: Sliding window traverses the entire data segment and obtains the average energy E(p,q) in each window; in, is the time domain window index, is the channel window index, is the floor function, W t The length of the sliding window for detecting the signal region initialized in step 1, W c The signal area initialized in step 1 detects the sliding window channel window length; S t Detect the sliding window time step for the signal region initialized in step 1; S c Detect the sliding window channel step for the signal region initialized in step 1; Step 4.2: Traverse the energy of each window and select the energy E(p,q)>E S The M windows are used as the active source signal detection area D m ,m=1,2,...,M,E S is the signal area energy threshold initialized in step 1; Step 4.3: Initialize the active source signal detection area index m=1.

6. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 5, the following method is used to extract the signal envelope in the current area, and estimate the signal start time t1 and end time t2 as well as the signal start channel C1 and end channel C2, and the signal source channel n s , specifically including the following steps: Step 5.1: Use smoothing filtering to extract the current signal detection area D m The preprocessed data x d (l,n) each channel signal envelope; Step 5.2: Perform peak detection on the extracted signal envelopes of each channel, estimate the start time t1, end time t2, start channel C1 and end channel C2 of the signal, and calculate the start time frame of the signal. End time frame And store the signal time range in Ω t , that is, Ω t =[L1,L2], the signal channel range is stored in Ω c , that is, Ω c =[C1,C2]; Step 5.3: Determine whether L1≠L2 and C1≠C2 are true. If so, proceed to step 5.4; otherwise, proceed to step 7. Step 5.4: Estimate the signal source channel n based on the signal spatiotemporal detection results s :

7. The method for detecting road holes using distributed acoustic sensing based on an active excitation source according to claim 1, characterized in that: In step 6, based on the spatiotemporal position of the signal extracted in step 5, the following method is used to extract the spatial correlation feature r(n), energy distribution kurtosis feature K, and energy symmetry feature S of the active source signal, specifically including the following steps: Step 6.1: Calculate the spatial correlation characteristic r(n) of the active source signal: in, Indicates that channel n is in Ω t The time domain average amplitude within the range is calculated as: Step 6.2: Extract Ω c The energy distribution kurtosis characteristic K in the range is: in, is the kurtosis of the energy distribution, calculated as: Among them, E ave (n) represents the time-domain average energy of the active source signal of channel n, represents the average energy of the active source signal, which is calculated as: Step 6.3: Extract signal source channel n s The energy distribution symmetry characteristic S nearby: where ξ is the mean square error about the symmetry of the energy distribution, calculated as: Step 6.4: For the signal channel range Ω c For each channel in , determine: Among them, soft is the soft threshold of spatial correlation feature, calculated as Ω c The average value of r(n) of all spatial correlation feature valley channels in ,th hard is the hard threshold of the spatial correlation feature initialized in step 1, th K The energy kurtosis feature determination threshold initialized in step 1, th S is the energy symmetry feature determination threshold initialized in step 1, and min(·) represents the smaller value of the two. If H(n) = 1, it is determined that there is a hole in channel n and it is recorded, otherwise it is not recorded.

8. The method for detecting road voids using distributed acoustic sensing based on an active excitation source according to claim 1, wherein: In step 7, the following method is used to determine whether the data has been processed, which specifically includes the following steps: Step 7.1: Determine whether the active source signal area index m=M. If so, DAS data processing ends and all channels with detected holes are output. Otherwise, set the active source signal area index m=m+1 and return to step 5.

Images