首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 31 毫秒
1.
The horizontal transversely isotropic model, with arbitrary symmetry axis orientation, is the simplest effective representative that explains the azimuthal behaviour of seismic data. Estimating the anisotropy parameters of this model is important in reservoir characterisation, specifically in terms of fracture delineation. We propose a travel‐time‐based approach to estimate the anellipticity parameter η and the symmetry axis azimuth ? of a horizontal transversely isotropic medium, given an inhomogeneous elliptic background model (which might be obtained from velocity analysis and well velocities). This is accomplished through a Taylor's series expansion of the travel‐time solution (of the eikonal equation) as a function of parameter η and azimuth angle ?. The accuracy of the travel time expansion is enhanced by the use of Shanks transform. This results in an accurate approximation of the solution of the non‐linear eikonal equation and provides a mechanism to scan simultaneously for the best fitting effective parameters η and ?, without the need for repetitive modelling of travel times. The analysis of the travel time sensitivity to parameters η and ? reveals that travel times are more sensitive to η than to the symmetry axis azimuth ?. Thus, η is better constrained from travel times than the azimuth. Moreover, the two‐parameter scan in the homogeneous case shows that errors in the background model affect the estimation of η and ? differently. While a gradual increase in errors in the background model leads to increasing errors in η, inaccuracies in ?, on the other hand, depend on the background model errors. We also propose a layer‐stripping method valid for a stack of arbitrary oriented symmetry axis horizontal transversely isotropic layers to convert the effective parameters to the interval layer values.  相似文献   

2.
While velocity contrasts are responsible for most of the events recorded in our data, the long wavelength behavior of the velocity model is responsible for the geometrical shape of these events. For isotropic acoustic materials, the wave dependency on the long (wave propagation) and short (scattering) wavelength velocity components is stationary with the propagation angle. On the other hand, in representing a transversely isotropic with a vertical symmetry axis medium with the normal moveout velocity, the anellepticity parameter η, the vertical scaling parameter δ, and the sensitivity of waves vary with the polar angle for both the long and short wavelength features of the anisotropic dimensionless medium parameters (δ and η). For horizontal reflectors at reasonable depths, the long wavelength features of the η model is reasonably constrained by the long offsets, whereas the short wavelength features produce very week reflections at even reasonable offsets. Thus, for surface acquired seismic data, we could mainly invert for smooth η responsible for the geometrical shape of reflections. On the other hand, while the δ long wavelength components mildly affects the recorded data, its short wavelength variations can produce reflections at even zero offset, with a behavior pattern synonymous to density. The lack of the long wavelength δ information will mildly effect focusing but will cause misplacement of events in depth. With low enough frequencies (very low), we may be able to recover the long wavelength δ using full waveform inversion. However, unlike velocity, the frequencies needed for that should be ultra‐low to produce long‐wavelength scattering‐based model information as δ perturbations do not exert scattering at large offsets. For a combination given by the horizontal velocity, η, and ε, the diving wave influence of η is absorbed by the horizontal velocity, severely limiting the η influence on the data and full waveform inversion. As a result, with a good smooth η estimation, for example, from tomography, we can focus the full waveform inversion to invert for only the horizontal velocity and maybe ε as a parameter to fit the amplitude. This is possibly the most practical parametrization for inversion of surface seismic data in transversely isotropic with vertical symmetry axis media.  相似文献   

3.
We address the issue of linearity and scale dependence in forward modelling of seismic data from well logs, for large ray parameters, wide angles or large offsets. We present a forward model, within the context of seismic‐to‐well matching, that is linearized in the elastic properties of the earth. This model preserves linearity at large ray parameters and can handle fine‐layering effects such as induced anisotropy. Starting from a low‐contrast small‐ray‐parameter model, we extend it to a large‐ray‐parameter model by fully linearizing the elastic‐property contrasts. Overall linearity of the forward model is extended by partitioning the compressional‐wave and shear‐wave velocity fields into two fundamental scales: a kinematic scale that governs wavefield propagation effects and a dynamic scale that governs wavefield scattering effects. This analysis reveals that the standard practice in forward modelling of strongly filtering the ratio of compressional‐wave velocity to shear‐wave velocity is well founded in the underlying physics. The partitioning of the velocity fields also leads naturally to forward modelling that accounts fully for stretch effects, to resolution of the angle‐of‐incidence versus ray‐parameter dichotomy in seismic‐amplitude analysis, and to full accounting for induced anisotropy and dispersion effects due to fine‐layering of isotropic media. With the onset of routine long‐offset acquisition and the compelling need to optimize asset management in order to maximize reserve recovery, this forward model recognizes the physics of seismic wave propagation and enables a more complete exploitation of amplitude information in pre‐critical seismic data.  相似文献   

