准噶尔盆地西北缘三叠系层序地层及有利勘探目标选择

Based on outcrop, borehole, seismic and regional geological data, the sequence stratigraphy, sedimentary facies of the Triassic in the western margin of the Zhugaer basin was studied, and favorable exploration target was forecasted. The major achievements include: (1) the Triassic in the western margin of the Zhugaer basin can be divided into 1 second-order sequence and 5 third-order sequences, which are, in ascending order, TSQ1, TSQ2, TSQ3, TSQ4, and TSQ5. TSQ1 is equivolent to Baikouquan formation, TSQ2 is equivolent to lower Kelamayi formation, TSQ3 is equivolent to upper Kelamai formation, TSQ4 is equivolent to lower and middle Baijiantan formation, and TSQ5 is equivolent to upper Baijiantan formation. Each sequence is divided into transgressive and regressive system tracts. Thus the sequence correlation framework is established. (2) The factors controlling development of sequences are analyzed, and it is believed that tectonic is the major controlling factor. Model of sequence development is summarized. (3)Through study on sedimentary facies, 6 types of facies are recognized: alluvial fan, fan delta, braided river, braided delta, delta and lake. Their microfacies are also recognized. In this study, it is proposed that the upper and lower Kelamayi formation（TSQ2、 TSQ3）is deposited by braided river instead of alluvial fan. This conclusion is of important theoretical and practical significance.(4) The sedimentary facies map of each sequence is compiled, and the sedimentary facies developed in each sequence is determined. In TSQ1, the sedimentary facies developed is alluvial fan and fan delta. In TSQ2, the sedimentary facies developed is mainly alluvial fan and fan delta in the north, and braided river and braided delta in the south. In TSQ3, the sedimentary facies developed is mainly braided river and braided delta. In TSQ4, the sedimentary facies developed is mainly braided delta in the north, and meandering delta in the south. In TSQ5, the sedimentary facies developed is mainly braided river and braided delta. (5) In the framework of sequence stratigrahpy, favorable areas for concealed traps are forecasted, and different types of traps are developed in different system tracts. (6) Favorable areas for future exploration are predicted.

砂砾岩油藏精细油藏描述–以准噶尔盆地西北缘八区下乌尔禾组油藏为例

Glutenite reservoir is one of the most important reservoir types in china. Because of its particularity of rock structure and pore structure, it is usually difficult in development, especially for its serious heterogeneity. On the basis of seismic, well logs, core data and production performance, the lower Wuerhe group can be divided into one second-order sequences, two third-order sequences and twenty two subsequences, corresponding to the five stages and twenty two minlayers. In addition, the fault systems are interpreted and the control action of fault systems to reservoir development is also described. The lower Wuerhe formation of 8th district belongs to fluvial-dominated fan delta sedimentation, according to the analysis of well logs, logging data and core data. It can be subdivided into two kinds of subfacies and nine kinds of microfacies. The fan delta plain subfacies mainly consist of braided channel, unconcentrated flow, mud flow and sieve deposit microfacies. The fan delta front subfacies include subaqueous distributary channel, subaqueous interdistributary channel, debris flow, subaqueous barrier and grain flow microfacies. Combined with the regional geological characteristics, the porosity model of lower Wuerhe formation is performed using core data. A permeability model based on the flow zone index is also formed according to the pore throat characteristics and flow property. Finally, the heterogeneity is analyzed. The result shows that the lower Wuerhe formation has a feature of middle-high heterogeneity, and it is controlled by material sources and sedimentary facies belt.

南海北部陆缘深水区构造－热演化与烃源岩成熟度史研究

