Chapter 4 Implementation and Experimental Results
4.5 Seismic Applications on Real Seismic Data
The system is also applied to detect direct wave and reflection wave in real seismic data. We obtain data from Seismic Unix System developed by Colorado School of Mine [5].
The real data with the size 3100×48 showed in Fig. 36 (a) is from Canadian Artic, which has 48 traces and 3100 samples per trace with sampling interval 0.002 seconds.
The horizontal axis is the trace number and the vertical axis is time t.
After envelope and threshold preprocessing [7], Fig. 36 (b) shows the result of envelope and Fig. 37 (a) shows the result of thresholding with the threshold 0.15. The result of peak detection is in Fig. 37 (b). We only choose points with y < 700 which includes points from direct wave, first reflection wave and second reflection wave as in Fig. 38 (a) where there are 88 points. Detected curves are plotted in Fig. 38 (b) with
K = 3 and we preset f
1 = 0. Here, the initial parameters mk,x and mk,yare random
number between 0 and 50, ak= −1, b
k = 1, and fk = 0. The cooling function is as (9) with a high enough temperature, Tmax = 1,000. There are Nt = 100 trials in a temperature, and the temperature decreases 500 times. Constants settings: αm = 1, αab= 0. 5, and αf = 10.
Since the second reflection wave is not a hyperbola in theory [3], we remove the points from the second-layer reflection wave. That is, we remove the points nearest to the bottom pattern in Fig. 38 (b). Remaining 65 points are plotted in Fig. 39 (a). We redo the experiment, and this time K = 2. Fig. 39 (b) shows the result. The detected parameters in Fig. 38 (b) and Fig. 39 (b) are listed in Table VI and Table VII.
Table VI
Detected parameters in Fig. 38 (b) with fixed f1 = 0 in image space 3100×48
mx my a b f
Direct wave 24.48 8.59 -25.69 0.038 0 (preset)
Reflection wave 24.83 28.83 -22.91 0.044 2,441.7
Second Reflection wave 24.74 49.59 -23.10 0.043 8,942.7
47
Table VII
Detected parameters in Fig. 39 (b) with fixed f1 = 0 in image space 3100×48
mx my a b f
Direct wave 24.53 11.20 -25.59 0.039 0 (preset)
Reflection wave 24.62 -2.81 -24.13 0.041 2,978.05
As mentioned in [3]-[5], two lines of direct wave is a pair of asymptotes of the hyperbola and a pair of asymptotes is a hyperbola with size zero. From the detected parameters in Table VIII, we obtain the equation of the direct wave
)
Rearrange (20) and take square root on both sides, the equations of two lines in image space are
48
(a)
(b)
Fig. 36. Experiment on real data -- (a): Real seismic data from Canadian Artic. (b):
Plot of envelope.
49
(a)
(b)
Fig. 37. Experiment on real data -- (a): Threshold 0.15. (b): Detect peak.
50
(a)
(b)
Fig. 38. Experiment on real data -- (a): Choose peak with y < 700. (b): Detection result of (a).
51
(a)
(b)
Fig. 39. Experiment on real data -- (a): Remove points nearest to the bottom pattern.
(b): Detection result.
52
(a)
(b)
Fig. 40. Plot detected curve on the original data -- (a): Detection result from Fig. 38 (b). (b): Detection result from Fig. 39 (b).
53
The other real data is Gulf of Cadiz’s seismic data. There are 48 traces and 2050 samples in a trace with sampling interval 0.004 seconds. Fig. 41 (a) shows the real data. The horizontal axis is the trace number and the vertical axis is time t.
After envelope and threshold preprocessing [7], Fig. 41 (b) shows the envelope and Fig. 42 shows the thresholding result with threshold 0.5. The detected peak in Fig.
42 are plotted in Fig. 43 (a) where there are 66 points and the size of image is 2050×48, where the horizontal axis is x and the vertical axis is y. The initial parameters mk,x and mk,y
are random number between 0 and 50, a
k = −1, bk = 1, and fk= 0. The cooling function is as (9) with a high enough temperature, Tmax = 1,000.
There are Nt = 100 trials in a temperature. The temperature decreases 500 times.
Constants settings: αm = 1, αab = 0.5, and αf = 10. Number of patterns is K = 2. The detection result is in Fig. 43 (b).
Since the points nearest to the pattern around y = 800 are from second-layer reflection wave. In theory, the second-layer reflection wave is not a hyperbola [3]. We remove those points and remaining 48 points are plotted in Fig. 44 (a). Fig. 44 (b) shows the detection result and Fig. 45 plots the detected curve in the original data.
Table VIII and Table IX list the detected parameters in image space in Fig. 43 (b) and Fig. 44 (b).
Table VIII
Detected parameters in Fig. 43 (b) in image space 2050×48
mx my a b f
Reflection wave 44.57 187.84 -6.93 0.144 2,116.8
Second-layer reflection wave
21.64 56.27 -25.13 0.040 12,537.7
Table IX
Detected parameters in Fig. 44 (b) in image space 2050×48
mx my a b f
Reflection wave 45.58 174.02 -7.00 0.143 2,519.30
54
(a)
(b)
Fig. 41. Experiment on real data -- (a): Real seismic data from Gulf of Cadiz. (b): Plot of envelope.
55
(a)
Fig. 42. Experiment on real data – Thresholding result of the envelope in Fig. 41 (b) with threshold 0.5.
56
(a)
(b)
Fig. 43. Experiment on real data -- (a): Detected peak from Fig. 42. (b): Detection result of (a), K = 2.
57
(a)
(b)
Fig. 44. Experiment on real data -- (a): Remove points nearest to the pattern around y
= 800 in Fig. 43 (b). (b): Detection result of (a), K = 1.
58
Fig. 45. Plot the detected curve in Fig. 44 (b) on the original data.
59
Chapter 5
Conclusions and Discussions
5.1 Conclusions
This paper is about a proposed system, which adopted the simulated annealing algorithm to detect patterns such as lines, circles, ellipses, and hyperbolas by finding their parameters in an unsupervised manner and global minimum fitting error related to points in an image. The iterative adjustment requires less memory space. Also, we define the distance from a point to a pattern and this makes the computation feasible, especially for hyperbola. Using four steps to adjust parameters from center, shape, angle to the size of the pattern can get fast convergence. Based on the average minimum distance from points to patterns, we have proposed a method to determine the number of patterns automatically. Experimental results on the detection of line, circle, ellipse, and hyperbola in images are successful. The detection results of line pattern of direct wave and hyperbolic pattern of reflection wave in one-shot seismogram are good, and can improve seismic interpretations and further seismic data processing.
5.2 Discussions
Parameter settings. In the cooling schedule, the value of T
max, and Nt, are set prior.For a trial which includes a change of center, a change of b and a, a change of θ, and a change of f for every pattern, there are three possible results to accept or reject the change determined by Metropolis criterion:
1. The new parameter has smaller error and it is accepted.
2. The new parameter has larger error and it is still accepted.
3. The new parameter has larger error and it is rejected.
The determination of Tmax, we considered the accept ratio of the larger-error trials. If the Tmax is not high enough, the trial with larger error will almost reject, that is, it always accept trial with smaller error, so it is possible to reach local minimum. Fig. 46 shows this situation, where Tmax = 1,500 iterations, initial center (0,0), a = 1, b = 1, θ
60
= 0, and f = 1. Fig. 47 shows the result when Tmax = 10. In Fig. 47 (b), the accept ratio of canter, angle, and size has increased and shows the good result for this simple example. Fig. 48 show the result when Tmax = 100,000. After 500 iterations, T = 4.1, but this temperature is not low enough and the high accept ratio of larger-error parameters results in the instability. To solve this problem, we can increase the number of iterations to 1,000 and show the result in Fig. 49. This still shows good result, but it takes more time. In conclusions, for temperature Tmax, we have to choose a high enough temperature that gives a high accept ratio of larger-error parameters.
Besides, we need many enough iterations to cool the temperature to ensure stability.
Also we find the setting of Tmax is proportional to the scale of input points. In Fig.
47, Tmax = 10 provides good result. In Fig. 50, we enlarge the scale of data by two.
Fig. 50 (a) with Tmax = 10 cannot give good result, but Fig. 50 (a) with Tmax = 100 can give good result. In our simulation experiments, we choose Tmax = 500 and 500 iterations to ensure high enough initial temperature and lower final temperature T ≈ 0.02.
As for Nt, if trials are not many enough, we cannot get good result. Larger Nt takes more time but gives more chances. So we can have as many trials as possible if the computational power is strong enough. Fig. 51 shows too few trials cannot provide good result and for this simple example we need only Nt = 10 to obtain a good result.
61
(a)
(b)
Fig. 46. Illustration of low initial temperature Tmax: (a) Detection result of Tmax = 1. (b) Accept ratio of larger-error parameters.
62
(a)
(b)
Fig. 47. Illustration of low initial temperature Tmax: (a) Detection result of Tmax = 10.
(b) Accept ratio of larger-error parameters.
63
(a)
(b)
Fig. 48. Illustration of high initial temperature Tmax: (a) Detection result of Tmax = 100,000. (b) Accept ratio of larger-error parameters.
64
(a)
(b)
Fig. 49. Illustration of high initial temperature Tmax with more iterations: (a) Detection result of Tmax = 100,000. (b) Accept ratio of larger-error parameters.
65
(a)
(b)
Fig. 50. Enlarge the scale of points by two: (a) Detection result of Tmax = 10. (b) Detection result of Tmax = 100.
66
(a)
(b)
Fig. 51. Relationship between Nt and detection result Tmax = 500: (a) Detection result of Nt= 1. (b) Detection result of Nt = 10.
Time consumption. As for the time consumption, Table III shows that the CPU
time is proportional to the number of patterns or the number of parameters. The larger number of parameters, the algorithm takes more time to obtain the solution.Memory requirement. For traditional HT, it needs an accumulation matrix. The
size of accumulation matrix grows as the number of parameters increases. Besides,67
the higher precision, the larger accumulator matrix is needed. On the other hand, SA algorithm for parameter detection needs only memories for the original parameters and the trials parameters. This depends on the number of patterns K. Furthermore, the parameters can be presented by SA algorithm with high precision since we do not need to quantize the parameter space as in the traditional HT.
Preprocessing. In seismic application, we have no constraint on the center.
However, for ideal case, the hyperbola has the center on x-axis, i.e. t = 0. In simulated seismic data, we can find that the center does not lie on the x-axis, since convolution produces a shift. So preprocessing is quite critical. Wavelet and deconvolution processing may be needed in the preprocessing to improve the detection result.
68
References
[1] P. V. C. Hough, “Method and means for recognizing complex patterns,” U.S.
Patent 3069654, 1962.
[2] R. O. Duda and P. E. Hart, “Use of Hough transform to detect lines and curves in pictures,” Communication ACM, vol. 15, 1972, pp. 11-15.
[3] M. M. Slotnick, Lseeons in Seismic Computing. The Society of Exploration Geophysicists, 1959.
[4] M. B. Dobrin and C. H. Savit, Introduction to Geophysical Prospecting.
McGraw-Hill Book Co., New York, 1988.
[5] O. Yilmaz, Seismic data processing. The Society of Exploration Geophysicists, 1987.
[6] P. Kearey and M. Brooks, An Introduction to Geophysical Exploration. Blackwall Science Publications, 1991.
[7] K. Y. Huang, K. S. Fu, T. H. Sheen, and S. W. Cheng, “Image processing of seismograms: (A) Hough transformation for the detection of seismic patterns. (B) Thinning processing in the seismogram,” Pattern Recognition, vol. 18, no.6, pp.
429-440, 1985.
[8] J. Basak and A. Das, “Hough transform networks: learning conoidal structures in a connectionist framework,” IEEE transaction on Neural Networks, vol. 13, 2002, pp. 381-392.
[9] K. Y. Huang, J. D. You, K. J. Chen, H. L. Lai, and A. J. Don, “Hough Transform neural network for seismic pattern detection,” International Joint Conference on
Neural Networks, July 16-21, 2006, pp. 4670-4675.
[10] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, “Optimization by simulated annealing,” Science, vol. 220, 1983, pp. 671-680.
[11] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, “Equation of state calculations by fast computing machines,” J Chem Phys., 1953, pp.1087-1092.
[12] R. W. Klein and R. C. Dubes, “Experiments in projection and clustering by simulated annealing,” Pattern Recognition, vol. 22, 1989, pp. 213-220.
[13] D. E. Brown and C. L. Huntley, “A practical application of simulated annealing to clustering,” Pattern Recognition, vol. 25, 1992, pp. 401-412.
[14] K. S. Al-Sultan, S. Z. Selim, “A global algorithm for the fuzzy clustering problem,” Pattern Recognition, vol. 26, 1993, pp. 1357-1361.
[15] S. Geman and D. Geman, “Stochastic relaxation, Gibbs distribution and Bayesian restoration of images,” IEEE Trans. Pattern Analysis and Machine Intelligence,
69
vol. 6, 1984, pp. 721-741.
[16] H. Szu, R. Hartley “Fast simulated annealing,” Physics Letters A, vol. 122, 1987, pp. 157-162.
[17] A. Corana, M. Marchesi, C. Martini, and S. Ridella, “Minimizing multimodal functions of continuous variables with the simulated annealing algorithm,”
Mathematical Software, ACM Transactions on, vol. 13, 1987, pp. 262-280.
[18] G. P. Barabino, G. S. Barabino, B. Bianco, and M. Marchesi, “A study on the performance of simplex methods for function minimization,” Circuits and
Computers IEEE International Conference on, ICCC 80, 1980, pp. 1150-1153.
[19] S. F. Masri, G. A. Bekey, F. B. Safford, “Global Optimization Algorithm Using Adaptive Random Search,” appl. Math. and Comput. vol. 7, 1980, pp. 353-376.