US 12,436,312 B2 · Method for measuring flow velocity of multilayer groundwater using distributed optical fiber with point-source active heating

Abstract

The present disclosure provides a method for measuring flow velocity of multilayer groundwater using distributed optical fiber with point-source active heating, comprising: setting and optimizing field test parameters; conducting a background temperature monitoring test; performing a point-source active heating distributed temperature measurement test, including conducting multiple rounds of heating tests; denoising the obtained multi-point source thermal plume attenuation signal data of the groundwater; searching for the peak value of temperature-permeability curve based on an automatic multi-scale peak search algorithm; performing secondary processing on the peak data; estimating the groundwater flow velocity. The method of the present disclosure solves the problem that the existing groundwater flow velocity measurement technology cannot simultaneously measure multiple points and multiple layers, and it can avoid groundwater contamination during the measurement process.

A method for measuring flow velocity of multilayer groundwater using distributed optical fiber with point-source active heating, comprising:

Step S1:
Setting and optimizing field test parameters;
Step S2:
Conducting a background temperature monitoring test, comprising performing multiple rounds of temperature measurement tests;
monitoring background temperature of a borehole along depth profile in each round of temperature measurement tests, with a monitoring time of more than 24 hours;
and performing a quality check on the temperature measurement test data after each round of temperature measurement tests;
Step S3:
Performing a point-source active heating distributed temperature measurement test, comprising performing multiple rounds of heating tests;
in each round of heating tests, controlling a heating device controller to turn on each electric heating unit for a predetermined time, so that the temperature of each groundwater monitoring point rises to a predetermined temperature value, then turning off each electric heating unit;
and stopping monitoring when the temperature of each groundwater monitoring point returns to the background temperature;
and performing a quality check on the temperature measurement test data after each round of heating tests;
Step S4:
Denoising multipoint-source thermal plume attenuation signal data of the groundwater obtained in Step 3;
S4.1 Original signal data trimming and conversion, comprising
S4.1.1 For the obtained multipoint-source thermal plume attenuation signal data of the groundwater, determine a depth range and depth interval of valid data based on a defined apex and end point, and, select an appropriate temperature data recording frequency and determine a time interval for data recording;
S4.1.2 Use the determined depth range to segment the original data, and select and obtain the required temperature data according to the determined time interval;
such that X is a given two-dimensional multi-point source thermal plume attenuation signal data of the groundwater, with a depth range L of [0, M], a time range T of [0, N], then the temperature data X is a two-dimensional array of shape N×M, X={xi,j}, i∈{1,2, . . . , N}, j∈{1,2, . . . , M};
S4.1.3 Perform a transposition operation on trimmed temperature data to obtain a two-dimensional array XT of shape M×N, XT={xi,j}, i∈{1,2, . . . , M}, j∈{1,2, . . . , N};
S4.2 Temperature data signal denoising, comprising
S4.2.1 Determine an appropriate cutoff frequency;
S4.2.2 Use a low-pass filter to perform low-pass filtering on the signal, wherein the low-pass filter retains low-frequency signals, that is, large-scale fluctuations in the data, and filters out high-frequency signals, that is, fast and small-scale fluctuations in the data;
wherein low-pass filtering methods comprise moving average filters, exponential smoothing filters, Butterworth low-pass filters, and Chebyshev filters;
for Butterworth filters, a relationship between amplitude and frequency of a nButterworth low-pass filter is expressed by the following formula:
❘
"\[LeftBracketingBar]"
H
⁡
(
j
⁢
ω
)
❘
"\[RightBracketingBar]"
=
1
1
+
(
ω
ω
c
)
2
⁢
n
,
where |H(jω)| is amplitude of frequency response;
ω is angular frequency;
ωc is cutoff frequency, which is a frequency point at which the filter begins to significantly attenuate the output signal;
n is the order of the filter, the higher the order, the faster the filter attenuates near the cutoff frequency;
Step S5:
Searching for a peak value of temperature-permeability curve based on an automatic multiscale-based peak detection algorithm, comprising
S5.1 Construct a local maximum spectrum:
assuming such that x is a given univariate uniformly sampled signal containing a periodic or quasi-periodic peak;
then first calculate the local maximum scan map;
S5.1.1 Remove linear trend data from the signal x;
S5.1.2 using a moving window method to find the local maximum value of signal x, the window size Wk varies as {Wk=2k|k=1,2, . . . , L}, where k is the k-th window of the signal, L=[N/2]−1, where [N/2] is the upper limit function, taking the smallest integer not less than N/2;
the expression of the local maximum mk,i is as follows:
m
k
,
i
=
{
0
,
x
i
-
1
>
x
i
-
k
-
1
^
x
i
-
1
>
x
i
+
k
-
1
r
+
α
,
otherwise
wherein xj-k-1∧xj-1 is the wedge product of xj-k-1 and xj-1;
such that, if x(t) is the local maximum value in the k window at time t, then mk,i=0;
otherwise mk,i=r+α;
where r is a uniformly distributed random number in a range [0,1], and α is a constant factor (α=1);
S5.1.3 Obtain a calculation result M of mk,i, the matrix M is the local maximum spectrum (LMS), and all elements of the matrix M of shape L×N are within the range of [0,1+r];
M
=
(
m
1
,
1
m
1
,
2
⋯
m
1
,
N
m
2
,
1
m
2
,
2
⋯
m
2
,
N
⋮
⋮
⋱
⋮
m
L
,
1
m
L
,
2
⋯
m
L
,
N
)
=
(
m
k
,
i
)
,
S5.2 Sum the local maximum matrix, comprising
S5.2.1 Sum each column of the matrix M, that is, for each time point t, calculate the sum of its local maximum marks at all K scales;
the global minimum λ=argmin(γk) of γ=[γ1, γ2, . . . , γl] represents the window with the most local maximum values;
it gives a vector of the same length as the time series, in which each element represents the total number of times each time point is marked as non-maximum at all scales;
γ
k
=
∑
i
-
1
N
m
k
,
i
,
for
⁢
k
∈
{
1
,
2
,
…
,
L
}
;
S5.2.2 Use the γ value obtained in the previous step to recalculate the LMS matrix M, delete all elements in mk,i that statisfy k>λ, and obtain a new λ×N matrix Mr={mk,i}, i∈{1,2, . . . , N}, k∈{1,2, . . . , λ};
S5.3 Calculate peak position, comprising
the location of the peak is determined by finding the location where the vector value is 0 in the previous step;
S5.3.1 Calculate the column-wise standard deviation of the matrix Mr;
σ
i
=
1
λ
-
1
⁢
∑
k
=
1
λ
[
(
m
k
,
i
-
1
λ
⁢
∑
k
=
1
λ
m
k
,
i
)
2
]
1
2
,
for
⁢
i
∈
{
1
,
2
,
…
,
N
}
;
S5.3.2 Find all indices i where σi=0;
store these values in a vector P=[P1, P2, . . . , P{circumflex over (N)}], where {circumflex over (N)} is the total number of detected peaks of signal x and P is the index of the detected peaks;
S5.4 Temperature-permeability curve peak data search, comprising
S5.4.1 For all the data obtained by S4, change direction along the time interval t, for the temperature and depth data at a certain time point ti, i∈{1,2, . . . , N}, use an Automatic multiscale-based peak detection (AMPD) algorithm to search for the peak value, and obtain a vector Pi=[P1, P2, . . . , P{circumflex over (M)}]i, where P is the index of the detected peaks, and the value of {circumflex over (M)} is the number of measurement points;
S5.4.2 According to the obtained peak index Pi, multiply the corresponding index by the depth interval defined in S4.1.1 to obtain the depth data ri=[r1, r2, . . . , r{circumflex over (M)}]i corresponding to the temperature peak point at a certain time point ti, and store it in a one-dimensional array;
S5.4.3 Change direction along the time interval, use the AMPD algorithm to all the data, to obtain the three-dimensional data of the depth change of the temperature peak point at different times, let the data be a three-dimensional matrix R[t, p], t is the time, p is the measurement point corresponding to the temperature peak, ri,j is the depth data corresponding to p, i∈{1,2, . . . , N}, j∈{1,2, . . . , {circumflex over (M)}};
R
=
[
r
11
⋯
r
1
⁢
M
^
⋮
⋱
⋮
r
N
⁢
1
⋯
r
N
⁢
M
^
]
=
{
r
i
,
j
}
;
Step S6:
Performing secondary processing on the peak data, comprising
based on the two-dimensional data of temperature variation over time at different depths obtained by S5, perform linear regression outlier analysis on the temperature data at different depths and outlier detection of linear data;
S6.1 Regression model fitting:
use ordinary least squares (OLS) to build a regression model for the temperature variation data at each depth over time;
the OLS calculation process is as follows:
wherein there are data points (x1, y1), (x2, y2), . . . , (xn, yn), where x1 is an independent variable and y1 is a dependent variable;
the linear model is expressed as:
y
=
β
0
+
β
1
⁢
x
;
wherein β0 and β1 are model parameters;
the goal of ordinary least squares is to minimize the sum of squared errors, that is:
S
=
∑
n
=
1
n
(
y
i
-
(
β
0
+
β
1
⁢
x
i
)
)
2
;
by taking partial derivatives with respect to β0 and β1 and setting them to 0, the optimal estimates of the two parameters are obtained;
S6.2 Calculate the studentized residuals for the established regression model:
r
i
=
e
i
s
e
⁢
1
-
h
ii
;
where:
ei is the residual of the i-th observation, that is, ei=yi−ŷi;
se is the estimated standard deviation of the residual,
∑
n
=
1
n
e
i
n
-
p
,
where n is the sample size and p is the number of parameters in the model;
hii is the diagonal element of the hat matrix, reflecting the degree of influence of the i-th observation on its own predicted value;
S6.3 Outlier Identification:
Values greater than 2 or 3 in the studentized residual calculation results are considered as outliers and the index of the data is recorded;
S6.4 Outlier processing:
According to the index of the outliers identified in the previous step, assign a null value to the value corresponding to the index;
then use a linear interpolation formula to process the null value, the linear interpolation formula is expressed as:
y
=
y
0
+
(
y
1
-
y
0
)
(
x
1
-
x
0
)
·
(
x
-
x
0
)
;
where:
x1, x0 are the x-coordinates of the known points, y1, y0 are the y-coordinates of the known points, and x is the x-coordinate of the corresponding y-value you want to estimate;
after secondary processing, the depth variation data of the temperature peak points at different measurement points over time are obtained, a three-dimensional matrix {circumflex over (R)}[t, p], t is the time, p is the measurement point corresponding to the temperature peak, and {circumflex over (r)} is the depth data corresponding to p, t∈{1,2, . . . , i, . . . , N}, p∈{1,2, . . . , j, . . . , {circumflex over (M)}};
R
^
=
[
r
^
11
⋯
r
^
1
⁢
M
^
⋮
⋱
⋮
r
^
N
⁢
1
⋯
r
^
N
⁢
M
^
]
=
{
r
^
i
,
j
}
;
Step S7:
Estimating the groundwater flow velocity, using the ordinary least squares method to establish a one-dimensional linear regression model for the temperature variation data at each depth over time, and obtaining the estimated value of the groundwater flow velocity at the depth from the fitted regression model equation;
for the matrix {circumflex over (R)}[t, p] obtained by S6.4, a one-dimensional linear regression model is established using the least squares method (OLS) for the two-dimensional data pj[t, {circumflex over (r)}] at a certain measurement point pj;
wherein, for data points (t1, {circumflex over (r)}1), (t2, {circumflex over (r)}2), . . . , (tn, {circumflex over (r)}n), where t is the independent variable and {circumflex over (r)} is the dependent variable;
the linear model is expressed as:
r
^
=
β
0
+
β
1
⁢
t
;
where β0 and β1 are model parameters;
the goal of ordinary least squares is to minimize the sum of squared errors, that is:
S
=
∑
n
=
1
n
(
r
^
i
-
(
β
0
+
β
1
⁢
t
i
)
)
2
;
by taking partial derivatives with respect to β0 and β1 and setting them to 0, the optimal estimates of the two parameters are obtained;
wherein β1 is the flow velocity at pj, wherein, if its value is negative, then the flow velocity is upward, and if its value is positive, then the flow velocity is downward;
Cycle j to obtain the vertical velocity distribution of groundwater at different measuring points;
wherein step S1 comprises:
analyzing physical parameters of the borehole and core, comprising lithology distribution, thermal conductivity, porosity, specific heat capacity and permeability;
setting and optimizing parameters of a distributed optical fiber temperature measurement host;
wherein the step of setting and optimizing parameters of a distributed optical fiber temperature measurement host comprises:
setting a monitoring length and heating point position of a temperature measuring optical fiber, and calibrating optical fiber refractive index parameters through a first water bath temperature calibration unit and a second water bath temperature calibration unit:
testing a monitoring temperature resolution and temperature accuracy of the distributed optical fiber temperature measurement host, and calibrating and optimizing optical fiber attenuation coefficient, temperature sensitivity coefficient and temperature compensation value of the distributed optical fiber temperature measurement host through the first water bath temperature calibration unit and the second water bath temperature calibration unit;
performing test in borehole groundwater, and setting maximum temperature, duration, and power of the heating based on expected temperature response and safety criteria.

480 more elements

Method for measuring flow velocity of multilayer groundwater using distributed optical fiber with point-source active heating

Life of this patent

Method for measuring flow velocity of multilayer groundwater using distributed optical fiber with point-source active heating

Life of this patent

Abstract

Claims7