A new correlation method for transit-time estimation

c'~S~eac in ~ l e a r

Energy.

Vol. I, pp. 219 to 229.

Pergamon Press 1977.

Printed in Great Britain.

A NEW CORRELATION METHOD FOR TRANSIT-TIME ESTIMATION

Hideaki Nishihara and Hide<) Konishi* Department of Nuclear Engineering, Kyoto University Yoshida, Sakyo-ku, Kyoto, 606 Japan

ABSTRACT A new correlation method for transit-time estimation is proposed. The transit-time of a propagating quantity between two points along the propagation path could be measured by noise analysis technique when this quantity conveys some measurable stochastic signal with it. This is usually performed by observing directly the peak in the cross-correlation function of the stochastic signal measured at two points in the propagation path, or indirectly by reading the slope of the phase-angular frequency curve. The proposed method presents the transit-time information directly in the same way as in the ordinary crosscorrelation technique, but with a better resolution when the conveyed noise signal is bandlimited. It also rejects extraneous noise contamination when such noise appears in a narrow frequency band. The usefulness of this method was tested in an air-water two-phase flow channel to estimate the transit-time of the bubbles. INTRODUCTION The estimation of the velocity of a propagating physical quantity or of the location of a signal transmitter is an important technology in nuclear reactor systems as well as in many other engineering fields. For such purposes use is made of a propagating stochastic signal conveyed by the propagating quantity. Correlating the stochastic signal measured by a pair of detectors placed in the propagation path of such a signal to obtain the transit-time can provide information of the propagation velocity if the distance between the two detectors is known or the distance between the detectors (or more often the difference in the distances from the transmitter to the detectors) if the propagation velocity is known. The detected fluctuation signals are usually processed either in cross-correlation functions in the time domain or in phase angle diagrams of CPSD's (cross-power spectral densities) in the frequency domain. People have used one method or the other that was readily available to them. In many aspects of noise analyses, time- and phase-domain analyses are found interchangeable, and information obtainable in one method is also obtained in the other with the same amount of information retained in it. No serious thought has hitherto been given to the choice of the two noise analysis methods for the estimation of transit-times. In many cases one of the two sufficed the purpose. The cross-correlation function produces a delta function like peak at the transit-time when the fluctuation is nearly white and little dispersed as it propagates along the path. The resolution of the peak will be poor for band-limited stochastic signals, and when the dispersion is large or background noise masks the signal the determination of the transit-time peak will be difficult because of the blurred peak. The measured stochastic signals are disturbed internally, for example, by resonance in the case of triangulating to the location of anomaly by acoustic means (Ref. i) or by the global noise in the case of BWR steam bubble velocity measurement (Ref. 2). A typical example of extraneous background noise is the 60 Hz noise (or the like) from power-lines. In such cases elaborate filtering is needed to eliminate these unwanted components. The other technique is to plot the phase angle of the CPSD vs. angular frequency and to read the transit-time from the slope of the plot (Ref. 3). For a simple delay time T the plot 0 should be a straight line with a negative slope of -T 0. This method, however, does not produce a direct presentation of the required quantity, and involves a guess work in drawing the straight line. Moreover, when narrow-band noises contaminate the signal the phase angle

*Present address:

Nippon Atomic Industry Group Co., Ltd., Kawasaki-shi, Kanagawa-ken, Japan 219

220

H. Nishihara and H. Konishi

plot becomes d i s c o n t i n u o u s a t such frequencies [also the phase angle is computed in the principal value of (-~/2, ~/2) making the plot discontinuous] and the estimation will be difficult. A noble method---called the method of cepstrum originally developed for seismic studies (Ref. 4)---was tried in a transit-time analysis of acoustic noise for boiling detection (Ref. 5). In short, the cepstrum is a Fourier transform of the logarithm of the power spectral density (PSD) of the sum of the signals, one of which delaying the other by a transit-time. As will be discussed below, this method involves an approximation; also, mixing the signal prior to the processing loses information on the direction of the signal propagation. In this paper a new technique of correlating stochastic signals is proposed for transit-time estimation and the various methods are compared. To demonstrate the usefulness of the proposed method, an air-water two-phase flow loop was constructed to simulate a BWR channel and the bubble velocity was estimated by detecting the light transmitted through the two-phase. VARIOUS CORRELATION SCHEMES The Fourier transforms of stochastic variables x(t) and y(t) are expressed

F[x(t)]

=

in polar forms as

IF[x]l ej@x,

(1)

and j@ Fly(t)] = ]F[y]J e where 0

x

