An Effective Algorithm to Correct the Orientation Error and Time Shift of Downhole Seismograms

Estimating the orientation error of a downhole seismometer, in general, is a two-parameter optimization problem. It requires determining both the time shift and orientation error at same time. Although this estimation in­ volves a tedious computation, the resolution of this optimization is still lim­ ited by the grid size selected for global searching. To improve the comput­ ing efficiency and resolution, we derive an analytic formulation for calcu­ lating the maximum cross-correlation function between synthetic and ro­ tated seismograms and apply it to an orientation-error-correction method developed by Chiu and Huang (2003). This modified method is capable of giving an accurate and efficient estimation of the orientation error in one step for a given time shift. This algorithm consists of two stages: (1) com­ puting the synthetic seismograms at the downhole station, and (2) search­ ing a time shift that yields a maximum cross-correlation between the ro­ tated seismograms (observations) and the synthetic seismograms at the downhole station. Results show that this modified algorithm allows for a more flexible selection of data for analysis, and gives a fast, stable and ac­ curate estimation of the orientation error. (


INTRODUCTION
For a downhole seismometer (velocity-type seismometer or accelerometer) without a built in compass, orientation errors must be corrected from recorded data after its first installation and upon any later reinstallation.In general, this is a two-parameter optimization problem.It requires a global search for the time shift due to unsynchronized recording, and the orientation error of the downhole seismometers due to the uncontrolled installation.Numerical optimiza tion can be done by grid search; the possible range of these two parameters (time shift and 1 Institute of Earth Sciences, Academia Sinica,Taipei, Taiwan, AOC 2 1nstitute of Seismology, National Chung-Cheng University, Chia-Yi, Taiwan, ROG orientation error) can be divided into discrete values with a small interval and the best solution can be found by selecting a pair of time shift and orientation error that gives the highest cross correlation coefficient.For example, the orientation error could be any value between 0 and 360° and the possible time shift can be any value between -2 and 2 s, if we assign the incre ment of the orientation error at 1 degree and time accuracy at 0.0 1 s, then we will have 1.44 X 105 grid points for numerical optimization with a resolution of one degree in orienta tion error and O.Ols in time shift.
Most correction algorithms based on the similarity of waveforms between the surface and the downhole waveforms (e.g., Seale and Archuletta 1989; Aster and Shearer 1991; Chiu et al. 1994).However, this similarity is low and limited in the narrow window of a seismogram.Chiu and Huang (2003) suggested selecting the SH wave on the free surface as the reference direction and to search the maximum cross-correlation coefficient between the downhole syn thetic and observed SH waves in various orientations.The corrected angle corresponding to the maximum cross-correlation coefficient then can be considered as the orientation error of a downhole seismometer.Their method includes more data points for analyses and thus reduces the estimation uncertainty; however, all these methods have common difficulties in increasing the computing efficiency and resolution.A rise in resolution requires reducing the grid size on both of the time shift and orientation error, thus causing a tremendous increase in the compu tation time.Although the CPU time is not a big issue, a reliable estimation in both the time shift and the orientation error is still a tedious process and its resolution is limit by the grid interval selected for both the time shift and the orientation error.Using an elegant optimization method such as the Genetic Algorithm (e.g., Goldberg 1989) can increase the searching effi ciency but it cannot improve the resolution.Developing an effective and robust algorithm that can reduce the computation time and increase the resolution is necessary.In this study, we make use of the advantages of the algorithm developed by Chiu and Huang (2003) and im prove its orientation error estimation to reduce the computing time and increase the resolution.

METHOD
The steps involved in the present algorithm are two-fold: first, computing the synthetic seismograms at a downhole station, and secondly, searching for orientation and time shift that give a maximum cross-correlation between the rotated seismograms (observations) and the synthetic SH waves at the downhole station.
The basic assumption behind this algorithm is that the synthetic (derived from the surface recording) and the corrected downhole seismograms (observation) have a maximum cross correlation coefficient (hereafter CCC) when both of their SH waves are in the same direction and the recordings are chronologically synchronized.Ideally, these two waveforms should be the same; however, deviations occur in an imperfect velocity model and an imperfect separa tion of the P, SH and SV waves.Nevertheless, compared with a surface recording, this syn thetic seismogram has much higher similarity with downhole recordings.
The reasons for selecting the SH waves as the reference were given by Chiu and Huang (2003 ), and one numerical method that can calculate the downhole synthetic seismogram in a layered medium using the surface recording as an input was given by (Huang and Chiu 1996).Here we will emphasize deriving a one-step algorithm for estimating the orientation error.To give a detailed description of this method, consider that a set of downhole seismograms in the where summation is taken over the discrete data within a window.
Let (3) The R x's in Equation (2) can be either a maximum or a minimum at 80 • To determine it, need to compute the second derivative of R , at e0• Le t A:::.)�2si , Equations (4) implies that Rx's is only a minimum at the e0 (we name it as em in hereafter) as given in Equation (3).Fortunately, we can determine the em ax from emin because the differ ence between these two 8s is always 180°.This relation can be obtained from equation (2) by substituting e by 8+180.After some manipulations we have Rx' s (8+180)=-Rx' / 8). (5) This equation further implies that the maximum and the minimum of Rx's always have an 180° difference in its location in 8space, and we can determine (}ma x from a given (}min .The relationship between e max and e mi n depends on the definition of e.If we define 0:::; e < 360° then we have e max = e min + 180° when e max = e min -180 ° when 0 :::; ; e min < 180°.180 :::; ; e min < 360°. (6a) (6b) In the following section, we will outline three tests to validate this method step by step.

