地震波谱元法数值模拟和偏移研究

With the development of both seismic theory and computer technology, numerical modeling technology of seismic wave has achieved great advancement during the past half century. The current methods under development include finite differentiation method (FDM), finite element method (FEM), pseudospectral method (PSM), integral equation method (IEM) and spectral element method (SEM). They exert their very important roles in every corner of seismology and seismic prospecting. Large quantity of researches towards spectral element method in the end of last century bring this method to a new era, which results in perfect solution of many difficult problems. However, parts of posterior works such as seismic migration and inversion which base on spectral element method have never been studied widely at least up to the present whereas are of importance to seismic imaging and seismic wave propagation. Based on previous work, this paper uses spectral element method to investigate the characteristics and laws of the seismic wave propagation in isotropic and anisotropic media. By thoroughly studying this high-accuracy method, we implement a kind of reverse-time pre- and post-stack migration based on SEM. In order to verify the validity of the SEM method, we have simulated the propagation of seismic wave in several different models. The simulation results show that: (1) spectral element method can be used to model any complex models and the computational results are comparable with the expected results and the analytic results; (2) the optimum accuracy can be achieved when the rank is between 4 and 9. When it is below 4, the dispersion may occur; and when it is above 9, the time step-length will be changed accordingly with the reducing space step-length in order to keep the computation stability. This will exponentially increase the computation time and at the same time the memory even if simulating the same media. This paper also applies explosive reflection surface imaging technology, time constancy principle of wave-filed extrapolation and least travetime raytracing technology of surface source to SEM pre- and post-stack migration of isotropic and anisotropic media. All imaging results derived by the above methods agree well with the real geological models and the position of interface and inflexions can also return to their right location well. This indicates that the method proposed in this paper is a kind of technology with high accuracy and robust stability. It can serve as an alternative method in real seismic data processing. All these work can boost the development of high-accuracy seismic imaging, and therefore have significant inference value.

电离层行进式扰动的GPS台网监测研究

In this dissertation, we investigated two types of traveling ionospheric disturbances (TIDs)/gravity waves (GWs) triggered separately by auroral energy input during super geomagnetic storms and solar terminator (ST) under quiet geomagnetic conditions (kp<3+) using TEC measurements from the global network of GPS receivers. Research into the generation and propagation of TIDs/GWs during storms greatly enhance our understandings on the evolution processes of energy transportation from the high-latitude’s magnetosphere to the low-latitude ionosphere and the conjugated effect of TIDs propagation between the northern and southern hemispheres. Our results revealed that the conjugacy of propagation direction between the northern and southern hemispheres was subject to the influence of Coriolis force. We also figure out the evolution processes of ionospheric disturbances at the global scale. These are important topics that had not been well addressed previously. In addition, we also obtained thee wave structures of medium scale TIDs excited by the solar terminator (ST) moving over the northern America and physical mechanisms involved. Our observations confirm that the ST is a stable and repetitive source of ionospheric wave disturbances and the evidence of solar terminator generated disturbances has been demonstrated experimentally via the GPS TEC measurement. The main researches and results of this dissertation are as follows. First, the global traveling ionospheric disturbances (TIDs) during the drastic magnetic storms of October 29–31, 2003 were analyzed using the Global Position System (GPS) total electron content (TEC) data observed in the Asian-Australian, European and North American sectors. We collected the most comprehensive set of the TEC data from more than 900 GPS stations on the International GNSS Services (IGS) website and introduce here a strategy that combines polynomial fitting and multi-channel maximum entropy spectral analysis to obtain TID parameters. Moreover, in collaboration with my thesis advisor, I have developed an imaging technique of 2-dimensional map of TIDs structures to obtain spatial and temporal maps of large scale traveling ionospheric disturbances (LSTIDs). The clear structures of TEC perturbations map during the passage of TIDs were displayed. The results of our study are summarized as follows: (1) Large-scale TIDs (LSTIDs) and medium-scale TIDs (MSTIDs) were detected in all three sectors after the sudden commencement (SC) of the magnetic storm, and their features showed longitudinal and latitudinal dependences. The duration of TIDs was longer at higher latitudes than at middle latitudes, with a maximum of about 16 h. The TEC variation amplitude of LSTIDs was larger in the North American sector than in the two other sectors. At the lower latitudes, the ionospheric perturbations were more complicated, and their duration and amplitude were relatively longer and larger. (2) The periods and phase speeds of TIDs were different in these three sectors. In Europe, the TIDs propagated southward; in North America and Asia, the TIDs propagated southwestward; in the near-equator region, the disturbances propagated with the azimuth (the angle of the propagation direction of the LSTIDs measured clockwise from due north with 0°) of 210° showing the influence of Coriolis force; in the Southern Hemisphere, the LSTIDs propagated conjugatedly northwestward. Both the southwestward and northeastward propagating LSTIDs are found in the equatorial region. These results mean that the Coriolis effect cannot be ignored for the wave propagation of LSTIDs and that the propagation direction is correlated with the polar magnetic activity. (3) The day (day of year: 301) before the SC (sudden commencement) of magnetic storm, we observed a sudden TEC skip disturbances (±10 TECU). It should be a response for the high flux of proton during the solar flare event, but not the magnetic storms. Next, the most comprehensive and dense GPS network’s data from North-America region were used in this paper to analyze the medium scale traveling ionospheric disturbances (MSTIDs) which were generated by the moving solar terminator during the quiet days in 2005. We applied the multi-channel maximum entropy spectral analysis to calculated TID parameters, and found that the occurrence of ST-MSTIDs depends on the seasonal variations. The results of our study are summarized as follows: (1) MSTIDs stimulated by the moving ST (ST-MSTIDs) are detected at mid-latitudes after the passage of the solar terminator with the life time of 2～3 hours and the variation amplitude of 0.2～0.8 TECU. Spectral analysis indicated that the horizontal wavelength, average period, horizontal phase velocity of the MSTIDs are around 300±150 km，150±80 m/s and 25±15 min, respectively. In addition, ST-MSTIDs have wave fronts elongating the moving ST direction and almost parallel to ST. (2) The statistical results demonstrate that the dusk MSTIDs stimulated by ST is more obvious than the dawn MSTIDs in summer. On the contrary, the more-pronounced dawn MSTIDs occurs in winter. (3) Further analysis indicates that the seasonal variations of ST-MSTIDs occurrence frequency are most probably related to the seasonal differences of the variations of EUV flux in the ionosphere region and recombination process during sunrise and sunset period at mid-latitudes. Statistical study of occurrence characteristics of TIDs using the GPS network in North-American and European during solar maximum, In conclusion, statistical studies of the propagation characteristics of TIDs, which excited by the two common origins including geomagnetic storms and moving solar terminator, were involved with global GPS TEC databasein this thesis. We employed the multichannel maximum entropy spectral analysis method to diagnose the characteristics of propagation and evolvement of ionospheric disturbances, also, the characteristics of their regional distribution and climatological variations were revealed by the statistic analysis. The results of these studies can improve our knowledge about the energy transfer in the solar-terrestrial system and the coupling process between upper and lower atmosphere (thermosphere-ionosphere-mesosphere). On the other hand, our results of the investigation on TIDs generated by particular linear origin such as ST are important for developing ionospheric irregularity physics and modeling the transionosphere radio wave propagation. Besides, the GPS TEC representation of the ST-generated ionospheric structure suggests a better possibility for investigating this phenomenon. Subsequently, there are scientific meaning of the result of this dissertation to deeply discuss the energy transfer and coupling in the ionosphere, as well as realistic value to space weather forecast in the ionosphere region.

