CN102866429A - Method for determining groundwater occurrence - Google Patents

Abstract

The invention discloses a method for determining groundwater occurrence. The method comprises the following steps of: removing interference signals from acquired geological radar data; carrying out computation of the low-frequency signal energy occupancy rate on the geological radar data subjected to interference signal removing, thus obtaining water-containing attribute data; and carrying out profile imaging treatment on the water-containing attribute data, thus displaying the condition of the groundwater occurrence. After the method for determining the groundwater occurrence provided by the invention is used, not only can the groundwater occurrence condition in a quaternary system stratum be obtained, but also the data obtained by the method is more exact and reliable compared with the data obtained by the direct explanation method in the prior art.

Description

Definite method that a kind of underground water distributes

Technical field

The present invention relates to the geology detecting field, particularly a kind of definite method of underground water distribution.

Background technology

Usually before and after coal mining, all will carry out geology detecting to stratal configuration, a main project of wherein surveying is exactly the distribution situation of determining underground water.

Geological radar is the major equipment that carries out geology detecting, and this equipment utilization wideband electromagnetic ripple is realized detection to underground latent objective body by the mechanism of transmission of electromagnetic wave in underground medium.Because it is simple to operate that geological radar has, detection accuracy is high, and the characteristics such as nondestructive detecting have widespread use at aspects such as engineering hydrogeologic survey and construction quality detections.

The main method of geologic radar detection technology on data is processed is to adopt the method for directly explaining at present, namely by after the source book of geologic radar detection is made the simple process such as some amplifications, filtering, stack, external appearance characteristic according to radar appearance, as reflect power, phase characteristic, variation characteristic of axle etc. information is directly made quantitative and qualitative analysis to reflected signal and is explained in the same way.

Because coal mining can cause ground subsidence, has destroyed the Quaternary Strata structure, causes loss and the variation of surface water, if adopt direct interpretation procedure, can't obtain the distribution situation of water in the Quaternary Strata, namely can't obtain underground water distribution situation.

Summary of the invention

In view of this, fundamental purpose of the present invention is to provide a kind of definite method of underground water distribution, and the method can obtain the distribution situation of water in the Quaternary Strata.

For achieving the above object, definite method that this underground water of the present invention distributes comprises the steps:

A, the Georadar Data of the Quaternary Strata structure that collects is removed undesired signal process;

B, to having removed the Georadar Data of undesired signal, calculate low frequency signal energy occupation rate, obtain moisture attribute data;

C, moisture attribute data is carried out the section imaging processing, explicitly descend the water distribution situation.

Preferably, the described removal undesired signal of steps A is processed and is comprised:

A1, the Georadar Data that collects is carried out offset correction, remove the signal drift noise of instrument self;

A2, the curve after the drift correction is carried out wavelet transformation, HF noise signal is suppressed;

A3, the signal behind the wavelet transformation hanged down cut filtering, the interference of excision signal DC component.

Preferably, the signal drift noise of described steps A 1 described removal instrument self comprises:

The drift parameter of each point in A11, the calculating Georadar Data;

A12, the drift parameter that original signal in the Georadar Data is had a few connect into the drift curve;

A13, original signal and drift curve are subtracted each other, form curve after the drift correction.

Preferably, the computing method of steps A 11 described drift parameters are: centered by the preset after being positioned at drift correction in the middle of the curve, set up a window that is formed by 101 sampling points, all original signal data in the window are averaged as the drift parameter of this sampling point.

Preferably, in the described steps A 2, choose the Moret wavelet function, scale parameter is 2, carries out wavelet transformation.

Preferably, in the described steps A 3, low to cut filtering parameter be 20MHz.

Preferably, described step B comprises:

B1, carry out the rolling time window and ask spectra calculation removing Georadar Data signal after the undesired signal; The rolling time window that B2, basis calculate is asked power spectrum, calculates low frequency signal energy occupation rate and calculates, and obtains moisture attribute data.

Preferably, described step B1 comprises:

B11, in geological radar image data time window, to limited discrete signal

Make auto-correlation:

$R_{\mathrm{ff}}\left(n\right)=\frac{1}{N}\underset{k=0}{\overset{N-1-n}{\mathrm{\Σ}}}f\left(k\right)f(k+n),$ 0≤n≤N-1

B12, to the autocorrelation certificate

Make FFT and calculate, try to achieve discrete power spectrum

Only get its positive frequency part

Enveloping curve as the power spectrum curve P of signal in the time window _f(ω)

