CN100383559C - Fixed-point decomposing method for three-dimensional seismic offset imaging - Google Patents

Fixed-point decomposing method for three-dimensional seismic offset imaging Download PDF

Info

Publication number
CN100383559C
CN100383559C CNB2006100472116A CN200610047211A CN100383559C CN 100383559 C CN100383559 C CN 100383559C CN B2006100472116 A CNB2006100472116 A CN B2006100472116A CN 200610047211 A CN200610047211 A CN 200610047211A CN 100383559 C CN100383559 C CN 100383559C
Authority
CN
China
Prior art keywords
matrix
imaging
diagonal
zero entry
fixed
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CNB2006100472116A
Other languages
Chinese (zh)
Other versions
CN1888935A (en
Inventor
石殿祥
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Individual
Original Assignee
Individual
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Individual filed Critical Individual
Priority to CNB2006100472116A priority Critical patent/CN100383559C/en
Publication of CN1888935A publication Critical patent/CN1888935A/en
Application granted granted Critical
Publication of CN100383559C publication Critical patent/CN100383559C/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

This is a research problem about the earthquake excursion imaging in prospecting physical geography field. It calculates the data of earthquake reflect wave by the handy arithmetic and puts into the disposal soft flat roof of earthquake data to form the imaging and confirm the best position for well or other underground resource. It adopts a new divide approaching method to achieve incomplete Cholesky decompose (because the coefficient matrix of the Cholesky decompose is denseness that can't be used in accounting). It prescribes some element is nonzero only a few case and calculate these nonzero cell according to the least error principle to complete the approach when the coefficient matrix is extraordinary rarefaction. It solves the problem about chaos phenomenon and badness boundary without avoiding in existing technology and avoids the direction aberration to assure the precision of excursion imaging for without dividing on direction.

Description

A kind of fixed-point decomposing method of three-dimensional seismic offset imaging
Technical field
The invention belongs to the problem that the seismic migration imaging aspect is studied in the exploration geophysics field, seismic exploration data is to be entered into magnetic medium and to form image in the seismic data processing software platform by field instrumentation equipment, the present invention proposes a kind of simple algorithm the earthquake reflected wave data is carried out migration imaging processing calculating, thereby is to understand the practicable reliable method that the underground medium structure is determined exploration well location or other subterranean resource present position exactly.
Background technology
Utilizing the earthquake reflected wave data information that stratal configuration is carried out to similarly is one of core technology of seismic prospecting.This technology that is called migration imaging (migration imaging), the spitting image of the CT imaging in the medical science, but its complexity but is to be difficult to analogy.The difficulty of seismic migration imaging is: the data amount is big especially; Signal source is uncertain; The dielectric structure complexity; Interference noise is powerful in the seismic wave propagation, and the unfavorable factor that diffraction, shielding, absorption etc. can't be grasped all produces bad influence to final imaging precision.Meanwhile, with regard to cost, the industrial input of a bite prospecting borehole generally reaches tens million of round Renminbi.Therefore the industrial migration imaging achievement fast and accurately that always requires is as disposing the exploration scheme, determining the foundation of well location.Therefore the three-dimensional seismic offset imaging development of technology is long-expected, in the geophysics field, the three-dimensional seismic offset imaging technology is to be subjected to common concern for a long time, to have attracted large quantities of outstanding scientists and engineering technical personnel to be engaged in the research of this respect from the theory to the commercial Application.
Over nearly two, 30 years, by means of rapid development of computer technology, handle on a large computer system or the network group of planes and set up the software platform that much differs from one another being specifically designed to seismic data, this is a requisite basic normal signal handling implement in the energy industry exploration.But any platform all presses for and improves speed and the precision that migration imaging is handled, and does not also have a kind of platform that gratifying three-dimensional seismic offset imaging technology can be provided at present, and large quantities of technical forces is being devoted to research and develop new method.Current, what be subjected to generally paying attention to is three-dimensional implicit schemes finite-difference migration method, and the advantage of this method is: adapt to the horizontal change of underground speed, the precision height; Shortcoming is that calculated amount is too big.
1, data volume problem: for the poststack data information of 10km * 10km, if line-spacing and CDP spacing are 25m, sampling interval is 4.00ms, data on each time point reach 400 * 400=160000 so, the equation exponent number that needs to find the solution on this time point is 160000 * 160000, is equation 3000-6000 time that 6 seconds data need be found the solution above-mentioned exponent number for record length.
2, counting yield problem:, calculate numerous and diverse known method for solving usually and can't be competent at because data volume is big.Existing mathematical theory of computation, general only research once or the algorithm of finding the solution several times.If, no longer be the problem of efficient and cost, whether the horizon problem of commercial Application has been arranged for thousands of so huge equations of exponent number of finding the solution of the 3D seismic data of small size very inferiorly.
3, direction distortion problem: because thousands of inferior second order wave equations of finding the solution in the wave field extrapolation process all use the data of having calculated at every turn, if can not control distortion in handling well, error transfer accumulation so is enough to make that end result has no practical significance.
So for a long time, people expectation be computing technique fast and accurately, if the problem of computing technique aspect has solved, what the problem of imaging aspect also will become is normal simple.To invert is relatively more commonly used in division in the existing known computing technique, roughly has alternating direction implicit (ADI) method to be also referred to as along X-Y direction division (Brown, 1983 etc.), to decompose (Rickett andClaerbout, 1998) based on the dimensionality reduction LU on spiral border.The latter is a kind of nonlinear method, can't avoid chaos phenomenon in the calculating, and boundary treatment also has problem.The former because circular symmetry is destroyed, can't overcomes the direction distortion end result is distorted in extrapolation process.Legend and application examples in the back will provide and specify.
Summary of the invention
1, goal of the invention: the present invention adopts a kind of new division approach method, realized that a kind of incomplete Cholesky decomposes (because the matrix of coefficients of Cholesky decomposition is that dense can not being used for calculated fully), element non-zero on the regulation only a few position, adopt the principle of error minimum to calculate these non-zero entry, finish under the sparse condition of assurance coefficient matrix height and approach, this and classic method not regulation non-zero entry position have essential difference.Owing to do not adopt along the direction division, avoided the direction distortion, guarantee the precision of migration imaging.
2, technical scheme: the present invention is achieved through the following technical solutions:
A kind of method that improves the geological data processing accuracy, it comprises and at first adopts conventional means to gather seismic exploration data, uses the fixed-point decomposing method of three-dimensional seismic offset imaging to handle to the seismic exploration data that collects then; The fixed-point decomposing method of described three-dimensional seismic offset imaging follows these steps to carry out:
A, the piece LU that at first matrix of coefficients of the wave field extrapolation equation of method of finite difference skew is done fixed point decompose:
Matrix is before decomposing
Figure C20061004721100061
T wherein k, E KBe sub-block matrix, T kBe triple diagonal matrix, E kBe diagonal matrix, make piece two diagonal angles and be decomposed into A=LL T+ R; Wherein L is the piece lower triangular matrix, and non-zero entry is arranged on principal diagonal, and show the non-zero entry except that adjacent in the left side of non-null matrix on the principal diagonal, and all the other positions are null matrix; L TTransposed matrix for L; R is the residual matrix that the incomplete Cholesky of piece decomposes; Null matrix does not need to carry out numerical evaluation, calculates the non-zero entry only relate on the principal diagonal and two negation elements of principal diagonal both sides, and the fixed point LU that obtains coefficient matrices A decomposes or claims fixed point BIC to decompose;
B, the residual error of then fixed point BIC being decomposed are calculated:
Residual matrix R=R L+ R UBeing one respectively has along the matrix of a row non-zero entry of diagonal, wherein R in upper and lower triangular portions LBe the non-zero entry matrix of last triangular portions, R fFor the non-zero entry matrix of following triangular portions, use x=L -TL -1B-Rx obtains a kind of iterative algorithm x (k+1)=L -TL -1B-Rx (k)
C, determine the pivot advantage according to seismic exploration data, the step size computation of extrapolating is selected suitable step-length;
D, the fixed point breakdown that uses the A that obtains in the above-mentioned steps and selected step-length obtain the image of the higher underground medium structure of precision;
E. by the image of the underground medium structure that obtains, determine the drilling well well location or reconnoitre other mineral resources reliable basis is provided for disposing the exploration scheme.
Preferably, the unit of R is in a small amount with respect to pivot in step a, and is controllable.
Preferably, if wish that in step c convergence is quicker, then gathering and pretreatment stage is taked trace interpolation or chosen suitable step-length and realize.
3, advantage and effect: the present invention is for solving the three-dimensional seismic offset imaging precision problem of long-term puzzlement exploration geophysics industry member, having adopted a kind of brand-new division approximation technique---the incomplete Cholesky that fixes a point decomposes, and has realized the wave field extrapolation fast and accurately of 3D seismic data is handled.Make Three dimensional finite difference method migration and imaging techniques realize the transformation from the theoretical research to the commercial Application.Before this because or to such an extent as to being that calculated amount is too huge cannot bear, or be to produce unacceptable distortion in the wave field extrapolation, make the method for finite difference three-D migration rest on the state of theoretical research for a long time.The circular symmetry that the present invention has kept wave field to propagate in data handling procedure is fully reliably guaranteed final imaging precision, because counting yield improves tens times than classic method, has guaranteed industrial practicality from assessing the cost again.The present invention has solved a technical barrier that cannot solve for a long time for exploration geophysics industry.Can confirm in the practical application that the computing method that adopted are calculated in institute's Data Processing production in the Liaohe Oil Field has tangible technique effect and remarkable economic efficiency in the enforcement.
Description of drawings
Accompanying drawing 1 is the classification chart of 3 D earthquake deflection method;
The various correction charts that accompanying drawing 2 is done for the distortion at wave field extrapolation operator circular symmetry;
Accompanying drawing 3 is the impulse response figure of CG software skew;
Accompanying drawing 4 is a wave field extrapolation circular symmetry design sketch;
Accompanying drawing 5 is another kind of wave field extrapolation circular symmetry design sketch;
Accompanying drawing 6 is a wave field extrapolation net symmetry design sketch of the present invention;
Accompanying drawing 7 is the preceding stacked profile map of Q271 line skew;
Accompanying drawing 8 is Q271 line offset effect figure;
Accompanying drawing 9 is another Q271 line offset effect figure;
Accompanying drawing 10 is Q271 line offset effect figure of the present invention;
Accompanying drawing 11 is three kinds of Q271 line offset effect comparison diagrams;
Accompanying drawing 12 is 0251 line offset effect figure;
Accompanying drawing 13 is another 0251 line offset effect figure;
Accompanying drawing 14 is 0251 line offset effect figure of the present invention;
Accompanying drawing 15 is three kind of 0251 line offset effect comparison diagram.
Embodiment
The present invention mainly is by computing technique fast and accurately, the earthquake reflected wave data information is calculated, the data result that obtains forms image comparatively accurately in the seismic data processing software platform, thereby obtains stratal configuration explored achievement comparatively accurately.Existing three-D migration formation method has a variety of, sees shown in the accompanying drawing 1.But the image that forms stratal configuration for seismic prospecting is difficult to obtain satisfied precision, mainly is because calculated amount is excessive, the accumulation of extrapolation process distortion error or produce chaos phenomenon in calculating, and the problem of aspect such as boundary treatment.The normal at present computing method that adopt have alternating direction implicit (ADI) method and spiral border Nonlinear Dimension Reduction LU factorization, spiral border Nonlinear Dimension Reduction LU factorization can not be avoided chaos phenomenon in calculating, 9 points, 13 points, 17 points, and sexangle McCllalen conversion are adopted in the various corrections that accompanying drawing 2 is done for the distortion at wave field extrapolation operator circular symmetry at 25 in the data.Accompanying drawing 3 is the impulse response of CG software skew, and as can be seen from the figure Pian Yi noise is serious.Accompanying drawing 4 is the impulse response dropping cut slice of foreign software CG skew, and as can be seen from the figure circular symmetry obviously distorts, and migration noise is serious.Accompanying drawing 5 is the impulse response dropping cut slice of software GR skew, and as can be seen from the figure circular symmetry obviously distorts, and migration noise is serious.The present invention will be further described below by embodiment:
Embodiment: for the poststack data information of 10km * 10km, adopting line-spacing and CDP spacing is 25m, and sampling interval is 4.00ms, and the matrix of coefficients that the data that obtain can be formed the wave field extrapolation equation is as follows:
Figure C20061004721100081
A, the LU that at first coefficient matrices A of wave field extrapolation equation is done fixed point decompose, wherein T k, E KBe sub-block matrix, T kBe triple diagonal matrix, E kMake piece two diagonal angles for diagonal matrix and be decomposed into A=LL T+ R; Wherein L is the piece lower triangular matrix, and only shows non-zero entry with principal diagonal left side adjacent on principal diagonal, and all the other positions are null matrix; L TTransposed matrix for L; R is the residual matrix that the incomplete Cholesky of piece decomposes; Null matrix does not need numerical evaluation, and this matrix computations only relates to two non-zero entry in principal diagonal and principal diagonal left side, so obtain the fixed point LU decomposed form of the coefficient matrices A of wave field extrapolation equation;
B, the residual error of then fixed point BIC being decomposed are calculated, and it is that a kind of sparse especially LU decomposes that the fixed point BIC that the present invention did decomposes.Its residual matrix R=R L+ R UBe one and row non-zero entry matrix, wherein a R along diagonal respectively arranged in upper and lower triangular portions LBe the non-zero entry matrix of last triangular portions, R fFor the non-zero entry matrix of following triangular portions, use x=L -TL -1B-Rx is easy to obtain a kind of iterative algorithm x (k+1)=L -TL -1B-Rx (k)
C, can determine the pivot advantage according to real data, the step size computation of extrapolating is selected suitable step-length;
In d, the software module in being embedded into the seismic data processing software platform, use fixed point breakdown, and the selected step-length of the A that obtains in the above-mentioned steps, can obtain the higher image of precision.
Image shown in the accompanying drawing 6 is a model data, can clearly find out its strict circular symmetry that keeps from figure, does not have migration noise.
The unit of R is in a small amount with respect to pivot in step a, and can control.
X in step b (k+1)=L -TL -1B-Rx (k)In order to improve the precision of algorithm.This is a kind of artificially method of control accuracy, can do balance between computer resource and computational accuracy during use.This is the characteristics that do not have in the existing business software.
If wish that in step c convergence is quicker, can gather and pretreatment stage by trace interpolation, or choose suitable step-length and accomplish the end in view.After as above handling, the data that obtain have formed comparatively desirable skew achievement section (as accompanying drawing 10) in the seismic data processing software platform.Accompanying drawing 7 is the preceding stacked section of Q271 line skew.Accompanying drawing 8 is CG software Q271 line offset effect figure, and it is relatively poor as can be seen from the figure to iris out the position imaging.Accompanying drawing 9 is GR software Q271 line offset effect figure, and it is undesirable as can be seen from the figure to iris out the position imaging.Accompanying drawing 10 is the fix a point Q271 line offset effect figure of decomposition method of the present invention, and the image imaging of as can be seen from the figure irising out the position is clear.Accompanying drawing 11 is three kinds of Q271 line offset effect comparison diagrams, the left side is GR software imaging effect figure, and middle is fixed point decomposition method imaging effect figure, and the right is foreign software CG imaging effect figure, by the intermediate image imaging is clear more as can be seen, effect is significantly better than the imaging effect of left and right two images.Accompanying drawing 12 is a GR software migration result, and the minor fold form that square frame is irised out the position in the image does not meet geologic rule, though the circled continuity is all right, level is abundant inadequately.Accompanying drawing 13 is a CG software migration result, and it is reliable that square frame is irised out the right flank form at position in the image, and left wing's reflection is lost, and the circled poor continuity, laterally is difficult to follow the trail of.Accompanying drawing 14 is the migration result of fixed point decomposition method of the present invention, position signal to noise ratio (S/N ratio) height in the circle as can be seen among the figure, and the position imaging is accurate in the square frame.Antiform is distinguishable, but continuity is not ideal, and circled continuity resolution is obviously improved, and the cap rock breakpoint is clear-cut.Accompanying drawing 15 is three kind of 0251 line offset effect comparison diagram, the left side is the image that software GR obtains, the centre is the image that the fixed point decomposition method obtains, the right is the image that foreign software CG shape obtains, as can be seen from the figure the symmetry of the interior intermediate image of rectangle frame obviously is better than the image on both sides, and imaging effect is significantly better than the image on left and right both sides.
By the enforcement of the technical program, the migration imaging achievement that obtains can reflect the structure of underground medium comparatively accurately, for disposing the exploration scheme, determine the exploration well location and reliable foundation being provided.

Claims (3)

1. the fixed-point decomposing method of a three-dimensional seismic offset imaging, be a kind of method that improves the geological data processing accuracy, it comprises and at first adopts conventional means to gather seismic exploration data, uses the fixed-point decomposing method of three-dimensional seismic offset imaging to handle to the seismic exploration data that collects then; The fixed-point decomposing method of described three-dimensional seismic offset imaging follows these steps to carry out:
A, the piece LU that at first matrix of coefficients of the wave field extrapolation equation of method of finite difference skew is done fixed point decompose:
Matrix is before decomposing
Figure C2006100472110002C1
T wherein k, E KBe sub-block matrix, T kBe triple diagonal matrix, E kBe diagonal matrix, make piece two diagonal angles and be decomposed into A=LL T+ R; Wherein L is the piece lower triangular matrix, and non-zero entry is arranged on principal diagonal, and show the non-zero entry except that adjacent in the left side of non-null matrix on the principal diagonal, and all the other positions are null matrix; L TTransposed matrix for L; R is the residual matrix that the incomplete Cholesky of piece decomposes; Null matrix does not need to carry out numerical evaluation, calculates the non-zero entry only relate on the principal diagonal and two negation elements of principal diagonal both sides, and the fixed point LU that obtains coefficient matrices A decomposes or claims fixed point BIC to decompose;
B, the residual error of then fixed point BIC being decomposed are calculated:
Residual matrix R=R L+ R UBeing one respectively has along the matrix of a row non-zero entry of diagonal, wherein R in upper and lower triangular portions LBe the non-zero entry matrix of last triangular portions, R TFor the non-zero entry matrix of following triangular portions, use x=L -TL -1B-Rx obtains a kind of iterative algorithm x (k+1)=L -TL -1B-Rx (k)
C, determine the pivot advantage according to seismic exploration data, the step size computation of extrapolating is selected suitable step-length;
D, the fixed point breakdown that uses the A that obtains in the above-mentioned steps and selected step-length obtain the image of the higher underground medium structure of precision;
E. the image by the underground medium structure that obtains is determined the drilling well well location or is surveyed for disposing the exploration scheme
2. a kind of method that improves the geological data processing accuracy according to claim 1 is characterized in that: the unit of R is in a small amount with respect to pivot in step a, and is controllable.
3. a kind of method that improves the geological data processing accuracy according to claim 1 is characterized in that: if wish that in step c convergence is quicker, then gathering and pretreatment stage is taked trace interpolation or chosen suitable step-length and realize.
CNB2006100472116A 2006-07-17 2006-07-17 Fixed-point decomposing method for three-dimensional seismic offset imaging Expired - Fee Related CN100383559C (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CNB2006100472116A CN100383559C (en) 2006-07-17 2006-07-17 Fixed-point decomposing method for three-dimensional seismic offset imaging

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CNB2006100472116A CN100383559C (en) 2006-07-17 2006-07-17 Fixed-point decomposing method for three-dimensional seismic offset imaging

Publications (2)

Publication Number Publication Date
CN1888935A CN1888935A (en) 2007-01-03
CN100383559C true CN100383559C (en) 2008-04-23

Family

ID=37578222

Family Applications (1)

Application Number Title Priority Date Filing Date
CNB2006100472116A Expired - Fee Related CN100383559C (en) 2006-07-17 2006-07-17 Fixed-point decomposing method for three-dimensional seismic offset imaging

Country Status (1)

Country Link
CN (1) CN100383559C (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108225274B (en) * 2017-12-26 2020-11-20 中国科学院电子学研究所 Conjugate gradient method bundle adjustment method based on incomplete decomposition pretreatment

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1073770A (en) * 1991-12-26 1993-06-30 切夫里昂研究和技术公司 Improve the method for architectonic seismic resolution
CN1111018A (en) * 1993-07-26 1995-11-01 埃克森生产研究公司 Migration velocity analysis using limited aperture migration
US6492811B1 (en) * 2000-11-22 2002-12-10 Koninklijke Philips Electronics N.V. Multishot echo planar imaging fluoroscopic technique

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1073770A (en) * 1991-12-26 1993-06-30 切夫里昂研究和技术公司 Improve the method for architectonic seismic resolution
CN1111018A (en) * 1993-07-26 1995-11-01 埃克森生产研究公司 Migration velocity analysis using limited aperture migration
US6492811B1 (en) * 2000-11-22 2002-12-10 Koninklijke Philips Electronics N.V. Multishot echo planar imaging fluoroscopic technique

Also Published As

Publication number Publication date
CN1888935A (en) 2007-01-03

Similar Documents

Publication Publication Date Title
US10795053B2 (en) Systems and methods of multi-scale meshing for geologic time modeling
CN106772583B (en) A kind of earthquake diffracted wave separation method and device
CN103122762B (en) Detection method and device for effective fractured intervals in unconventional shale oil and gas reservoir
EP3371629B1 (en) Representing structural uncertainty in a mesh representing a geological environment
CN103926619B (en) Reverse time migration method of three-dimensional VSP data
CN102692644B (en) Depth domain common-image gather generation method
WO2014078355A1 (en) Unstructured grids for modeling reservoirs
CN103487835A (en) Multi-resolution wave impedance inversion method based on model constraints
CN103245970A (en) Pre-stack seismic wide angle retrieval method
Cainelli et al. On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations
CN111638556A (en) Magnetotelluric forward modeling method and device based on geospatial solution strategy and storage medium
Wang et al. Seismic velocity inversion transformer
CN103376463B (en) A kind of Inverse modeling method based on faults control
Bai et al. 3-D non-linear travel-time tomography: imaging high contrast velocity anomalies
US6324478B1 (en) Second-and higher-order traveltimes for seismic imaging
CN100383559C (en) Fixed-point decomposing method for three-dimensional seismic offset imaging
CN106662665B (en) The interpolation and convolution of rearrangement for the processing of faster staggered-mesh
CN103217715B (en) Multiple dimensioned regular grid Static Correction of Tomographic Inversion method
CN115375867B (en) Method, system, equipment and medium for calculating geothermal resource quantity by using grid model
Zhou et al. An iterative factored topography-dependent eikonal solver for anisotropic media
CN104750954A (en) Method and device for simulating earthquake waves in complex anisotropic media
CN100383558C (en) Block pursuing method for three-dimensional seismic offset imaging
CN100547434C (en) A kind of inverse operator approach method of three-dimensional seismic offset imaging
YU et al. Joint inversion of gravity and seismic data based on common gridded model with random density and velocity distributions
Gong et al. Application of multi-level and high-resolution fracture modeling in field-scale reservoir simulation study

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
C17 Cessation of patent right
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20080423

Termination date: 20110717