论文部分内容阅读
Abstract As 3D high-precision seismic exploration is more and more widely used in seismic data acquisition， traveltime tomographic inversion based on first arrivals is developed from 2D to 3D. However， magnanimity data of 3D traveltime inversion brings about the problem of data storage； the absence of first arrivals with near offset reduces the precision of shallow layer； the utilization of prior information， such as small refraction and micrologging data， can improve the precision of 3D traveltime inversion. Therefore， we make some improvements in 3D traveltime inversion method. We take compression storage for large and sparse matrix， propose virtual receivers technology， and add prior information to tomographic inversion linear equations. The application in 3D real data indicates that the improvements can effectively improve 3D traveltime tomographic inversion.
Key words：3D seismic exploration； 3D traveltime inversion method； 3D traveltime tomographic inversion
INTRODUCTION
The incompleteness of data projection angles in geophysical field gets tomography method not to precisely image on the velocity distribution of subsurface stratigrapher. But compared with other imaging methods， tomography still is a good method of studying subsurface physical properties change and geological structure. In order to apply traveltime tomography well， people have already studied tomography algorithms， such as ray tracing method and solution method of tomographic inversion linear equations， and make some achievements. However， particularly for 3D traveltime tomographic inversion， the coefficient matrix of tomographic inversion linear equations needs a lot of memory space， first arrivals with near offset in data acquisition are inadequate and even absent， and the utilization of nearsurface velocity survey data can reduce inversion’s multisolution. Regardless of the problems above， the original tomographic inversion algorithm can not obtain a good inversion result. Therefore， we attempt to adopt different methods to solve the problems respectively.
1. THEORY AND THE METHOD
1.1 Basic Principle
The mathematical theory of traveltime tomographic inversion is Radon Transform. Under the hypothesis of linear condition， according to the traveltime equation
This paper uses SIRT to resolve the linear equations above. Because traveltime tomography is a very complex problem， various constraint conditions need to be added in the process of tomographic inversion to ensure the convergence and stability of inversion process. The conventional constraint conditions mainly include the minimum and maximum constraint of velocity， correction quantity constraint of slowness， inner iteration number of times， velocity extrapolation， velocity smoothness and so on. Based on the conventional constraint conditions， we further make some improvements. 1.2 The Improvements of Tomographic Inversion 1.2.1 The Compression Storage and Solution of Large and Sparse Matrix
The paper adopts row compression storage technology to storage large and sparse matrix. We use three vectors which are row， col and val to represent the large and sparse matrix. Vector row stores total number of nonzero elements of all rows before some row， its data type is integer and its dimension equals row number of the matrix plus one. Vector col stores colum index of all non-zero elements in original matrix， its data type is integer and its dimension equals the number of all non-zero elements. Vector val stores value of all non-zero elements in original matrix， its data type is float or double and its dimension equals the number of all non-zero elements.
SIRT algorithm doesn’t involve operation of an element in the matrix， but product of between a matrix or its transpose and a vector. So realization of the product of between compression storage matrix or its transpose and a vector can easily resolve compression storage linear equations.
The product of compression storage matrix and a vector is followed below
Where the row number of original matrix is m， the column number is n， the length of vector s is n， vector t is the product of compression storage matrix and vector and its length is m， row（i） is the number of all non-zero elements before row number i， row（i+1） is the number of all non-zero elements before row number i+1， val（j） and col（j） are the value and column index of jth non-zero element respectively.
The product between the transpose of compression storage matrix and a vector is followed below
W}here the row number of original matrix is m， the column number is n， λ（j） is the product between jth column of original matrix and vector s， the length of vector s is m， row（i） is the number of all non-zero elements before row number i， row（i+1） is the number of all non-zero elements before row number i+1， val（j） and col（j） are the value and column index of jth non-zero element respectively.
1.2.2 The Virtual Receiver Technology
In order to enhance the accuracy of shallow layer in near-surface velocity model， we add one or more virtual receivers between the shot point and the nearest receiver away from it. We can calculate the traveltime from shot point to virtual receiver according to average velocity obtained by near-surface survey and the distance between the shot point and virtual receiver. Then the calculated traveltime of virtual receivers and picked first break time of actual receivers are used to inverse near-surface velocity. So this can increase the number of ray in the shallow layer and further can get the velocity of more partition units in the shallow layer to be corrected. The detailed principle is followed as
Where point （X1，0） and point （X2，0） are actual receivers’ position， point （x1，0） and point （x2，0） are virtual receivers’ position. The travel time of virtual receivers（x1，0） and （x2，0） is calculated by t1=x1/v1 and t2=x2/ v2respectively， where velocity v1 and v2 are obtained by near-surface survey data.
We use an ideal model with inclined fault to test the effectiveness of this method. The size of the model is 240m×100m×60m， and the size of partitioned unit is 2m×2m×1m. The velocity of the layer above inclined fault and under one is 500 m/s and 1300 m/s respectively. The number of shot points is 20， and the number of seismic traces is 7200. Initial model’s velocity is gradient， and the number of iteration time in the whole inversion process is 20. From the comparison of Figure 5 and Figure 7， we can see that adding virtual receivers in the inversion process gets the shape of inclined fault to be more close to the ideal model， in the meantime the velocity of the layer above inclined fault also gets more close to ideal velocity value.
1.2.3 Prior Information Constraint Tomographic
Inversion According to traveltime equation （1）， we can obtain the inversion linear equations expressed by M seismic traces and N unknown numbers， that is When target function reaches minimum， ？ζ（？s）=0 exists， equation （9） is solved， and the computed ？s is optimal solution. It can be seen from the target function that the constraint conditions don’t work if β=0， in addition， the bigger the β is， the stronger the constraint is.
We use an ideal model to test the effectiveness of this method. The ideal model has an undulating surface and a lens inside it. The size of the model is 240m×100m×60m， the size of partitioned unit is 2m×2m×1m， the number of shot points is 20， the number of seismic traces is 7200， and four wells for micro-logging are added in it. The four wells’ position is （11，23）， （41，31）， （71，39）， （101，47） respectively and their depth is 30m， 31m， 29m， 27m respectively. Initial model’s velocity is gradient， and the number of iteration time in the whole inversion process is 20.
From the comparison of Figure 11 and Figure 13， we can see that adding prior microlog information constraint in the inversion process gets the shape of lens to become more clear， in the meantime the velocity of lens also gets close to ideal velocity value.
2. EXAMPLES
In order to test the application effectiveness of the improvement methods， we use 3D real seismic data of some exploration area in Shengli Oilfield. The exploration area is flat and its surface strucure is divided into three layers which are low velocity layer， weathering layer and high velocity layer. Initial velocity model is gradient and its size is 1400m×7225m×30m， iteration time is 20. The interpretation result of small refraction and micrologging data shows that the thickness of low velocity layer and weathering velocity layer ranges from 4m to 12m. From Figure 15 and Figure 17， we can see that the thickness of both low velocity layer and weathering velocity one obtained by improved inversion method is more close to known interpretation result than original method due to the control of virtual receivers. At the same time， velocity value obtained by using improved inversion method is also more accurate because of the constraint of small fraction and micro-logging data and the control of virtual receivers.
CONCLUSIONS
The row compression storage technology can save a lot of memory space， the virtual receiver technology can effectively control the inversed depth of shallow layer if low velocity layer and weathering velocity one are more shallow， and adding prior information constraint can enhance the inversion accuracy. Ideal model test and real data test both show that the improvement of 3D traveltime tomographic inversion method can obtain a better inversion result than original one.
REFERENCES
[1] Williamson， P. R. （1990）. Tomographic Inversion in Reflection Seismology. geophys. J. int.， 100， 255-274.
[2] Toshiki， W.， et al. （1995）. Improvement in Accuracy of Seismic Traveltime Tomography by Adding a Priori Information as Constraints on Velocity Reconstruction. J. Soc. Mat. Sci.， 44， 874-879.
[3] Chen， B. F.， Xiong， D. Y.， Ren， X. Q.， & Cheng， C. H. （2008）. The First-Arrival Tomographic Inversion and Its Application to Identify Thick Near-Surface Structures. SEg Expanded Abstracts， 27， 3255-3259.
[4] Cheng， G.， Zhang， B. J. （2008）. Compression Storage and Solution of Large and Sparse Matrix in Traveltime Tomography of Reflection Seismic Data. Progress in geophysics， 23（3）， 674-680.
Key words：3D seismic exploration； 3D traveltime inversion method； 3D traveltime tomographic inversion
INTRODUCTION
The incompleteness of data projection angles in geophysical field gets tomography method not to precisely image on the velocity distribution of subsurface stratigrapher. But compared with other imaging methods， tomography still is a good method of studying subsurface physical properties change and geological structure. In order to apply traveltime tomography well， people have already studied tomography algorithms， such as ray tracing method and solution method of tomographic inversion linear equations， and make some achievements. However， particularly for 3D traveltime tomographic inversion， the coefficient matrix of tomographic inversion linear equations needs a lot of memory space， first arrivals with near offset in data acquisition are inadequate and even absent， and the utilization of nearsurface velocity survey data can reduce inversion’s multisolution. Regardless of the problems above， the original tomographic inversion algorithm can not obtain a good inversion result. Therefore， we attempt to adopt different methods to solve the problems respectively.
1. THEORY AND THE METHOD
1.1 Basic Principle
The mathematical theory of traveltime tomographic inversion is Radon Transform. Under the hypothesis of linear condition， according to the traveltime equation
This paper uses SIRT to resolve the linear equations above. Because traveltime tomography is a very complex problem， various constraint conditions need to be added in the process of tomographic inversion to ensure the convergence and stability of inversion process. The conventional constraint conditions mainly include the minimum and maximum constraint of velocity， correction quantity constraint of slowness， inner iteration number of times， velocity extrapolation， velocity smoothness and so on. Based on the conventional constraint conditions， we further make some improvements. 1.2 The Improvements of Tomographic Inversion 1.2.1 The Compression Storage and Solution of Large and Sparse Matrix
The paper adopts row compression storage technology to storage large and sparse matrix. We use three vectors which are row， col and val to represent the large and sparse matrix. Vector row stores total number of nonzero elements of all rows before some row， its data type is integer and its dimension equals row number of the matrix plus one. Vector col stores colum index of all non-zero elements in original matrix， its data type is integer and its dimension equals the number of all non-zero elements. Vector val stores value of all non-zero elements in original matrix， its data type is float or double and its dimension equals the number of all non-zero elements.
SIRT algorithm doesn’t involve operation of an element in the matrix， but product of between a matrix or its transpose and a vector. So realization of the product of between compression storage matrix or its transpose and a vector can easily resolve compression storage linear equations.
The product of compression storage matrix and a vector is followed below
Where the row number of original matrix is m， the column number is n， the length of vector s is n， vector t is the product of compression storage matrix and vector and its length is m， row（i） is the number of all non-zero elements before row number i， row（i+1） is the number of all non-zero elements before row number i+1， val（j） and col（j） are the value and column index of jth non-zero element respectively.
The product between the transpose of compression storage matrix and a vector is followed below
W}here the row number of original matrix is m， the column number is n， λ（j） is the product between jth column of original matrix and vector s， the length of vector s is m， row（i） is the number of all non-zero elements before row number i， row（i+1） is the number of all non-zero elements before row number i+1， val（j） and col（j） are the value and column index of jth non-zero element respectively.
1.2.2 The Virtual Receiver Technology
In order to enhance the accuracy of shallow layer in near-surface velocity model， we add one or more virtual receivers between the shot point and the nearest receiver away from it. We can calculate the traveltime from shot point to virtual receiver according to average velocity obtained by near-surface survey and the distance between the shot point and virtual receiver. Then the calculated traveltime of virtual receivers and picked first break time of actual receivers are used to inverse near-surface velocity. So this can increase the number of ray in the shallow layer and further can get the velocity of more partition units in the shallow layer to be corrected. The detailed principle is followed as
Where point （X1，0） and point （X2，0） are actual receivers’ position， point （x1，0） and point （x2，0） are virtual receivers’ position. The travel time of virtual receivers（x1，0） and （x2，0） is calculated by t1=x1/v1 and t2=x2/ v2respectively， where velocity v1 and v2 are obtained by near-surface survey data.
We use an ideal model with inclined fault to test the effectiveness of this method. The size of the model is 240m×100m×60m， and the size of partitioned unit is 2m×2m×1m. The velocity of the layer above inclined fault and under one is 500 m/s and 1300 m/s respectively. The number of shot points is 20， and the number of seismic traces is 7200. Initial model’s velocity is gradient， and the number of iteration time in the whole inversion process is 20. From the comparison of Figure 5 and Figure 7， we can see that adding virtual receivers in the inversion process gets the shape of inclined fault to be more close to the ideal model， in the meantime the velocity of the layer above inclined fault also gets more close to ideal velocity value.
1.2.3 Prior Information Constraint Tomographic
Inversion According to traveltime equation （1）， we can obtain the inversion linear equations expressed by M seismic traces and N unknown numbers， that is When target function reaches minimum， ？ζ（？s）=0 exists， equation （9） is solved， and the computed ？s is optimal solution. It can be seen from the target function that the constraint conditions don’t work if β=0， in addition， the bigger the β is， the stronger the constraint is.
We use an ideal model to test the effectiveness of this method. The ideal model has an undulating surface and a lens inside it. The size of the model is 240m×100m×60m， the size of partitioned unit is 2m×2m×1m， the number of shot points is 20， the number of seismic traces is 7200， and four wells for micro-logging are added in it. The four wells’ position is （11，23）， （41，31）， （71，39）， （101，47） respectively and their depth is 30m， 31m， 29m， 27m respectively. Initial model’s velocity is gradient， and the number of iteration time in the whole inversion process is 20.
From the comparison of Figure 11 and Figure 13， we can see that adding prior microlog information constraint in the inversion process gets the shape of lens to become more clear， in the meantime the velocity of lens also gets close to ideal velocity value.
2. EXAMPLES
In order to test the application effectiveness of the improvement methods， we use 3D real seismic data of some exploration area in Shengli Oilfield. The exploration area is flat and its surface strucure is divided into three layers which are low velocity layer， weathering layer and high velocity layer. Initial velocity model is gradient and its size is 1400m×7225m×30m， iteration time is 20. The interpretation result of small refraction and micrologging data shows that the thickness of low velocity layer and weathering velocity layer ranges from 4m to 12m. From Figure 15 and Figure 17， we can see that the thickness of both low velocity layer and weathering velocity one obtained by improved inversion method is more close to known interpretation result than original method due to the control of virtual receivers. At the same time， velocity value obtained by using improved inversion method is also more accurate because of the constraint of small fraction and micro-logging data and the control of virtual receivers.
CONCLUSIONS
The row compression storage technology can save a lot of memory space， the virtual receiver technology can effectively control the inversed depth of shallow layer if low velocity layer and weathering velocity one are more shallow， and adding prior information constraint can enhance the inversion accuracy. Ideal model test and real data test both show that the improvement of 3D traveltime tomographic inversion method can obtain a better inversion result than original one.
REFERENCES
[1] Williamson， P. R. （1990）. Tomographic Inversion in Reflection Seismology. geophys. J. int.， 100， 255-274.
[2] Toshiki， W.， et al. （1995）. Improvement in Accuracy of Seismic Traveltime Tomography by Adding a Priori Information as Constraints on Velocity Reconstruction. J. Soc. Mat. Sci.， 44， 874-879.
[3] Chen， B. F.， Xiong， D. Y.， Ren， X. Q.， & Cheng， C. H. （2008）. The First-Arrival Tomographic Inversion and Its Application to Identify Thick Near-Surface Structures. SEg Expanded Abstracts， 27， 3255-3259.
[4] Cheng， G.， Zhang， B. J. （2008）. Compression Storage and Solution of Large and Sparse Matrix in Traveltime Tomography of Reflection Seismic Data. Progress in geophysics， 23（3）， 674-680.