Abstract In order to provide basic data for evaluation of the petroleum potential in the deep water area of the northern margin of the South China Sea (SCS), present-day thermal regime and basin tectonothermal evolution are reconstructed and the maturation history of the Cenozoic major source rocks in the study area is derived. The present-day geothermal regime in the deep water area of the northern margin of SCS is defined according to the geothermal gradient, thermal properties and heat flow data. Tectonic subsidence history is reconstructed based on borehole and seismic data, and accordingly the stretching episodes are determined from the subsidence pattern. Heat flow history in the deep water area of the northern margin of SCS is estimated on a finite time, laterally non-uniform and multi-episode stretching model. Maturation history of the main source rocks in the study area is estimated through EASYRo% kinetic model and thermal history, and the potential of petroleum in the deep water area of the northern margin of SCS is evaluated based on the data above. The results show that the present-day geothermal regime in the deep water area of the northern margin of SCS is characterized by “hot basin” with high geothermal gradient (39.1±7.4℃/km) and high heat flow (77.5±14.8 mW/m2), and that the Qiongdongnan Basin (QDNB) underwent three stretching episodes and consequently suffered three heating episodes (Eocene, Oligocene and Pliocene time) with highest paleo-heat flow of 65～90 mW/m2 at the end of the Pliocene, that the Pearl River Mouth Basin (PRMB) two stretching and two heating episodes (Eocene, Oligocene time) with highest paleo-heat flow of 60～70 mW/m2 at the end of the Oligocene, and that the source rocks matured drastically responding to the heating episodes. There are four hydrocarbon generation kitchens in the deep water area of the northern margin of SCS which are favor of its bright petroleum perspective. Tectonothermal analysis indicates that the present-day geothermal regime which is characterized with “hot basin” in the deep water area of the PRMB resulted mainly from the Cenozoic stretching as well as faulting and magmatic activities during the Neotectonic period, and that the Pliocene heating episode of the QDNB is coupled with the transition from sinistral to dextral gliding of the Red Rive fault, and that the deep water basins in the northern margin of SCS are typical of multiple rifting which caused multi-episode heating process.

边坡工程的极限分析上限法

Theory of limit analysis include upper bound theorem and lower bound theorem. To deal with slope stability analysis by limit analysis is to approximate the real solution from upper limit and lower limit. The most used method of limit analysis is upper bound theorem, therefore it is often applied to slope engineering in many cases. Although upper bound approach of limit analysis can keep away from vague constitutive relation and complex stress analyses, it also can obtain rigorous result. Assuming the critical surface is circular slip surface, two kinematically admissible velocity fields for perpendicular slice method and radial slice method can be established according to the limit analysis of upper bound theorem. By means of virtual work rate equation and strength reduction method, the upper-bound solution of limit analysis for homogeneous soil slope can be obtained. A log-spiral rotational failure mechanism for homogeneous slope is discussed from two different conditions which represent the position of shear crack passing the toe and below the toe. In the dissertition, the author also establishes a rotational failure mechanics with combination of different logarithmic spiral arcs. Furthermore, the calculation formula of upper bound solution for inhomogeneous soil slope stability problem can be deduced based on the upper bound approach of rigid elements. Through calculating the external work rate caused by soil nail, anti-slide pile, geotechnological grid and retaining wall, the upper bound solution of safety factor of soil nail structure slope, slip resistance of anti-slide pile, critical height of reinforced soil slope and active earth pressure of retaining wall can be obtained by upper bound limit analysis method. Taking accumulated body slope as subject investigated, with study on the limit analysis method to calculate slope safety factor, the kinematically admissible velocity fields of perpendicular slice method for slope with broken slip surface is proposed. Through calculating not only the energy dissipation rate produced in the broken slip surfaces and the vertical velocity discontinuity, but also the work rate produced by self-weight and external load, the upper bound solution of slope with broken slip surface is deduced. As a case study, the slope stability of the Sanmashan landslide in the area of the Three Gorges reservoir is analyzed. Based on the theory of limit analysis, the upper bound solution for rock slope with planar failure surface is obtained. By means of virtual work-rate equation, energy dissipation caused by dislocation of thin-layer and terrane can be calculated; furthermore, the formulas of safety factor for upper bound approach of limit analysis can be deduced. In the end, a new computational model of stability analysis for anchored rock slope is presented after taking into consideration the supporting effect of rock-bolts, the action of seismic force and fissure water pressure. By using the model, not only the external woke-rate done by self-weight, seismic force, fissure water pressure and anchorage force but also the internal energy dissipation produced in the slip surface and structural planes can be totally calculated. According to the condition of virtual work rate equation in limit state, the formula of safety factor for upper bound limit analysis can be deduced.

大型反倾岩体变形破坏机制及稳定性分析———以金沙江宽谷坝址龙蟠边坡为例