NUMERICAL TESTS
In the first test, we consider a noise-free experiment that assumes the synthetic and the observed waveforms at downhole are identical when orientation error has been corrected.The waveform selected for testing and the results of this test are shown in Fig. 1.Two orthogonal components of displacement waveforms x and y are considered as the downhole seismograms as shown in (a) and (b) in Fig. 1.If the downhole seismometer has an orientation error of 37.3°, the synthetic seismogram can be obtained by making a counter-clockwise rotation on wave forms x and y by a 8=37.3°using Equation (1).The resulting waveform is shown in (c).Based on the waveforms (a), (b) and (c), we can estimate the e m in by applying Equation(3) which gives 217.3°.This emin corresponding to a emax =37.3° and the corrected waveform as shown in (d).The CCC for various orientations is given at the bottom of Fig. 1.This result gives an exact estimation in the orientation error ( Bmax =37.3° and CCC=l at 8 m a x ).The CCC plot also indicates that e max and emin have a difference of 180°.
The shape of the CCC depends upon the seismograms x, y and s.In the second test, we select another set of waveforms [(a) and (b) in Fig. 2] for testing and let the orientation error be 37.621°.In addition, we let the synthetic seismogram have a time shift of 1.54 s.The final seismogram is shown in Fig. 2(c).This waveform set has a richer signal than that of the dis placement waveform in the previous test.Furthermore, we can increase the decimal shift of the orientation error to 3 digits.In order to search for the optimized orientation error and time shift, we set the possible time shift in a range between -2.5 and 2.5 s with an increment of 0. 005s.We should emphasize here that this grid point increment is the same with the sampling interval of these of waveforms.For each time shift, we can calculate its maximum CCC using equations (3) and ( 5).Comparing maximum CCC among all the possible time shifts, we found that the global maximum CCC is 1, which occurred at 8max=37.62 1° and t= -1.54 s.Since the positive time shift is defined as an advance in the time shift, the negative time shift repre sents the time delay of the observed seismogram relative to the synthetic seismograms or the surface recording.The corrected waveform is shown in (d), which is identical to the wave form in (c).Again, in this noise-free case, we can obtain an exact solution for both the time shift and orientation error.
In the previous two tests, we calculated Bmax indirectly from Bm in using Equation (5).To show the detail of the CCC and the relationship between e max and e m in, we now calculate the i�tc ) ��-3 � I 0 10 20 30 40 50 60   CCC over a range of the time shift between -2.5 and 2.5 s and orientation error in a range between 0° and 360° and present the result in a 3D plot as shown in Fig. 3.As expected, this two-parameter optimization has several local maxima and minima.The global maximum and minimum of these CCC are marked by solid dots in the figure and their corresponding time shift, orientation error and CCC are labeled beside the dots.For this noise-free case, the maximum CCC is 1 while the value for the minimum CCC is -1.These two extremes have the same tsri while the 8 has a difference of 180° .This numerical example also shows that @max and emin' have a difference of 180°.Equation ( 3) is an exact solution in numerical form.It doesn't have any limitation on resolution.However, the numerical calculation may introduce an estimation error.The pos sible errors would be caused by the numerical error due to the rounding calculations of a computer, data noise, and the error caused by the grid point selected for finding the time shift.The rounding error of a computer can be ignored because it is relatively small.The amount of the data noise is case by case and will not be discussed here.In the next test, we will discuss the effect of the grid size of the time shift on the estimation of the orientation error and time shift.In this test, we still use the same waveforms as those in the Fig. 2. If the grid point coincides with the data sampling points, we always can obtain an exact time shift and orienta tion error and the corresponding CCC is 1 (as shown in Fig. 2).In the following test, we will examine the estimation error when the grid points do not fall on the sampling points.To design