4.
5.
Most sedimentary rocks are anisotropic, yet it is often difficult to accurately incorporate anisotropy into seismic workflows because analysis of anisotropy requires knowledge of a number of parameters that are difficult to estimate from standard seismic data. In this study, we provide a methodology to infer azimuthal P‐wave anisotropy from S‐wave anisotropy calculated from log or vertical seismic profile data. This methodology involves a number of steps. First, we compute the azimuthal P‐wave anisotropy in the dry medium as a function of the azimuthal S‐wave anisotropy using a rock physics model, which accounts for the stress dependency of seismic wave velocities in dry isotropic elastic media subjected to triaxial compression. Once the P‐wave anisotropy in the dry medium is known, we use the anisotropic Gassmann equations to estimate the anisotropy of the saturated medium. We test this workflow on the log data acquired in the North West Shelf of Australia, where azimuthal anisotropy is likely caused by large differences between minimum and maximum horizontal stresses. The obtained results are compared to azimuthal P‐wave anisotropy obtained via orthorhombic tomography in the same area. In the clean sandstone layers, anisotropy parameters obtained by both methods are fairly consistent. In the shale and shaly sandstone layers, however, there is a significant discrepancy between results since the stress‐induced anisotropy model we use is not applicable to rocks exhibiting intrinsic anisotropy. This methodology could be useful for building the initial anisotropic velocity model for imaging, which is to be refined through migration velocity analysis.  相似文献   

6.
This paper presents a new explicit method for the estimation of layered vertical transverse isotropic (VTI) anisotropic parameters from walkaway VSP data. This method is based on Dix‐type normal moveout (NMO) inversion. To estimate interval anisotropic parameters above a receiver array, the method uses time arrivals of surface‐related double‐reflected downgoing waves. A three‐term NMO approximation function is used to estimate NMO velocity and a non‐hyperbolic parameter. Assuming the vertical velocity is known from zero‐offset VSP data, Dix‐type inversion is applied to estimate the layered Thomsen anisotropic parameters ?, δ above the receivers array. Model results show reasonable accuracy for estimates through Dix‐type inversion. Results also show that in many cases we can neglect the influence of the velocity gradient on anisotropy estimates. First breaks are used to estimate anisotropic parameters within the walkaway receiver interval. Analytical uncertainty analysis is performed to NMO parameter estimates. Its conclusions are confirmed by modelling.  相似文献   

7.
To better understand (and correct for) the factors affecting the estimation of attenuation (Q), we simulate subsurface wave propagation with the Weyl/Sommerfeld integral. The complete spherical wavefield emanating from a P‐wave point source surrounded by a homogeneous, isotropic and attenuative medium is thus computed. In a resulting synthetic vertical seismic profile, we observe near‐field and far‐field responses and a 90° phase rotation between them. Depth dependence of the magnitude spectra in these two depth regions is distinctly different. The logarithm of the magnitude spectra shows a linear dependence on frequency in the far‐field but not in those depth regions where the near‐field becomes significant. Near‐field effects are one possible explanation for large positive and even negative Q‐factors in the shallow section that may be estimated from real vertical seismic profile data when applying the spectral ratio method. We outline a near‐field compensation technique that can reduce errors in the resultant Q estimates.  相似文献   