地震波模拟：非规则网格离散的数值解法、吸收边界及其实际应用

This dissertation presents a series of irregular-grid based numerical technique for modeling seismic wave propagation in heterogeneous media. The study involves the generation of the irregular numerical mesh corresponding to the irregular grid scheme, the discretized version of motion equations under the unstructured mesh, and irregular-grid absorbing boundary conditions. The resulting numerical technique has been used in generating the synthetic data sets on the realistic complex geologic models that can examine the migration schemes. The motion equation discretization and modeling are based on Grid Method. The key idea is to use the integral equilibrium principle to replace the operator at each grid in Finite Difference scheme and variational formulation in Finite Element Method. The irregular grids of complex geologic model is generated by the Paving Method, which allow varying grid spacing according to meshing constraints. The grids have great quality at domain boundaries and contain equal quantities of nodes at interfaces, which avoids the interpolation of parameters and variables. The irregular grid absorbing boundary conditions is developed by extending the Perfectly Matched Layer method to the rotated local coordinates. The splitted PML equations of the first-order system is derived by using integral equilibrium principle. The proposed scheme can build PML boundary of arbitrary geometry in the computational domain, avoiding the special treatment at corners in a standard PML method and saving considerable memory and computation cost. The numerical implementation demonstrates the desired qualities of irregular grid based modeling technique. In particular, (1) smaller memory requirements and computational time are needed by changing the grid spacing according to local velocity; (2) Arbitrary surfaces and interface topographies are described accurately, thus removing the artificial reflection resulting from the stair approximation of the curved or dipping interfaces; (3) computational domain is significantly reduced by flexibly building the curved artificial boundaries using the irregular-grid absorbing boundary conditions. The proposed irregular grid approach is apply to reverse time migration as the extrapolation algorithm. It can discretize the smoothed velocity model by irregular grid of variable scale, which contributes to reduce the computation cost. The topography. It can also handle data set of arbitrary topography and no field correction is needed.

复杂介质中地震波场数值模拟研究——基于广义正交多项式迭积微分算子法

