METHOD FOR DETERMINING PHYSICAL CHARACTERISTICS OF A HOMOGENEOUS MEDIUM AND ITS BOUNDARIES

Abstract

The harmonic wave, which is oscillation of a physical value along one direction of propagation in a homogeneous medium, is recorded by means of sensors along the direction of propagation of the oscillation at least at five points equally spaced from each other. The output signals of the sensors are converted into the corresponding complex spectral amplitudes corresponding to the frequency decomposition of the output signals. A model of harmonic wave propagation in a homogeneous medium is created, in which for any oscillation frequency the wave is represented as descending and ascending exponentially decaying harmonic waves propagating in opposite directions. The absolute values of the complex spectral amplitudes of the output signals of the sensors at each frequency are used as input data for equations comparing the absolute values of the complex amplitudes with the created model of wave propagation. By solving the obtained equations, the total complex amplitudes of the descending and ascending waves and the complex propagation constant of oscillations at each frequency are determined and the characteristics of the boundaries of the homogeneous medium are determined basing on the ratio of the complex amplitudes of the descending and ascending waves, and the characteristics of the homogeneous medium are determined basing on the phase velocity and attenuation coefficient of the wave.

Claims

1. A method for determining physical characteristics of a homogeneous medium and its boundaries, comprising: recording a harmonic wave propagating in a homogeneous medium and representing an oscillation of a physical value along one direction of propagation in the homogeneous medium, wherein recording comprises recording said physical value by means of sensors along the direction of propagation of the oscillation at least five at points equally spaced from each other, converting the output signals of the sensors, with the help of a computer system, using a spectral analysis method, into the corresponding complex spectral amplitudes corresponding to the frequency decomposition of the output signals, creating a model of harmonic wave propagation in the homogeneous medium, in which for any oscillation frequency the wave is represented as descending and ascending exponentially decaying harmonic waves propagating in opposite directions, wherein the model depends on the complex amplitudes of the descending and ascending waves and the complex constant of oscillation propagation, using, with the help of the computer system, the absolute values of the complex spectral amplitudes of the output signals of the sensors at each frequency as input data for equations comparing the absolute values of the complex amplitudes with the created model of wave propagation, determining, by solving the obtained equations, the total complex amplitudes of the descending and ascending waves and the complex propagation constant of oscillations at each frequency, and determining the physical characteristics of the boundaries of the homogeneous medium basing on the ratio of the complex amplitudes of the descending and ascending waves, and the physical characteristics of the homogeneous medium basing on the phase velocity and attenuation coefficient of the wave.

2. Method of claim 1, comprising creating fluctuations of the physical value along the direction of propagation in the medium by artificial means.

3. Method of claim 2, wherein the homogeneous medium is a segment of a fluid-filled well drilled in the formation, the fluctuations of the physical value are produced by a pump connected to the well or arranged inside it, the physical characteristics of the homogeneous medium are compressibility and viscosity of the fluid filling the well, and the physical characteristics of the boundaries of the homogeneous medium are the coefficient of productivity of the formation adjacent to this segment.

4. Method of claim 2, wherein the homogeneous medium is a carrier of electromagnetic waves, and the oscillations of the physical value are produced by an emitter of electromagnetic oscillations.

5. Method of claim 1, wherein the conversion of the output signals of the sensors is carried out using a discrete Fourier transformation.

6. Method of claim 1, wherein the recording of said physical value at least at five points is carried out simultaneously by means of sensors, each being installed at a corresponding point.

7. Method of claim 1, wherein the recording of said physical value at least at five points is carried out sequentially by successively moving at least one sensor in the direction of propagation of fluctuation of the physical value.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0013] The invention is illustrated by drawings, where

[0014] FIG. 1 is a graph of Pol.sub.8(y) polynomial;

[0015] FIG. 2 is a graph of Pol.sub.8(y) polynomial in the presence of an additional sensor;

[0016] FIG. 3 shows the general structure of the solutions;

[0017] FIG. 4 shows an example implementation of the method.

DETAILED DESCRIPTION OF THE INVENTION

[0018] The following description discloses the essence of the present invention, which permits determination of all the parameters of a harmonic wave, relying only on the absolute values of the complex amplitudes.

[0019] For a group of 2N+1 equidistant sensors (N=0, 1, 2, . . . ) we have, for each frequency v:

p.sub.lp(z.sub.l)=e.sup.2iv.sup.l(AG.sup.l+BG.sup.l)

l=N,N+1, . . . ,1,0,1, . . . N 1,N(1)

where A, B are the complex amplitudes of the descending and ascending waves, G is the complex number, the antenna transfer function, p.sub.lthe values of the measured value at the point number l, z.sub.l is the coordinate of the point number l. Here, unknown time delays .sub.l between the starting time of the sensors are also introduced. These equations are valid for describing wave propagation in any time-independent homogeneous medium.

[0020] By calculating the absolute values of the complex values p.sub.i, the phase dependence is eliminated:

|p.sub.l|.sup.2=|AG.sup.l+BG.sup.l|.sup.2=(a.sup.2g.sup.2l+b.sup.2g.sup.2l+2ab cos(2l))(2)

where:

Ge.sup.zge.sup.i, g=e.sup.Re()z, =m()z(3)

a=|A|, b=|B|, =arg(A)arg(B)

[0021] and zdistance between adjacent sensors.

[0022] With N=2 there are 5 real equations for 5 real values: a, b, g, . Algebraic transformations make it possible to reduce these equations to a single polynomial equation of the 8th order for

[0023] y=cos(2). The latter is solved numerically.

[0024] The equations have three discrete symmetries:

.fwdarw.+(3.1)

.fwdarw.,.fwdarw.(3.2)

a.fwdarw.b, b.fwdarw.a, g.fwdarw.1/g,(3.3)

[0025] Thus, in all cases, at least 4 solutions are obtained for y from one for y=cos(2):

=.sub.0, =.sub.0, .sub.0=a cos(y)(4)

[0026] Let's describe the procedure for constructing a solution. Determine

[00001]$\begin{array}{cc}{D}_0{.Math.p_0.Math.}^2& \left(4.1\right)\\ D_1{.Math.p_1.Math.}^2+{.Math.p_{-1}.Math.}^2& \left(4.2\right)\\ D_2{.Math.p_2.Math.}^2+{.Math.p_{-2}.Math.}^2& \left(4.3\right)\\ Q_1{.Math.p_1.Math.}^2-{.Math.p_{-1}.Math.}^2& \left(4.4\right)\\ Q_2{.Math.p_2.Math.}^2-{.Math.p_{-2}.Math.}^2& \left(4.5\right)\\ M_1{D}_2+2.Math..Math.D_0-2.Math..Math.D_1.Math.y,M_2{D}_1-2.Math..Math.D_0.Math.y,w\left(y\right)=\frac{M_1}{M_2},& \left(5\right)\end{array}$

[0027] Then the equation of the 8th order for y is as follows:

P.sub.8(y)=0

P.sub.8(y)=(Q.sub.1M.sub.1Q.sub.2M.sub.2).sup.2(M.sub.1.sup.24M.sub.2.sup.2)4(1y.sup.2)[D.sub.0(2M.sub.2.sup.2D.sub.0(M.sub.12yM.sub.2))(M.sub.12yM.sub.2)(M.sub.1.sup.24M.sub.2.sup.2)(Q.sub.22Q.sub.1y).sup.2M.sub.2.sup.4](6)

[0028] Among its roots, one must choose those that meet the conditions

m(y)=0, 1y1, w(y)2(7)

[0029] Since the argument of the transfer function of the antenna is restored only up to discrete symmetries (3.1-2), we obtain four from one y. Additional considerations are required to eliminate discrete uncertainty, for example, if the estimated phase velocity c, is known, one can write

[00002]$\begin{array}{cc}\frac{2.Math..Math..Math..Math.v.Math..Math..Math..Math.z}{c}& \left(8\right)\end{array}$

and choose the one closest to (8) from the four .

[0030] Other quantities are expressed via y as follows.

Calculate g as

[0031] [00003]$\begin{array}{cc}{g}_{}^2=\frac{1}{2}.Math.\left(w\sqrt{w^2-4}\right),g_{+}.Math.g_{-}=1& \left(9\right)\end{array}$

Then, a and b are found to be