8.
A major complication caused by anisotropy in velocity analysis and imaging is the uncertainty in estimating the vertical velocity and depth scale of the model from surface data. For laterally homogeneous VTI (transversely isotropic with a vertical symmetry axis) media above the target reflector, P‐wave moveout has to be combined with other information (e.g. borehole data or converted waves) to build velocity models for depth imaging. The presence of lateral heterogeneity in the overburden creates the dependence of P‐wave reflection data on all three relevant parameters (the vertical velocity VP0 and the Thomsen coefficients ε and δ) and, therefore, may help to determine the depth scale of the velocity field. Here, we propose a tomographic algorithm designed to invert NMO ellipses (obtained from azimuthally varying stacking velocities) and zero‐offset traveltimes of P‐waves for the parameters of homogeneous VTI layers separated by either plane dipping or curved interfaces. For plane non‐intersecting layer boundaries, the interval parameters cannot be recovered from P‐wave moveout in a unique way. Nonetheless, if the reflectors have sufficiently different azimuths, a priori knowledge of any single interval parameter makes it possible to reconstruct the whole model in depth. For example, the parameter estimation becomes unique if the subsurface layer is known to be isotropic. In the case of 2D inversion on the dip line of co‐orientated reflectors, it is necessary to specify one parameter (e.g. the vertical velocity) per layer. Despite the higher complexity of models with curved interfaces, the increased angle coverage of reflected rays helps to resolve the trade‐offs between the medium parameters. Singular value decomposition (SVD) shows that in the presence of sufficient interface curvature all parameters needed for anisotropic depth processing can be obtained solely from conventional‐spread P‐wave moveout. By performing tests on noise‐contaminated data we demonstrate that the tomographic inversion procedure reconstructs both the interfaces and the VTI parameters with high accuracy. Both SVD analysis and moveout inversion are implemented using an efficient modelling technique based on the theory of NMO‐velocity surfaces generalized for wave propagation through curved interfaces.  相似文献   

9.
Reflection full waveform inversion can update subsurface velocity structure of the deeper part, but tends to get stuck in the local minima associated with the waveform misfit function. These local minima cause cycle skipping if the initial background velocity model is far from the true model. Since conventional reflection full waveform inversion using two‐way wave equation in time domain is computationally expensive and consumes a large amount of memory, we implement a correlation‐based reflection waveform inversion using one‐way wave equations to retrieve the background velocity. In this method, one‐way wave equations are used for the seismic wave forward modelling, migration/de‐migration and the gradient computation of objective function in frequency domain. Compared with the method using two‐way wave equation, the proposed method benefits from the lower computational cost of one‐way wave equations without significant accuracy reduction in the cases without steep dips. It also largely reduces the memory requirement by an order of magnitude than implementation using two‐way wave equation both for two‐ and three‐dimensional situations. Through numerical analysis, we also find that one‐way wave equations can better construct the low wavenumber reflection wavepath without producing high‐amplitude short‐wavelength components near the image points in the reflection full waveform inversion gradient. Synthetic test and real data application show that the proposed method efficiently updates the background velocity model.  相似文献   

10.
Surface waves in seismic data are often dominant in a land or shallow‐water environment. Separating them from primaries is of great importance either for removing them as noise for reservoir imaging and characterization or for extracting them as signal for near‐surface characterization. However, their complex properties make the surface‐wave separation significantly challenging in seismic processing. To address the challenges, we propose a method of three‐dimensional surface‐wave estimation and separation using an iterative closed‐loop approach. The closed loop contains a relatively simple forward model of surface waves and adaptive subtraction of the forward‐modelled surface waves from the observed surface waves, making it possible to evaluate the residual between them. In this approach, the surface‐wave model is parameterized by the frequency‐dependent slowness and source properties for each surface‐wave mode. The optimal parameters are estimated in such a way that the residual is minimized and, consequently, this approach solves the inverse problem. Through real data examples, we demonstrate that the proposed method successfully estimates the surface waves and separates them out from the seismic data. In addition, it is demonstrated that our method can also be applied to undersampled, irregularly sampled, and blended seismic data.  相似文献   

11.
Refracted arrivals are analysed to estimate the near‐surface anisotropy of marine sediments using a vertical‐cable (VC) configuration. In the presence of dip, the horizontal and vertical ray‐slownesses are obtained from the observed apparent slownesses in the up‐ and downdip directions using a sum or difference at each azimuth. The multiple azimuths generated by a VC geometry permit the ray‐slowness distribution of the marine sediments to be determined. An inversion procedure is developed to provide dip and anisotropy parameters for refractive layers from the measured refraction traveltimes in multilayered azimuthally isotropic and anisotropic media. Two sets of transversely isotropic models are used to analyse the azimuthal variations of apparent and ray slownesses. In the first set, we fix the anisotropic parameters of the models but vary the dip (0°, 5° and 10°) to test the effects of the presence of dip. In the second set, we vary the P‐wave anisotropy strength (5.2%, 10.3%, 15.8% and 22.0%) to examine the sensitivity and accuracy of ray‐slowness approximations which are independent of dip. We test this inversion procedure on synthetic P‐wave VC data calculated for six different models by a finite‐difference method. The results of applications to real VC data acquired from the North Sea are also presented.  相似文献   