In this paper, we propose a new numerical modeling method – Convolutional Forsyte Polynomial Differentiator (CFPD), aimed at simulating seismic wave propagation in complex media with high efficiency and accuracy individually owned by short-scheme finite differentiator and general convolutional polynomial method. By adjusting the operator length and optimizing the operator coefficient, both global and local informations can be easily incorporated into the wavefield which is important to invert the undersurface geological structure. The key issue in this paper is to introduce the convolutional differentiator based on Forsyte generalized orthogonal polynomial in mathematics into the spatial differentiation of the first velocity-stress equation. To match the high accuracy of the spatial differentiator, this method in the time coordinate adopts staggered grid finite difference instead of conventional finite difference to model seismic wave propagation in heterogeneous media. To attenuate the reflection artifacts caused by artificial boundary, Perfectly Matched Layer (PML) absorbing boundary is also being considered in the method to deal with boundary problem due to its advantage of automatically handling large-angle emission. The PML formula for acoustic equation and first-order velocity-stress equation are also derived in this paper. There is little difference to implement the PML boundary condition in all kind of wave equations, but in Biot media, special attenuation factors should be taken. Numerical results demonstrate that the PML boundary condition is better than Cerjan absorbing boundary condition which makes it more suitable to hand the artificial boundary reflection. Based on the theories of anisotropy, Biot two-phase media and viscous-elasticity, this paper constructs the constitutive relationship for viscous-elastic and two-phase media, and further derives the first-order velocity-stress equation for 3D viscous-elastic and two-phase media. Numerical modeling using CFPD method is carried out in the above-mentioned media. The results modeled in the viscous-elastic media and the anisotropic pore elastic media can better explain wave phenomena of the true earth media, and can also prove that CFPD is a useful numerical tool to study the wave propagation in complex media.

中高层大气对低层大气声重波响应的分析

Acoustic Gravity waves (AGW) play an important role in balancing the atmospheric energy and momentum budget. Propagation of gravity wave in the atmosphere is one of the important factors of changing middle and upper atmosphere and ionosphere. The purpose of this dissertation is to study the propagation of gravity wave in a compression atmosphere whit means of numerical simulation and to analyze the response of middle and upper atmosphere to pulse disturbance from lower atmosphere. This work begins with the establishment of 2-D fully nonlinear compressible atmospheric dynamic model in polar coordinate, which is used ton numerically study gravity wave propagation. Then the propagation characteristics of acoustic gravity wave packets are investigated and discussed. We also simulate the response of middle and upper atmosphere to pulse disturbance of lower atmosphere in background winds or without background winds by using this model and analyze the data we obtained by using Fourier Transform (FT), Short-time Fourier Transform (STFT) and Empirical Mode Decomposition (EMD) method which is an important part of Hilbert-Huang Transform (HHT). The research content is summarized in the following: 1. By using a two-dimensional full-implicit-continuous-Eulerian (FICE) scheme and taking the atmospheric basic motion equations as the governing equations, a numerical model for nonlinear propagation of acoustic gravity wave disturbance in two-dimensional polar coordinates is solved. 2. Then the propagation characteristics of acoustic gravity wave packets are investigated and discussed. Results of numerical simulation show that the acoustic gravity wave packets propagate steadily upward and keep its shape well after several periods. 3. We simulate the response of middle and upper atmosphere to pulse disturbance of lower atmosphere in background winds or without background winds by using this model, and obtain the distribution of a certain physical quantity in time and space from earth’s surface to 300km above. The results reveal that the response of ionosphere occurs at a large horizontal distance from the source and the disturbance becomes greater with increasing of height. The situation when the direction of the background wind is opposite to or the same as the direction of disturbed velocity of gravity-wave is studied. The results show that gravity wave propagating against winds is easier than those propagating along winds and the background wind can accelerate gravity wave propagation. Just upon the source, an acoustic wave component with period of 6 min can be found. These images of simulation are similar to observations of the total electron content (TEC) disturbances caused by the great Sumatra-Andaman earthquake on December 26 in 2004. 4. Using the EMD method the disturbed velocity data of a certain physical quantity in time and space can be decomposed into a series of intrinsic mode function (IMF) and a trend mode respectively. The results of EMD reveal impact of the gravity wave frequency under the background winds.

离散裂缝诱导各向异性粘弹性地震波模拟与特征分析

The real media always attenuate and distort seismic waves as they propagate in the earth. This behavior can be modeled with a viscoelastic and anisotropic wave equation. The real media can be described as fractured media. In this thesis, we present a high-order staggered grid finite-difference scheme for 2-D viscoelastic wave propagation in a medium containing a large number of small finite length fractures. We use the effective medium approach to compute the anisotropic parameters in each grid cell. By comparing our synthetic seismogram by staggered-grid finite-difference with that by complex-ray parameter ray tracing method, we conclude that the high-order staggered-grid finite-difference technique can effectively used to simulate seismic propagation in viscoelastic-anisotropic media. Synthetic seismograms demonstrate that strong attenuation and significant frequency dispersion due to viscosity are important factors of reducing amplitude and delaying arrival time varying with incidence angle or offset. On the other hand, the amount of scattered energy not only provides an indicator of orientation of fracture sets, but can also provide information about the fracture spacing. Analysis of synthetic seismograms from dry- and fluid-filled fractures indicates that dry-filled fractures show more significant scattering on seismic wavefields than fluid-filled ones, and offset-variations in P-wave amplitude are observable. We also analyze seismic response of an anticlinal trap model that includes a gas-filled fractured reservoir with high attenuation, which attenuates and distorts the so-called bright spot.