Wherein N is number of samples; Described time window rolls downwards from initial sampling point.

Preferably, the described calculating low frequency signal of step B2 energy occupation rate is: calculate 1/4th following power spectrum energy of Ground Penetrating Radar Signal antenna dominant frequency and the ratio of the whole power spectrum energy of antenna.

Preferably, described step B further comprises B3: moisture attribute data is carried out two-dimentional running mean, eliminate the wherein interference of geology catastrophe point.

Preferably, describedly moisture attribute data carried out two-dimentional running mean be:

$x_{\mathrm{\Σ}}\left(t\right)=\frac{1}{N2-N1+1+M2-M1+1}\underset{i=N1}{\overset{N2}{\mathrm{\Σ}}}\underset{j=M1}{\overset{M2}{\mathrm{\Σ}}}{x}_{\mathrm{ij}}\left(t\right)$ N1＜N2?M1＜M2

Wherein: N1 is the start channel number, and N2 stops the road number; M1 is initial number of samples; M2 stops number of samples.

Preferably, described step C comprises:

C1, moisture attribute data is carried out normalized;

C2, carry out colored section imaging, the moisture attribute data after the normalized is shown in the bitmap mode.

Preferably, described step C2 comprises:

C21, carry out color range modulation, construct 8 scale-of-two color tables, have 256 color ranges, the moisture attribute data after the interior normalized of the corresponding preset range of each color;

C22, carry out colored section imaging according to the color range of C21, the moisture attribute data after the normalized is shown in the bitmap mode.

As seen from the above technical solutions, definite method that this underground water of the present invention distributes is processed owing to first the Quaternary Strata structure Georadar Data that collects is removed undesired signal; To having removed the Georadar Data of undesired signal, calculate low frequency signal energy occupation rate again, obtain moisture attribute data; At last moisture attribute data is carried out the section imaging processing, explicitly descend the water distribution situation.Therefore, not only can obtain the distribution situation of water in the Quaternary Strata, and the data that obtain than the method for the direct explanation of prior art more accurately, reliable.

Description of drawings

Fig. 1 is the processing flow chart of definite method one preferred embodiment of underground water distribution of the present invention;

Fig. 2 a is middle geologic radar detection raw data curve embodiment illustrated in fig. 1 and the drift curve map that calculates;

Fig. 2 b is the data and curves figure after the curve of Fig. 2 a passes through offset correction;

Fig. 2 c is the data and curves figure after the curve of Fig. 2 b passes through wavelet transformation;

Fig. 2 d is that the curve of Fig. 2 c is cut filtered data and curves figure through hanging down;

Fig. 2 e is the stratum low frequency signal energy occupation rate change curve according to the raw data acquisition of Fig. 2 a;

Fig. 3 is the original section image of geologic radar detection corresponding to Fig. 2 a;

Fig. 4 is the pretreated profile image of original section image process shown in Figure 3;

Fig. 5 is the profile image of the bitmap mode of Fig. 4 after normalization;

Fig. 6 is that Fig. 5 is through the water percentage profile image after the two-dimentional running mean.

Embodiment

Developing simultaneously referring to accompanying drawing, the present invention is described in detail for specific embodiment.

The invention provides definite method that a kind of underground water distributes, the method can obtain the distribution situation of water in the Quaternary Strata.Below lifting a specific embodiment is elaborated.

As shown in Figure 1, definite method that the underground water of a preferred embodiment of the present invention distributes adopts the geological radar image data, specifically comprises the steps:

Step 101, the Quaternary Strata structure Georadar Data that collects is carried out offset correction, remove the signal drift noise of instrument self.

In this step, the original signal that geological radar collects is shown in Fig. 2 a, and the profile image of its raw data as shown in Figure 3.In the present embodiment, the method for the profile image of geological radar image data and demonstration raw data belongs to prior art, repeats no more here.

In this step, the method for removing the signal drift noise of instrument self comprises three steps:

The first step: the drift parameter that calculates each point in the Georadar Data.

Particularly, the preset N point in the middle of the original signal of Fig. 2 a centered by the N point, is set up a window W who is comprised of 101 sampling points as example, and all original signal data in the window are averaged as the drift parameter of this sampling point.

Second step: the drift parameter that original signal is had a few connects into the drift curve, shown in Fig. 2 a.

The 3rd step: original signal and drift curve are subtracted each other, and curve after the formation drift correction is shown in Fig. 2 b.

Step 102, the curve after the drift correction is carried out wavelet transformation, HF noise signal is suppressed.

