Unformatted text preview: SPE 121393
A New Method for Continual Forecasting of Interwell Connectivity in
Waterfloods Using an Extended Kalman Filter
Daoyuan Zhai, Jerry M. Mendel, Feilong Liu, University of Southern California Copyright 2009, Society of Petroleum Engineers
This paper was prepared for presentation at the 2009 SPE Western Regional Meeting held in San Jose, California, USA, 24–26 March 2009.
This paper was selected for presentation by an SPE program committee following review of information contained in an abstract submitted by the author(s). Contents of the paper have not been
reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material does not necessarily reflect any position of the Society of Petroleum Engineers, its
officers, or members. Electronic reproduction, distribution, or storage of any part of this paper without the written consent of the Society of Petroleum Engineers is prohibited. Permission to
reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of SPE copyright. Abstract
This paper is based on a relatively simple parametric
model that characterizes the system function between a
specific producer and each of its contributing injectors.
The model has only two parameters for each producerinjector pair; so, if N injectors are assumed to contribute
to a producer, there will be 2N unknown parameters. An
adaptive strategy, using an Extended Kalman Filter
(EKF), is used to estimate the 2N parameters, which are
then used to generate N numeric InjectorProducerRelationship (IPR) values for the N producerinjector
pairs. The IPR values allow one to assess how well an
injector influences the producer.
This same model and an EKF were first used in Liu, et
al [5]. The modified EKF used in this paper avoids
problems that can arise when processing real data and
provides additional information that is useful for future
research. Our modified EKF is applied to real data from a
section of an oil field. A validation strategy for the
estimated IPR values is developed in terms of “prediction
errors.” A strategy is also presented for choosing an
optimal set of injectors that affect a producer. Finally, a
simple method is presented for converting producercentric IPR values to injectorcentric IPR values so that
reservoir engineers can easily see which producers are
being affected by a specific injector.
I. Introduction
One of the most important issues in waterflood
management is the ability to detect preferential subsurface
flow trends, i.e., where the water flows. Recently, one
approach for doing this was to estimate the interwell
relationship of each injector and its contributing
producers using only the production and injection rates.
More specifically, the interwell relationship, referred
here as the “InjectorProducerRelationship (IPR)”, was
evaluated by estimating parameters of a model that
characterizes the reservoir. After obtaining the IPR
values, one can ascertain whether an injector is
contributing to a producer and, if so, by how much. In addition, one can often infer other features about the
reservoir from the IPR values, e.g., the directional sweep
efficiency of a given pattern and presence of directional
fractures.
[1], [2], [4], and [7][11] introduce many different
methods, yet all aim to achieve the above goal. The most
recent works are Albertoni and Lake [1] and Yousef, et al.
[11], both of which model the reservoir as a continuous
impulse response that converts input signals (injection
rates) to output signals (production rates). Liu, et al [5]
have pointed out the following difficulties that [1] and
[11] have encountered: (a) the parameters in the model
and the IPRs are assumed stationary over the window of
measurements for which processing occurs, and so when
the IPRs change the processing needs to be redone for the
new situation, and this may not be practical because the
reservoir is dynamic and it may be very difficult to
recognize when a change has occurred; and, (b) the CM
model is somewhat complex, and although it is
characterized by two parameters, to use it, the primary
bottom hole pressure impact also needs to be determined.
Additionally, Albertoni and Lake [1] model the
production rate as a weighted sum of delayed injection
rates from the injectors assumed to contribute to a
producer; but, in general, how many delayed injection
rates to include in the model is unknown.
To overcome the above difficulties, Liu, et al [5]
proposed to use a relatively simple twoparameter autoregressive model to characterize the impulse response
between a single injector and a single producer. The area
under the impulse response is used as a measure of the
IPR. If N injectors are assumed to influence a producer,
this producercentric model contains exactly 2N
parameters. This impulse response model is then used to
establish an equivalent State Variable Model (SVM),
which is then used by an Extended Kalman Filter (EKF)
to forecast and track the IPR values. This paper continues
to use the same impulse response model as in [5];
however, we modify the SVM in some important ways so
that the (modified) EKF can provide some new 2 SPE 121393 information, and so that some practical difficulties can be
overcome when the modified EKF is used to process real
data.
The rest of this paper is organized as follows: Section II
provides a brief recapitulation of the SVM that is in [5];
Section III describes our modified SVM; Section IV
describes the modified EKF that is associated with the
modified SVM; Section V provides a strategy for using
our EKF to process data for an entire section; Section VI
sets forth a validation strategy for our modified EKF;
Section VII describes a procedure for choosing the initial
set of candidate injectors that may (or may not) influence
a specific injector; Section VIII explains how to convert
producercentric IPR values into injectorcentric IPR
values; and, Section IX includes our conclusions as well
as some suggested topics for further research. A. Reservoir Model
Our reservoir is modeled as a collection of continuoustime impulse responses that convert injection rates into a
production rate. A producercentric reservoir model with
N independent injectors and one producer is depicted in
Fig.1, where i1 (t ), i2 (t ),..., iN (t ) , n1 (t ), n2 (t ),..., nN (t ) and im,1 (t ), im,2 (t ), ..., im, N (t ) are the actual injection rates that
flow into the reservoir, the corresponding injection rate
measurement noise, and the measured injection rates,
respectively; p(t ) , n p (t ) and pm (t ) are the actual
production rate, the corresponding production rate
measurement noise, and the measured production rate,
respectively; p cj (t ) represents the amount of production
p(t ) (2) where α = e − aT and γ = bα T are two parameters that
determine the model and T is the sampling period. The
total production rate at a specific producer is the sum of
the individual production rates contributed by injectors
that influence this producer. Assume N injectors influence
a producer; then, for this producer:
N N j =1 j =1 P ( z ) = ∑ Pjc ( z ) = ∑ H j ( z ) f (rj , k j ) I j ( z ) (3) γ j f (rj , k j ) z −1
=∑
I j ( z)
−1 2
j =1 (1 − α j z )
N where P ( z ) , Pjc ( z ) and I j ( z ) are the Z transforms of the II. Liu and Mendel’s Original SVM rate in γ z −1
(1 − α z −1 )2 H ( z) = caused by the jth injector ; and, f (rj , k j ) ( j = 1,..., N ) are scale functions which show
how much of each injection rate flows in the direction
toward the producer, where f (rj , k j ) is viewed as a
linear or nonlinear scalar function of the distance, rj , and
the permeability, k j , between the jth injector and the
producer. Because noisefree data, i1 (t ), i2 (t ),..., iN (t ) and
p(t ) , are not directly available, we use their measured total production rate, the production rate contributed by
the jth injector, and the injection rate of the jth injector,
respectively.
From (2), the model of the subsystem between the
producer and the jth injector can be expressed as:
(1 − α j z −1 ) 2 Pjc ( z ) = γ j f (rj , k j ) z −1 I j ( z ) (4) which, when transformed back into time domain, is:
p c (k + 1) − 2α j p c (k ) + α j 2 p c (k − 1) = γ j f (rj , k j )i j (k ) (5)
j
j
j Liu, et al. [5] defined the IPR as the area under the
sampled impulse response for each injector, and showed
that it can be calculated as: IPR j = γ j f ( rj , k j )
(1 − α j ) 2 (6) Observe that three parameters— α j , γ j and f (rj , k j ) —
need to be estimated to compute this IPR j ; however,
because γ j and f (rj , k j ) are multiplicative in (5) and
(6), the number of parameters can be reduced from three
to two by letting γ ′ ≡ γ j f (rj , k j ) . Note, also, that only
j values, im,1 (t ), im,2 (t ),..., im, N (t ) and pm (t ) , for our data the measured injection rates, im, j (k ) , are available for processing.
B. InjectorProducer Model for Each Subsystem processing; hence, (5) and (6) can be reexpressed, as: pc (k + 1) − 2α j pc (k ) + α j 2 pc (k − 1)
j
j
j After comparing different models and getting advice
from petroleum engineers, the following two parameter
autoregressive (AR) model was chosen in [5] to represent
the impulse response between a single producer and a
single injector:
h(t ) = bte − at = γ ′ (im, j (k ) − n j (k )) = γ ′im, j (k ) + n pc (k )
j
j
IPR j = (1) Because only sampled injection and production rates are
available for processing, a discretetime version of (1) is
used; its Ztransform is: (7) j γ ′j
(1 − α j ) 2 (8) Note that in (7), we set n pc (k ) = γ ′j n j (k ) so that we will
j be able to use im, j (k ) instead of i j (k ) as our input and
n pc ( k ) is the standard additive state equation noise.
j SPE 121393 3 ⎧
x11 (k )
⎤
⎡
⎪
⎥
⎢
x12 (k )
⎪
⎥
⎢
⎪
⎥
⎢
x14 (k )
⎪
⎥
⎢
2
2 x11 (k ) x14 (k ) − x11 (k ) x13 (k ) + x12 (k )im ,1 (k ) ⎥
⎪
⎢
⎪
⎥
⎢
x21 (k )
⎪
⎥
⎢
x22 (k )
⎪
⎥
⎢
⎪ x(k + 1) = ⎢
⎥ + n (k )
x24 (k )
⎪
⎥ x
⎢
2
⎪
⎢ 2 x21 (k ) x24 (k ) − x21 (k ) x23 (k ) + x22 (k )im ,2 ( k ) ⎥
⎪
⎥
⎢
⎨
M
⎥
⎢
⎪
⎢
⎥
xN 1 ( k )
⎪
⎢
⎥
⎪
xN 2 ( k )
⎢
⎥
⎪
⎢
⎥
⎪
xN 4 ( k )
⎢
⎥
⎪
2
⎢ 2 xN 1 (k ) xN 4 (k ) − xN 1 (k ) xN 3 (k ) + xN 2 ( k )im , N ( k ) ⎦
⎥
⎣
⎪
⎪
⎪
c
c
⎪ p(k + 1) = p1c (k + 1) + p2 (k + 1) + L + pN (k + 1) + n p (k + 1)
⎪
= [ 0 0 0 1 0 0 0 1 L 0 0 0 1] x(k + 1) + n p (k + 1)
⎪
⎩ C. StateVariable Model
To use an EKF, one must first construct an SVM [6].
To establish the SVM, one starts with the secondorder
finitedifference equation (7) that is described by two
state variables p c (k − 1) and p cj ( k ) . In addition, the
j unknown parameters α j and γ ′ are also treated as state
j
variables so that they can be estimated by an EKF; hence,
(7) is actually described by the following 4 × 1 state
vector:
x j (k ) = ⎡ x j1 (k ), x j 2 ( k ), x j 3 ( k ), x j 4 (k ) ⎤ '
⎣
⎦
= ⎡α j ( k ), γ ′j (k ), p c (k − 1), p c (k ) ⎤ '
j
j
⎣
⎦ (9) Aggregating the state vectors of all N injectors, one
obtains a complete 4N × 1 state vector, as:
x(k ) = [ x1 (k ) ', x2 (k ) ',L , xN (k ) '] ' III. New Modified SVM
The SVM in (11) has been modified by us in two
important ways:
(1) IPR is treated as a state variable, and is estimated
directly. Its error (pseudo) variance can be used to
provide upper and lower bounds for the IPR estimates.
(2) Because α and IPR must be positive numbers, and
sometimes their estimates become negative due to strong
noises and other uncertainty factors, α and IPR are
used as state variables instead of α and γ ′ . Regardless of the sign of the estimated values of α and IPR ,
their squared values always give positive estimates for α
and IPR.
Consequently,
in
this
paper
c
c
⎡ IPR j (k )
α j ( k ) p j (k − 1) p j (k ) ⎤ ' is the state
⎣
⎦ vector,
(10) Using standard techniques [6], the SVM for the entire
producercentric system is given in (11) where
nx (k ) = ⎡ nx1 (k ) ' nx2 (k ) ' K nxN (k ) '⎤ ' and n p (k + 1)
⎣
⎦
are additive zeromean white noises with covariance
matrix Qx and variance rk +1 , respectively. Each
component nx j (k ) is ⎡ nα j (k ) nγ ′j (k ) 0 n pc (k + 1) ⎤ ' .
⎢
⎥
j
⎣
⎦
The explicit form of Qx is: x j (k ) , for the jth injector. As in the previous section,
the
complete
state
x(k ) = [ x1 (k ) ', x2 (k ) ',L , xN ( k ) '] ' . L , rnα N , rnγ ′ , 0, rnpc ⎤
N
N ⎦ vector is Using (8), γ ′ (k ) can be expressed in terms of IPR j (k )
j
and α j (k ) as: γ ′j (k ) = IPR j (k )(1 − α j (k )) 2 (13) which allows us to reexpress (7) in terms of the new state
variables as:
2 Qx = diag ⎡ rnα 1 , rnγ ′ , 0, rnpc , rnα 2 , rnγ ′ , 0, rnpc ,L
1
2
1
2
⎣ (11) 4 p c (k + 1) − 2 α j p c (k ) + α j p c (k − 1)
j
j
j
2 (12) 2 = IPR j (1 − α j ) 2 im, j (k ) + n pc (k )
j (14) 4 SPE 121393 x11 (k )
⎧
⎡
⎤
⎪
⎢
⎥
x12 (k )
⎪
⎢
⎥
⎪
⎢
⎥
x14 (k )
⎪
⎢
⎥
2
4
2
2
2 x12 (k ) x14 (k ) − x12 (k ) x13 (k ) + x11 ( k )(1 + x12 ( k )) 2 im ,1 ( k ) ⎥
⎪
⎢
⎪
⎢
⎥
x21 (k )
⎪
⎢
⎥
x22 (k )
⎪
⎢
⎥
⎪
⎢
⎥ + n (k )
x24 (k )
⎪ x(k + 1) = ⎢
⎥ x
2
4
2
2
2
⎪
⎢ 2 x22 (k ) x24 (k ) − x22 (k ) x23 (k ) + x21 (k )(1 + x22 (k )) im ,2 (k ) ⎥
⎪
⎢
⎥
⎨
M
⎢
⎥
⎪
⎢
⎥
xN 1 ( k )
⎪
⎢
⎥
⎪
xN 2 ( k )
⎢
⎥
⎪
⎢
⎥
⎪
xN 4 ( k )
⎢ 2
⎥
⎪
2
2
2
4
⎢ 2 xN 2 (k ) xN 4 (k ) − xN 2 (k ) xN 3 (k ) + xN 1 (k )(1 + xN 2 (k )) im , N (k ) ⎥
⎣
⎦
⎪
⎪
⎪
c
c
⎪ p(k + 1) = p1c (k + 1) + p2 (k + 1) + L + pN (k + 1) + n p (k + 1)
⎪
= [ 0 0 0 1 0 0 0 1 L 0 0 0 1] x (k + 1) + n p (k + 1)
⎪
⎩ 1
0
0
0 ⎤
⎡
⎢
0
1
0
0 ⎥
⎢
⎥
0
0
0
1 ⎥
Aj = ⎢
⎢
⎥
4 x j 2 ( k ) x j 4 ( k ) − 4 x 32 ( k ) x j 3 ( k ) +
⎢
⎥
j
2
2
4
2
− x j 2 (k ) 2 x j 2 (k ) ⎥
⎢ 2 x j1 (k )(1 + x j 2 (k )) im , j (k )
2
2
4 x j1 (k )(1 + x j 2 (k )) x j 2 (k )im , j (k )
⎢
⎥
⎣
⎦ Consequently, our modified SVM is given (15) where,
nx ( k ) and n p (k + 1) are as described in Section II.C. matrix, the explicit forms of which are given below for
our NL state equation in (15):
⎛ A1
⎜
0
Fx = ⎜
⎜ M
⎜
⎝0 (16) (18) What distinguishes our SVM in (15) from the more
general SVM in (16) is that in (15) h [•] is linear, and can
be expressed as: y(k + 1) = Hx(k + 1) + ny (k + 1) of the state vector, and nx ( k ) and ny (k ) correspond to
the additive zeromean white noises for the state and
measurement equations, respectively. In the EKF, f [•] (20) where
H = [ 0 0 0 1 0 0 0 1 L 0 0 0 1] (21) ˆ(
is linearized about x k  k ) , i.e., ˆ
[ x ( k ) − x ( k  k ) ] + nx ( k ) 0 ⎞
⎟
M ⎟
O O 0 ⎟
⎟
L 0 AN ⎠
0 L
A2 O where the explicit form of Aj is given in (19). where, in general, f [•] and h [•] are nonlinear functions ˆ
ˆ
x(k + 1) ≈ f [ x(k  k ), k ] + Fx [ x(k  k ), k ] × (19) where Fx = ∂f [ x(k ), k ] ∂x(k ) is a 4N × 4N Jacobian IV. EKF Processing
Details about the EKF can be found in [3] and [6].
Because the EKF for our modified SVM is quite similar
to that of the original SVM (see [5]) we briefly state its
structure.
The nonlinear (NL) SVM for an EKF has the following
general form:
⎧ x(k + 1) = f [ x(k ), k ] + nx (k )
⎪
⎨
⎪ y (k + 1) = h [ x(k + 1), k + 1] + n y (k )
⎩ (15) (17) is a 1 × 4N vector.
The EKF has two stages, Predictor and Corrector, and
is summarized as follows:
ˆ(0
1. Initialize the EKF with x  0) , P (0  0) , Qk and rk . SPE 121393 5 Predictor ( k = 0,1,L ): 2. ˆ
ˆ
x (k + 1  k ) = f [ x ( k  k ), k ] (22) ˆ
ˆ
P (k + 1 k ) = Fx [ x(k  k ), k ] P (k  k ) Fx′ [ x( k  k ), k ] + Qk (23) Corrector ( k = 0,1,L ): 3. ˆ
ˆ
x(k + 1 k + 1) = x(k + 1 k ) + K (k + 1) ×
ˆ
{ y (k + 1) − Hx(k + 1 k )} K (k + 1) = P(k + 1 k ) H ′
HP(k + 1 k ) H ′ + rk +1 P (k + 1  k + 1) = [ I − K (k + 1) H ] P( k + 1 k ) (24) (25) (26) V. Real Data Processing
The data used in this paper is from a section of a real oil
field. Most of the injectors in this section have two
completions; and because the injection rates of different
completions for the same injector are separate, each
completion is treated independently, i.e., each completion
is treated as an independent “injector” in our SVM. The
total number of producers in this field is denoted by N and
the total number of injector completions is denoted by M;
and the well names have been removed and relabeled.
The data starts from Jan. 1st, 2005, which is labeled as
day 1 in our figures, and ends on Jul. 31st, 2008, which is
labeled as day 1308. Both injection and production rates
are welltest data 1 . The sampling rate of injection rate is
once per day and the sampling rate of production rate is,
on average, 15 days. For days when production rate
measurements were not available, we used the value from
the last available data until the next available data.
Our strategy for processing the entire section has been
to apply our producercentric EKF to each producer
separately. We choose a producer labeled as P130 to give
an example of this process.
A very important first step is to choose an initial set of
injectors that are possibly influencing a producer. A
relatively simple and commonly used way for doing this
is to use expert knowledge provided by petroleum
engineers that are familiar with this field. For this section
it is known that there are parallel fractures along each
well that align 45o to NE. It is believed that injectors
located along this fracture alignment are more likely to
contribute to a producer. Knowing this, one strategy for
choosing the injectors that contribute to a producer is to
draw an ellipse centered about the producer whose major
axis is along the fracture alignment, and to assume the
injectors inside the ellipse are influencing that producer.
Additional expert knowledge was provided about the size
of the ellipse, i.e., in general the major axis of the ellipse
1 Pump Off Controller (POC) data for production rates are
presently not available to us. should be 700 feet long and its minor axis should 500 feet
long. Because these lengths are subjective, they may vary
from producer to producer; hence, a strategy for choosing
an optimal ellipse size is presented in Section VII.
We began by using a 700′ × 500′ ellipse for a producer.
Completions inside this ellipse are included in our EKF
model. To give the reader a clearer idea of this, P130 and
injectors inside the ellipse are depicted in Fig. 2(a) by
circle and triangles, respectively. Note that the upper and
lower triangles at the same location represent short and
long completion of an injector, respectively. And there are
46 completions in this model.
After applying our EKF to this single producer46
completions model, we obtained the IPR curves that are
depicted in Fig. 3(a).
46 completions is not a small number, so it is very
likely that some of these completions may be irrelevant to
P130. Observe, from Fig. 3(a), that there are many IPR
curves that only have very small values. We assumed,
therefore, that the completions that have very low IPR
values do not influence P130 and should be eliminated
from the model. More specifically, we computed the
mean of the IPR values for the most recent month (i.e.,
days 12781308) for each completion and compared it
with a chosen threshold. A completion was kept in the
model only if its IPR value was greater than or equal to
the threshold; otherwise, that completion was eliminated
from the model.
Our strategy for choosing a good threshold is
dependent on the number of completions inside the
ellipse, and was inspired by the fact that an existing
common, simple and practical way for deciding the
impacts of injectors in waterflood management is to
assign equal weights to them. For example, if N injectors
are thought to be impacting a producer, then one assumes
each of them has an impact weight of 1/N. We modified
this by using ρ / N as the threshold, where we found, by
trial and error, that 80% of 1/N does a good job as a
threshold.
The IPR curves in Fig. 3(a) are in absolute values, but,
to use a percentage threshold they need to be normalized.
This is done by dividing each IPR curve by the sum of all
the IPR curves. These IPR curves in percentages are
depicted in Fig. 3(b). Also shown on that figure is the
80%/N =1.74% threshold line (dotted) when N = 46.
After eliminating all of the completions whose IPRs
fell below 1.74%, only 17 completions were left. The
remaining completions and eliminated completions are
depicted in Fig. 2(b) by regular size triangles and smaller
triangles, respectively.
Our EKF was then reapplied to this single producer,
but for only the 17 remaining completions. The resulting
IPR values are depicted in Fig. 4. Observe that, although
some portion of the three lowest IPR curves are still
below the threshold, we are only considering the mean
IPRs of the most recent month, which occurs at the rightend of the data, and none of the mean IPR values for that
month fall below this threshold. Consequently, none of
the 17 injectors were eliminated, and therefore, our 6 SPE 121393 processing for P130 is completed. If, perchance, some of
the 17 injectors had been eliminated, we would have then
repeated this procedure until a situation was reached
where no more of the injectors are eliminated 2 .
The EKF procedure that we have described is
performed for each of the producers in the section. If
parallel processing is available then this processing can be
done in parallel because of its producercentric nature.
VI. Validation Using History Matching
How does one validate our EKF method? Unless a field
test is performed, there are no truth data available as there
would be when synthetic data or reservoir simulation data
are used (as in [5]). Although field tests can be performed,
they are disruptive and expensive. Another approach, the
one we have used is history matching, i.e., we go back in
history and locate a time point at which a significant step
change in an injection rate has occurred (as suggested to
us by petroleum engineers), and use the estimated IPR at
that point to “forecast” the production rate at some future
time. This forecasted production rate is then compared
with the real production rate, since the latter is available
to us as a part of the historical data record we began with.
Forecasting of the production rate is accomplished by
using the predictor equation (21) of our EKF. We
illustrate this process next by an example.
By carefully inspecting the historical injection rates for
all the remaining 17 completions for P130, we noticed
that most of them had a significant rise in their injection
rates around day 1045; therefore, day M=1045 was
chosen as the starting day for our forecasts. Starting from
day M, we used the estimated parameter values obtained
at that day, i.e.: ˆ
ˆ
ˆ
x j ( M  M ) = ⎡ IPR j ( M  M ), α j ( M  M ),
⎢
⎣
(27)
ˆ c ( M − 1 M ), p cj ( M  M ) ⎤ ', j = 1,...N
ˆ
pj
⎦
and the planned injection rates, i.e., i j ( M ),L, i j ( M + 30)
( j = 1,..., N ) , to forecast the next month’s production
rates, daybyday, i.e.: ˆ
⎧ x( M + 1 M ) = f [ x( M  M ) ]
⎪ˆ
⎨
ˆ
⎪ p ( M + 1 M ) = Hx( M + 1 M )
⎩ˆ
ˆ
⎧ x( M + 2  M ) = f [ x( M + 1  M ) ]
⎪ˆ
⎨
ˆ
⎪ p ( M + 2  M ) = Hx( M + 2  M )
⎩ˆ
M (28) ˆ
⎧ x ( M + 30  M ) = f [ x( M + 29  M ) ]
⎪ˆ
⎨
ˆ
⎪ p ( M + 30  M ) = Hx( M + 30  M )
⎩ˆ
where H is the same as in (21).
2 During such an iteration process, the threshold would be kept
at 80%/N, where N is the initial number of injectors, in this case,
46. N is not chosen to be the surviving number of injectors
because the resulting threshold would become so large that it is
possible that too many injectors would be eliminated. After one month, new injection and production rate
measurements [ im, j (k ) and pm (k ) ( j = 1,..., N ;
k = M + 1,..., M + 30) ] are available, so we update our
Section IV EKF using those new measurements, to obtain
ˆ(
the updated state vector, i.e., x M + 30  M + 30) , that are
then used to forecast another month’s production rates by
using (28) again. This process is repeated until we have
run out of data history.
We observed that IPR curves are not flat during an
entire month (e.g., see Fig. 4); hence, we computed their
gradients at the day the forecast starts, and assumed the
IPRs would follow the gradients during the following one
month’s forecasts, e.g., the first month’s forecasts start at
ˆ
M=1045, and the gradient of
IPR ( M  M ) and
j ˆ
α j (M  M ) , g IPR j ( M ) and g αj ( M ) ( j = 1,..., N ) , can be approximated as:
⎧
ˆ
ˆ
IPR j ( M + 1  M ) − IPR j ( M  M )
⎪g
(M ) =
⎪ IPR j
( M + 1) − M
(29)
⎨
ˆ
ˆ
α j (M + 1  M ) − α j (M  M )
⎪
⎪ g α j (M ) =
( M + 1) − M
⎩
It follows, therefore, that:
⎧ IPR ( M + 1  M ) = g
ˆ
ˆ
( M ) + IPR j ( M  M )
j
IPR j
⎪
(30)
⎨
ˆ
ˆ
⎪ α j ( M + 1 M ) = g α j ( M ) + α j ( M  M )
⎩ Because we assume IPR j and α j follow the gradient for the entire month, (30) provides the following daily
interpolation formulas that were used by us ( d = 1,...,30 ):
⎧ IPR ( M + d  M ) = g
ˆ
(M ) +
j
IPR j
⎪
⎪
ˆ
IPR j ( M + d − 1 M )
⎪
⎨
ˆ
⎪ α j (M + d  M ) = g α (M ) +
j
⎪
⎪
ˆ
α j (M + d − 1  M )
⎩ (31) Results have shown that forecasts incorporating these
interpolations outperformed the forecasts that did not use
them, so they have been used in all of our results that are
described below.
The overall historical production rate for P130 is
plotted as the solid line in Fig. 5, whereas the dotted line
shows the monthly forecasts, starting at day 1045.
Observe that the production starts to increase after day
1045 and that the forecasted production has detected this
increase, although there are some forecasting errors. For
viewing convenience, we zoom in to the prediction
interval and this is plotted in Fig. 6(a), and, the prediction
error is plotted in Fig. 6(b). SPE 121393 Some summaries of the monthly errors are given in
Table 1. It contains monthly error summaries for different
values of a scalar parameter, s, which is used in Section
VII to determine the size of the ellipse for establishing the
initial set of injectors. Our present case corresponds to the
row in Table 1 for which s=1.0.
The monthly forecasting began on day 1045 and
continued for 8 months. In Table 1, the average daily
prediction errors for each month, and the average daily
error over the entire eightmonth period, have been
computed. Additionally, ErrorProductionRatio (EPR),
which is the ratio of the average daily prediction error to
the average daily production rate over 8 months, has been
computed. EPR is useful because petroleum engineers
have told us that the measurementnoise level for this
section is around 1020%, e.g. if the measured production
is 300 barrels, then about 3060 barrels is noise. The EPR
for s = 1 is 10.23% which is a very good value.
VII. Optimization of Initial Set of Injectors Using
History Matching
An issue that plays a very critical role in our EKF
processing is how to choose a set of initial injectors for
the producercentric model. Although it has been
previously mentioned that an ellipse with 700 feet major
axis and 500 feet minor axis is used to select the initial set
of injectors, in general, an optimal size for the ellipse may
be very different for different producers. If the ellipse is
too small, one would miss some influential injectors from
the very beginning. On the other hand, if the ellipse is too
large it could include many irrelevant injectors. Such
injectors may significantly bias the first round of EKF
processing, because they can affect the threshold
elimination process that has been described in Section V.
Inspired by techniques used in standard optimization
problems, and the history matching results described
Section V, our strategy for finding an optimal ellipse size
(for each producer) was to use the average daily
prediction error as an objective function and to minimize
it with respect to ellipse size. To do this, we used an
ellipse with 700 feet major axis and 500 feet minor axis as
a standard ellipse and then scaled its major and minor
axes by the same scalar, s; hence, as different values of s
are used one obtains ellipses of different sizes but with the
same ratio of major to minoraxes.
In theory, s should be discretized very finely and also
range from a small enough value to a large enough value
so that the computed minimal point of the objective
function is its global minimum. Unfortunately, using too
many values of s is computationally very costly; hence, in
this study we chose s = {0.5, 0.6, 0.7, 0.8, 0.9, 1.0}. We
have observed that when s = 1.0 (which is the case
discussed in Section V) the ellipse is already large enough
and usually includes more than 20 injectors, and (as we
will see below), when s = 0.5 the ellipse only includes a
few neighboring injectors. Additionally, we observed,
after several tests, that s does not have to be discretized
very finely, e.g., when s changes from 0.95 to 0.9 the
ellipse only shrinks a little and no injectors from the s =
0.95 ellipse are left out of the s = 0.9 ellipse. 7 EKF processing for different values of s proceeds in
exactly the same manner as described in Section IV for s
= 1. Results for s = 0.9, 0.6 and 0.5 are given in Figs. 715. Summaries of the prediction errors for s = {0.5, 0.6,
0.7, 0.8, 0.9, 1.0} are also given in Table 1. Observe that
the objective function is minimized when s = 0.6. Thus,
for P130, the optimal ellipse used to choose the initial
injector set is one with 0.6 × 700 = 420 feet major axis
and 0.6 × 500 = 300 feet minor axis.
VIII. IPR Table: Producer Centric to Injector
Centric Conversion
The previous sections have presented a complete
procedure to process only one producercentric model. Of
course, this procedure has to be applied to every producer
in the entire section (or in the entire field, if such data are
available) The IPR results of doing this can then be
summarized in a table, which in our case has N rows
(producers) and M columns (injectors). A small portion of
this table is depicted in Table 2. The nonzero values in
this table are for those injectors that have survived our
EKF processingprocedure. Zero values are for the
remaining injectors. This table allows one to quickly look
up the IPR value between any producerinjector pair.
Reservoir engineers are more interested in an injectorcentric viewpoint than in a producercentric viewpoint,
because they have control over water allocation at each
injector; hence, for them it would be better if the results
were injectorcentric. Observe that Table 2 can be viewed
in two different ways. Data are entered into it one row at a
time (producercentric), but it can then be viewed
columnwise (injectorcentric). Unfortunately, the
numerical IPR values are still difficult to interpret. One
solution to this is to normalize the IPR values into
percentage values by dividing the values in each column
by the sum of values in the column. This allows us to
observe what percent of the total water that is allocated to
an injector goes to a specific producer. We show such
results in Table 3 for two injectors. Observe that very
large amounts of water from injectors I322 and I344 go
to producers P124 and P146, respectively.
Fig. 16 depicts a novel graphical way to represent these
injectorcentric results. Observe that we have connected
an injector to all of the producers that have nonzero IPRs
using lines of different widths. The higher the %IPR value
is, the thicker the line is, and the lower the %IPR value is
the thinner the line is.
IX. Conclusion and Future Research Topics
This paper has continued the work of [5]. After extensive
testing of the EKF on real data, important modifications
have been made by us to the SVM, allowing the square
root of IPR values to be estimated directly, which avoids
estimating negative results for IPR values. We have also:
described a way to integrate important expert knowledge
about the oil field during the initialization of the EKF
processing; presented a method to eliminate unimportant
injectors from the model; and shown how to validate the
accuracy of our EKF results by using productionrate
historymatching. Our validation process has also led us 8 SPE 121393 [2] A. N. AraqyeMartinez, “Estimation of autocorrelation and
its use in sweep efficiency calculation,” MS thesis, U. of
Texas at Austin, Austin, Texas, 1993.
[3] S. Haykin, Kalman Filtering and Neural Networks, John
Wiley & Sons, New Jersey, 2001.
[4] K. Heffer, K. J., Fox, R. J., and McGill, C. A., “Novel
Techniques Show Links Between Reservoir Flow
Directionality, Earth Stress, Fault Structure and
Geomechanical Changes in Mature Waterfloods,” Paper
SPE30711 Presented at SPE Annual Technical Conference
and Exhibition, Dallas, Oct., 1995.
[5] F. Liu, J. M. Mendel and A. M. Nejad “Forecasting InjectorProducer Relationships from Production and Injection Rates
Using an Extended Kalman Filter”, SPE 110520, Presented
at SPE Annual Technical Conference and Exhibition,
Anaheim, CA, Nov. 2007.
[6] J. M. Mendel, Lessons in Estimation Theory for Signal
Processing, Communications, and Control. Prentice Hall
PTR, Upper Saddle River, NJ, 1995.
[7] M. N. Panda and A. K. Chopra, “An integrated approach to
estimate well interaction, ” Paper SPE 39563, presented at
1998 SPE India Oil and Gas Conference and Exhibition,
New Delhi, Feb., 1719, 1998.
[8] De Sant’Anna Pizarro, J. O., “Estimating injectivity and
lateral autocorrelation in heterogeneous media,” Ph.D
dissertation, U. of Texas Austin, Texas, 1998.
[9] B. T. Refunjol, “Reservoir characterization of North Buck
Draw
field
based
on
tracer
response
and
production/injection analysis,” MS Thesis, Univ. of Texas
Austin, Texas, 1996.
[10] T. Soeriawinata and M. Kelkar, “Reservoir management
using production data,” Paper SPE 52224, presented at
1999 SPE MidContinent Operation Symposium,
Oklahoma City, Oklahoma, March 1999.
[11] A. A. Yousef, P. Gentil, J. L. Jensen and L. W. Lake, “A
Capacity Model to Infer Interwell Connectivity From
Production and Injection Rate Fluctuatuions,” SPE 95322,
Presented at SPE ATCE 2005, Dallas, Oct., 2005. to develop a method for selecting an optimum ellipse size
(for each producer) for choosing the initial set of injectors
by minimizing historymatching errors. Finally, we
showed how to convert a table of producercentric IPR
values to a table of injector centric IPR values.
Although the EKF gives good results, there still is room
for improvements. Using more frequently sampled POC
data will very likely improve the estimation. By using
POC data, it may be possible to group several producers
together, in which case our reservoir model changes from
a singleproducer–multipleinjector model to a multiple
producer–multiple injector model. The EKF is still
applicable to such a model. Finally, it may be possible to
treat each actual injection rate as a state variable so that
noisy injection rates do not have to be used in the
predictor equations. The downside to doing this is an
increase in the dimension of the EKF.
Acknowledgments
This study was funded by the Center of Excellence for
Research and Academic Training on Interactive Smart
Oilfield Technologies (CiSoft). CiSoft is a joint
University of Southern CaliforniaChevron initiative. The
authors would like to thank Prof. Iraj Ershaghi, Dave Tuk,
Jim Brink and John Houghton for their valuable
discussions. The authors would also like to thank
Hyokeyong Lee for generously providing computer
programs that preprocess the raw data, and other
members from Prof. Cyrus Shahabi’s group for their
useful discussions on data issues.
References
[1] A. Albertoni and L. W. Lake, “Inferring Interwell
Connectivity only from WellRate Fluctuations in
Waterfloods,” SPE Reservoir Evaluation & Engineering,
Vol. 6, Number 1, Feb., 2003. Table 1. Measures of average daily prediction error of production rates for eight months, including:
overall average daily prediction error (barrels/day), and average daily error to average daily production
rate ratio (EPR), all for different size ellipses (s).
st S=1.0
S=0.9
S=0.8
S=0.7
S=0.6
S=0.5 1 month
7.3067
3.7455
7.2293
4.4729
5.8975
9.7345 nd rd th th th th th 2 month 3 month 4 month 5 month 6 month 7 month 8 month
16.259
66.977
34.464
34.588
23.108
75.625
19.180
15.506
38.252
26.247
38.229
22.287
43.321
13.253
15.273
33.060
16.235
43.111
13.786
47.965
13.034
14.829
26.305
17.878
43.224
14.721
35.270
13.643
14.161
26.667
18.717
43.360
14.663
51.121
20.796
16.310
26.735
18.209
48.188
10.872
56.535
11.865 average
34.688
28.779
26.583
13.268
13.234
13.833 EPR
10.23%
8.49%
7.84%
3.91%
3.90%
4.08% SPE 121393 9 Table2. Upper left portion of the IPR table for the entire section of 191 producers (rows) and 374 injectors
(columns) 3 . (Well names have been relabled.)
P1
P2
P3
P4
P5
P6
P7
P8
P9
P10
P11
P12
P13
P14
P15
P16
P17
P18
P19
P20
P21
P22
P23
P24
P25
P26
P27 I1
0.0211
0
0
0
0
0.0683
0
0
0
0.0946
0.1062
0.0678
0.0512
0
0
0
0
0
0.0360
0
0
0
0
0
0
0
0 I2
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0.0654 I3
0.1571
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 I4
0
0
0
0
0
0
0
0
0
0.0852
0.0782
0.0541
0
0
0
0
0
0
0
0
0
0
0
0
0
0.0536
0 I5
0
0
0
0
0
0
0
0
0.1429
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 I6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 I7
0.0692
0
0
0
0
0
0.0339
0
0
0.1477
0
0.1014
0
0
0
0
0
0
0
0
0
0
0
0
0
0.0504
0 I8
0
0
0
0
0
0
0
0
0.1334
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 I9
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0 I10
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0.0402
0.0329
0.0543
0
0
0
0
0
0.0911
0
0
0
0 Table3. Normalized IPRs for injectors I322 and I344. Entries are in %. (Well names have relabled.)
P124
I322 P170 P5 P116 P167 P191 41.77 21.42 20.51 16.16 0.54 0.20 P146
I344 P139 P180 P39 P135 41.77 17.70 14.32 14.14 12.07 Figure 1: The Reservoir Model for reservoir with N injectors and single producer [5]. 3 All the entries in the IPR table have not been normalized yet. 10 SPE 121393 Figure 2: When s = 1, (a) P130’s local area and 46 initial completions included in our model; and (b) completions after our
elimination process. (Well names have been relabeled.) Figure 3: When s = 1, (a) IPR curves for the 46 initial completions; and (b) normalized IPR curves (solid lines) and the elimination
threshold (dotted line). Figure 4: When s = 1, (a) IPR curves for the 17 remaining completions after the elimination process; and, (b) normalized IPR curves
(solid lines) and the elimination threshold (dotted line). SPE 121393 11 Figure 5: Historical production rate for P130 (solid line) and predicted production rate (dotted line) when s = 1. Predictions begin at
day 1045. Figure 6: When s = 1, (a) Historical production rate forP130 (solid line), and predicted production rate (dotted line); and, (b) error
between the real production rate and the predicted production rate. Figure 7: When s = 0.9, (a) P130’s local area and initial completions included in our model; and (b) completions after our elimination
process. (Well names have been relabeled.) 12 SPE 121393 Figure 8: When s = 0.9, (a) IPR curves for the 10 remaining completions after the elimination process; and, (b) normalized IPR curves
(solid lines) and the elimination threshold (dotted line). Figure 9: When s = 0.9, (a) Historical production rate for P130 (solid line), and predicted production rate (dotted line); and (b) error
between the real and the predicted production rate. Figure 10: When s = 0.6, (a) P130’s local area and initial completions included in our model; and (b) completions after our
elimination process. (Well names have been relabeled.) SPE 121393 13 Figure 11: When s = 0.6, (a) IPR curves for the 9 remaining completions after the elimination process; and, (b) normalized IPR curves
(solid lines) and the elimination threshold (dotted line). Figure 12: When s = 0.6, (a) Historical production rate for P130 (solid line), and predicted production rate (dotted line); and (b) error
between the real and the predicted production rate. Figure 13: When s = 0.5, (a) P130’s local area and initial completions included in our model; and (b) completions after our
elimination process. (Well names have been relabeled.) 14 SPE 121393 Figure 14: When s = 0.5, (a) IPR curves for the 8 remaining completions after the elimination process; and, (b) normalized IPR curves
(solid lines) and the elimination threshold (dotted line). Figure 15: When s = 0.5, (a) Historical production rate for P130 (solid line), and predicted production rate (dotted line); and (b) error
between the real and the predicted production rate. Figure 16: Graphical representations of normalized IPRs for (a) injector I322 and (b) injector I344 (see Table 3). (Well names have
been relabeled.) ...
View
Full Document
 Spring '11
 tan
 producer, IPR, injector, The Producer

Click to edit the document details