旋转轴对称介质中的qP波一阶射线追踪

Elastic anisotropy is a very common phenomenon in the Earth’s interior, especial for sedimentary rock as important gas and oil reservoirs. But in the processing and interpretation of seismic data, it is assumption that the media in the Earth’s interior is completely elastic and isotropic, and then the methods based on isotropy are used to deal with anisotropic seismic data, so it makes the seismic resolution lower and the error on images is caused. The research on seismic wave simulation technology can improve our understanding on the rules of seismic wave propagation in anisotropic media, and it can help us to resolve problems caused by anisotropy of media in the processing and interpretation of seismic data. So researching on weakly anisotropic media with rotated axis of symmetry, we study systematically the rules of seismic wave propagation in this kind of media, simulate the process with numerical calculation, and get the better research results. The first-order ray tracing (FORT) formulas of qP wave derived can adapt to every anisotropic media with arbitrary symmetry. The equations are considerably simpler than the exact ray tracing equations. The equations allow qP waves to be treated independently from qS waves, just as in isotropic media. They simplify considerably in media with higher symmetry anisotropy. In isotropic media, they reduce to the exact ray tracing equations. In contrast to other perturbation techniques used to trace rays in weakly anisotropic media, our approach does not require calculation of reference rays in a reference isotropic medium. The FORT-method rays are obtained directly. They are computationally more effective than standard ray tracing equations. Moreover the second-order travel time corrections formula derived can be used to reduce effectively the travel time error, and improve the accuracy of travel time calculation. The tensor transformation equations of weak-anisotropy parameters in media with rotated axis of symmetry derived from the Bond transformation equations resolve effectively the problems of coordinate transformation caused by the difference between global system of coordinate and local system of coordinate. The calculated weak-anisotropy parameters are completely suitable to the first-order ray tracing used in this paper, and their forms are simpler than those from the Bond transformation. In the numerical simulation on ray tracing, we use the travel time table calculation method that the locations of the grids in the ray beam are determined, then the travel times of the grids are obtained by the reversed distance interpolation. We get better calculation efficiency and accuracy by this method. Finally we verify the validity and adaptability of this method used in this paper with numerical simulations for the rotated TI model with anisotropy of about 8% and the rotated ORTHO model with anisotropy of about 20%. The results indicate that this method has better accuracy for both media with different types and different anisotropic strength. Keywords: weak-anisotropy, numerical simulation, ray tracing equation, travel time, inhomogeneity

复杂介质边界元法数值模拟研究

A major impetus to study the rough surface and complex structure in near surface model is because accuracy of seismic observation and geophysical prospecting can be improved. Wave theory study about fluid-satuated porous media has important significance for some scientific problems, such as explore underground resources, study of earth's internal structure, and structure response of multi-phase porous soil under dynamic and seismic effect. Seismic wave numerical modeling is one of the effective methods which understand seismic propagation rules in complex media. As a numerical simulation method, boundary element methods had been widely used in seismic wave field study. This paper mainly studies randomly rough surface scattering which used some approximation solutions based on boundary element method. In addition, I developed a boundary element solution for fluid saturated porous media. In this paper, we used boundary element methods which based on integral expression of wave equation to study the free rough surface scattering effects of Kirchhoff approximation method, Perturbation approximation method, Rytov approximation method and Born series approximation method. Gaussian spectrum model of randomly rough surfaces was chosen as the benchmark model. The approximation methods result were compared with exact results which obtained by boundary element methods, we study that the above approximation methods were applicable how rough surfaces and it is founded that this depends on and ( here is the wavenumber of the incident field, is the RMS height and is the surface correlation length ). In general, Kirchhoff approximation which ignores multiple scatterings between any two surface points has been considered valid for the large-scale roughness components. Perturbation theory based on Taylor series expansion is valid for the small-scale roughness components, as and are .Tests with the Gaussian topographies show that the Rytov approximation methods improves the Kirchhoff approximation in both amplitude and phase but at the cost of an extra treatment of transformation for the wave fields. The realistic methods for the multiscale surfaces come with the Born series approximation and the second-order Born series approximation might be sufficient to guarantee the accuracy of randomly rough surfaces. It could be an appropriate choice that a complex rough surface can be divided into large-, medium-, and small-scale roughness components with their scattering features be studied by the Kirchhoff or Rytov phase approximations, the Born series approximation, and the perturbation theory, respectively. For this purpose, it is important to select appropriate parameters that separate these different scale roughness components to guarantee the divided surfaces satisfy the physical assumptions of the used approximations, respectively. In addition, in this paper, the boundary element methods are used for solving the porous elastic wave propagation and carry out the numerical simulation. Based on the fluid-saturated porous model, this paper analyses and presents the dynamic equation of elastic wave propagation and boundary integral equation formulation of fluid saturated porous media in frequency domain. The fundamental solutions of the elastic wave equations are obtained according to the similarity between thermoelasticity and poroelasticity. At last, the numerical simulation of the elastic wave propagation in the two-phase isotropic media is carried out by using the boundary element method. The results show that a slow quasi P-wave can be seen in both solid and fluid wave-field synthetic seismograms. The boundary element method is effective and feasible.