12.
River water temperature is an important water quality parameter that also influences most aquatic life. Physical processes influencing water temperature in rivers are highly complex. This is especially true for the estimation of river heat exchange processes that are highly dependent on good estimates of radiation fluxes. Furthermore, very few studies were found within the stream temperature dynamic literature where the different radiation components have been measured and compared at the stream level (at microclimate conditions). Therefore, this study presents results on hydrometeorological conditions for a small tributary within Catamaran Brook (part of the Miramichi River system, New Brunswick, Canada) with the following specific objectives: (1) to compare between stream microclimate and remote meteorological conditions, (2) to compare measured long‐wave radiation data with those calculated from an analytical model, and (3), to calculate the corresponding river heat fluxes. The most salient findings of this study are (1) solar radiation and wind speed are parameters that are highly site specific within the river environment and play an important role in the estimation of river heat fluxes; (2) the incoming, outgoing, and net long‐wave radiation within the stream environment (under the forest canopy) can be effectively calculated using empirical formula; (3) at the study site more than 80% of the incoming long‐wave radiation was coming from the forest; (4) total energy gains were dominated by solar radiation flux (for all the study periods) followed by the net long‐wave radiation (during some periods) whereas energy losses were coming from both the net long‐wave radiation and evaporation. Conductive heat fluxes have a minor contribution from the overall heat budget (<3·5%); (5) the reflected short‐wave radiation at the water surface was calculated on average as 3·2%, which is consistent with literature values. Results of this study contribute towards a better understanding of river heat fluxes and water temperature models as well as for more effective aquatic resources and fisheries management. Copyright © 2011 John Wiley & Sons, Ltd.  相似文献   

13.
The double‐square‐root equation is commonly used to image data by downward continuation using one‐way depth extrapolation methods. A two‐way time extrapolation of the double‐square‐root‐derived phase operator allows for up and downgoing wavefields but suffers from an essential singularity for horizontally travelling waves. This singularity is also associated with an anisotropic version of the double‐square‐root extrapolator. Perturbation theory allows us to separate the isotropic contribution, as well as the singularity, from the anisotropic contribution to the operator. As a result, the anisotropic residual operator is free from such singularities and can be applied as a stand alone operator to correct for anisotropy. We can apply the residual anisotropy operator even if the original prestack wavefield was obtained using, for example, reverse‐time migration. The residual correction is also useful for anisotropic parameter estimation. Applications to synthetic data demonstrate the accuracy of the new prestack modelling and migration approach. It also proves useful in approximately imaging the Vertical Transverse Isotropic Marmousi model.  相似文献   

14.
In areas of complex geology such as the Canadian Foothills, the effects of anisotropy are apparent in seismic data and estimation of anisotropic parameters for use in seismic imaging is not a trivial task. Here we explore the applicability of common‐focus point (CFP)‐based velocity analysis to estimate anisotropic parameters for the variably tilted shale thrust sheet in the Canadian Foothills model. To avoid the inherent velocity‐depth ambiguity, we assume that the elastic properties of thrust‐sheet with respect to transverse isotropy symmetry axis are homogeneous, the reflector below the thrust‐sheet is flat, and that the anisotropy is weak. In our CFP approach to velocity analysis, for a poorly imaged reflection point, a traveltime residual is obtained as the time difference between the focusing operator for an assumed subsurface velocity model and the corresponding CFP response obtained from the reflection data. We assume that this residual is due to unknown values for anisotropy, and we perform an iterative linear inversion to obtain new model parameters that minimize the residuals. Migration of the data using parameters obtained from our inversion results in a correctly positioned and better focused reflector below the thrust sheet. For traveltime computation we use a brute force mapping scheme that takes into account weakly tilted transverse isotropy media. For inversion, the problem is set up as a generalized Newton's equation where traveltime error (differential time shift) is linearly dependent on the parameter updates. The iterative updates of parameters are obtained by a least‐squares solution of Newton's equations. The significance of this work lies in its applicability to areas where transverse isotropy layers are heterogeneous laterally, and where transverse isotropy layers are overlain by complex structures that preclude a moveout curve fitting.  相似文献   

