首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 15 毫秒
1.
A seismic variant of the distorted Born iterative inversion method, which is commonly used in electromagnetic and acoustic (medical) imaging, has been recently developed on the basis of the T‐matrix approach of multiple scattering theory. The distorted Born iterative method is consistent with the Gauss–Newton method, but its implementation is different, and there are potentially significant computational advantages of using the T‐matrix approach in this context. It has been shown that the computational cost associated with the updating of the background medium Green functions after each iteration can be reduced via the use of various linearisation or quasi‐linearisation techniques. However, these techniques for reducing the computational cost may not work well in the presence of strong contrasts. To deal with this, we have now developed a domain decomposition method, which allows one to decompose the seismic velocity model into an arbitrary number of heterogeneous domains that can be treated separately and in parallel. The new domain decomposition method is based on the concept of a scattering‐path matrix, which is well known in solid‐state physics. If the seismic model consists of different domains that are well separated (e.g., different reservoirs within a sedimentary basin), then the scattering‐path matrix formulation can be used to derive approximations that are sufficiently accurate but far more speedy and much less memory demanding because they ignore the interaction between different domains. However, we show here that one can also use the scattering‐path matrix formulation to calculate the overall T‐matrix for a large model exactly without any approximations at a computational cost that is significantly smaller than the cost associated with an exact formal matrix inversion solution. This is because we have derived exact analytical results for the special case of two interacting domains and combined them with Strassen's formulas for fast recursive matrix inversion. To illustrate the fact that we have accelerated the T‐matrix approach to full‐waveform inversion by domain decomposition, we perform a series of numerical experiments based on synthetic data associated with a complex salt model and a simpler two‐dimensional model that can be naturally decomposed into separate upper and lower domains. If the domain decomposition method is combined with an additional layer of multi‐scale regularisation (based on spatial smoothing of the sensitivity matrix and the data residual vector along the receiver line) beyond standard sequential frequency inversion, then one apparently can also obtain stable inversion results in the absence of ultra‐low frequencies and reduced computation times.  相似文献   

2.
量子遗传算法作为一种高效的优化算法,仍存在容易陷入局部极值的缺点.为提高算法的高效性,并探讨将算法应用于大地电磁二维反演的可行性和有效性,本文对算法进行了改进,并通过一维两层D型和四层HK型模型数值试验验证了改进的有效性.然后将改进后的算法引入二维大地电磁反演,在引入滑动子空间思想,同时只考虑最简化反演条件的前提下,对一个简单的二维低阻地电模型进行了传统量子遗传算法和改进量子遗传算法的反演,结果说明了将量子遗传算法引入基于子空间的二维大地电磁反演是可行的和有效的,而且改进的算法效果要优于传统算法.最后对实测资料进行了反演,得到了比较好的结果.  相似文献   

3.
Velocity model building and impedance inversion generally suffer from a lack of intermediate wavenumber content in seismic data. Intermediate wavenumbers may be retrieved directly from seismic data sets if enough low frequencies are recorded. Over the past years, improvements in acquisition have allowed us to obtain seismic data with a broader frequency spectrum. To illustrate the benefits of broadband acquisition, notably the recording of low frequencies, we discuss the inversion of land seismic data acquired in Inner Mongolia, China. This data set contains frequencies from 1.5–80 Hz. We show that the velocity estimate based on an acoustic full‐waveform inversion approach is superior to one obtained from reflection traveltime inversion because after full‐waveform inversion the background velocity conforms to geology. We also illustrate the added value of low frequencies in an impedance estimate.  相似文献   

