WO2006130010A2

WO2006130010A2 - Method for processing sampled data

Info

Publication number: WO2006130010A2
Application number: PCT/NO2006/000077
Authority: WO
Inventors: Jan Erik Nordtvedt; Stig Olsen
Original assignee: Epsis As
Priority date: 2005-06-03
Filing date: 2006-02-28
Publication date: 2006-12-07
Also published as: NO20052674D0; NO322656B1; NO20052674A

Description

METHOD FOR PROCESSING SAMPLED DATA

The present invention concerns a method for processing sampled data for the purpose of data transmittal and/or data storage.

Prior art

Real time data is playing an increasingly important role in today's energy and petroleum industry. A distinct development trend is that more and more realtime data is becoming available for operators and engineers to aid analysis and decision-making within production optimization, operations and maintenance as well as within drilling. Although acquired at field location or in wells, the data may be transmitted to off-site locations and available in real-time for analysis and decision-making. For most onshore or offshore assets, typically several thousand data sources may be measured, many of these sources are polled on a second-by-second basis.

Although this development offers a series of opportunities with respect to making sound operational decisions, it also poses quite a few challenges. For example, with a high sampling rate, any sensor may give millions of data points each year. To store and administer this data to make them easily retrievable whenever needed becomes a practical hurdle in the utilization of the data. Also using the data points directly in analysis applications may slow down the calculation significantly. In addition, data outliers may cause unwarranted alarms to go off, resulting in additional and unnecessary work being initiated.

This invention aids the practical use of production data in the oil and gas industry, by providing a system that can automatically de-noise and compress the data from any sensor or system - without any human intervention. In addition, there is no need to set up the invented system specifically for any data source, e.g., the standard deviation of the signal to de-noise is detected automatically as a part of the de-noising procedure. This means that the filter system can be deployed on the site installation, office location or anywhere in between without time-consuming set-up or maintenance routines. Data filtering and data compression becomes important in order to reduce the amount of data and to extract the important information from said data. This lead to the development of de-noising methods in order distinguishes the important trends from the noise, and compression in order to reduce the amount of data without losing the important information the data carries. The advantages of de-noising and compression of signals are numerous and may be summarized as follows:

• Speed up of analysis by data compression. • Remove outliers that seriously can affect regression results.

• Compress large amount of data for efficient storage or transmission.

• Detect underlying trends or changes in signals.

• Obtain visual better-looking signals (makes decision-making easier).

In the last decade the literature has been dominated by a technique called waveshrink or wavelet filtering. Wavelet filtering has proved to be an efficient way of finding structure in a signal contaminated with noise without assuming any particular parametric regression model on the function f eL²(m).

y,^. = f, +oε, ' = 1.2,3, n. (1 )

where ε is assumed to be a normal distributed variables with mean μ and standard deviation σ. The noise standard deviation might or might not be known.

The reasons for the popularity of wavelet filtering is due to its spatial adaptively, computational effectively and asymptotic optimal properties. Most signals found are non-stationary signals and the wavelets spatial adaptively make them ideal for analyzing such signals.

A large number of wavelet filtering methods have been proposed, and some of these methods have been considered in petroleum industry. It is however mainly outside the petroleum industry the main applications of wavelet filtering can be found. Among the numerous disciplines where wavelet filtering is extensively used are seismic and geophysical signal processing, medical and biomedical imaging, computer vision and communications. Several studies comparing various wavelet filtering methods are also published. Other studies containing limited amounts of wavelet filtering techniques are also published, often in relation with recent proposed methods.

Although numerous wavelet-based techniques have been presented, several problems still remain before an effective tool exists. Some wavelet-based techniques have proven to give good performance in some cases; however, no comprehensive comparison exists guiding the selection of wavelet techniques for use with production data. For most of the proposed techniques several parameters such as noise level, primary resolution level, compression ratio, wavelet type as well as threshold values and wavelet shrinkage rule need to be determined. This makes these techniques time consuming and requires expert users for effective use and sound performance. Solving all these problems holds the promise of automatic data filtering and compression, being of crucial importance for smart production where analysis modules need to access data with minimum human intervention between the data source and the analysis modules.

One object of the present invention is to introduce an improved method for wavelet filtering. This will include proposing new improved noise, level, threshold and compression ratio estimators that perform well compared to others and substantial better than some of the most common used wavelet filtering methods.

Another object with the present invention is to provide a method for determining the standard deviation of a signal automatically for optimal noise reduction, i.e. to remove as much noise as possible from the signal without over smoothing the signal.