In this step, can choose the Moret wavelet function, scale parameter is 2, and transformation results is shown in Fig. 2 c.

In this step, the detailed process of carrying out wavelet transformation is same as the prior art, repeats no more here.

Step 103, hang down and cut filtering, the interference of excision signal DC component.

In this step, low to cut filtering parameter be 20MHz, and the result is shown in Fig. 2 d after the filtering.

In this step, hang down the detailed process of cutting filtering same as the prior art, repeat no more here.

Above-mentioned steps 101-103 is actual to be the removal of carrying out various undesired signals, can guarantee to obtain like this reliability of the distribution situation of water in the Quaternary Strata.

Step 104, carry out the rolling time window and ask spectra calculation removing Georadar Data signal after the undesired signal.

Particularly, the spectral amplitude of removing the Georadar Data signal after the undesired signal is non-negative and be even symmetry.Its power spectrum | F (ω) | ²Can give prominence to those | F (ω) |＞1 principal ingredient, suppress those | F (ω) |＜1 submember, so power spectrum can be given prominence to the main frequency composition of signal.

In geological radar image data time window, to limited discrete signal

(wherein N is number of samples, and the present invention is 128) makes first auto-correlation:

$R_{\mathrm{ff}}\left(n\right)=\frac{1}{N}\underset{k=0}{\overset{N-1-n}{\mathrm{\Σ}}}f\left(k\right)f(k+n),$ 0≤n≤N-1

Then to the autocorrelation certificate

Make FFT and calculate, try to achieve discrete power spectrum Only get its positive frequency part

Enveloping curve as the power spectrum curve P of signal in the time window _f(ω).Time window rolls downwards from initial sampling point like this, has just consisted of the rolling time window.

Step 105, carry out low frequency signal energy occupation rate and calculate, obtain moisture attribute data.

The Quaternary Strata water cut is very large on the dynamic characteristic impact of Electromagnetic Wave Propagation, especially power spectrum characteristic.Because water-bearing media electric conductivity increases, so the increasing of medium electro-magnetic wave absorption coefficient, the especially absorption of high-frequency signal, so stratum power spectrum LF-response feature reflects the moisture degree of saturation in stratum.That is to say, low frequency signal energy occupation rate has embodied moisture attribute, can be called moisture attribute data here.

The low frequency signal energy occupation rate in other words computing method of moisture attribute data realizes in the following ways: calculate 1/4th following power spectrum energy of Ground Penetrating Radar Signal antenna dominant frequency and the ratio of the whole power spectrum energy of antenna, its result is shown in Fig. 2 e.

Step 106, moisture attribute data is carried out two-dimentional running mean, eliminate the wherein interference of geology catastrophe point.

Computing method are as follows:

$x_{\mathrm{\Σ}}\left(t\right)=\frac{1}{N2-N1+1+M2-M1+1}\underset{i=N1}{\overset{N2}{\mathrm{\Σ}}}\underset{j=M1}{\overset{M2}{\mathrm{\Σ}}}{x}_{\mathrm{ij}}\left(t\right)$ N1＜N2?M1＜M2

Wherein: N1 is the start channel number, and N2 stops the road number; M1 is initial number of samples; M2 stops number of samples.

Step 107, moisture attribute data is carried out normalized.

Because moisture attribute data is all below 1, in order to make things convenient for image to show, these data are carried out normalized, all normalize in the 0-255 data area, implementation is as follows: the maximal value of finding out all data, all data be multiply by 255, will take advantage of again the result divided by above-mentioned definite maximal value.

Step 108, carry out colored section imaging, the moisture attribute data after the normalized is shown in the bitmap mode.

In actual applications, in order to see intuitively the distribution situation of underground water, the present embodiment has carried out colored section imaging processing.

At first, carry out the color range modulation, construct 8 scale-of-two color tables, have 256 color ranges, the moisture attribute data after the normalized in the corresponding preset range of each color.

Carry out at last colored section imaging, the moisture attribute data after the normalized is shown in the bitmap mode.

In fact, in order to see intuitively the effect of the interference front and back of eliminating the geology catastrophe point.After above-mentioned steps 105, can carry out first step 107 and step 108 one time, to eliminate the colored profile image (as shown in Figure 5) before the interference of geology catastrophe point, in order to compare with colored profile image (as shown in Figure 6) after the interference of eliminating the geology catastrophe point.