A large number of catastrophic accidents were aroused by the instability and destruction of anti-dip rock masses in the worldwide engineering projects, such as hydropower station, mine, railways and so on. Problems in relation to deformation and failure about anti-dip rock slopes are significant for engineering geology research. This dissertation takes the Longpan slope in the Jinsha River as a case to study the deformation mechanism of large-scale anti-dip rock masses and the slope stability analysis method. The primary conclusions are as follows. The Dale Reach of Jinsha River, from Longpan to the debouchment of Chongjiang tributary, is located in the southeastern margin of the Qinghai-Tibet Plateau. Longpan slope is the right embankment of Dale dam, it is only 26 km to the Shigu and 18 km to Tiger Leaping Gorge. The areal geology tectonic structures here area are complicated and blurry. Base on the information of geophysical exploration (CSAMT and seismology) and engineering geological investigation, the perdue tectonic pattern of Dale Reach is put forward for the first time in this paper. Due to the reverse slip of Longpan fault and normal left-rotation of Baihanchang fault, the old faulted valley came into being. The thick riverbed sediments have layered characters of different components and corresponding causes, which attribute to the sedimentary environments according with the new tectonic movements such as periodic mountain uplifting in middle Pleistocene. Longpan slope consists of anti-dip alternate sandstone and slate stratums, and the deformable volume is 6.5×107m3 approximately. It was taken for an ancient landslide or toppling failure in the past so that Dale dam became a vexed question. Through the latest field surveying, displacement monitoring and rock masses deforming characters analyses, the geological mechanism is actually a deep-seated gravitational bending deformation. And then the discrete element method is used to simulate the deforming evolution process, the conclusion accords very well with the geo-mechanical patterns analyses. In addition strength reduction method based on DEM is introduced to evaluate the factor of safety of anti-dip rock slope, and in accordance with the expansion way of the shear yielding zones, the progressive shear failure mechanism of large-scale anti-dip rock masses is proposed for the first time. As an embankment or a close reservoir bank to the lower dam, the stability of Longpan slope especially whether or not resulting in sliding with high velocity and activating water waves is a key question for engineering design. In fact it is difficult to decide the unified slip surface of anti-dip rock slope for traditional methods. The author takes the shear yielding zones acquired form the discrete element strength reduction calculation as the potential sliding surface and then evaluates the change of excess pore pressure and factor of stability of the slope generated by rapid drawdown of ponded water. At the same time the dynamic response of the slope under seismic loading is simulated through DEM numerical modeling, the following results are obtained. Firstly the effective effect of seismic inertia force is resulting in accumulation of shear stresses. Secondly the discontinuous structures are crucial to wave transmission. Thirdly the ultimate dynamic response of slope system takes place at the initial period of seismic loading. Lastly but essentially the effect of earthquake load to bringing on deformation and failure of rock slope is the coupling effect of shear stresses and excess pore water pressure accumulation. In view of limitations in searching the critical slip surface of rock slope of the existing domestic and international software for limit equilibrium slope stability analyses, this article proposes a new method named GA-Sarma Algorithm for rock slope stability analyses. Just as its name implies, GA-Sarma Algorithm bases on Genetic Algorithm and Sarma method. GA-Sarma Algorithm assumes the morphology of slip surface to be a broken line with traceability to extend along the discontinuous surface structures, and the slice boundaries is consistent with rock mass discontinuities such as rock layers, faults, cracks, and so on. GA-Sarma Algorithm is revolutionary method that is suitable for global optimization of the critical slip surface for rock slopes. The topics and contents including in this dissertation are closely related to the difficulties in practice, the main conclusions have been authorized by the engineering design institute. The research work is very meaningful and useful for the engineering construction of Longpan hydropower station.

波动方程模拟与偏移的实用化方法研究