15.
In 2005, a multicomponent ocean bottom node data set was collected by BP and BHP Billiton in the Atlantis field in the Gulf of Mexico. Our results are based on data from a few sparse nodes with millions of shots that were analysed as common receiver azimuthal gathers. A first‐order look at P‐wave arrivals on a common receiver gather at a constant offset reveals variation of P‐wave arrival time as a function of azimuth indicating the presence of azimuthal anisotropy at the top few layers. This prompted us to investigate shear arrivals on the horizontal component data. After preliminary processing, including a static correction, the data were optimally rotated to radial (R) and transverse (T) components. The R component shows azimuthal variation of traveltime indicating variation of velocity with azimuth; the corresponding T component shows azimuthal variation of amplitude and phase (polarity reversal). The observed shear‐wave (S‐wave) splitting, previously observed azimuthal P‐wave velocity variation and azimuthal P‐wave amplitude variation, all indicate the occurrence of anisotropy in the shallow (just below the seafloor) subsea sediment in the area. From the radial component azimuthal gather, we analysed the PP‐ and PS‐wave amplitude variation for the first few layers and determined corresponding anisotropy parameter and VP/VS values. Since fracture at this depth is not likely to occur, we attribute the observed azimuthal anisotropy to the presence of microcracks and grain boundary orientation due to stress. The evidence of anisotropy is ubiquitous in this data set and thus it argues strongly in favour of considering anisotropy in depth imaging for obtaining realistic subsurface images, at the least.  相似文献   

16.
We estimate velocity anisotropy factors from seismic traveltime tomographic data and apply a correction for anisotropy in the inversion procedure to test possible improvements on the traveltime fit and the quality of the resulting tomographic images. We applied the anisotropy correction on a traveltime data set obtained from the investigation of the foundation structure of a monumental building: a Byzantine church from the 11th century AD, in Athens, Greece. Vertical transverse isotropy is represented by one axis of symmetry and one anisotropy magnitude for the entire tomographic inversion grid. We choose the vertical direction for the symmetry axis by analysing the available data set and taking into account information on the character of the foundations of the church from the literature and past excavations. The anisotropy magnitude is determined by testing a series of values of anisotropy and examining their effect on the tomographic inversion results. The best traveltime fit and image quality are obtained with an anisotropy value (Vmax/Vmin) of 1.6, restricted to the high velocity structures in the subsurface. We believe that this anisotropy value, which is significantly higher than the usual values reported for near‐surface geological material, is related to the fabric of the church foundations, due to the shape of the individual stone blocks and the layout of the stonework. Inversion results obtained with the correction for anisotropy indicate that both the traveltime fit and the image quality are improved, providing an enhanced reconstruction of the velocity field, especially for the high‐velocity features. Based on this enhanced and more reliable reconstruction of velocity distribution, an improved image of the subsurface material character was made possible. In particular, the pattern and state of the church foundations and possible weak ground material areas were revealed more clearly. This improved subsurface knowledge may assist in a better design of restoration measures for monumental buildings such as Byzantine churches.  相似文献   

17.
Seismic attenuation in Faroe Islands basalts   总被引:2,自引:1,他引:1  
We analysed vertical seismic profiling (VSP) data from two boreholes at Glyvursnes and Vestmanna on the island of Streymoy, Faroe Islands, to determine the magnitude and causes of seismic attenuation in sequences of basalt flows. The work is part of SeiFaBa, a major project integrating data from vertical and offset VSP, surface seismic surveys, core samples and wireline log data from the two boreholes. Values of effective seismic quality factor (Q) obtained at Glyvursnes and Vestmanna are sufficiently low to significantly degrade the quality of a surface reflection seismic image. This observation is consistent with results from other VSP experiments in the North Atlantic region. We demonstrate that the most likely cause of the low values of effective Q at Glyvursnes and Vestmanna is a combination of 1D scattering and intrinsic attenuation due to seismic wave‐induced fluid flow within pores and micro‐cracks. Tests involving 3D elastic wave numerical modelling with a hypothetical basalt model based on field observations, indicate that little scattering attenuation is caused by lateral variations in basalt structure.  相似文献   

