GB2265710A

GB2265710A - Data compression of FTIR Spectra

Info

Publication number: GB2265710A
Application number: GB9206965A
Authority: GB
Inventors: Anthony Kevin Booer
Original assignee: Services Petroliers Schlumberger SA
Current assignee: Services Petroliers Schlumberger SA
Priority date: 1992-03-31
Filing date: 1992-03-31
Publication date: 1993-10-06
Anticipated expiration: 2012-03-31
Also published as: GB9206965D0; GB2265710B

Abstract

A method for processing IR spectra, e.g. of drilling fluids or cement slurries, to achieve data compression while retaining important features, comprises the following steps: generate spectrum, obtain autocorrelation sequence for spectrum, select a desired order for an auto regressive model, calculate auto regressive coefficient from the autocorrelation sequence and the selected order so as to construct an autoregressive model, calculate the poles of the autoregressive model, and solve the characteristic equation. <IMAGE>

Description

METHOD FOR PROCESSING INFRARED SPECTRAL DATA The present invention relates to a method for processing spectroscopic data. In particular, the invention relates to a method of processing Fourier Transform Infra Red spectral data so as to compress the data contained therein while retaining the important characteristic. However, it will be appreciated that the invention is not restricted to this alone.

The technique of FTIR analysis is well known. An absorption or transmission spectrum, typically in the wavenumber (frequency) domain, is generated from an interferogram produced by the analyser. Each spectrum might typically contain 2000 datapoints and so it is desirable to reduce the number of datapoints for a spectrum when a large number of spectra need to be analysed, especially if individual features must be identified. Peak picking procedures for data extraction have problems in estimating absorption in closely spaced peaks or dealing with noise.

It is an object of the present invention to provide a method for processing infrared spectral data so as to achieve data compression while retaining important spectral information.

The method according to the present invention comprises the following steps: 1 Generate spectrum 2 Obtain autocorrelation sequence for spectrum 3 Select a desired order (p) for an autoregressive model 4 Calculate auto regressive coefficient from the autocorrelation sequence and the selected order so as to construct an autoregressive model [n] 5 Calculate the poles (Zk)of the autoregressive model