Another object of the present invention is to provide a method for automated wavelet filtering and data compression to obtain improved results for data preparation.

Summary of the invention

The present invention introduces a comprehensive analysis of the noise reduction and compression performance of different types of wavelet based techniques. New methods for estimation of noise level, primary resolution level, threshold value as well as wavelet and wavelet shrinkage methods are presented. In addition it is also presented a methodology that automatically selects the optimal compression ratio for any given signal.

The present invention provides a method for removing noise from pressure and rate signals. The following procedure is then used to perform real-time data cleansing:

1. Perform a real-time wavelet transform

2. Perform real-time noise estimation

3. Perform outlier removal 4. Perform signal de-noising

5. Perform data compression

The objects mentioned herein and other objects are obtained by providing a method for handling sampled data for the purpose of data transmittal and data storage, characterized by the features of claim 1.

Preferred embodiments of the method according to the present invention are specified by the characteristic features in the following dependent claims 2-12.

An application of the method according to the invention is described in claim 13.

It appears that much of the focus in wavelet filtering has been to remove all of the noise in the signal resulting in over smoothing the signal since a substantial amount of the signal itself is removed along with the noise. Using the approach according to the present invention a large portion of the noise is removed without destroying or over smoothing the underlying signal.

One advantage of the present invention is that the new filtering techniques and recommendations according to the invention can produce up to 90% noise reduction. In situations with few transients, more than 95% removal of noise can be observed.

The present invention will now be described with reference to the accompanying figures, wherein Figure 1 shows an area of application of the method according to the present invention.

Figure 2 shows calculation of scaling coefficients according to the state of the art. Figure 3 shows calculation of scaling coefficients according to the state of the art.

Figure 4 shows the process of new points according to the present invention added in figure 2. Figure 5 shows update of the median according to the method in the present invention.

Figure 6 shows a graph of collected data before and after filtering according to the present invention.

Figure 7 shows a graph of the error estimate versus filter level according to the present invention.

Figure 8 shows a graph of the filter level versus number of data points according to the present invention.

Figure 9 shows a change of wavelet coefficient according to the present invention. Figure 10a-c shows a compression method according to the present invention.

Figure 11 shows a graph for determining optimal compression ratio according to the invention.

Figure 1 shows a field where the method in the present invention can be used. The automatic filtering and compression system can be put in between the sensor 4 and the control system 2 on the site location (onshore/offshore) filtering and compression data before entering into the control system 2 itself (thus, it becomes a part of the data acquisition for the control system). Alternatively, the filtering and compression system may be placed on the servers 6/historians at the site before data storage and/or transmission. As an alternative or supplement, the system can be implemented on the office location in association with historians or application servers/systems filtering and compressing data before or in association with analysis. Also, the system may be implemented on off-site networks as shown in Figure 1.

The methods according to the present invention are used to remove outliers, reduce noise and compress real-time data. As an initial step of the procedure, the wavelet for the specific data type is selected and the window size for the outlier removal procedure is selected. Various types of wavelet adapt well to different types of signals, and due to the different properties of rate and pressure signals, different wavelets should be used. Various types of wavelets adapt well to different types of signals, and due to the different properties of rate and pressure signals, different wavelets should be used. Based on analyses of production data, the following wavelets have been selected to exemplify: • Pressure, temperature, level: spline39

• Flow rates, choke sizes: Spline4246

Spline39 indicates that the low pass and high pass filter used are of a different filter size, spline39meaning that the low pass filter has a filter size of 3 and the high pass filter has a size of 9.

It should be noted that the wavelets used in the present invention are not limited to the wavelets according to the example, and that various wavelets can be used for obtaining optimal results according to the input signal. The inventors of the present invention recognize that the selection of good wavelets as well as filter settings can be crucial. With this motivation, the performance of several known de-noising methods was investigated and new methods proposed were feasible.

Real time wavelet transform Normally the wavelet transform is performed on a fixed signal, and modifications are required if the wavelet transform should be utilized on continuous data streams. Using the normal wavelet transform in the case of real-time data would require great computational effort, since the wavelet transform rapidly will be recalculated as new points entered. In addition, boundary effects in a wavelet transform will disturb the last points. These points are the ones normally monitored by an operator, and they might in addition be connected to alarms to detect abrupt signal changes. Consequently, it will be vastly desirable to transform these points without interfering boundary effects.