4.
We have previously applied three‐dimensional acoustic, anisotropic, full‐waveform inversion to a shallow‐water, wide‐angle, ocean‐bottom‐cable dataset to obtain a high‐resolution velocity model. This velocity model produced an improved match between synthetic and field data, better flattening of common‐image gathers, a closer fit to well logs, and an improvement in the pre‐stack depth‐migrated image. Nevertheless, close examination reveals that there is a systematic mismatch between the observed and predicted data from this full‐waveform inversion model, with the predicted data being consistently delayed in time. We demonstrate that this mismatch cannot be produced by systematic errors in the starting model, by errors in the assumed source wavelet, by incomplete convergence, or by the use of an insufficiently fine finite‐difference mesh. Throughout these tests, the mismatch is remarkably robust with the significant exception that we do not see an analogous mismatch when inverting synthetic acoustic data. We suspect therefore that the mismatch arises because of inadequacies in the physics that are used during inversion. For ocean‐bottom‐cable data in shallow water at low frequency, apparent observed arrival times, in wide‐angle turning‐ray data, result from the characteristics of the detailed interference pattern between primary refractions, surface ghosts, and a large suite of wide‐angle multiple reflected and/or multiple refracted arrivals. In these circumstances, the dynamics of individual arrivals can strongly influence the apparent arrival times of the resultant compound waveforms. In acoustic full‐waveform inversion, we do not normally know the density of the seabed, and we do not properly account for finite shear velocity, finite attenuation, and fine‐scale anisotropy variation, all of which can influence the relative amplitudes of different interfering arrivals, which in their turn influence the apparent kinematics. Here, we demonstrate that the introduction of a non‐physical offset‐variable water density during acoustic full‐waveform inversion of this ocean‐bottom‐cable field dataset can compensate efficiently and heuristically for these inaccuracies. This approach improves the travel‐time match and consequently increases both the accuracy and resolution of the final velocity model that is obtained using purely acoustic full‐waveform inversion at minimal additional cost.  相似文献   

5.
Mud volcanism is commonly observed in Azerbaijan and the surrounding South Caspian Basin. This natural phenomenon is very similar to magmatic volcanoes but differs in one considerable aspect: Magmatic volcanoes are generally the result of ascending molten rock within the Earth's crust, whereas mud volcanoes are characterised by expelling mixtures of water, mud, and gas. The majority of mud volcanoes have been observed on ocean floors or in deep sedimentary basins, such as those found in Azerbaijan. Furthermore, their occurrences in Azerbaijan are generally closely associated with hydrocarbon reservoirs and are therefore of immense economic and geological interest. The broadside long‐offset transient electromagnetic method and the central‐loop transient electromagnetic method were applied to study the inner structure of such mud volcanoes and to determine the depth of a resistive geological formation that is predicted to contain the majority of the hydrocarbon reservoirs in the survey area. One‐dimensional joint inversion of central‐loop and long‐offset transient electromagnetic data was performed using the inversion schemes of Occam and Marquardt. By using the joint inversion models, a subsurface resistivity structure ranging from the surface to a depth of approximately 7 km was determined. Along a profile running perpendicular to the assumed strike direction, lateral resistivity variations could only be determined in the shallow depth range using the transient electromagnetic data. An attempt to resolve further two‐dimensional/three‐dimensional resistivity structures, representing possible mud migration paths at large depths using the long‐offset transient electromagnetic data, failed. Moreover, the joint inversion models led to ambiguous results regarding the depth and resistivity of the hydrocarbon target formation due to poor resolution at great depths (>5 km). Thus, 1D/2D modelling studies were subsequently performed to investigate the influence of the resistive terminating half‐space on the measured long‐offset transient electromagnetic data. The 1D joint inversion models were utilised as starting models for both the 1D and 2D modelling studies. The results tend to show that a resistive terminating half‐space, implying the presence of the target formation, is the favourable geological setting. Furthermore, the 2D modelling study aimed to fit all measured long‐offset transient electromagnetic Ex transients along the profile simultaneously. Consequently, 3125 2D forward calculations were necessary to determine the best‐fit resistivity model. The results are consistent with the 1D inversion, indicating that the data are best described by a resistive terminating half‐space, although the resistivity and depth cannot be determined clearly.  相似文献   

6.
In this paper, we present the uncertainty analysis of the 2D electrical tomography inverse problem using model reduction and performing the sampling via an explorative member of the Particle Swarm Optimization family, called the Regressive‐Regressive Particle Swarm Optimization. The procedure begins with a local inversion to find a good resistivity model located in the nonlinear equivalence region of the set of plausible solutions. The dimension of this geophysical model is then reduced using spectral decomposition, and the uncertainty space is explored via Particle Swarm Optimization. Using this approach, we show that it is possible to sample the uncertainty space of the electrical tomography inverse problem. We illustrate this methodology with the application to a synthetic and a real dataset coming from a karstic geological set‐up. By computing the uncertainty of the inverse solution, it is possible to perform the segmentation of the resistivity images issued from inversion. This segmentation is based on the set of equivalent models that have been sampled, and makes it possible to answer geophysical questions in a probabilistic way, performing risk analysis.  相似文献   