wherein hk = Ak exp(j6 Zk = exp [(ak +j 2n: fk)?l and Ak = amplitude Ok = phase ak = decay rate fk = frequency = < T = sampling interval The present method is applicable to HIR spectra.The values of such spectra are not the original interference pattern but are the result of subsequent FT analysis.

However, the autocorrelation sequence for a spectrum can be obtained by using an inverse Fourier transform on this data. The autocorrelation sequence obtained can then be used to formulate an autoregressive (AR) model of order p. The order p of the model is selected such that it is greater than the number of peaks expected in the spectrum (since the exact number are not known and the model would not work if the order is less than the number of peaks). The spectrum will also contain noise which should also be modelled which also required a larger order than this expected number. In order to avoid problems with the model assigning more than one peak where only one is appropriate, it is desirable to ensure that the model order is typically less than twice the number of expected peaks.

The AR model can be used to express the spectrum H(eis0) as

and the coefficients a[M] (M = l...p) determined by calculating the solution of the AR Yule Walker normal equations

It is assumed, from the present invention, that the underlying mechanism which results in a measured FTIR spectrum is that of a discrete number of exponentially damped resonances.The model x for the Time domain data x[l]...x[Nl (sampling interval T) is

where the complex constants hk and Zk are defined as hk =Aexp(jO (jOk) Zk = exp [(αk + j2# fx)T] [V1 where amplitude Ak, phase ok. decay rate ak, frequency fk and real, j = < and ak < 0 so IZkl < l.

Equation [El can also be written as

The matrix of elements Zk is Vandermonde (t = Zi-1j) and if Zk are known, the equation is simply a set of simultaneous equations (in a least squares sense) for the complex amplitudes hk. It can be shown that Zk are the roots of the characteristic equation associated with a linear difference equation

for p+l < n < N where a[Ml are chosen to minimise

Once the coefficients a[m] have been determined from [11], the linear difference equation [VII1 can be used to construct a time domain representation x' of the AR spectrum thus

where #[n]= #e #[n]= 6e if n = 0 0 if n # 0 The values Zk are the poles of the AR model and can be calculated by finding the roots of the coefficient polynomial a[m]. The values calculated from Zk and x'[m] (substituted for x[n]) can be used in equation [V'] to calculate the values of hk.

The present invention will now be described in further detail with reference to the accompanying drawings, in which: - Figures 1 (a) and (b) comprise examples of ETIR spectra of a cement slurry and a dried mud powder respectively; - Figures 2(a) and (b) comprise the autoregressive model spectrum for the cement slurry and the residual error between the model and the raw spectrum; - Figures 3(a) and (b) comprise corresponding spectra and error for.a dried mud sample; - Figures 4(a) and (b) comprise wavenumber/amplitude and decay timel wavenumber parameters for the AR model of cement data; - Figures 5(a) and (b) comprise reduced parameters for Figures 4(a) and (b); - Figures 6(a) and (b) and 7(a) and (b) show corresponding parameters for dried mud data;; - Figure 8 shows a reconstructed cement spectrum; and - Figure 9 shows a reconstructed mud spectrum.

The invention finds particular use in the decomposition of FTIR spectra for quantitative analysis of oilfield cements and drilling muds. Such spectra are particularly complex but contain useful chemical and physical information. The automatic analysis of a number of spectra is desirable but problems may occur if the number of datapoints in each spectrum is large since the time taken for analysis will increase greatly. Matlab codes for FTIR spectrum analysis and spectrum synthesis are shown in Appendices A and B below. Typical examples of cement and dried mud FTlR spectra are shown in Figures l(a) and (b) and as can be seen, both show a large amount of structure, at least some of which is due to noise or interference. In each case, it is expected that the spectrum will include about 25 peaks as an AR model order of 50 has been chosen.

Other similar spectra might contain 10-20 significant peaks so that the mode order can be chosen accordingly. Once the model order has been selected hk and Zk can be calculated for each of the 50 parameters in accordance with the method described generally above. The amplitude Ak can be obtained from hk Ak = Ihkl [X] and the frequency fk and decay time Xk from Zk. The decay time xk is obtained from the decay exponent ak thus

A complex set of parameters (Ak, Fk, Tk) can be obtained for a given period model order p.

Figures 4 and 6 show plots in the frequency (wavenumber) domain of Ak and xk.

As the model order is 50, there are 50 lines on each plot. In order to reduce the number of lines on each plot, only 25 significant points are taken. These are graded according to the energy of each "peak" which is related to the product Ak Tk. The values of the reduced parameters are shown in the table below and plotted in Figures 5 and 7. Table 1: Estimated wavenumbers, amplitudes, and decay times for cement and dried mud spectra. Values in parentheses are wavenumbers of some known components.

cement Dried Mud wavenumber amplitude decay time wavenumber amplitude decay time fk Xkx104 Xk fk Xk x 104 Xk 90 132 26 439 93 11 131 100 24 506 106 11 262 96 58 572 99 10 338 86 25 634 106 11 411 73 22 705 96 12 519 83 20 783 103 13 613 368 24 1 850 68 10 675 593 31 918 111 12 742 627 24 994 126 10 809 448 17 1041 100 9 886 367 14 1110 (1107) 123 10 949 371 14 1171 102 10 1108 (1113) 272 18 1234 78 10 1632 (1635) 502 22 1445 (1445) 128 10 1678 274 16 1642 (1637) 146 12 3024 159 13 2854 77 11 3096 307 15 2928 (2927) 111 12 3165 630 19 3132 87 7 3219 426 16 3202 103 7 3304 556 26 3260 110 7 3387 585 22 3327 122 7 3451 389 14 3393 132 8 3510 396 15 3460 118 8 3573 269 16 3523 101 7 3639 160 20 3635 (3621) 179 13 The values obtained from the 25 peaks can be used to reconstruct spectra which are shown in Figures 2(a) and 3(a). The residual error between the reconstructed spectra and the measured spectra are shown in Figures 2(b) and 3(b) respectively and are obtained by subtraction of spectra.

Claims

A method of processing Fourier Transform infra-red (FTIR) data comprising: a) generating a FnR spectrum of a sample being analysed; b) obtaining an auto-correlation sequence for the spectrum; c) selecting a desired order (p) for an auto-regressive model; d) calculating an auto-regressive coefficient from the autocorrelation sequence and the selected order (p) so as to construct an auto-regressive model x[n]; e) calculating the poles (Zk) of the auto-regressive model; and

wherein hk = Ak exp(jEk) Zk = exp [(αk + j 2# fk)T] and Ak = amplitude ok = phase ak = decay rate fk = frequency = < T = sampling interval A method as claimed in claim 1, wherein the autocorrelation sequence is obtained by using an inverse Fourier Transform on the original spectral data.

A method as claimed in claim 1 or 2, wherein the order (p) is selected so as to be greater than the expected number of peaks in the spectrum but less than twice this number.

A method as claimed in any preceding claim for compositional analysis of a drilling fluid on cement slurry.

A Matlab code for FTIR spectrum analysis function (fr, tau, am, X, R] = ftir2(s,f, p) %FTIR2 % ftir2(spectrum, frequency, order) - analyze FTIR spectrum for peak frequencies, heights, and decay rates n = length(s); fmin = min(f); fmax = max(f); % acs = (ifft(s.*s)); % auto-correl'n EA,v) = levinson(acs(l:p+1)); % AR parameters A = 1 ; A]; x = filter(l,A,tsqrt(v);zeros(n - 1,1)]); % impulse response X = abs(fft(x)); % spectrum R = roots(A); % pole locations an = angle(R)/(2*pi); fr = (an + (an < O)) * (fmax-fmin) + fmin; % unvrap frequency tau = - 1 ./ log(abs(R)); % decrement % V = vandermonde(R,0:(2*p-1)); am = V \ x(1:2*p); % solve L2 problem B Mat lab code for spectrum synthesis A synthetic spectrum may be reconstructed from model parameters using the following code: function X = recon(rt,am, l,p); %RECON % % X = recon(rt,am,l,p) % % reconstruct spectrum from roots and amplitudes % 1 is length of xtn] sequence, p is spectrum length.

%.

X = vandermonde(rt,O:(p-1)); x = V * am; X = abs(fft((x ; zeros(l-p,1)]));