Along with the widespread and in-depth applications in petroleum prospecting and development, the seismic modeling and migration technologies are proposed with a higher requirement by oil industrial, and the related practical demand is getting more and more urgent. Based on theories of modeling and migration methods for wave equation, both related with velocity model, I thoroughly research and develop some methods for the goal of highly effective and practical in this dissertation. In the first part, this dissertation probes into the layout designing by wave equations modeling, focusing on the target-oriented layout designing method guided by wave equation modeling in complicated structure areas. It is implemented by using the fourth order staggered grid finite difference (FD) method in velocity-stress 2D acoustic wave equations plus perfectly matched layer (PML) absorbing boundary condition. To design target-oriented layout: (a) match the synthetic record on the surface with events of subsurface structures by analyzing the snapshots of theoretical model; (b) determine the shot-gather distance by tracking the events of target areas and measuring the receiving range when it reaches the surface; (c) restrict the range of valid shot-gather distance by drawing seismic windows in single shot records; (d) choose the best trace distance by comparing the resolution of prospecting targets from the simulated records with different trace distance. Eventually, we obtained the observation system parameters, which achieve the design requirements. In the second part, this dissertation presents the practical method to improve the 3D Fourier Finite Difference (FFD) migration, and carefully analyzes all the factors which influence 3D FFD migration’s efficiency. In which, one of the most important parameters of migration is the extrapolating step. This dissertation presents an efficient 3D FFD migration algorithm, which use FFD propagator to extrapolate wavefields over big layers, and use Born-Kirchhoff interpolator to image wavefields over small layers between the big ones. Finally, I show the effectiveness of this hybrid migration method by comparing migration results from 3D SEG/EAGE model with different methods.

华北及周边地区层析成像研究

Our study deals with the high resolution body wave tomography in North china and adjacent areas(30°N-43°N，100°E-130°E), where earthquakes occurred many times in history and has a very complicated geological structure. 6870 events recorded at 273 digital seismic stations from CDSN during 1996－2002 and stations settled by Seislab of IGCAS in Bohai Bay area, including 1382 local earthquakes and 5488 teleseismic earthquakes are used in this study. In the data we used, the average number of received stations is greater than 5, the error of picking up direct arrival time is 0.1－0.5s. Before the inversion, we use Checkerboard method to confirm the reliability of result of Local events; use Restoring Resolution Test to confirm the reliability of result of teleseismic events. We also analyzed the effect of different parameters in the inversion. Based the analysis above, the model used in this paper is divided into small blocks with a dimension of 0.33°in the latitude and longitude directions and 5km、15km、30km in depth, and initial velocity model. Using pseudobending method to calculate the ray traveling path, LSQR algorithm to inverse, finally, we got the body velocity images below 25km and above 480km in this area using Joint- inversion with local events and teleseismic events. We made the conclusion at last: (1)at top zone of the south of Sichuan Basin , there exits low velocity anomalies, below 40km is the high velocity zone extend to 300km; (2) Above the 40km of Ordos block exits low velocity zone, while below 40km until 240km, the high velocity anomalies are interlaced by low velocity anomalies. Below 300km, the anomalies are unclear any more; (3) On the whole, the velocity structure below 400km on the mantle transition zone of Eastern China area shows its changes from low velocity to high velocity.

基于广义Lippmann-Schwinger方程的地震波场模拟算法与应用

In the last several decades, due to the fast development of computer, numerical simulation has been an indispensable tool in scientific research. Numerical simulation methods which based on partial difference operators such as Finite Difference Method (FDM) and Finite Element Method (FEM) have been widely used. However, in the realm of seismology and seismic prospecting, one usually meets with geological models which have piece-wise heterogeneous structures as well as volume heterogeneities between layers, the continuity of displacement and stress across the irregular layers and seismic wave scattering induced by the perturbation of the volume usually bring in error when using conventional methods based on difference operators. The method discussed in this paper is based on elastic theory and integral theory. Seismic wave equation in the frequency domain is transformed into a generalized Lippmann-Schwinger equation, in which the seismic wavefield contributed by the background is expressed by the boundary integral equation and the scattering by the volume heterogeneities is considered. Boundary element-volume integral method based on this equation has advantages of Boundary Element Method (BEM), such as reducing one dimension of the model, explicit use the displacement and stress continuity across irregular interfaces, high precision, satisfying the boundary at infinite, etc. Also, this method could accurately simulate the seismic scattering by the volume heterogeneities. In this paper, the concrete Lippmann-Schwinger equation is specifically given according to the real geological models. Also, the complete coefficients of the non-smooth point for the integral equation are introduced. Because Boundary Element-Volume integral equation method uses fundamental solutions which are singular when the source point and the field are very close，both in the two dimensional and the three dimensional case, the treatment of the singular kernel affects the precision of this method. The method based on integral transform and integration by parts could treat the points on the boundary and inside the domain. It could transform the singular integral into an analytical one both in two dimensional and in three dimensional cases and thus it could eliminate the singularity. In order to analyze the elastic seismic wave scattering due to regional irregular topographies, the analytical solution for problems of this type is discussed and the analytical solution of P waves by multiple canyons is given. For the boundary reflection, the method used here is infinite boundary element absorbing boundary developed by a pervious researcher. The comparison between the analytical solutions and concrete numerical examples validate the efficiency of this method. We thoroughly discussed the sampling frequency in elastic wave simulation and find that, for a general case, three elements per wavelength is sufficient, however, when the problem is too complex, more elements per wavelength are necessary. Also, the seismic response in the frequency domain of the canyons with different types of random heterogeneities is illustrated. We analyzed the model of the random media, the horizontal and vertical correlation length, the standard deviation, and the dimensionless frequency how to affect the seismic wave amplification on the ground, and thus provide a basis for the choice of the parameter of random media during numerical simulation.