地震波传播模拟的边界元―体积元积分方程法及其近似求解

To improve the efficiency of boundary-volume integral equation technique, this paper is involved in the approximate solutions of boundary-volume integral equation technique. Firstly, based on different interpretations of the self-interaction and extrapolation operators of the resulting boundary integral equation matrix, two different hybrid BEM+Born series modeling schemes are formulated and validated through comparisons with the full-waveform BE numerical solutions for wave propagation simulation in a semicircular alluvial valley and a complex fault model respectively. Numerical experiments indicate that both the BEM+Born series modeling schemes are suitable for complex geological structures and significantly improve computational efficiency especially for the cases of high frequencies and multisource seismic survey. Then boundary-volume integral equation technique is illuminated in detail and verified by modeling wave propagation in complex media. Furthermore, the first-order and second-order Born approximate solutions for the volume-scattering waves are studied and quantified by numerical simulation in different random medium models. Finally, preconditioning generalized minimal residual method is applied to solve boundary-volume integral equation and compared with Gaussian elimination method. Numerical experiments indicate this method makes the calculations more efficient.

粘弹性介质中的地震波传播与波形反演研究

The real earth is far away from an ideal elastic ball. The movement of structures or fluid and scattering of thin-layer would inevitably affect seismic wave propagation, which is demonstrated mainly as energy nongeometrical attenuation. Today, most of theoretical researches and applications take the assumption that all media studied are fully elastic. Ignoring the viscoelastic property would, in some circumstances, lead to amplitude and phase distortion, which will indirectly affect extraction of traveltime and waveform we use in imaging and inversion. In order to investigate the response of seismic wave propagation and improve the imaging and inversion quality in complex media, we need not only consider into attenuation of the real media but also implement it by means of efficient numerical methods and imaging techniques. As for numerical modeling, most widely used methods, such as finite difference, finite element and pseudospectral algorithms, have difficulty in dealing with problem of simultaneously improving accuracy and efficiency in computation. To partially overcome this difficulty, this paper devises a matrix differentiator method and an optimal convolutional differentiator method based on staggered-grid Fourier pseudospectral differentiation, and a staggered-grid optimal Shannon singular kernel convolutional differentiator by function distribution theory, which then are used to study seismic wave propagation in viscoelastic media. Results through comparisons and accuracy analysis demonstrate that optimal convolutional differentiator methods can solve well the incompatibility between accuracy and efficiency, and are almost twice more accurate than the same-length finite difference. They can efficiently reduce dispersion and provide high-precision waveform data. On the basis of frequency-domain wavefield modeling, we discuss how to directly solve linear equations and point out that when compared to the time-domain methods, frequency-domain methods would be more convenient to handle the multi-source problem and be much easier to incorporate medium attenuation. We also prove the equivalence of the time- and frequency-domain methods by using numerical tests when assumptions with non-relaxation modulus and quality factor are made, and analyze the reason that causes waveform difference. In frequency-domain waveform inversion, experiments have been conducted with transmission, crosshole and reflection data. By using the relation between media scales and characteristic frequencies, we analyze the capacity of the frequency-domain sequential inversion method in anti-noising and dealing with non-uniqueness of nonlinear optimization. In crosshole experiments, we find the main sources of inversion error and figure out how incorrect quality factor would affect inverted results. When dealing with surface reflection data, several frequencies have been chosen with optimal frequency selection strategy, with which we use to carry out sequential and simultaneous inversions to verify how important low frequency data are to the inverted results and the functionality of simultaneous inversion in anti-noising. Finally, I come with some conclusions about the whole work I have done in this dissertation and discuss detailly the existing and would-be problems in it. I also point out the possible directions and theories we should go and deepen, which, to some extent, would provide a helpful reference to researchers who are interested in seismic wave propagation and imaging in complex media.

纵波方位属性裂缝检测技术综合研究