7.
This paper presents a comparison between subsurface impedance models derived from different deterministic and geostatistical seismic inversion methodologies applied to a challenging synthetic dataset. Geostatistical seismic inversion methodologies nowadays are common place in both industry and academia, contrasting with traditional deterministic seismic inversion methodologies that are becoming less used as part of the geo‐modelling workflow. While the first set of techniques allows the simultaneous inference of the best‐fit inverse model along with the spatial uncertainty of the subsurface elastic property of interest, the second family of inverse methodology has proven results in correctly predicting the subsurface elastic properties of interest with comparatively less computational cost. We present herein the results of a benchmark study performed over a realistic three‐dimensional non‐stationary synthetic dataset in order to assess the performance and convergence of different deterministic and geostatistical seismic inverse methodologies. We also compare and discuss the impact of the inversion parameterisation over the exploration of the model parameter space. The results show that the chosen seismic inversion methodology should always be dependent on the type and quantity of the available data, both seismic and well‐log, and the complexity of the geological environment versus the assumptions behind each inversion technique. The assessment of the model parameter space shows that the initial guess of traditional deterministic seismic inversion methodologies is of high importance since it will determine the location of the best‐fit inverse solution.  相似文献   

8.
Seismic conditioning of static reservoir model properties such as porosity and lithology has traditionally been faced as a solution of an inverse problem. Dynamic reservoir model properties have been constrained by time‐lapse seismic data. Here, we propose a methodology to jointly estimate rock properties (such as porosity) and dynamic property changes (such as pressure and saturation changes) from time‐lapse seismic data. The methodology is based on a full Bayesian approach to seismic inversion and can be divided into two steps. First we estimate the conditional probability of elastic properties and their relative changes; then we estimate the posterior probability of rock properties and dynamic property changes. We apply the proposed methodology to a synthetic reservoir study where we have created a synthetic seismic survey for a real dynamic reservoir model including pre‐production and production scenarios. The final result is a set of point‐wise probability distributions that allow us to predict the most probable reservoir models at each time step and to evaluate the associated uncertainty. Finally we also show an application to real field data from the Norwegian Sea, where we estimate changes in gas saturation and pressure from time‐lapse seismic amplitude differences. The inverted results show the hydrocarbon displacement at the times of two repeated seismic surveys.  相似文献   

9.
We present a novel methodology for 3D gravity/magnetic data inversion. It combines two algorithms for preliminary separation of sources and an original approach to 3D inverse problem solution. The first algorithm is designed to separate sources in depth and to remove the shallow ones. It is based on subsequent upward and downward data continuation. For separation in the lateral sense, we approximate the given observed data by the field of several 3D line segments. For potential field data inversion we apply a new method of local corrections. The method is efficient and does not require trial-and-error forward modeling. It allows retrieving unknown 3D geometry of anomalous objects in terms of restricted bodies of arbitrary shape and contact surfaces. For restricted objects, we apply new integral equations of gravity and magnetic inverse problems. All steps of our methodology are demonstrated on the Kolarovo gravity anomaly in the Danube Basin of Slovakia.  相似文献   

10.
Time‐lapse 3D seismic reflection data, covering the CO2 storage operation at the Snøhvit gas field in the Barents Sea, show clear amplitude and time‐delay differences following injection. The nature and extent of these changes suggest that increased pore fluid pressure contributes to the observed seismic response, in addition to a saturation effect. Spectral decomposition using the smoothed pseudo‐Wigner–Ville distribution has been used to derive discrete‐frequency reflection amplitudes from around the base of the CO2 storage reservoir. These are utilized to determine the lateral variation in peak tuning frequency across the seismic anomaly as this provides a direct proxy for the thickness of the causative feature. Under the assumption that the lateral and vertical extents of the respective saturation and pressure changes following CO2 injection will be significantly different, discrete spectral amplitudes are used to distinguish between the two effects. A clear spatial separation is observed in the distribution of low‐ and high‐frequency tuning. This is used to discriminate between direct fluid substitution of CO2, as a thin layer, and pressure changes that are distributed across a greater thickness of the storage reservoir. The results reveal a striking correlation with findings derived from pressure and saturation discrimination algorithms based on amplitude versus offset analysis.  相似文献   