To achieve this, the present method uses a moving wavelet window approach. Based on asymptotical consideration, this window should be as large as possible, and based on computational effort, as small as possible. Also the filter size and filter level will place limitation on the sample size. If a filter with size 4 is used for filtering to level 5, a sample size of 4⁵=1024 is required in order to calculate the approximate signal at level 5. In general the number of signal points required is (filter size)^level . Figure 2 shows determination of the scaling coefficients using two levels and filter size 2. Boundary points 2 can be seen to the right, affecting more and more points with higher filter level.

Figure 2 shows determination of the scaling coefficients for the wavelet analysis after eight data 10 points using two levels and a filter size 2. Leaving the boundary points 12 in fig. 2 out of the calculations will result in a wavelet transform 20 as shown in fig. 3. If the window size in this example is chosen to eight, and only two filter level is considered, fig. 4 illustrates what will happen as a new point 28 appears. The new point 28 will require the calculation of two new scaling coefficients 26. The same applies to the wavelet coefficient where two new coefficients are added. In addition to the calculation of four new coefficients (two wavelet + two scaling), four points 30 need to be deleted once the window size reaches its limit, and they are. Thus, the transform will only result in a small number of calculations compared to the normal wavelet transform, and it will eliminate the influence of boundary effects.

There will however be a delay in the approximate signal and the detailed signal depending on the filter size and the approximate level, and the approximate and detail signals at these levels will only be close to real time.

Real time noise estimation

The noise estimate can be performed by first calculating the median value of the detailed wavelet coefficients at level one in the wavelet transform. This is denoted MAD. The noise estimate σ can then be calculated:

MAD σ = (2)

0.6745

In this invention, MAD is normalizing to the type of wavelet used and the data stream in question and by excluding all zero coefficients. As a result, a robust noise estimation is obtained. For real time data, the moving median noise estimation is given as:

MAD_mmtng__wMm ^) ^{" a} where the constant a is selected depending on various signals. The update of the median is fast and quickly accomplished with removal and insertion in a continuously sorted list as shown in fig. 5.

The present invention introduces new values for the noise normalization value a, producing an accurate estimate of the standard deviation, i.e. the noise estimation σ . RMAD, Robust Median value, is determined by determining the median of the wavelet coefficients at level one and removing all zero wavelet coefficients. An optimal noise normalization value a through regression analysis, using a series of data sets for pressure and production rate with known standard deviation. Two new estimators named Medium and High Robust Median Estimator are provided for the earlier mentioned two types of wavelets used in production data:

Medium RMAD estimator for use with rate data using spline 4246 is proposed as;

_ RMAD (4) ^{σ ~} 0.5164 High RMAD estimator for use with pressure data using spline 39 is proposed as;

_ RAdAD (5J ^σ ~ 0.4193 This set-up requires no user input and can hence be done automatically for any data source. The values for the constant a are mostly preferential in connection with the mentioned wavelet, but a deviation of at least 50% can also be acceptable.

The new estimators combined with the correct wavelet is not sensitive to the number of rate changes, number of zero rate, number of points (unless very small) or the magnitude of the rate changes. For other wavelets then the two mentioned above, clear dependency of for instance number of rate changes was detected.

It should be noted that the above-mentioned coefficients in formulas (4) and (5) are valid with the proposed spline wavelet, and that the scope of the invention is not limited to the use of these wavelets, and that other coefficient values can be determined through regression analysis based on the input signal and the selected wavelet.

Outlier removal Outliers can be removed using a median filter, a window filter and requires an input of the window size and an input of the estimate of the noise. The window size allows us to specify the largest number of successive outliers that can be removed. The value of this parameter should be chosen higher than the largest number of successive outliers to be removed and less than the number of data points in the smallest transient. In general, it should not be difficult to set a quite general value for the window size, depending on the data source resolution.

The method works such that every point that has larger deviation from the trend than the estimated noise level is defined as outliers. In order to avoid removing any points that are solely affected by noise, three times the noise estimate is used for the outliers. Thus, if the investigated data point has a value three times higher than the median value in the current window, the value is replaced by the median value.

A data set with three standard deviation white noise and up to seven successive outliers was generated in a random manner and with random outlier amplitude. By the use of the median filter and noise estimation all outliers were removed, single outliers as well as successive outliers. A zoom in on the data can be seen in figure 6.

Real time signal de-noising

The process above has described how to obtain the wavelet and the scaling coefficient. It is known that a signal can be filtered by removing or shrinking the wavelet coefficients in the following process;

Z = D(Y,λ) ^W

S = W^A (Z) where S is the signal, W is the wavelet transform, D is the threshold operator, A is the threshold value and W¹ is the inverse wavelet transform.

The idea behind wavelet shrinkage is that only a few large wavelet coefficients are due to the signal while the rest and small ones are due to noise, meaning that the signal can be de-noised by shrinking all small coefficients toward zero and keeping large ones.

For utilizing wave shrink, several parameters need to be determined. This includes noise estimate (described earlier), filter level, wavelet (described earlier), shrinkage rule and threshold value. The following optimal settings have been identified according to the present invention:

Primary resolution level Some work as described earlier, has been done regarding determination of the optimal resolution level. The common idea is to filter up to a fixed or estimated resolution level and ignoring the coarsest levels. However, in terms of minimizing the difference or square difference error, filtering up to maximum possible level defined in the wavelet transform appears in most cases to be optimal. Figure 7 shows error decrease in filtering as filter level is increasing for a pressure and a rate signal containing 32.000 points. The error shows a monotonic decrease as the primary resolution level is increasing in the filtering.

This, however, assumes that the optimal thresholds are chosen at each level, and it would be erroneously to apply this approach together with a threshold that would seriously over smooth the signal.

However, the gain of filtering the last levels is usually small. A more practical level for filtering would be to filter at least one level and then increase the level as the number of points increase until the number of levels reach nine, above which small enhancement is found. Figure 8 shows optimal filter level for various signal sizes in the sense that only minor improvement in filtering is possible above these levels.

One should also stop filtering when there is small improvement in adding additional levels, as shown in fig. 8. This gives the formula for the level estimation as

where n is the number of data points.

Shrinkage Based on a large number of filtering using various wavelets and shrinkage functions the following shrinking rule, named SCAD, is preferred to perform shrinkage for production data:

(8)

where λ is the threshold value, S is the signal, and a =3.7, based on a Bayesian argument. This is a keep, shrink or kill rule.

It should be noted that the present invention is not limited to using the SCAD shrinkage function, but other shrinkage functions can be used if desirable, such as Hard and Soft Shrinkage, Garrote Shrinkage, or Firm Shrinkage.

Threshold The thresholds will typically decrease exponential as the level increase. The threshold will in addition to the level, depend on the noise and the number of data points. It appears that an empirical threshold of the type

RED = axσx{b x log(«))^feve/ (9)

can be applicable. The parameters a and b are determined through regression analysis, and n is the number of data points. For pressure data using spline39 and high RMAD the following formula for determining threshold is obtained:

RED_hlgh = 2.2σ(0.381og(«)y^CTe/ (10)

or for rate data using spline4246 and medium RMAD the following formula is obtained:

The use of the above-mentioned formulas for thresholds gave excellent results for pressure data and rate data, and they performed close to optimal for the various shrinkage rules. It should be noted that the values for a and b can vary since they are determined by regression analysis, and will therefore be affected by use of other models in connection therewith. Also, the values for a and jb in connection with the spline wavelet are mostly preferred, but a deviation of at least 50% can be acceptable. Other models for threshold, such as first and second order polynomials were also attempted, but these models were not successful. It should be noted that various models for threshold can be selected depending on the input signal, and that method described herein is only an example for clarification.

Figure 9 shows that a change in a wavelet coefficient can cause many updates in the signal. The filtering and inverse transform is slightly more complicated than the forward transform. The reason for this is that changing wavelet coefficients at high levels will cause changes in many signal elements and therefore require many updates each time a change occurs.

In real time filtering, the forward transform is first performed as described earlier in figure when new points enter. This will only cause a limited amount of new calculations. The new wavelet coefficients are then tested whether or not they should be kept, reduced or killed depending on the shrinkage rule and threshold values. The threshold value often comprises a noise estimate. If some of the wavelet coefficients are changed, the calculation as shown in fig. 8 must be performed in order to update the signal.

This means that a point will be more and more filtered as higher-level detail coefficients are included in the filtering, and the operator will see a gradual signal smoothing.

De-noising in sum

For pressure, temperature and level data de-noising can be accomplished by using noise estimator high RMAD (5), level estimation (7), spline39, SCAD shrinkage (8) and high RED threshold (10).

For rate and choke size this should be accomplished by using noise medium RMAD (4), spline4246, level estimation (7), SCAD shrinkage (8) and low RED threshold (11).

Data compression After wavelet filtering, one should apply a compression filter. This will in particular be important for rate data. The main purpose of the compression will be to get rid of rate changes due to noise more than to simply reduce the share number of data points. The rate changes will complicate the calculations. For accurately compressing rate data, outliers are removed firstly following de- noising of data with an appropriate wavelet filter as described in previous section. If data is not de-noised prior to compression, the largest noise peaks will typically be preserved, actually meaning that the noise in the signal can be emphasised.

The method for compression can start by eliminating points. Figure 10 shows three figures a)-c). Figure 10a) shows the signal to be compressed. Figure 10 b) has compressed the data 2:1 while figure 10 c) has compressed the data 6:1. The points prior to the largest changes are always kept and are marked with referral numbering 36 in figure 10 b). After removal of points, they are interpolated back to the original spacing using the compressed data set. These points shaped as ellipses. In figure 11 c) it can be seen that too many data points are removed. This means that the square difference between the interpolated subset and the de-noised subset will be high.