[00004]$\begin{array}{cc}{a}_{}^2=\frac{1}{2}.Math.\left(u{v}_{+}\right)& \left(10\right)\\ b_{}^2=\frac{1}{2}.Math.\left(u{v}_{+}\right)& \left(11\right)\\ v_{}=\frac{1}{\sqrt{w^2-4}}.Math.\frac{Q_2-2.Math..Math.Q_1.Math.y}{w-2.Math..Math.y}& \left(12\right)\end{array}$

[0032] The last value, cos(), is equal to

[00005]$\begin{array}{cc}\cos \left(\right)=\frac{D_0-u}{2.Math..Math.\mathrm{ab}}& \left(13\right)\end{array}$

and it is not sensitive to the uncertainty of the sign. However, there is another simple uncertainty to obtain ,

.Math.(14)

[0033] There are 8 different solutions for G, since the symmetries with respect to g, which transform g into g1, exist along with four from one y. There are 4 different solutions for A and B due to the symmetry, along with the solution .Math.. With the symmetry the reflection coefficient

[00006]$\begin{array}{cc}R\frac{B}{A}=\frac{b}{a}.Math.e^{-i.Math..Math.}& \left(15\right)\end{array}$

is changed as follows

[00007]$\begin{array}{cc}a.Math.b.Math.R\frac{B}{A}=\frac{b}{a}.Math.e^{-i.Math..Math.}.Math.\frac{a}{b}.Math.e^{-i.Math..Math.}=\frac{1}{\stackrel{\_}{R}}& \left(16\right)\end{array}$

With the .Math. symmetry, we obtain

R.Math.R(17)

[0034] Thus, it is possible to determine R with a 4-fold discrete uncertainty. The true value of R can be found with the help of additional considerations, such as calculations by formula (8).

[0035] Usually there is more than one root of the equation of the 8th order, which meet the necessary constraints (7):

m(y)=0, 1y1, w(y)2(18)

[0036] FIG. 1 is a graph of the Pol.sub.8(y) polynomial for A=1.5, B=1+0.1i, G=1.0462+0.33992i. Crosses indicate solutions Pol.sub.8(y)=0. Solid circles indicate a subset of solutions that satisfy the condition w(y)2. The true solution is shown in the square.

[0037] Addition of one more sensor makes it possible to eliminate the above uncertainty and obtain the true root for y. Let's demonstrate this in the following example. FIG. 2 is a graph of the Pol.sub.8(y) polynomial for the same input data as for FIG. 1. Crosses indicate solutions Pol.sub.8(y)=0. Solid circles indicate a subset of solutions that satisfy the condition w(y)2. The true solution is shown in the square. The triangle denotes the solution remaining when taking into account the data of the additional sixth sensor.

[0038] The general structure of the solutions is shown in FIG. 3. FIG. 3a shows the solutions for G, and FIG. 3b shows the solutions for A, B. 3 data sets are shown: 1) G 2) A, B. For each set, the true solution is shown, along with the solutions for methods based on the use of 5 sensors, 5s, and 6 sensors, 6s. The graph on the left: G. A cross is a true value, circles are: solutions for 5 sensors, triangles: solutions for 6 sensors. The graph on the right: A, B. An arrow+line on the right: true A, an arrow+line on the right: true B, filled circles: A from the solutions for 5 sensors, empty circles: B from the solutions for 5 sensors, triangles pointing down: A from the solutions for 6 sensors, triangles pointing up: B from solutions for 6 sensors.

[0039] Solutions for G are divided into 8-fold sets, each corresponding to the solution P8(y)=0. In all cases there is one solution for method 5s or 6s, which coincides with the true solution, for all values A, B, G. However, all 8 solutions for method 6s will always correspond to the true value y=cos (2), therefore method 6s is preferred. The uncertainty of the total time shift associated with a shift in the reference time of the true records by a constant value is used to set

m(A)=0.

[0040] As a practical application of the proposed method, we consider an example given in FIG. 4, where 1 is a pump, 2 is a pipe or pipe system connecting the pump and a well, 3 is a well, 4 is the earth surface, 5 is a harmonic pressure wave profile in the well, at a fixed point in time and on one of the dominant frequencies, 6 is a formation, 7 is a profile of a harmonic pressure wave in the formation, at a fixed point in time and on one of the dominant frequencies, 8 is a system of six equidistant sensors. This example in no way limits the application of the method and is given by way of illustration.