井间地震速度层析成像方法与软件系统

Crosshole Seismic tomography has been broadly studied and applied in the fields of resource exploration and engineering exploration because of its special observing manner and better resolution than normal seismic exploration. This thesis will state the theory and method of Crosshole Seismic tomography. Basing on the previous studies，the thesis studied the initial velocity model，ray-tracing method, and developed the three-dimension tomography software. All the cells that a ray passes through are of the same velocities if the paths from transmitters to receivers are straight. The cells that the each ray passes through are recorded, and rays that pass through each cell are calculated. The ray average velocity which passes through a cell is set as the cell velocity. Analogously we can make a initial node velocity model because the velocity sum is calculated on the all cells which own to a certain node, and the cell number is summed about each nodes，the ratio of the velocity sum to the all cells number is set as the node velocity. The inversion result from the initial node velocity model is better than that of the average velocity model. Ray-bending and Shortest Path for Rays (SPR) have shortcomings and limitations respectively. Using crooked rays obtained from SPR rather than straight lines as the starting point can not only avoid ray bending converging to the local minimum travel time path, but also settle the no smooth ray problem obtained by SPR. The hybrid method costs much computation time, which is roughly equal to the time that SPR expends. The Delphi development tool based on the Object Pascal language standard has an advantage of object-oriented. TDTOM (Three Dimensions Tomography) was developed by using Delphi from the DOS version. Improvement on the part of inversion was made, which bring faster convergence velocity. TDTOM can be used to do velocity tomography from the first arrival travel time of the seismic wave, and it has the good qualities of friendly user interface and convenient operation. TDTOM is used to reconstruct the velocity image for a set of crosshole data from Karamay Oil Field. The geological explanation is then given by comparing the inversion effects of different ray-tracing methods. High velocity zones mean the cover of oil reservoir, and low velocity zones correspond to the reservoir or the steam flooding layer.

旋转轴对称弱各向异性介质模型的射线追踪方法

Ray tracing is a rapid and effective method for wave field calculation. Not only in the field of seismic-wave theory, but also in the field of seismic inversion and migration imaging,the seismic ray tracing method has become one of the most important methods. In anisotropic media, group velocity and phase velocity have different propagation directions. The seismic wave propagates along the direction of group velocity , it does not depend on the direction of phase velocity. Ray angle is a complex function with respect to phase angle, it is difficult to measure and calculate. But most rocks are weak anisotropic, so the expression of phase velocity can be simplified greatly. Based on the approximate expression of phase velocity this thesis for rotating axisymmetric weak anisotropic media deduces an expression of the partial derivative of phase velocity and an expression of group velocity with the method of linear approximation. This paper uses the fourth order Runge-Kutta method together with the two-dimensional interpolation and linear interpolation to obtain the parameters of the physical locations. At last the paths of seismic wave in rotating axisymmetric weak anisotropic media are computed. According to the analysis of the computational results, it indicates that the method developed in this paper has strong adaptability, high computational efficiency and high accuracy for rotating axisymmetric weak anisotropic media.

最优烃类指标方法的研究与应用