and 8

y

Y,

(2)

represent the phase of x and y, respectively. xy

The CPSD of x and y is

(~] = F*[x]-F[y]

IFex] I'IFey] j

e J (@Y-OX) ,

and the cross-correlation function is given as the inverse Fourier transform of Eq. ~xy(t) = F -l[¢xy (~)]

(3) (3). (4)

When y(t) is a repetition of x(t) with a delay time of % 0, i.e. y(t) = x(t-T O) ,

(5)

the Fourier transform becomes

-j~t o F[y]

=

F[x] e

(6)

Hence

~xy (~)

=

iF]X] j2

e

-J~0

,

(7)

of which inverse transform is Cxy (t) = %xx(t)*6(t-tO ) =

~xx (t-~O) •

(8) (9)

Thus the cross-correlation function will have a peak at t . Also the phase component of Eq. (3) becomes exp(-j~%0), and the phase-angular frequency w011 be a straight line with a slope of -t 0 . In the cepstrum technique the two signals are added together with an attenuation in intensity of the delayed signal by a factor e(
(iO)

A new correlation method for transit-time estimation

sz

(m) = #

where # (m) is the PSD of x. mate relationship holds:

xx

(m) (i + 2 a cos~T 0 + 2 ) ,

If we take the logarithm of Eq.

lOg#zz(m) = lOg%xx(m) + 2 e cosmT 0.

221

(ii)

(ll), the following approxi-

(12)

The desired information on r 0 is contained in the second term. If the first term is a slowly varying function of ~, it could be removed by appropriate lifterin@+, leaving the ripple component of the second term. The Fourier transform of the liftered log spectrum would be e[6(t-T0)+~(t+T0)], giving two peaks at ~ and -T 0. In the CPSD, Eq. (7), or in the cross-correlation function, Eq. (8), the transit-time information is contained only in the second term. The gain term of Eq. (7), which is the PSD of x, folds into the second term in the course of transformation to produce Eq. (9). If the first term of Eq. (7) is removed before the inverse transformation, it would give a delta function which has an unambiguous peak at the transit-time T 0. This would be equivalent to the cepstrum, but does not involve approximation. Inverse transformation Eq. (4) to obtain the crosscorrelation function can be considered as the process of transforming the needed simple delay time information by weighting it with the PSD which is not necessary in obtaining the transittime. The weighting also blurrs the peak. Since IF[x]l and IF[y]l are the square roots of PSD's of x and y, i.e., /~ (m) and /~ (~) respectively, Eq. (3) can be rewritten: xx yy J (8y-8 x) ~xy (~) = /~xx(m)~yy(~) We define the complex coherence function

r

(m) xy

F

xy

e

(13)

(m) by%%

xY(~)

(14)

/~xx (m) #yy (m)

Then J (@y-@x) r

xy

(w) = e

(15)

We denote the inverse Fourier transform of Fxy(m) by Yxy(T), Yxy(T) = F-l[Fxy(m) ] . When Eq.

(16)

(5) holds, we have Yxy(T) = 6(T-T0) ,

(17)

which produces a well defined peak at T 0. In practice the upper limit of the inverse transformation (which is not infinite but a finite fixed number) should be chosen carefully, since this correlation scheme rejects the weighting by the amplitude as in the case of calculating the cross-correlation function but carries out the integration with equal weight for all frequencies even at a high frequency where the signal to noise ratio is generally poor. Comparing the y-function and cepstrum, we observe that the two processes are similar but the former does not involve approximation nor needs a complicated liftering process.

%In the cepstrum terminology, paraphrased words are used for operations in the frequency domain that are generally performed in the time domain. The paraphrased words are made by the interchange of early consonant groups. Thus cepstrum is spectrum paraphrased, and lif__tering means a "filtering" operation on a frequency-dependent function as we normally do on a time-dependent function. %fCoherence function is usually defined as F (<~)2, a real function which does not contain the transit-time information, xy

222

H. Nishihara and H. Konishi

APPLICATION TO TWO-PHASE FLOW VELOCITY MEASUREMENT Experimental

set-up

2

I TEST SECTION 2 TUNGUSTEN LAMPS 3 PHOTOCOI~IVECELLS

-=I I¢::I =3

4 RESISTANCE PROBE 5 MIXING CHAMBER

6 7 8 9

II

SEPARATION TANK RESERVOIR PUMP ROTAMETER(WATER)

I0 AIR TANK II ROTAMETER~IR)

Fig. i.

Experimental

set-up.

> X-Y PLOTTER PREPROCESSI NG