Fractured oil and gas reservoir is an important type of oil and gas reservoir, which is taking a growing part of current oil and gas production in the whole world. Thus these technologies targeted at exploration of fractured oil and gas reservoirs are drawing vast attentions. It is difficult to accurately predict the fracture development orientation and intensity in oil and gas exploration. Focused on this problem, this paper systematically conducted series study of seismic data processing and P-wave attributes fracture detection based on the structure of ZX buried mountain, and obtained good results. This paper firstly stimulated the propagation of P-wave in weak anisotropic media caused by vertical aligned cracks, and analyzed the rule of P-wave attributes’ variation associated with observed azimuth, such as travel-time, amplitude and AVO gradient and so on, and quantitatively described the sensitive degree of these attributes to anisotropy of fracture medium. In order to further study the sensitive degree of these attributes to anisotropy of fractures, meanwhile, this paper stimulated P-wave propagation through different types and different intensity anisotropic medium respectively and summarized the rule of these attributes’ variation associated with observed azimuth in different anisotropic medium. The results of these studies provided reliable references for predicting orientation, extensity and size of actual complicated cracked medium by P-wave azimuth attributes responses. In the paper, amounts of seismic data processing methods are used to keep and recover all kinds of attributes applied for fracture detection, which guarantee the high accurate of these attributes, thus then improve the accurate of fracture detection. During seismic data processing, the paper adopted the three dimensional F-Kx-Ky field cone filter technique to attenuate ground roll waves and multiple waves, then enhances the S/N ratio of pre-stack seismic data; comprehensively applying geometrical spread compensation, surface consistent amplitude compensation, residual amplitude compensation to recover amplitude; common azimuth processing method effectively preserves the azimuthal characteristics of P-wave attributes; the technique of bend ray adaptive aperture pre-stack time migration insures to obtain the best image in each azimuth. Application of these processing methods guaranteed these attributes’ accuracy, and then improved the accuracy of fracture detection. After comparing and analyzing a variety of attributes, relative wave impedance (relative amplitude) attribute is selected to inverse the orientation of fracture medium; attenuation gradient and corresponding frequency of 85% energy are selected to inverse the intensity of fracture medium; then obtained the fracture distribution characteristics of lower Paleozoic and Precambrian in ZX ancient buried mountains. The results are good accord with the characteristics of faults system and well information in this area.

特大地震激发地球自由振荡的数值模拟研究

Large earthquakes, such as the Chile earthquake in 1960 and the Sumatra-Andaman earthquake on Dec 26, 2004 in Indonesia, have generated the Earth’s free oscillations. The eigenfrequencies of the Earth’s free oscillations are closely related to the Earth’s internal structures. The conventional methods, which mainly focus on calculating the eigenfrequecies by analytical ways, and the analysis on observations can not easily study the whole processes from earthquake occurrence to the Earth’s free oscillation inspired. Therefore, we try to use numerical method incorporated with large-scale parallel computing to study on the Earth’s free oscillations excited by giant earthquakes. We first give a review of researches and developments of the Earth’s free oscillation, and basical theories under spherical coordinate system. We then give a review of the numerical simulation of seismic wave propagation and basical theories of spectral element method to simulate global seismic wave propagation. As a first step to study the Earth’s free oscillations, we use a finite element method to simulate the propagation of elastic waves and the generation of oscillations of the chime bell of Marquis Yi of Zeng, by striking different parts of the bell, which possesses the oval crosssection. The bronze chime bells of Marquis Yi of Zeng are precious cultural relics of China. The bells have a two-tone acoustic characteristic, i.e., striking different parts of the bell generates different tones. By analysis of the vibration in the bell and the spectrum analysis, we further help the understanding of the mechanism of two-tone acoustic characteristics of the chime bell of Marquis Yi of Zeng. The preliminary calculations have clearly shown that two different modes of oscillation can be generated by striking different parts of the bell, and indicate that finite element numerical simulation of the processes of wave propagation and two-tone generation of the chime bell of Marquis Yi of Zeng is feasible. These analyses provide a new quantitative and visual way to explain the mystery of the two-tone acoustic characteristics. The method suggested by this study can be applied to simulate free oscillations excited by great earthquakes with complex Earth structure. Taking into account of such large-scale structure of the Earth, small-scale low-precision numerical simulation can not simply meet the requirement. The increasing capacity in high-performance parallel computing and progress on fully numerical solutions for seismic wave fields in realistic three-dimensional spherical models, Spectral element method and high-performance parallel computing were incorporated to simulate the seismic wave propagation processes in the Earth’s interior, without the effects of the Earth’s gravitational potential. The numerical simulation shows that, the results of the toroidal modes of our calculation agree well with the theoretical values, although the accuracy of our results is much limited, the calculated peaks are little distorted due to three-dimensional effects. There exist much great differences between our calculated values of spheroidal modes and theoretical values, because we don’t consider the effect the Earth’ gravitation in numerical model, which leads our values are smaller than the theoretical values. When , is much smaller, the effect of the Earth’s gravitation make the periods of spheroidal modes become shorter. However, we now can not consider effects of the Earth’s gravitational potential into the numerical model to simulate the spheroidal oscillations, but those results still demonstrate that, the numerical simulation of the Earth’s free oscillation is very feasible. We make the numerical simulation on processes of the Earth’s free oscillations under spherically symmetric Earth model using different special source mechanisms. The results quantitatively show that Earth’s free oscillations excited by different earthquakes are different, and oscillations at different locations are different for free oscillation excited by the same earthquake. We also explore how the Earth’s medium attenuation will take effects on the Earth’s free oscillations, and take comparisons with the observations. The medium attenuation can make influences on the Earth’s free oscillations, though the effects on lower-frequency fundamental oscillations are weak. At last, taking 2008 Wenchuan earthquake for example, we employ spectral element method incorporated with large-scale parallel computing technology to investigate the characteristics of seismic wave propagation excited by Wenchuan earthquake. We calculate synthetic seismograms with one-point source model and three-point source model respectively. Full 3-D visualization of the numerical results displays the profile of the seismic wave propagation with respect to time. The three-point source, which was proposed by the latest investigations through field observation and reverse estimation, can better demonstrate the spatial and temporal characteristics of the source rupture processes than one-point source. Primary results show that those synthetic signals calculated from three-point source agree well with the observations. This can further reveal that the source rupturing process of Wenchuan earthquake is a multi-rupture process, which is composed by at least three or more stages of rupture processes. In conclusion, the numerical simulation can not only solve some problems concluding the Earth’s ellipticity and anisotropy, which can be easily solved by conventional methods, but also finally solve the problems concluding topography model and lateral heterogeneity. We will try to find a way to fully implement self-gravitation in spectral element method in future, and do our best to continue researching the Earth’s free oscillations using the numerical simulations to see how the Earth’ lateral heterogeneous will affect the Earth’s free oscillations. These will make it possible to bring modal spectral data increasingly to bear on furthering our understanding of the Earth’s three-dimensional structure.