If the data set is de-noised, the square difference between the filtered data set and the compressed data set interpolated back to the filtered data set spacing will be below the noise level if the insignificant points are removed, and above if significant points are removed. This is shown in fig. 11. In the point (40), where the squared noise estimate crosses the removed energy from the signal, the optimal compression ratio can be found, as indicated in figure 11.

If too many small rate changes are present, the procedure above might fail, but the error introduced in these cases is believed to be small. It appears that this procedure in most cases will work satisfactory.

The procedure used for rate data assumes that the underlying signal is constant between transients, but the invention is not limited to this assumption, and the technique can easily be extended to other situations.

Claims

Patent claims:

1. Method for processing sampled data for the purpose of data transmittal and/or data storage, characterized in that the method comprises one or more of the following steps:

- perform real time discrete wavelet transformation of said sampled data using a moving wavelet window method,

- perform a real time noise estimation:

_ RMAD ^{σ = a} where RMAD, Robust Median value, is determined by determining the median of the wavelet coefficients at level one and removing all zero wavelet coefficients, and determine an optimal noise normalization value a through regression analysis, - perform outlier removal of the signal using a median filter, a window filter, in addition to an input signal of window size and said noise estimation σ ,

- filter the signal by removing and/or shrinking said wavelet coefficients by using said noise estimation and wavelet, and determining filter level, shrinkage rule and threshold value, and - performing data compression by determining an optimal compression ratio.

2. Method according to claim 1 , characterized in that said sampled data comprises production data, such as pressure, temperature and level data, flow rates and choke sizes.

3. Method according to claim 2, characterized in that for said pressure, temperature and level data, flow rates and choke sizes, a spline wavelet is selected, preferably a spline39 wavelet for said temperature and level data, and preferably a spline 4246 wavelet for said flow rated and choke sizes.

4. Method according to claim 1-3, characterized in that said optimal noise normalization value a preferably is in the range of 0,26 < a ≤ 0,78 , and most preferably a « 0.5164 for use with a data rate using spline4246 wavelet.

5. Method according to claim 1-3, characterized in that said optimal noise normalization value a preferably is in the range of 0,21 < a ≤ 0,62 , and most preferably a « 0.4193 for use with a data rate using spline39 wavelet.

6. Method according to one of the preceding claims, characterized in that an outlier is defined to be three times the noise estimate.

7. Method according to claim 6, characterized in that said filter level estimation is calculated by

l

where n is the number of data points.

8. Method according to claim 1 , characterized in that a SCAD shrinkage rule is used in connection with production data:

_gs

where λ is the threshold value, S is the signal, and a =3.7, based on a Bayesian argument.

9. Method according to one of the preceding claims, characterized in that an empirical threshold RED can be determined by

RED = a x σ x (b x log{n)f^veI

where a and b are determined by regression analysis.

10. Method according to claim 9, characterized in that 1.1 < a < 3.3 , and most preferably a » 2.2 , and preferably 0.19 < b < 0.57 , and most preferably 5 « 0,38 for use with data rate using spline39 and high RMAD.

11. Method according to claim 9, characterized in that 2.0 ≤ a ≤ 6.0 , and most preferably a « 4.0 , and preferably 0.13 ≤ b ≤ 0.39 , and most preferably b « 0,26 for use with data rate using spline4246 and medium RMAD.

12. Method according to one of the preceding claims, characterized in that the optimal compression ratio is determined by finding the point where the squared noise estimation crosses the removed energy from the signal.

13. Use of the method according to one of the preceding claims, for collection of production data within the petroleum industry, for improved data processing, such as automatic filtering and data compression, and speedy data preparation, for monitoring and storing said data.