11.
交错模型指的是由几组走向不同的线性构造(近似二维)在空间(纵向、横向)上组合而成的模型.对于电性结构,交错模型表现为纵向或横向上的构造存在电性主轴方位的变化.从构造维性的角度来看,交错模型是一种由二维模型组合而成的特殊的三维模型,直接对其进行大地电磁二维反演,不易获得可靠的反演结果.本文针对交错模型的特点,提出分频段-分区段反演方案.该方案首先需要借助于阻抗张量成像技术,在频率域确定组成交错模型的各线性构造的电性主轴,然后,针对不同的频段、区段,选择对应电性主轴的数据进行反演,通过初始模型的构建将不同电性主轴方位的反演结果对接起来.本文通过三维理论模型的研究,系统展示了分频段-分区段反演的全过程,归纳得到:在分频段反演时先做低频段反演,在分区段反演时先做二维性更不显著一侧的反演.最后,本文将这一技术用于郯庐断裂带中南段一条实测剖面的反演中,其结果与常规二维反演结果相比较,深部的信息更为丰富,且与其他已有地质、地球物理结果的可对比性更好,表明在构造复杂地区,大地电磁分频段-分区段二维反演具有较高的模型分辨率和可靠性.  相似文献   

12.
起伏地形下CSAMT二维正反演研究与应用   总被引:8,自引:5,他引:8       下载免费PDF全文
雷达 《地球物理学报》2010,53(4):982-993
CSAMT在山区金属矿勘查中,采用各种滤波和相位积分之类的处理方法,校正因地形起伏和局部电性不均匀引起的静态效应,往往难保奏效,开发消除静态效应的新方法是提高CSAMT资料处理与解释水平和方法应用效果的重要研究课题. 本文以如何消除地形影响为重点,对起伏地形下CSAMT二维大地三维源地电模型,采用加权余弦数值积分法,进行波数域电磁场二维有限单元法正演. 为模拟复杂地形地电模型,选取交叉对称网格三角形剖分法,实现了在国内常用赤道电偶极装置的CSAMT二维正演计算;在二维正演的基础上,开发了基于奥克姆反演法的CSAMT二维反演技术,研制出一套起伏地形下CSAMT二维正反演处理与解释方法技术系统. 通过理论模型试算和实测数据处理证实,本系统能有效地削减起伏地形影响. 在找矿应用中,该系统反演的电阻率断面,极大地消除了起伏地形影响和静态效应,突显出清晰的控矿构造和矿体的异常,取得了重要成效.  相似文献   

13.
我国琼东南盆地深水区勘探程度低,工区内只有一口刚打到研究层段黄流组顶部的探井,针对该情况,介绍了一种在无井或少井情况下利用速度谱资料进行地震反演的方法.该方法首先利用离工区较近的井进行有井反演,然后结合有井反演过程的认识,利用速度谱资料提取伪井信息进行伪波阻抗反演,并充分利用已有的地质信息来完善无井反演的结果,最后利用频谱成像技术对反演结果进行验证.结果表明,为了获取准确的伪波阻抗数据,需要准确的速度谱数据,且在构造比较平缓、能量集中、质量较好的道集提取伪井资料,并充分利用已知地震、地质信息和频谱成像分析结果,可以减小无井反演结果的多解性.  相似文献   

14.
Time‐lapse refraction can provide complementary seismic solutions for monitoring subtle subsurface changes that are challenging for conventional P‐wave reflection methods. The utilization of refraction time lapse has lagged behind in the past partly due to the lack of robust techniques that allow extracting easy‐to‐interpret reservoir information. However, with the recent emergence of the full‐waveform inversion technique as a more standard tool, we find it to be a promising platform for incorporating head waves and diving waves into the time‐lapse framework. Here we investigate the sensitivity of 2D acoustic, time‐domain, full‐waveform inversion for monitoring a shallow, weak velocity change (?30 m/s, or ?1.6%). The sensitivity tests are designed to address questions related to the feasibility and accuracy of full‐waveform inversion results for monitoring the field case of an underground gas blowout that occurred in the North Sea. The blowout caused the gas to migrate both vertically and horizontally into several shallow sand layers. Some of the shallow gas anomalies were not clearly detected by conventional 4D reflection methods (i.e., time shifts and amplitude difference) due to low 4D signal‐to‐noise ratio and weak velocity change. On the other hand, full‐waveform inversion sensitivity analysis showed that it is possible to detect the weak velocity change with the non‐optimal seismic input. Detectability was qualitative with variable degrees of accuracy depending on different inversion parameters. We inverted, the real 2D seismic data from the North Sea with a greater emphasis on refracted and diving waves’ energy (i.e., most of the reflected energy was removed for the shallow zone of interest after removing traces with offset less than 300 m). The full‐waveform inversion results provided more superior detectability compared with the conventional 4D stacked reflection difference method for a weak shallow gas anomaly (320 m deep).  相似文献   