Cross-Correlation Coefficient
•: : .•::. . ., �: :c�_ '. _ : .��: Again, the worst case is still at the center of two sampling points.For this case, the maximum error in orientation estimation is 0.03°.Of course, the estimation errors shown in this test will change as the data or B changes.Overall, the estimation error due to the selecting the grid point is proportional to the time difference between the sampling point and grid point.This error is small and a decrease in grid interval will cause a decrease in the estimation error.
The data selected for the last test is from the downhole data set recoded at the Dahan Downhole Array, Hu alien, Tai wan.This earthquake occurred on 21 September 1999 in the southwest of the Array.The local magnitude of this earthquake was 5.4 and its Epicenter was about 30 km from the array.Before proceeding farther, a window must be selected for the data f o r analysis.To select a proper window, both the nonlinear effect and the effect of the imper fect separation of SH waves from a three-component recording need to be considered.The non-linear effects only occur when ground shaking is strong while the imperfect separation of SH waves is always found in the later arrivals (the later portion of a waveform).Since the non linear behavior was not found in this data set, the selection of a 12-s window is mainly for avoiding in imperfect separation of the SH waves.The original seismogram of this data was recorded on acceleration; however, we selected the velocity waveforms for estimating the orientation errors because it is better than an acceleration or displacement waveform for such study (Chiu and Huang 2003).This data set was evaluated by Chiu and Huang (2003) and proved to be one of the best recordings suitable for estimating the orientation errors.The orientation errors of these downhole accelerometers at depths of 50, 100 and 200 m estimated by an independent approach (Chiu and Huang 2003) were 40°, 114° and 285°, respectively.In this case, all downhole recordings were recorded in the same time system and no time shift correction was necessary.The B mi n estimated using equation (3) were 219.451°, 289.7310° and 104.8816° which corresponding to Bmax were 39.459 1°, 113.7310° and 284.8816°.The comparisons of the synthetic and the corrected waveforms are shown in Fig. 5. Three pairs of waveforms from the top to the bottom are synthetic-observed SH waveform pairs at depths of 50, 100 and 200 m.It is clear these waveform pairs are very similar.
Further comparisons of these velocity waveforms before and after orientation correction are shown in Fig. 6.The amplitudes of these waveforms are plotted on the same scale and the vertical distance between two traces is proportional to the depth of the recording accelerometers.Before the correction (top), no consistent phases (arrivals) could be found among these waveforms.After the correction, all the phases (bottom) displayed a very good alignment after the orientation correction.From the corrected waveforms, very clear up-going and down going primary S waves were obvious.Both of these results suggest that the orientation errors of these downhole accelerometers have been corrected accurately.

DISCUSSION AND CONCLUSIONS
This algorithm has its disadvantages and advantages when compared with other methods.The disadvantage is that we must know the velocity model in advance.In general, this doesn't become a serious problem because all the downhole array sites must have a velocity well logging before installation.On the other hand, the main advantage of this algorithm is in its   ------'------'-------'------'------'------' 0 2 4 6 8 10 12 Qi Synthetic > -0.4 1-------------�" "--� �.tl"-...  A good match in phase and the underestimation of the amplitude as shown in Figs. 5 and 6 imply that the average velocity model is sufficient.But, a further improvement on waveform matching requires a refined velocity model.
In most cases, the shallow structure can be a good approximation for a horizontal layered medium.For such cases, we can apply the Haskell-Thomson matrix to calculate the downhole seismograms.This study did not discuss the case of a non-horizontal layered medium.However, if the synthetic seismograms for such a velocity model can be obtained, we can extend this method for use in a non-horizontal layered medium.
SH-SV plan has two perpendicular components x,. and Y r If x; and Y ; are rotated counter-clock w, ise by an angle e to correct the orientation error, then we have two new components x � and Yi• i.e.

Fig. 1 .
Fig. 1.Displacement waveforms in (a) and (b) are x and y components in the Cartesian coordinate system.(c) is the projection of (a) and (b) in a direction of new x axis (x') after a 37.3° counter-clockwise rotation.In the Test Model 1, waveforms in (a) and (b) are considered as observation, while (c) is considered as the noise free synthetic seismogram at the downhole station.The corrected waveform is given in (d).The cross-correlation coefficient between various rotated compo nents of the observation (x in equation (2)) and the synthetic seismogram s is shown in the bottom of the Figure.A dotted line marks the orientation error of 37 .3°and a star marks the maximum cross-correlation coefficient.

Fig. 2 .
Fig. 2. The selected waveforms and results for the Test Model 2. The format of this figure is similar to that of the Fig.I except that the waveform in (c) has a time shift 1.535 s and the cross-correlation coefficients are plotted for various time shifts.

Fig. 3 .Fig. 4 .
Fig. 3.The 3D plot of the cross-correlation coefficients based on the waveforms in Test Model 2. The range of the time shift is between -2.5 and 2.5 s and the orientation range is between 0° and 360°.Two dots in this figure mark the global maximum and minimum.Three numbers in the labels next to the dots are the corresponding time shift, orientation error and cross correlation coefficient respectively.

Fig. 5 .
Fig. 5. Comparisons between the corrected and synthetic waveforms at the depths of 50 (top), 100 (middle) and 200 m (bottom).The synthetic SH waves are calculated base on the ID layered model.

Fig. 6 .
Fig. 6.The uncorrected (top) and the corrected (bottom) velocity waveforms for orientation error.The amplitude of these SH waves are plotted on the same scale and the distance between two traces is proportional to the depth of the accelerometers.