Nowadays, with the development of reservoir exploration, the method of exploration is growing. Together with well information and laboratory results, seismic information with high quality can predict reservoir successfully. Hydrocarbon Indicator is a method, which picks the most sensitive rock properties of hydrocarbons, scans the aim area with rock physics tools, and then indicates the area of reservoir. Obviously, the more is the difference between brine and oil/gas, the better this method works. Which parameter can be used as the Optimal Hydrocarbon Indicators is still in discussion. The author introduced several kinds of Hydrocarbon Indicators in this thesis. After analyzing the response of different parameter to reservoir, together with seismic information, the reservoir can be predicted. In this paper, the reservoir of Zhunge’er is studied to prove this kind of method is suitable for real exploration in China. Besides, the author chose Haila’er reservoir to testify whether this method could be used in metamorphic reservoir other than sandstone reservoir. The results highlighted the meaning of Optimal Hydrocarbon Indicators on reservoir identification. The author also mentioned some thoughts for the development of hydrocarbon indicators in the future.

贝尔凹陷苏德尔特构造带古潜山构造与储层研究

With the develop ment of oil and gas exploration, the conventional struc ture exploration era has gradually been substituted by the concealed reser voir exploration technology. Hill poll becomes one of the most important areas in the future exploration. This paper is based on the three-dimensional seismic interpretation of Sudeerte structure. In terms of the overall character istics of Sudeerte structure, we use the coherent cube and the time slice to interpret the fault in plane. For the interpretation，we especially used the well to adjust the interpretation results. The results of seismic attribution analysis, spectrum decomposition and post- stack seismic inversion forecast that hill pools reservoir are dist ributed in several bands along the north-northeast to northeast-east. Xing'anling Group shows that the potential reser voirs are mainly distributed along Bei 14 –Bei 40 and De (99-212) - Bei (16-1) and Budate Group distributed along Bei 14 –Bei 40 in northwestern direction and De (99-212) - Bei (16-1) in north eastern direction. At the same time, by analyzing the structure and the reservoir, and combining with other data, three models are built. The characteristics of reservoirs dist ribution are concl uded，and potent ional favorable exploration dire ctions are predi cted.

数据驱动型多次波衰减方法的研究

In the practical seismic profile multiple reflections tend to impede the task of even the experienced interpreter in deducing information from the reflection data. Surface multiples are usually much stronger, more broadband, and more of a problem than internal multiples because the reflection coefficient at the water surface is much larger than the reflection coefficients found in the subsurface. For this reason most attempts to remove multiples from marine data focus on surface multiples, as will I. A surface-related multiple attenuation method can be formulated as an iterative procedure. In this essay a fully data-driven approach which is called MPI —multiple prediction through inversion (Wang, 2003) is applied to a real marine seismic data example. This is a pretty promising scheme for predicting a relative accurate multiple model by updating the multiple model iteratively, as we usually do in a linearized inverse problem. The prominent characteristic of MPI method lie in that it eliminate the need for an explicit surface operator which means it can model the multiple wavefield without any knowledge of surface and subsurface structures even a source signature. Another key feature of this scheme is that it can predict multiples not only in time but also in phase and in amplitude domain. According to the real data experiments it is shown that this scheme for multiple prediction can be made very efficient if a good initial estimate of the multiple-free data set can be provided in the first iteration. In the other core step which is multiple subtraction we use an expanded multi-channel matching filter to fulfil this aim. Compared to a normal multichannel matching filter where an original seismic trace is matched by a group of multiple-model traces, in EMCM filter a seismic trace is matched by not only a group of the ordinary multiple-model traces but also their adjoints generated mathematically. The adjoints of a multiple-model trace include its first derivative, its Hilbert transform and the derivative of the Hilbert transform. The third chapter of the thesis is the application for the real data using the previous methods we put forward from which we can obviously find the effectivity and prospect of the value in use. For this specific case I have done three group experiments to test the effectiveness of MPI method, compare different subtraction results with fixed filter length but different window length, invest the influence of the initial subtraction result for MPI method. In terms of the real data application, we do fine that the initial demultiple estimate take on a great deal of influence for the MPI method. Then two approaches are introduced to refine the intial demultiple estimate which are first arrival and masking filter respectively. In the last part some conclusions are drawn in terms of the previous results I have got.

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

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.