CN102866429B - 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, definite method that particularly a kind of underground water distributes.

Background technology

Usually before and after coal mining, all will carry out geology detecting to stratal configuration, wherein a main project of Tan Ceing is exactly to determine phreatic distribution situation.

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, the detection accuracy height, 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 Processing is to adopt the method for directly explaining at present, just by after the source book of geologic radar detection is made simple process such as some amplifications, filtering, stack, external appearance characteristic according to radar appearance, as reflect power, phase characteristic, Zhou variation characteristic or the like information is directly made qualitative and quantitative interpretation to reflected signal in the same way.

Because coal mining can cause ground subsidence, has destroyed the Quaternary system stratal configuration, causes the loss and the variation of surface water, if adopt direct interpretation procedure, can't obtain the distribution situation of water in the Quaternary system stratum, promptly 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 this method can obtain the distribution situation of water in the Quaternary system stratum.

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

A, the geological radar data of the Quaternary system stratal configuration that collects are removed undesired signal handle;

B, to having removed the geological radar 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, show underground water distribution situation.

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

A1, the geological radar data that collect are 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 geological radar data;

A12, the drift parameter that original signal in the geological radar 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: a preset in the middle of the curve is the center after the drift correction to be positioned at, set up a window of forming 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 power spectrum to calculate removing geological radar 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 a 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, the moisture attribute data after the normalized is shown in the bitmap mode according to the color range of C21.

As seen from the above technical solutions, definite method that this underground water of the present invention distributes is handled owing to earlier the Quaternary system stratal configuration geological radar data that collect are removed undesired signal; To having removed the geological radar 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, show underground water distribution situation.Therefore, not only can obtain the distribution situation of water in the Quaternary system stratum, 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 be embodiment illustrated in fig. 1 in geologic radar detection raw data curve 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 cuts filtered data and curves figure for the curve of Fig. 2 c 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 the geologic radar detection of Fig. 2 a correspondence;

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 for Fig. 5 through the water percentage profile image after the two-dimentional running mean.

Embodiment

Followingly develop simultaneously with reference to accompanying drawing that the present invention is described in detail for specific embodiment.

The invention provides definite method that a kind of underground water distributes, this method can obtain the distribution situation of water in the Quaternary system stratum.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 system stratal configuration geological radar data that collect are 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 geological radar data.

Particularly, be example with the preset N point in the middle of the original signal of Fig. 2 a, be the center with the N point, set up a window W who forms by 101 sampling points, 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, it is same as the prior art to hang down the detailed process of cutting filtering, repeats no more here.

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

Step 104, carry out the rolling time window and ask power spectrum to calculate removing geological radar data-signal after the undesired signal.

Particularly, the spectral amplitude of removing the geological radar 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 a number of samples, and the present invention is 128) makes auto-correlation earlier:

$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 constituted the rolling time window.

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

Quaternary system stratum water cut is very big to the dynamic characteristic influence 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 that low frequency signal energy occupation rate has embodied moisture attribute, can be called moisture attribute data here.

The low frequency signal energy occupation rate computing method of moisture attribute data in other words 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, image shows for convenience, these data are carried out normalized, all normalize in the 0-255 data area, concrete enforcement is as follows: the maximal value of finding out all data, all data be multiply by 255, will take advantage of the result again 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 phreatic distribution situation intuitively, 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 colored section imaging at last, the moisture attribute data after the normalized is shown in the bitmap mode.

In fact, in order to see the effect of the interference front and back of eliminating the geology catastrophe point intuitively.After above-mentioned steps 105, can carry out step 107 and step 108 earlier one time, to eliminate the colored profile image (as shown in Figure 5) before the interference of geology catastrophe point, so that 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 system stratum, and the data that obtain than the method for the direct explanation of prior art more accurately, reliable.

Claims (11)

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

A, the geological radar data of the Quaternary system stratal configuration that collects are removed undesired signal handle;

B, to having removed the geological radar 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, show underground water distribution situation;

Wherein, described step B comprises:

B1, carry out the rolling time window and ask power spectrum to calculate removing geological radar 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;

Wherein, 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 a number of samples; Described time window rolls downwards from initial sampling point.

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

A1, the geological radar data that collect are carried out drift 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 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 geological radar data;

A12, the drift parameter that original signal in the geological radar 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 is characterized in that: the computing method of steps A 11 described drift parameters are:

A preset in the middle of the curve is the center after the drift correction to be positioned at, and sets up a window of being made up of 101 sampling points, and 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 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 1 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.

8. definite method as claimed in claim 1 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.

9. definite method as claimed in claim 8 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＜N2M1＜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.

10. definite method as claimed in claim 8 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.

11. definite method as claimed in claim 10 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, the moisture attribute data after the normalized is shown in the bitmap mode according to the color range of C21.

Images