Improved Correlation Identification of Subsurface Using All Phase FFT Algorithm

Abstract

The correlation identification of the subsurface is a novel electrical prospecting method which could suppress stochastic noise. This method is increasingly being utilized by geophysicists. It achieves the frequency response of the underground media through division of the cross spectrum of the input & output signal and the auto spectrum of the input signal. This is subject to the spectral leakage when the cross spectrum and the auto spectrum are computed from cross correlation and autocorrelation function by Discrete Fourier Transformation (DFT, "To obtain an accurate frequency response of the earth system, we propose an improved correlation identification method which uses all phase Fast Fourier Transform (APFFT) to acquire the cross spectrum and the auto spectrum. Simulation and engineering application results show that compared to existing correlation identification algorithm the new approach demonstrates more precise frequency response, especially the phase response of the system under identification.

1.Introduction

The correlation identification is widely applied in many fields, such as automation, biomedicine and so on. For instance, a fast visual evoked potential measurement technique-system identification with pseudo-random sequences is widely used in clinical diagnosis and research of central nervous system [1]. It uses Pseudo Random Binary Signal (PRBS) as input signals and obtains the impulse response and frequency response of the system under identification [2]. PRBS is often used as identification signals for the system, because it has properties similar to white noise, can guarantee a good identification accuracy and easy to achieve in the project. The electrical prospecting is also a process of system identification where the subsurface is an unknown linear time-invariant system. The correlation identification method is a universal system identification method, which can be applied to any active source electrical prospecting methods. Besides the complex resistivity, we can get DC resistivity, IP and the electromagnetic parameters [3]. Meanwhile, due to the low degree of correlation between PRBS and random noise signals[4], random signals are widely used in the field of communications, to suppress echo effects in digital video broadcasting [5] [6] [7]. Correlation identification combines this feature with system identification to suppress random noise and obtain accurate system identification results. Duncan [8] and Quincy [9] applied correlation identification in an electromagnetic sounding system, and Li [3] and Prasad [10] used it to obtain the complex resistivity spectrum. This is a type of spectral induced polarization (SIP, “Zhao [11] designed and manufactured an instrument based on correlation identification into practice, but the reduced correlation identification effect inhibits its popularization and application.  

Fourier transform is one of the most important fundamental and mathematical tools . According to different types of original signals, the Fourier transform is divided into four types: Discrete aperiodic Fourier series, Discrete periodic Discrete Fourier transform, Continuous aperiodic Fourier integral and Continuous periodic discrete-time Fourier transform[12], which are represented in the following table, Table 1. FFT is an efficient algorithm for the DFT, called fast Fourier transform. When an accurate single frequency signal is used, the FFT is remarkably stable [13]. Actually, when dealing with complex signals, accurate amplitude and phase response of the subsurface always cannot be obtained because discrete Fourier transformation (DFT) brings about side lobe effect. This significantly reduces the correlation identification effect. 

Table 1. Four classes of Fourier transform 

All phase fast Fourier transform (APFFT) is an effective spectrum analysis algorithm. The APFFT reduces spectral leakage and acquires precise amplitude and phase spectrum by modifying the signal truncating method of traditional DFT. It is widely used in discrete spectrum correction and filtering but in this paper we apply it in correlation identification which expands its application area.  

The proposed method uses APFFT in correlation identification instead of fast Fourier transform (FFT) to obtain the cross spectrum and the auto spectrum. The new method demonstrates obtaining an increased accuracy in the amplitude and phase spectrum of the underground media. The method is valid for theoretical and real SIP data. 

2. Methods

The model of correlation identification of subsurface is illustrated in Fig. 1, in which the impulse response of the subsurface to be identified h(t), the input PRBS signal v(t), the output signal r(t), the noise n(t), and the final output voltage signal acquired by the receiver is u(t).  

Fig. 1. The model of correlation identification 

The impulse response of the subsurface (referring to the underground medium between the field source and the measuring point) can be deduced from the cross correlation function of the transmitted signal and the received signal. Correlation identification method is often dealt with in frequency domain but it is also impacted by the spectral leakage. The APFFT can mitigate the problem and improves the practicability. 

2.1 Correlation Identification Method 

Correlation identification method uses PRBS as the input signal. After achieving the results of correlation operations between the input and output signals, the frequency response as well as the impulse response can be finally derived [14]. It is an effective way to eliminate random noise which takes advantage of the low correlation between the PRBS and the random noise signal. This is also the core algorithm of the article. That is the reason why we use correlation identification method to calculate system response. 

The relationship between the input and output signals is given by the following equation: 

u (t) r=(t)+n(t) v=(t)*h(t)+n(t)        (1)

Correlation operation is carried out with the input signal v(t)on both sides of the equation (1):  

v(t)⊗u (t)=v(t)⊗v(t)*h(t) v(t)+⊗n(t)        (2)

The noise n(t), which refers to random noise here, v (t)⊗ n(t)can be omitted because of the independence between the PRBS signals and the noise. Rewrite the formula (2) after using R_vu(t) and R_vv(t) to denote the cross correlation and autocorrelation respectively.

R_vu(t)=R_vv(t)*h(t)       (3)

We can easily get h(t)since the autocorrelation of the input signal is a Dirac delta function as formula (4), where c is a constant.  

We can easily get h(t)since the autocorrelation of the input signal is a Dirac delta function as formula (4), where c is a constant.  

R_vu(t)=R_vv(t)*h(t)==cδ (t)*h (t)= ch(t)       (4)

In frequency domain we can easily get the frequency response H (ω)of the system by division of the cross spectrum and auto spectrum. 

H(ω)=R_vu(ω)/R_vv(ω)       (5)

We compute the cross spectrum and auto spectrum via DFT which needs to truncate the cross correlation and autocorrelation, but truncating produces spectral leakage, causing frequency, amplitude and phase to be inaccurate. Therefore, it is incapable to get precise H (ω ) , and we cannot achieve the desired result of correlation identification.  

2.2 APFFT Algorithm 

FFT algorithm will result in the frequency leakage because of truncating and signal discontinuity after periodization. The proposal of APFFT algorithm is fundamentally proposed to solve the truncation error caused by the truncation of signal data as well as to diminish the discontinuity by preprocessing. To further illustrate the advantages of APFFT, we will compare the performance of FFT and APFFT.  

We assume that there is a single-frequency complex exponential sequence: 

x(n)=e^{j(2Πf₀nΔt+φ₀)}       (6)

where f₀ is the signal frequency and ω₀ is the phase. The N points FFT of x(n) is as formula (7): 

\(X(k)=\frac{1}{N} \sum_{n=0}^{N-1} x(n) e^{-j \frac{2 \pi k n}{N}}=\frac{1}{N} \frac{\sin \left(\pi\left(k-k_{0}\right)\right)}{\sin \left(\pi\left(k-k_{0}\right) / N\right)} \cdot e^{j \varphi_{0}} \cdot e^{j \frac{(N-1)}{N}\left(k-k_{0}\right) \pi}, k \in[0, N-1]\)       (7)

The amplitude spectrum of FFT is 

The phase spectrum of FFT is 

The APFFT algorithm consists of an all phase preprocessing unit which constructs a new signal sequence by a continuously cyclic shift of the original signal sequence and sums up at the same position, and then performs FFT on the new sequence. The preprocessing procedure is illustrated as Fig. 2. 

Fig. 2. The diagram of all phase preprocessing unit 

For baseband signal x(n), sampling 2N-1 points in one cycle to form a time sequence. For a certain sample point x(N) in the time series, there are N vectors of N dimension including this sample point: 

x₀=[x(N) , x(N+1), ···, x(2N-1]^T

x1=[x(N-1) , x(N), ···, x(2N-2]^T

 ·····

x_N-1=[x(1) , x(2), ···, x(N)]^T

Extend all vectors periodically, then truncate every vector with an N-length rectangular window and make the x(N) be the first element of every truncated vector, 

y₀=[x(N) , x(N+1), ···, x(2N-1]^T

y₁=[x(N) , x(N+1), ···, x(N-1]^T

 ·····

Y_N-1=[x(N) , x(1), ···, x(N-1]^T

We can get all phase data vector by adding the same position elements of all vectors

Suppose Y_i(k) and X_i(k)are the FFT of y_i (i=0,1 ···,N-1) and x_i (i=0,1 ···,N-1), respectively, the relationship of Y_i(k) and X_i(k) according to the shift property of FFT [15] are shown as formula (8)

\(Y_{i}(k)=X_{i}(k) e^{j \frac{2 \pi k i}{N}}\)       (8)

All phase spectrum is made up of the sum of Y_i(k) which is as formula (9).

\(\begin{aligned} X_{A P}(k) &=\frac{1}{N} \sum_{i=0}^{N-1} Y_{i}(k) \\ &=\frac{1}{N} \sum_{i=0}^{N-1} X_{i}(k) e^{j \frac{2 \pi k i}{N}} \\ &=\frac{1}{N^{2}} \sum_{i=0}^{N-1} \sum_{n=0}^{N-1} x(n-i) \cdot e^{-j \frac{2 \pi}{N} k n} \cdot e^{j \frac{2 \pi k i}{N}} \\ &=\frac{1}{N^{2}} e^{\varphi_{0}} \frac{\sin ^{2}\left(\pi\left(m-k_{0}\right)\right)}{\sin ^{2}\left(\pi\left(m-k_{0}\right) / N\right)} \end{aligned}\)       (9)

The m is the normalized frequency offset. The amplitude spectrum of APFFT is        \(\frac{1}{N^{2}} \cdot\left(\frac{\sin \left(\pi\left(m-k_{0}\right)\right)}{\sin \left(\pi\left(m-k_{0}\right) / N\right)}\right)^{2}\). The phase spectrum of APFFT is ϕ₀ .

Comparing equation (7) and (9), the result of equation (9) is actually the square of that of equation (7, “This means the change will seem more apparent from the APFFT than the normal FFT. The amplitude of X_AP(k) underlines that APFFT algorithm can reduce spectral leakage because it suppresses the side lobe and the major lobe is relatively reinforced.

Another crucial characteristic of APFFT algorithm is that its phase is constant and isn’t influenced by the frequency shift. Therefore the phase spectrum of the APFFT is independent of the frequency which means the phase estimation can be directly obtained regardless of any changes of frequency [16]. As the traditional FFT is not free of frequency influence, APFFT gets almost zero error in the presence of frequency offset. In that way, the real phase of signal can be acquired by APFFT algorithm when the signal is non-integer truncated. In conclusion, it is quite suitable to use the APFFT to realize the conversion between the frequency and time domain in the correlation identification algorithm.

2.3 APFFT based Correlation Identification Method

Based on the above principle, we propose an APFFT correlation identification method to obtain more precise identification result, which is represented as Fig. 3

Fig. 3. The block diagram of the improved correlation identification method

The specific process of improved correlation identification using APFFT algorithm is as follows:

(1) Define the discrete sampling input as v(n) and output signal as u(n).

(2) Do correlation, get the autocorrelation R_vv(n) of v(n) as well as the cross correlation R_vu(n )  of v (n) and u(n) .

(3) Using the principle described in Fig. 2. Truncated R_vv(n) into N vectors of N dimension, and every vector includes the same point R_vv(N)  . Do the same to R_vu(n) and R_vu(N). Take R_vv(n) as an example. Fig. 4 shows the truncated vectors after aligning the same points.

Fig. 4. The truncated vectors of R_vv(n)

(4) For all truncated subvectors, R_vv(n) , cyclic shift at their truncated position, add the same position values of each vector, then select its main value interval and let R_vv(N) be the first position of the sequences. Consequently, all phase data vector R_vv-AP(n)  (length is N) is obtained. R_vu-AP(n) will also be obtained through the same method. Fig. 5 shows the process of cyclic shift and summary of R_vv-AP(n) 

Fig. 5. Cyclic shift and summary of R_vv-AP(n) 

(5) Suppose that R_vv-AP(k)  is the DFT of R_vv-AP(n)  , R_vu-AP(n)  should be R_vu-AP(n)  . Divide R_vu-AP(k) by R_vv-AP(k) and H_AP(k) can be obtained, which is the frequency response of subsurface, All are shown as formula (17, “

\(H_{\mathrm{AP}}(k)=\frac{R_{u u-A P}(k)}{R_{w-A P}(k)}=\frac{\sum_{n=0}^{N-1} e^{-i \frac{2 \pi}{N} n k} R_{w-A P}(n)}{\sum_{n=0}^{N-1} e^{-i \frac{2 \pi}{N} n k} R_{v u-A P}(n)}, k=0,1, \cdots, N-1\)       (17)

(6) We can obtain the amplitude and phase response from the frequency response.

3. Simulation and Comparison

We will verify by simulation that the improved correlation identification is more accurate compared to traditional correlation identification.

The Cole-Cole relaxation model illustrated in Fig. 6 is widely used to fit experimental SIP data which simulates the subsurface.

Fig. 6. Cole-Cole relaxation model

The expression of the Cole-Cole relaxation model can be written as

\(Z(\omega)=R_{0}\left[1-m\left(1-\frac{1}{1+(i \omega \tau)^{c}}\right)\right]\)       (18)

In this formula, R0 is Direct Current (DC) resistivity, m the chargeability, τ the time constant, and c the frequency dependence [17].

The parameters in the formula can be typical assigned as Z(0)=1000,m=0.333, r=0.01, c=1 which are the values of the Cole-Cole relaxation model. In this way the impulse response of the Cole-Cole relaxation model is

\(h(t)=\frac{1}{4} \delta(t)+\frac{250}{8} e^{-375 t}\)       (19)

The typical values of the input PRBS are N=8192,q=5,a=1, under which better identification results can be produced. In the formula, N means the cycle length, q the number of periods and a the amplitude of PRBS. The sampling rate is 1 kHz. A Hanning window is used.

\(w(n)=0.5\left(1-\cos \left(2 \pi \frac{n}{M}\right)\right), 0 \leq n \leq M\)       (20)

Fig. 7 and Fig. 8 respectively show the amplitude and phase spectrum comparison between correlation identification and proposed correlation identification.

Fig. 7. The comparison of amplitude spectrum

Fig. 8. The comparison of phase spectrum

According to Fig. 7, both traditional and proposed correlation identification methods obtain accurate amplitude spectrum, but the traditional identification result in Fig. 7(b) is obviously less accurate. Fig. 8 also shows that the phase spectrum of proposed correlation identification is more accurate than that of traditional correlation identification. Both amplitude and phase spectrum of proposed correlation identification methods can rebuild a very similar spectrum.

Related to the expected noise under field conditions, we can add a 10V random Gaussian noise signal with an average of 0, a standard deviation of 1 to the output and signal-to-noise ratio is 0. The comparison of amplitude and phase spectrum shows in Fig. 9 and Fig. 10.

Fig. 9(b) shows that amplitude spectrum of traditional correlation identification has some large peaks much higher than the amplitude of the system, definitely affected by 10V random noise therefore it fails to identify the amplitude of unknown system. Fig. 9(c) is the result of proposed correlation identification method. Although the amplitude spectrum that improved method obtains is also affected by noise and indeed there is some fluctuation at each frequency point compared with Fig. 9(a), its error is smaller but unavoidable.

Fig. 9. The comparison of amplitude spectrum (with 10V random noise)

In addition, we can see that there are several large peaks in amplitude spectrum Fig. 7(b) and 9(b, “The reason is that when FFT(R_vv(t)) is approximate to 0, the denominator in the formula \(H(w)=\frac{F F T\left(R_{v u}(t)\right)}{F F T\left(R_{w}(t)\right)}\),should be 0. When doing division at this kind of frequency, there will be an abnormal phenomenon showing a high peak.

Fig. 10 represents the phase that two methods identify. We find that phase is sensitive to noise. Fig. 10(b), the phase result of traditional method, can be seen that there are some unknown phase jumps at some frequency points affected by noise. This is very different from the theoretical value. After simulating multiple times, we found the results shown are different each time. This is because the noise is random and the phase is very sensitive to noise, making the result of traditional method unstable and different each time. The phase that improved method gives shows as Fig. 10(c) is much better than Fig. 10(b) and does not change with different random noise. The improved method better identifies the phase of unknown systems

Fig. 10. The comparison of phase spectrum (with 10V random noise)

Additional comparative experiments were performed changing the parameters of input signal and calculating the errors respectively. Below are equations, in which Q₀ and Q₁ represent values of error.

\(e(k)=\hat{g}(k)-g(k)\)       (21)

\(Q_{0}=10 \log \left\{\left[\sum_{k=0}^{N-1} e^{2}(k)\right] /\left[\sum_{k=0}^{N-1} g^{2}(k)\right]\right\}\)       (22)

\(Q_{1}=10 \log \left\{\left[\sum_{k=0}^{N-1} e^{2}(k) \cdot g^{2}(k)\right]^{\top}\left[\sum_{k=0}^{N-1} g^{2}(k)\right]\right\}\)       (23)

g(k)is theoretical impulse response value of Cole-Cole relaxation model, while g(k) is actual impulse response value obtained by identification. The smaller Q₀ and Q₁ resulted in more accurate identification results

N, q , a and f_s in two tables are variables in PRBS. By controlling variables in this way, Table 2 indicates the identification result of the improved method is significantly more accurate especially when the number of input signal cycles is greater than one. Table 3 shows that when cycle length N increases, the result of the presented method continues to improve. Therefore, we can conclude that the present algorithm is able to acquire more accurate amplitude and phase spectrum of unknown system.

Table 2. The error of correlation identification and improved method when q is variable​​​​​​​

Table 3. The error of correlation identification and improved method when N is variable​​​​​​​

4.Application

In addition to mathematical formula derivation and simulation we also apply the improved method to field experiment. In the experiment there are sending signals and receiving signals obtained from field experiments. The receiving electrode distance is 0.2 m, the sending electrode distance is 52 m. The number of cycles of sending PRBS signal is 35 and every cycle has 1024 points. The length of receiving signal is the same as the sending signal. The sampling rate of v(n) and u(n) is 64Hz.

At Sifang Mountain we arranged experimental profile of frequency-domain induced polarization (IP) method at about 95m to the east of ZF0-1. Instrument model is 85321-pip, which is the frequency domain exciter developed by Champion Geophysical Technology Ltd. The profile length is 2000m, the dot pitch is 20m, use 25 four-channel receivers end to end and arrange 100 parallel measurement channels at a time. Power was supplied from the middle of the measuring line and set 34 power-supply polar distances totally. We made 100 parallel measurement channels to measure all the power supplies, and collected the signal of sending current at the same time.

After we compared the results of conventional correlation identification method with that of improved correlation identification method. Fig. 11 and Fig. 12 were made, which respectively represent amplitude and phase spectrum of identifying using traditional and improved correlation identification method.

Fig. 11. Amplitude of field experiments using traditional and proposed correlation identification method

Fig. 12. Phase of field experiments using traditional and proposed correlation identification method​​​​​​​

In terms of the amplitude, we can find that there are many abnormal values in Fig. 11(a), which is the identifying result of traditional method, while the curve of the improved method is smoother in Fig. 11(b, “It suppresses the distinct peaks of amplitude, which are close to theoretical values.

The range of phase jump amplitude using traditional method is -20 to 45 degrees, while using the improved method is -20 to 10 degrees. Therefore, the curve of Fig. 12(b) is stabler and better. In summary, the improved correlation identification method can obtain a more accurate phase spectrum when the number of the input signal is greater than one.

The identifying result of field experiments proves that improved correlation identification method will acquire stable and accurate amplitude and phase spectrum of unknown systems, when the number of input signal is unequal to one.

5. Conclusion

An improved correlation identification method, which combines correlation identification and APFFT algorithm is presented in this paper. Correlation identification can be used to obtain the spectrum of unknown systems, but its precision is not satisfied especially when the input signal has many cycles. Therefore, we can use APFFT algorithm to improve correlation identification method. The APFFT algorithm can be able to reduce spectrum leakage and obtain accurate phase value. In this paper, a series of simulation images and data results are used for error analysis. The result of the improved algorithm is compared with the results of the unmodified algorithm with the assumption system model. This is demonstrated from both image and data. Both have obtained a stable and satisfactory result. Meanwhile, the results of field experiments are also in line with expectations. Therefore, it is effective and practical to use APFFT to improve the correlation identification algorithm. Simulation and application prove that the improved method is able to acquire a more precise amplitude and phase spectrum. This could be applied to SIP, magnetotelluric method (MT) and other active prospecting methods.

However there exists one weakness that the improved method is affected by very strong noise and the identification process may fail in that case.