I

i i

FAST FOURIER TRANSFORM i PSD I CPSD | , COHERENCE J.

~

X-Y PLOTTER

I INVERSE FAST FOURIERTRANSFORM i AUTOCORRELATION FUNCTION CROSS-CORRELATION FUNCTIONF ¥-FUNCTION ,J

Fig. 2.

Data processing

X-Y PLOTTER

scheme.

A new correlation method for transit-time estimation

223

A typical example of correlation technique for the transit-time measurement is found in BWR application to measure the steam velocity. TO prove the usefulness of the proposed new correlation method, or the method of y-function, an air-water two-phase flow loop was constructed. Figure 1 shows the schematic diagram of the experimental set-up. Air and water are mixed in the mixing chamber to produce bubble flow. The test section, made of transparent acrylic plate, had a square cross-section of 1 × 1 cm 2 and a length of 60 cm. The two-phase mixture is separated in the separation tank, and the water goes back to the reservoir. The flow rates of air and water are measured by rotameters, and a resistance probe placed at i0 cm upstream of the channel exit monitors the bubble velocity and void fraction. Two tungsten lamps are placed along one side of the test section, and two photoconductive cells are placed on the other side each facing the lamp. One pair of the lamp and the cell is set at the center of the test section and the other at a location chosen in each experiment. For the illustrative runs reported here it was set at 42 mm upstream of the channel mid-point.

0,12

0.i0

0.08

0.06