15.
The index flood method is widely used in regional flood frequency analysis (RFFA) but explicitly relies on the identification of ‘acceptable homogeneous regions’. This paper presents an alternative RFFA method, which is particularly useful when ‘acceptably homogeneous regions’ cannot be identified. The new RFFA method is based on the region of influence (ROI) approach where a ‘local region’ can be formed to estimate statistics at the site of interest. The new method is applied here to regionalize the parameters of the log‐Pearson 3 (LP3) flood probability model using Bayesian generalized least squares (GLS) regression. The ROI approach is used to reduce model error arising from the heterogeneity unaccounted for by the predictor variables in the traditional fixed‐region GLS analysis. A case study was undertaken for 55 catchments located in eastern New South Wales, Australia. The selection of predictor variables was guided by minimizing model error. Using an approach similar to stepwise regression, the best model for the LP3 mean was found to use catchment area and 50‐year, 12‐h rainfall intensity as explanatory variables, whereas the models for the LP3 standard deviation and skewness only had a constant term for the derived ROIs. Diagnostics based on leave‐one‐out cross validation show that the regression model assumptions were not inconsistent with the data and, importantly, no genuine outlier sites were identified. Significantly, the ROI GLS approach produced more accurate and consistent results than a fixed‐region GLS model, highlighting the superior ability of the ROI approach to deal with heterogeneity. This method is particularly applicable to regions that show a high degree of regional heterogeneity. Copyright © 2014 John Wiley & Sons, Ltd.  相似文献   

16.
Stochastic optimization methods, such as genetic algorithms, search for the global minimum of the misfit function within a given parameter range and do not require any calculation of the gradients of the misfit surfaces. More importantly, these methods collect a series of models and associated likelihoods that can be used to estimate the posterior probability distribution. However, because genetic algorithms are not a Markov chain Monte Carlo method, the direct use of the genetic‐algorithm‐sampled models and their associated likelihoods produce a biased estimation of the posterior probability distribution. In contrast, Markov chain Monte Carlo methods, such as the Metropolis–Hastings and Gibbs sampler, provide accurate posterior probability distributions but at considerable computational cost. In this paper, we use a hybrid method that combines the speed of a genetic algorithm to find an optimal solution and the accuracy of a Gibbs sampler to obtain a reliable estimation of the posterior probability distributions. First, we test this method on an analytical function and show that the genetic algorithm method cannot recover the true probability distributions and that it tends to underestimate the true uncertainties. Conversely, combining the genetic algorithm optimization with a Gibbs sampler step enables us to recover the true posterior probability distributions. Then, we demonstrate the applicability of this hybrid method by performing one‐dimensional elastic full‐waveform inversions on synthetic and field data. We also discuss how an appropriate genetic algorithm implementation is essential to attenuate the “genetic drift” effect and to maximize the exploration of the model space. In fact, a wide and efficient exploration of the model space is important not only to avoid entrapment in local minima during the genetic algorithm optimization but also to ensure a reliable estimation of the posterior probability distributions in the subsequent Gibbs sampler step.  相似文献   