18.
Seismic anisotropy provides important constraints on deformation patterns of Earth's material. Rayleigh wave dispersion data with azimuthal anisotropy can be used to invert for depth-dependent shear wavespeed azimuthal anisotropy, therefore reflecting depth-varying deformation patterns in the crust and upper mantle. In this study, we propose a two-step method that uses the Neighborhood Algorithm(NA) for the point-wise inversion of depth-dependent shear wavespeeds and azimuthal anisotropy from Rayleigh wave azimuthally anisotropic dispersion data. The first step employs the NA to estimate depthdependent VSV(or the elastic parameter L) as well as their uncertainties from the isotropic part Rayleigh wave dispersion data. In the second step, we first adopt a difference scheme to compute approximate Rayleigh-wave phase velocity sensitivity kernels to azimuthally anisotropic parameters with respect to the velocity model obtained in the first step. Then we perform the NA to estimate the azimuthally anisotropic parameters Gc/L and Gs/L at depths separately from the corresponding cosine and sine terms of the azimuthally anisotropic dispersion data. Finally, we compute the depth-dependent magnitude and fast polarization azimuth of shear wavespeed azimuthal anisotropy. The use of the global search NA and Bayesian analysis allows for more reliable estimates of depth-dependent shear wavespeeds and azimuthal anisotropy as well as their uncertainties.We illustrate the inversion method using the azimuthally anisotropic dispersion data in SE Tibet, where we find apparent changes of fast axes of shear wavespeed azimuthal anisotropy between the crust and uppermost mantle.  相似文献   

19.
SIM‐France is a large connected atmosphere/land surface/river/groundwater modelling system that simulates the water cycle throughout metropolitan France. The work presented in this study investigates the replacement of the river routing scheme in SIM‐France by a river network model called RAPID to enhance the capacity to relate simulated flows to river gauges and to take advantage of the automated parameter estimation procedure of RAPID. RAPID was run with SIM‐France over a 10‐year period and results compared with those of the previous river routing scheme. We found that while the formulation of RAPID enhanced the functionality of SIM‐France, the flow simulations are comparable in accuracy to those previously obtained by SIM‐France. Sub‐basin optimization of RAPID parameters was found to increase model efficiency. A single criterion for quantifying the quality of river flow simulations using several river gauges globally in a river network is developed that normalizes the square error of modelled flow to allow equal treatment of all gauging stations regardless of the magnitude of flow. The use of this criterion as the cost function for parameter estimation in RAPID allows better results than by increasing the degree of spatial variability in optimization of model parameters. Likewise, increased spatial variability of RAPID parameters through accounting for topography is shown to enhance model performance. Copyright © 2011 John Wiley & Sons, Ltd.  相似文献   

20.
This research incorporates the generalized likelihood uncertainty estimation (GLUE) methodology in a high‐resolution Environmental Protection Agency Storm Water Management Model (SWMM), which we developed for a highly urbanized sewershed in Syracuse, NY, to assess SWMM modelling uncertainties and estimate parameters. We addressed two issues that have long been suggested having a great impact on the GLUE uncertainty estimation: the observations used to construct the likelihood measure and the sampling approach to obtain the posterior samples of the input parameters and prediction bounds of the model output. First, on the basis of the Bayes' theorem, we compared the prediction bounds generated from the same Gaussian distribution likelihood measure conditioned on flow observations of varying magnitude. Second, we employed two sampling techniques, the sampling importance resampling (SIR) and the threshold sampling methods, to generate posterior parameter distributions and prediction bounds, based on which the sampling efficiency was compared. In addition, for a better understanding of the hydrological responses of different pervious land covers in urban areas, we developed new parameter sets in SWMM representing the hydrological properties of trees and lawns, which were estimated through the GLUE procedure. The results showed that SIR was a more effective alternative to the conventional threshold sampling method. The combined total flow and peak flow data were an efficient alternative to the intensive 5‐min flow data for reducing SWMM parameter and output uncertainties. Several runoff control parameters were found to have a great effect on peak flows, including the newly introduced parameters for trees. Copyright © 2013 John Wiley & Sons, Ltd.  相似文献   

设为首页 | 免责声明 | 关于勤云 | 加入收藏

Copyright©北京勤云科技发展有限公司  京ICP备09084417号