0,(2{I

0.012

0.00

-0.~ 0.0

0.i

0.2

0,3

0.4

0,5

0.6

0.7

0,8

0.9

1,0

DELAy T[MIE (SEC) (a)

0.48

0.40

0.32

*= 0,24

_~ 0.16 0.08 0.00 -0.08 0.0

I

I

i

0.1

0.2

0.3

I

1

1

0.4 0.5 0,6 DE~Y TI~ (S£C)

[

I

I

0.7

0.8

0.9

1.0

(b) Fig. 3.

Analysis of experimental data. (U) Cross-correlation function. (b) y-function. Water flow rate = 5 Z/min, void fraction = 5.67 ~.

224

H. N i s h i h a r a

and H. Konishi

The light from the lamps are s c a t t e r e d as it passes t h r o u g h the two-phase m i x t u r e and the f l u c t u a t i o n s of the light intensity are m e a s u r e d as a change in the r e s i s t i v i t y of the cells. E x t r a n e o u s noise can be added to the system either at the l i g h t - s o u r c e or at the cell ends to simulate v a r i o u s situations. The cells, H a m a m a t s u T V CdS cells P-440, have r e s i s t i v i t y of 55 k~ at I0 lux and 7.8 k~ at i00 lux. A c c o r d i n g to the manufacturer's specification the rise and the decay times are 20 and i0 msec at i0 lux, respectively. A resistor is c o n n e c t e d to each cell in series and a DC v o l t a g e of ±15 volts is a p p l i e d to read the f l u c t u a t i n g v o l t a g e t h r o u g h the cell. The f l u c t u a t i o n s are r e c o r d e d on a tape r e c o r d e r for s u b s e q u e n t processing. Data p r o c e s s i n g The m e a s u r e d signals are recorded r e p r o d u c e d signals are low-passed

on the tape recorder for a p p r o x i m a t e l y one minute. The through a n t i - a l i a s i n g filters set at 40 Hz and AD c o n v e r t -

0,20 0.16 -'\

O,12 -

",,

0.08 /

\

\ ~n

0,04

0.00

-0.04

-0,08 0,0

I 0.1

L

I

t

I

I

I

i

I

0,2

0,3

0.4

0,5

0.6

0,7

0,8

0,9

DELAY

TIME

1.0

(SEC)

(a)

0.20 0,16 0.12

=~ 0.08 0,04

o.oo

-0.04 --

-0.08

0,0

I

J

I

I

A

I

I

I

I

0,I

0,2

0,3

0.4

0,5

0,6

0,7

0.8

0,9

1,0

DELAY TIME (SEe)

(b) Fig.

4.

A n a l y s i s of e x p e r i m e n t a l data. (~) C r o s s - c o r r e l a t i o n function. (b) y-function. W a t e r flow rate = 3 ~/min, void f r a c t i o n = 9.49 ~.

A new correlation method for transit-time estimation

ed in 12 shown in duration magnetic scheme.

bits this time tape

225

on a small computer, FACOM U-200, with a sampling time of i0 msec for the cases paper. The data size was 1024 for each of the two channels. Thus the signal actually used for the analysis was about i0 sec. The digital data stored on a are fed to a larger computer, FACOM 230-75. Figure 2 shows the data processing

The raw data are removed of linear trends, normalized by the standard deviation and cosinetaper windowed in the preconditioning stage. They are then directly Fourier transformed by an FFT (fast Fourier transform) subroutine (Ref. 6), and analyzed in PSD's, CPSD, coherence function and frequency response function (both amplitude and phase). Fourier inverse transformation of these quantities gives the autocorrelation functions, cross-correlation function and the newly proposed y-function. The various computed functions are output in the form of figures by an X-Y plotter together with the original data.

0.32

0.24

0.16

j

~= 0.08

0,00

-0.08

-0.16

-0.241 0.0

0.I

I

I

I

I

I

I

l

l

0,2

0,3

0.4

0.5

0.6

0,7

0,8

0.9

1,0

DELAY TIME (SEE) (a)

0.40

0.32

0.24

=z 0.16

0.08

0,00

-0,08

-0.16 0.0

I

L

I

I

0,1

0,2

0,3

0.4

L 0.5

I

I

I

I

0.6

0.7

0.8

0.9

1.0

DELAY TIME (SEC) (b)

Fig. 5.

Analysis of experimental data. (a) Cross-correlation function. (b) y-function. Water flow rate = 4.4 //min, void fraction = 7.97 %. The light source is modulated by a sine wave.

226

H. Nishihara and H. Konishi

Experimental results Three groups of runs were performed. The first group is the experiments with no extraneous noise added to the system. This corresponds to the BWR experiments in which the global noise is eliminated and only the local noise is registered by the detectors. Figure 3 compares the cross-correlation function and y-function when the flow velocity is comparatively large (water flow rate = 5 ~/min, and void fraction = 5.67 ~). Both of the curves show sharp peaks at the transit-time. When the flow velocity is reduced, the fluctuation becomes slower or tends to contain more of the low frequency components. Figure 4 shows the results when the water flow rate was 3 //min, and the void fraction 9.49 ~. The cross-correlation peak becomes less sharp while the y-function retains a sharp peak. Next set of runs are the cases when sinusoidal perturbations are added to the power supplied

0,44

0,36

0,28

0,20

0.12

0,04

-0.04

-0.12 0.0

I

I

I

I

I

l

I

I

I

0.1

0.2

0.3

0.4

0,5

0.6

0.7

0,8

0.9

1,0

DELAYT]ME (SEC) (a)

0.24

0.20

0.16 F.-

0,12

~

0,08

0,04

I

0,00

-0,04 0.0

I

I

I

I

I

I

[

I

I

0.i

0,2

0.3

0.q

0.5

0.6

0.7

0.8

0,9

1,0

DELAY TIME (SEC) (b)

Fig. 6.

Analysis of experimental data. (a) Cross-correlation function. (b) y-function. Water flow rate = 2.4 //min, void fraction = 16.5 ~. The light source is modulated by a low-pass white noise, with a cutoff set at 1 Hz.

A new correlation method

for transit-time

estimation

227

to the lamps. The detected signals become more like sine waves on which perturbations induced by the b~bbles are superposed. Figure 5 gives an example when a sine wave of 5 Hz is added to a constant (non-zero) input voltage to the lamps. The power to the lamps will then contain two sinusoidal components, one at 5 Hz and the other at I0 Hz. The cross-correlation function is a composite of two cosine curves, while the y-function has a peak at the correct transit-time. The water flow rate was 4.4 Z/min and the void fraction 7.97 ~ for this run. Other runs were made for sinusoidal input of different frequencies. This set of experiments shows that when the bandwidth of the extraneous background noise is narrow the y-function peak is not affected by the background. This is understandable since the y-function is the Fourier inverse transform of the phase diagram; the phase diagram is disturbed by such noise only at discrete frequencies but the disturbance smears out through the integration in the course of inverse transformation.

0,64 0.48

-

0.32

0,16 0.00 -0.16 -0,32

~

-0,48

0.0

i

i

I

0.i

0.2

0.3

I

I

I

0.4 0.5 0.6 DELAYTIME (SEC)

I

I

I

0.7

0.8

0,9

1.0

(a) 0.24 0.20 I 0,16 --

z

0,12 -

0,08 -

0.~

-

0~ -0.04 ~ 0.0

0,1

I

I

0,2

0,3

1

I

]

0.4 0,5 0.6 DELAY TIME (5EC)

I

I

]

0,7

0,8

0,9

1,0

(b) Fig. 7.

Analysis of experimental data. (a) Cross-correlation function. (b) y-function. Water flow rate = 2.5 Z/min, void fraction = 15.8 ~. The light source is modulated by a low-pass white noise of a medium intensity.

228

H. Nishihara

and H. Konishi

Lastly, a low-pass white noise is added to the source. The transit-time information is contained in the the fluctuation of the bubbly flow which covers a fairly large frequency range, whereas the low-pass noise appears in a low frequency range with a zero time delay for the two detected signals. Thus two peaks are expected to occur in the cross-correlation function; one is the transit-time delay peak and the other the common mode peak with a zero time delay. The extraneously added noise can be considered as simulating the global noise encountered in BWR noise studies. Since the noise generator is built to repeat signals with an interval of 0.5 Hz, the generated signal was tape-recorded and played back at a speed 1/20 of the recording speed to eliminate this effect. Figure 6 shows the case when the cut-off frequency of the added noise was 1 Hz. AS expected, the cross-correlation function shows two peaks. The large low-frequency zerolag peak masks the other peak at the transit delay time. On the other hand, the y-function retains a sharp peak at the transit-time.

0.22

0.18

0.14 z

=

<

\ 0.10 \

0,06 -

\

\ \

0.02 -

-0,02 --

-0.06 0,0

I

I

I

I

I

L

I

I

I

0.1

0,2

0,3

0.4

0.5

0,6

0.7

0,8

0.9

0,8

0.9

1,0

DELAY TIME (SEC)

(a)

0,20 1 0,16

0.12

F

0,08

0.04

0.00

-0.04

-0.08 0,0

I

I

I

L

I

I

I

0.1

0,2

0,3

0,4

0,5

0.6

0.7

I 1.0

DELAY TIME (S£C)

(b)

Fig. 8.

Analysis of experimental data. (a) Cross-correlation function. (b) y-function. Water flow rate = 2.2 Z/min, void fraction = 17.6 ~. The light source is modulated by a low-pass white noise of a large intensity

A new correlation method for transit-time estimation

229

Figures 7 and 8 are the cases when the amplitude of the extraneous noise was increased. The cut-off frequency was the same, and the water flow rate and the void fraction were similar to the cases shown in Fig. 6. With the increase in the amplitude of the extraneous low-pass noise, the transit-time information becomes progressively buried in the cross-correlation function, while in the y-function the transit-time delay peak is still recognizable although the corm~on mode peak becomes larger at the zero delay time. CONCLUSIONS The proposed new method for estimating the transit-time of a propagating stochastic information was compared with the ordinary cross-correlation technique and was found more effective in producing a well defined time-delay peak; with a due attention to the choice of the frequency range in the inverse-transformation integral, the new method automatically produces the time-delay peak even when the measured signals are contaminated with a narrow-band limited extraneous noise. The relationship between the two methods often used for the transit-time estimation, i.e., the method of cross-correlation function in the time domain and that of CPSD phase angle diagram in the frequency domain, was made clear and found not interchangeable. The new method is essentially a variation of the latter method, but it presents the result directly in a time domain diagram. It is expected that the newly proposed method will find applications in many engineering problems, nuclear or otherwise. ACKNOWLEDGMENTS The use of digital computers at the Kyoto University Computation Center is acknowledged. also thank Mr. S. Asai for his help in performing the experiment.

We

REFERENCES (i)

Hideaki Nishihara, Resonant acoustic noise spectra of nucleate coolant boiling, J. Nucl. Sci. Technol. ii, 1 (1974).

(2)

See, for example, T. Nomura, et al., BWR noise spectra and application of noise analysis to FBR, Paper presented at S M O ~ - I (1974).

(3)

See, for example, D. Wach, The analysis of at-power neutron flux noise in the frequency range of vibrating reactor structures, Paper presented at SMORN-I (1974).

(3)

B. P. Bogert, M. J. R. Healy and J. W. Tukey, The quefrency alanysis of time series for echoes: cepstrum, psuedo-autocovariance, cross-cepstrum and saphe cracking, M. Rosenblatt, ed. (1963) Proceedings of the Symposium on Time Series Analysis, Chapter 15, Wiley, New York.

(5)

S. Ozaki and S. Iida, Application of cepstrum for the detection of location of boilingnoise source in reactor coolant, Paper C-3, 1975 Fall Meeting on Reactor Physics and Engineering, At. Energy Soc. of Japan (1975).

(6)

J. S. Bendat and A. G. P. Piersol dures, Wiley, New York.

(1971) Random Data Analysis and Measurement Proce-