地震波场模拟及叠前深度偏移的李群方法

The processes of seismic wave propagation in phase space and one way wave extrapolation in frequency-space domain, if without dissipation, are essentially transformation under the action of one parameter Lie groups. Consequently, the numerical calculation methods of the propagation ought to be Lie group transformation too, which is known as Lie group method. After a fruitful study on the fast methods in matrix inversion, some of the Lie group methods in seismic numerical modeling and depth migration are presented here. Firstly the Lie group description and method of seismic wave propagation in phase space is proposed, which is, in other words, symplectic group description and method for seismic wave propagation, since symplectic group is a Lie subgroup and symplectic method is a special Lie group method. Under the frame of Hamiltonian, the propagation of seismic wave is a symplectic group transformation with one parameter and consequently, the numerical calculation methods of the propagation ought to be symplectic method. After discrete the wave field in time and phase space, many explicit, implicit and leap-frog symplectic schemes are deduced for numerical modeling. Compared to symplectic schemes, Finite difference (FD) method is an approximate of symplectic method. Consequently, explicit, implicit and leap-frog symplectic schemes and FD method are applied in the same conditions to get a wave field in constant velocity model, a synthetic model and Marmousi model. The result illustrates the potential power of the symplectic methods. As an application, symplectic method is employed to give synthetic seismic record of Qinghai foothills model. Another application is the development of Ray+symplectic reverse-time migration method. To make a reasonable balance between the computational efficiency and accuracy, we combine the multi-valued wave field & Green function algorithm with symplectic reverse time migration and thus develop a new ray+wave equation prestack depth migration method. Marmousi model data and Qinghai foothills model data are processed here. The result shows that our method is a better alternative to ray migration for complex structure imaging. Similarly, the extrapolation of one way wave in frequency-space domain is a Lie group transformation with one parameter Z and consequently, the numerical calculation methods of the extrapolation ought to be Lie group methods. After discrete the wave field in depth and space, the Lie group transformation has the form of matrix exponential and each approximation of it gives a Lie group algorithm. Though Pade symmetrical series approximation of matrix exponential gives a extrapolation method which is traditionally regarded as implicit FD migration, it benefits the theoretic and applying study of seismic imaging for it represent the depth extrapolation and migration method in a entirely different way. While, the technique of coordinates of second kind for the approximation of the matrix exponential begins a new way to develop migration operator. The inversion of matrix plays a vital role in the numerical migration method given by Pade symmetrical series approximation. The matrix has a Toepelitz structure with a helical boundary condition and is easy to inverse with LU decomposition. A efficient LU decomposition method is spectral factorization. That is, after the minimum phase correlative function of each array of matrix had be given by a spectral factorization method, all of the functions are arranged in a position according to its former location to get a lower triangular matrix. The major merit of LU decomposition with spectral factorization (SF Decomposition) is its efficiency in dealing with a large number of matrixes. After the setup of a table of the spectral factorization results of each array of matrix, the SF decomposition can give the lower triangular matrix by reading the table. However, the relationship among arrays is ignored in this method, which brings errors in decomposition method. Especially for numerical calculation in complex model, the errors is fatal. Direct elimination method can give the exact LU decomposition But even it is simplified in our case, the large number of decomposition cost unendurable computer time. A hybrid method is proposed here, which combines spectral factorization with direct elimination. Its decomposition errors is 10 times little than that of spectral factorization, and its decomposition speed is quite faster than that of direct elimination, especially in dealing with a large number of matrix. With the hybrid method, the 3D implicit migration can be expected to apply on real seismic data. Finally, the impulse response of 3D implicit migration operator is presented.

华南地壳岩石弹性和变形性质的高温高压实验研究