17.
Z‐axis tipper electromagnetic and broadband magnetotelluric data were used to determine three‐dimensional electrical resistivity models of the Morrison porphyry Cu–Au–Mo deposit in British Columbia. Z‐axis tipper electromagnetic data are collected with a helicopter, thus allowing rapid surveys with uniform spatial sampling. Ground‐based magnetotelluric surveys can achieve a greater exploration depth than Z‐axis tipper electromagnetic surveys, but data collection is slower and can be limited by difficult terrain. The airborne Z‐axis tipper electromagnetic tipper data and the ground magnetotelluric tipper data show good agreement at the Morrison deposit despite differences in the data collection method, spatial sampling, and collection date. Resistivity models derived from individual inversions of the Z‐axis tipper electromagnetic tipper data and magnetotelluric impedance data contain some similar features, but the Z‐axis tipper electromagnetic model appears to lack resolution below a depth of 1 km, and the magnetotelluric model suffers from non‐uniform and relatively sparse spatial sampling. The joint Z‐axis tipper electromagnetic inversion solves these issues by combining the dense spatial sampling of the airborne Z‐axis tipper electromagnetic technique and the deeper penetration of the lower frequency magnetotelluric data. The resulting joint resistivity model correlates well with the known geology and distribution of alteration at the Morrison deposit. Higher resistivity is associated with the potassic alteration zone and volcanic country rocks, whereas areas of lower resistivity agree with known faults and sedimentary units. The pyrite halo and ≥0.3% Cu zone have the moderate resistivity that is expected of disseminated sulphides. The joint Z‐axis tipper electromagnetic inversion provides an improved resistivity model by enhancing the lateral and depth resolution of resistivity features compared with the individual Z‐axis tipper electromagnetic and magnetotelluric inversions. This case study shows that a joint Z‐axis tipper electromagnetic–magnetotelluric approach effectively images the interpreted mineralised zone at the Morrison deposit and could be beneficial in exploration for disseminated sulphides at other porphyry deposits.  相似文献   

18.
Remote sensing can provide multi-spatial resolution, multi-temporal resolution multi-spectral band and multi-angular data for the observation of land surface. At present, one of research focuses is how to make the best of these data to retrieve geophysical parameters in conjunction with their a priori knowledge and simul-taneously consider the influence of data uncertainties on inversion results[1-5]. The essence of remote sensing lies in inversion. It is difficult to precisely retrieve parame…  相似文献   

19.
二连盆地由众多小凹陷组成,伊和凹陷是其中之一.重力资料在盆地"探边摸底"、"定凹选带"等方面发挥着重要作用.本文基于布格重力异常资料,利用2.5D建模反演技术对其中横跨主洼槽的两条测线进行可视化建模及正演拟合.其中,初始模型的几何参数取自二维地震解释结果;物性参数中的纵向密度取自测井声波波速与密度的相关公式.正演拟合为建模反演技术的关键环节,密度横向变化,是在拟合过程中通过人机交互试错调整实现.通过正演拟合,揭示了主洼槽沉积层横向和纵向均存在较大的密度不均匀性,主洼槽有利于油气生成.与地震解释结果对比,二者在赛汉组-第四纪、腾格尔组深度相差不大,但凹陷底部拟合结果与地震解释结果有明显差异.通过钻井结果验证,正演拟合各沉积层深度与钻井揭示的各层深度结果一致,二维地震解释结果偏浅,其可能原因是时深转换选取的速度值偏小造成.  相似文献   

20.
Parametric method of flood frequency analysis (FFA) involves fitting of a probability distribution to the observed flood data at the site of interest. When record length at a given site is relatively longer and flood data exhibits skewness, a distribution having more than three parameters is often used in FFA such as log‐Pearson type 3 distribution. This paper examines the suitability of a five‐parameter Wakeby distribution for the annual maximum flood data in eastern Australia. We adopt a Monte Carlo simulation technique to select an appropriate plotting position formula and to derive a probability plot correlation coefficient (PPCC) test statistic for Wakeby distribution. The Weibull plotting position formula has been found to be the most appropriate for the Wakeby distribution. Regression equations for the PPCC tests statistics associated with the Wakeby distribution for different levels of significance have been derived. Furthermore, a power study to estimate the rejection rate associated with the derived PPCC test statistics has been undertaken. Finally, an application using annual maximum flood series data from 91 catchments in eastern Australia has been presented. Results show that the developed regression equations can be used with a high degree of confidence to test whether the Wakeby distribution fits the annual maximum flood series data at a given station. The methodology developed in this paper can be adapted to other probability distributions and to other study areas. Copyright © 2014 John Wiley & Sons, Ltd.  相似文献   

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

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