By the above embodiments as seen, definite method that this underground water of the present invention distributes not only can obtain the distribution situation of water in the Quaternary Strata, and the data that obtain than the method for the direct explanation of prior art more accurately, reliable.

Claims (13)

1. the definite method that underground water distributes is characterized in that, comprises the steps:

A, the Georadar Data of the Quaternary Strata structure that collects is removed undesired signal process;

B, to having removed the Georadar Data of undesired signal, calculate low frequency signal energy occupation rate, obtain moisture attribute data;

C, moisture attribute data is carried out the section imaging processing, explicitly descend the water distribution situation.

2. definite method as claimed in claim 1 is characterized in that: the described removal undesired signal of steps A is processed and is comprised:

A1, the Georadar Data that collects is carried out offset correction, remove the signal drift noise of instrument self;

A2, the curve after the drift correction is carried out wavelet transformation, HF noise signal is suppressed;

A3, the signal behind the wavelet transformation hanged down cut filtering, the interference of excision signal DC component.

3. definite method as claimed in claim 2, it is characterized in that: the signal drift noise of described steps A 1 described removal instrument self comprises:

The drift parameter of each point in A11, the calculating Georadar Data;

A12, the drift parameter that original signal in the Georadar Data is had a few connect into the drift curve;

A13, original signal and drift curve are subtracted each other, form curve after the drift correction.

4. definite method as claimed in claim 3, it is characterized in that: the computing method of steps A 11 described drift parameters are:

Centered by the preset after being positioned at drift correction in the middle of the curve, set up a window that is formed by 101 sampling points, all original signal data in the window are averaged as the drift parameter of this sampling point.

5. definite method as claimed in claim 2, it is characterized in that: in the described steps A 2, choose the Moret wavelet function, scale parameter is 2, carries out wavelet transformation.

6. definite method as claimed in claim 2 is characterized in that: in the described steps A 3, low to cut filtering parameter be 20MHz.

7. definite method as claimed in claim 2, it is characterized in that: described step B comprises:

B1, carry out the rolling time window and ask spectra calculation removing Georadar Data signal after the undesired signal;

The rolling time window that B2, basis calculate is asked power spectrum, calculates low frequency signal energy occupation rate and calculates, and obtains moisture attribute data.

8. definite method as claimed in claim 7, it is characterized in that: described step B1 comprises:

B11, in geological radar image data time window, to limited discrete signal

Make auto-correlation:

$R_{\mathrm{ff}}\left(n\right)=\frac{1}{N}\underset{k=0}{\overset{N-1-n}{\mathrm{\Σ}}}f\left(k\right)f(k+n),$ 0≤n≤N-1

B12, to the autocorrelation certificate

Make FFT and calculate, try to achieve discrete power spectrum Only get its positive frequency part

Enveloping curve as the power spectrum curve P of signal in the time window _f(ω)

Wherein N is number of samples; Described time window rolls downwards from initial sampling point.

9. definite method as claimed in claim 7, it is characterized in that: the described calculating low frequency signal of step B2 energy occupation rate is: calculate 1/4th following power spectrum energy of Ground Penetrating Radar Signal antenna dominant frequency and the ratio of the whole power spectrum energy of antenna.

10. definite method as claimed in claim 7, it is characterized in that: described step B further comprises B3: moisture attribute data is carried out two-dimentional running mean, eliminate the wherein interference of geology catastrophe point.

11. definite method as claimed in claim 10 is characterized in that: describedly moisture attribute data is carried out two-dimentional running mean be:

$x_{\mathrm{\Σ}}\left(t\right)=\frac{1}{N2-N1+1+M2-M1+1}\underset{i=N1}{\overset{N2}{\mathrm{\Σ}}}\underset{j=M1}{\overset{M2}{\mathrm{\Σ}}}{x}_{\mathrm{ij}}\left(t\right)$ N1＜N2?M1＜M2

Wherein: N1 is the start channel number, and N2 stops the road number; M1 is initial number of samples; M2 stops number of samples.

12. definite method as claimed in claim 10 is characterized in that: described step C comprises:

C1, moisture attribute data is carried out normalized;

C2, carry out colored section imaging, the moisture attribute data after the normalized is shown in the bitmap mode.

13. definite method as claimed in claim 12 is characterized in that: described step C2 comprises:

C21, carry out color range modulation, construct 8 scale-of-two color tables, have 256 color ranges, the moisture attribute data after the interior normalized of the corresponding preset range of each color;

C22, carry out colored section imaging according to the color range of C21, the moisture attribute data after the normalized is shown in the bitmap mode.

Images