[0041] Suppose there is a well 3, drilled in the ground, and filled with fluid, for example, water, or oil. The well 3 may be in communication with at least one permeable formation 6 that intersects it, for example, an oil-bearing formation. Suppose there is a pump 1, installed on surface 4, and either pumping fluid into the formation 6 through the pipe 2, or pumping fluid from the formation 6 through the well 3. As a rule, the pump, regardless of the specifics of its technical implementation, produces periodic fluctuations in the pressure in the fluid filling the well, in addition to the main quasistationary change of pressure in the well. Thus, the well 3 is filled with pressure waves 7, propagating up and down in it. The well is used to place the system 8 of six equidistant pressure sensors, such as high-speed pressure gauges or hydrophones, and periodic pressure fluctuations in the well are recorded that correspond to the pump operation. Sensors can record data in memory, and can transfer them to the surface immediately after recording, if they are connected to a suitable data transmission system, such as, for example, a geophysical cable. In the first case, data analysis involves removing sensors to the surface and uploading data to a computer, in the second case data can be analyzed without removing sensors from the well. Thus, the data set p.sub.l(t) is obtained, where t is time, and l=1, 2, . . . , N is the sensor number. The data obtained are analyzed by a computer program as follows. A discrete Fourier transformation is performed and the complex amplitudes are obtained P.sub.l(v.sub.i), where v.sub.i is a set of frequencies. The dominant frequencies v.sub.iD are determined, for which the modules of complex amplitudes are maximal, equations (1) are solved for each of v.sub.iD, where P.sub.l(v.sub.iD) are used as input data, and complex amplitudes A.sub.iD, B.sub.iD and the transfer function of the antenna G.sub.iD are determined. The G.sub.iD contains information about the phase velocity and attenuation coefficient of pressure waves at dominant frequencies and, thus, it can be used to determine the rheology of the fluid, in particular its compressibility and viscosity, while the ratio A.sub.iD/B.sub.iD can be used to determine the coefficient of productivity of a permeable formation.

[0042] In another modification of the method of the present invention, the same situation is considered as described above, however, instead of a set of six sensors, a single sensor is used, which sequentially records pressure at six equally distant depths. The obtained data p.sub.l(t), l=1 . . . 6 is then used in the same way as in the method described above, while the time shift between the measurements made by one sensor at different depths is insignificant due to the mathematical structure of the method. The only limitation in this case is the condition of stability of the pump operation and the stationarity of all main parameters of the well and the formation for the entire time of data collection.

[0043] It is implied that a model of wave propagation in said medium was previously created, in which for each oscillation frequency a wave is represented as descending and ascending exponentially decaying harmonic waves propagating in opposite directions, and the model depends on the complex amplitudes of the descending and ascending waves and on the complex constant of propagation of fluctuations. Since the model depends on a set of geometrical and physical parameters of the medium, and the values characterizing the wave propagation are functions of these parameters and frequencies, it is possible to choose the parameters of the medium so that these values coincide with the measured ones. An example is the propagation of a weak pressure pulse in a rigid pipe filled with fluid of a constant diameter, in this case the physical parameters defining the phase velocity and attenuation coefficient at a certain frequency are density, bulk modulus, viscosity, and diameter of the pipe, while an example of parameters affecting the reflection coefficient is the ratio of diameters of the adjacent pipe sections.

[0044] As it is indicated above, the obtained results can be used to determine the physical characteristics of the boundaries of a homogeneous medium basing on the interpretation of the ratio of the complex amplitudes of the descending and ascending waves (reflection coefficient), as well as to determine the physical characteristics of the medium itself basing on the interpretation of the components of the complex propagation constant, namely phase velocity and wave attenuation coefficient. Here is an example.

[0045] In the case of waves propagating in a hard pipe of radius r filled with a viscous fluid having density p with the phase velocity of bulk waves c and viscosity , the following relations are satisfied in the low-frequency approximation

[00008]$=i.Math.\frac{}{c}+\frac{4.Math..Math.}{.Math..Math.c.Math..Math.r^2}$

Where i is the imaginary unit, is the circular frequency, and the reflection coefficient in the zone of the joint of pipes of radii r.sub.1 and r.sub.2

[00009]$R=\frac{r_2^2.Math..Math.r_1^2}{r_2^2+r_1^2}$

Thus, knowing r.sub.1 and measuring y and R by the method described above, one can determine

[00010]$c,\frac{}{},\frac{r_2}{r_1}.$