The South continent of China lies to southeast of Eurasia block. It is an active area from the view of crust growth and continent spread and is a transition zone between continental crust and oceanic crust. The compressional wave velocities and anisotropies of typical crustal metamorphic rocks were determined at high temperature (up to 1000 ℃) and high pressure(up to 800MPa). The experimental results show that the velocities generally increase with pressure, and is unaffected by temperature up to around 550 ℃. But the velocities of all experimental samples start to drop above a temperature point. For an example, this greatly reduce the speed of wave propagation in amphibolite and serpentinite above 760 ℃ and above 550 ℃ respectively, which may be due to dehydrate of amphibole and serpentine. P-wave anisotropy coefficients of those rocks range from 2% to 10% at 800MPa and 500 ℃. The anisotropies decrease with increasing pressure at room temperature, but hardly change as function of temperature at constant 800MPa or 600MPa pressure. The average velocity of the six crustal rocks is 6.28km/s under the condition of 800MPa and 550 ℃, which is consistent with the result of deep seismic sounding data. Based on this experimental result, we deduce there may exist a lot of felsic granulites and amphibolites at the depth of 15-25km underground. With increasing temperature and pressure, the deformation behavior of the rocks undergoes from localized brittle fracture, semi-brittle deformation (cataclastic flow or semi-brittle faulting, semi-brittle flow) to homogeneous crystal-plastic flow. This transition is associated with mechanical behavior and micro-mechanism. It is very important to understanding earthquake source mechanics, the strength of the lithosphere and the style of deformation. The experiments were conducted at temperature of 600-1000 ℃, confining pressure of 500MPa, and stain rates of 10~(-4)-10~(-6) S~(-1). For fine-grained natural amphibolite, the results of experiments show that brittle faulting is major failure mode at temperature <600 ℃, but crystal-plastic deformation is dominate at temperature >800 ℃, and there is a transition with increasing temperature from sembrittle faulting to cataclastic flow and sembrittle flow at temperature of 670-750 ℃. For medium-grained natural Felsic granulite, the results of experiments show that brittle faulting is major failure mode at temperature <500 ℃, but crystal-plastic deformation is dominate at temperature >700 ℃, and there is a transition with increasing temperature from semibrittle faulting to cataclastic flow and sembrittle flow at temperature of 500-600 ℃.

泥岩裂缝储层地震识别的理论和方法

During the exploration of fractured reservoirs, worldwide difficult problems will be encountered: how to locate the fractured zones, how to quantitatively determine the azimuth, density, and distribution of the fractures, and how to compute the permeability and porosity of the fractures. In an endeavor to solve these problems, the fractured shale reservoir in SiKou area of ShengLi oil field was chosen as a study area. A study of seismic predictive theory and methods for solving problems encountered in fractured reservoir exploration are examined herein. Building on widely used current fractured reservoir exploration techniques, new seismic theories and methods focusing on wave propagation principles in anisotropic medium are proposed. Additionally, integrated new seismic data acquisition and processing methods are proposed. Based on research and application of RVA and WA methods from earlier research, a new method of acoustic impedance varying with azimuth (IPVA) creatively is put forth. Lastly combining drilling data, well log data, and geologic data, an integrated seismic predictive method for cracked reservoir bed was formed. A summary of the six parts of research work of this paper is outlined below. In part one, conventional geologic and geophysical prediction methods etc. for cracked reservoir exploration are examined, and the weaknesses of these approaches discussed. In part two, seismic wave propagation principles in cracked reservoirs are studied. The wave equation of seismic velocity and attenuation factor in three kinds of fracture mediums is induced, and the azimuth anisotropy of velocity and attenuation in fracture mediums is determined. In part three, building on the research and application of AVA and WA methods by a former researcher, a new method of acoustic impedance creatively varying with azimuth (IPVA) is introduced. A practical software package utilizing this technique is also introduced. In part four, Base on previously discussed theory, first a large full azimuth 3d seismic data (70km~2) was designed and acquired. Next, the volume was processed with conventional processing sequence. Then AVA, WA, and IPVA processing was applied, and finally the azimuth and density of the fractures were quantitatively determined by an integrated method. Predictions were supported by well data that indicate the approach is highly reliable. in part five, geological conditions contributing to cracked reservoir bed formation are analyzed in the LuoJia area resulting in the discovery that the main fractured zones are related to fault distribution in the basin, that also control the accumulation of the oil and gas, the generation mechanisms and types of fractured shale reservoirs are studied. Lastly, by using full 3D seismic attributes, azimuth and density of cracked reservoir zones are successfully quantitative predicted. Using an integrated approach that incorporates seismic, geologic and well log data, the best two fractured oil prospects in LouJia area are proposed. These results herein represent a break through in seismic technology, integrated seismic predictive theory, and production technology for fractured reservoirs. The approach fills a void that can be applied both inside China, and internationally. Importantly, this technique opens a new exploration play in the ShengLi oil field that while difficult has substantial potential. Properly applied, this approach could play an important role toward stabilizing the oil field' production. In addition, this technique could be extended fracture exploration in other oil fields producing substantial economic reward.