首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 78 毫秒
1.
When performing forward modelling and inversion of Magnetic Resonance Sounding (MRS) data, the water-content distribution is typically assumed to be horizontal (1D case). This assumption is fully justified because MRS is often used for characterizing continuous aquifers in a nearly flat environment. However, MRS can also be used in areas with sharp topographical variations. Following a review of the standard MRS equations when using a coincident transmitter/receiver loop, the mathematical terms potentially affected by tilting of the loop are discussed. We present the results of a numerical modelling exercise, studying a case where the surface is not horizontal and the loop cannot be considered to be parallel to the top of the aquifer. This shows that maximum variations in the MRS-signal amplitude are caused mainly by north- or south-dipping slopes. Slope effects depend on the loop size (a larger loop produces a larger error) especially in the presence of shallow water. With a geomagnetic-field inclination of 65° and a slope angle ≤ 10°, the topography causes a maximum variation in amplitude of less than 10%. Near magnetic poles and equator, the slope effect is lower and undetectable in most cases. It was found that within a 10% range of variation in the amplitude, errors introduced into inversions are within the typical uncertainty for MRS inversion and hence no topographic corrections are necessary. Thus, a significant effect from non-horizontal topography might be expected only when data uncertainty is lower than the slope effect (the slope effect is lower than equivalence when data quality is poor). Today, most field data sets are inverted using the modulus of the MRS signal, but some new developments consider the complex signal (both modulus and phase). However, inversion of complex MRS signals, which would provide a higher sensitivity to groundwater distribution, may be affected by slope effect. Thus, the slope orientation and dip angle should be accurately measured in the field when the phase of MRS signals is inverted too.  相似文献   

2.
Geoelectrical profiling with multi-electrode systems has become an important tool for monitoring dike embankments bordering rivers. Profiles running perpendicular to the dike axis are affected by the dike topography, with the amplitude of this effect dependent on the surface geometry and the choice of the electrode configuration. Investigations using seven different electrode configurations have shown that some configurations are less sensitive to the topography than others.The topography correction method (TCM) is an important tool for processing data from measurements at river dikes. This method is generally recommended for flank angles steeper than 10°. The topography effect is calculated by two-dimensional finite element modelling. The resulting synthetic data of a homogeneous dike body are used to apply a topographic correction for each measurement.The topographic effect and correction procedure is demonstrated for synthetic dike data and for a data set from a river dike in Thai Binh province (Vietnam). The topography can be ignored for flank angles less than 25° if an averaged Half-Wenner electrode configuration is used. This configuration has proved to be less affected by undulated topography and the focusing effect of averaging the two data sets provides reliable structural information without the need for time-consuming data inversion.  相似文献   

3.
Time‐domain marine controlled source electromagnetic methods have been used successfully for the detection of resistive targets such as hydrocarbons, gas hydrate, or marine groundwater aquifers. As the application of time‐domain marine controlled source electromagnetic methods increases, surveys in areas with a strong seabed topography are inevitable. In these cases, an important question is whether bathymetry information should be included in the interpretation of the measured electromagnetic field or not. Since multi‐dimensional inversion is still not common in time‐domain marine controlled source electromagnetic methods, bathymetry effects on the 1D inversion of single‐offset and multi‐offset joint inversions of time‐domain controlled source electromagnetic methods data are investigated. We firstly used an adaptive finite element algorithm to calculate the time‐domain controlled source electromagnetic methods responses of 2D resistivity models with seafloor topography. Then, 1D inversions are applied on the synthetic data derived from marine resistivity models, including the topography in order to study the possible topography effects on the 1D interpretation. To evaluate the effects of topography with various steepness, the slope angle of the seabed topography is varied in the synthetic modelling studies for deep water (air interaction is absent or very weak) and shallow water (air interaction is dominant), respectively. Several different patterns of measuring configurations are considered, such as the systems adopting nodal receivers and the bottom‐towed system. According to the modelling results for deep water when air interaction is absent, the 2D topography can distort the measured electric field. The distortion of the data increases gradually with the enlarging of the topography's slope angle. In our test, depending on the configuration, the seabed topography does not affect the 1D interpretation significantly if the slope angle is less or around 10°. However, if the slope angle increases to 30° or more, it is possible that significant artificial layers occur in inversion results and lead to a wrong interpretation. In a shallow water environment with seabed topography, where the air interaction dominates, it is possible to uncover the true subsurface resistivity structure if the water depth for the 1D inversion is properly chosen. In our synthetic modelling, this scheme can always present a satisfactory data fit in the 1D inversion if only one offset is used in the inversion process. However, the determination of the optimal water depth for a multi‐offset joint inversion is challenging due to the various air interaction for different offsets.  相似文献   

4.
We use a total of 839,369 PcP, PKPab, PKPbc, PKPdf, PKKPab, and PKKPbc residual travel times from [Bull. Seism. Soc. Am. 88 (1998) 722] grouped in 29,837 summary rays to constrain lateral variation in the depth to the core-mantle boundary (CMB). We assumed a homogeneous outer core, and the data were corrected for mantle structure and inner-core anisotropy. Inversions of separate data sets yield amplitude variations of up to 5 km for PcP, PKPab, PKPbc, and PKKP and 13 km for PKPdf. This is larger than the CMB undulations inferred in geodetic studies and, moreover, the PcP results are not readily consistent with the inferences from PKP and PKKP. Although the source-receiver ambiguity for the core-refracted phases can explain some of it, this discrepancy suggest that the travel-time residuals cannot be explained by topography alone. The wavespeed perturbations in the tomographic model used for the mantle corrections might be too small to fully account for the trade off between volumetric heterogeneity and CMB topography. In a second experiment we therefore re-applied corrections for mantle structure outside a basal 290 km-thick layer and inverted all data jointly for both CMB topography and volumetric heterogeneity within this layer. The resultant CMB model can explain PcP, PKP, and PKKP residuals and has approximately 0.2 km excess core ellipticity, which is in good agreement with inferences from free core nutation observations. Joint inversion yields a peak-to-peak amplitude of CMB topography of about 3 km, and the inversion yields velocity variations of ±5% in the basal layer. The latter suggests a strong trade-off between topography and volumetric heterogeneity, but uncertainty analyses suggest that the variation in core radius can be resolved. The spherical averages of all inverted topographic models suggest that the data are best fit if the actual CMB radius is 1.5 km less than in the Earth reference model used (i.e. the average outer core radius would be 3478 km).  相似文献   

5.
The magnitudes of the initial amplitude of the magnetic resonance sounding (MRS) signals from an aquifer located in a layered electrically conductive earth, are nonlinear functions of water content distribution. Occam's inversion method is adapted to the nonlinear inversion problem. In the case of an electrically conductive medium, the Jacobian matrix is analytically evaluated at the beginning of the inversion. And the uniqueness of the inversion can be partially solved by imposing the flattest and smoothest model constraints on the optimization problem. Synthetic MRS signals from resistive and conductive earth, as well as field data, have been inverted by Occam's method. The results indicate that with the help of Occam's inversion, a true model can be obtained from an initial model of homogeneous water content. Furthermore, for noise-free MRS signals, both the flattest and smoothest models reveal correct water content distributions. When signals are contaminated by noises, the case is different; and the smoothest model might have a lower water content distributing in a larger range than that of the true model, while which might be obtained by utilizing the flattest model Occam's inversion.  相似文献   

6.
Based on the line integral (LI) and maximum difference reduction (MDR) methods, an automated iterative forward modelling scheme (LI‐MDR algorithm) is developed for the inversion of 2D bedrock topography from a gravity anomaly profile for heterogeneous sedimentary basins. The unknown basin topography can be smooth as for intracratonic basins or discontinuous as for rift and strike‐slip basins. In case studies using synthetic data, the new algorithm can invert the sedimentary basins bedrock depth within a mean accuracy better than 5% when the gravity anomaly data have an accuracy of better than 0.5 mGal. The main characteristics of the inversion algorithm include: (1) the density contrast of sedimentary basins can be constant or vary horizontally and/or vertically in a very broad but a priori known manner; (2) three inputs are required: the measured gravity anomaly, accuracy level and the density contrast function, (3) the simplification that each gravity station has only one bedrock depth leads to an approach to perform rapid inversions using the forward modelling calculated by LI. The inversion process stops when the residual anomalies (the observed minus the calculated) falls within an ‘error envelope’ whose amplitude is the input accuracy level. The inversion algorithm offers in many cases the possibility of performing an agile 2D gravity inversion on basins with heterogeneous sediments. Both smooth and discontinuous bedrock topography with steep spatial gradients can be well recovered. Limitations include: (1) for each station position, there is only one corresponding point vertically down at the basement; and (2) the largest error in inverting bedrock topography occurs at the deepest points.  相似文献   

7.
Introduction Since the middle of the century, gravitational isostasy has been a fundamental hypothesis for inverting the gravity data to find the crust thickness. Geophysicists have done a lot of researches on using gravity data to investigate the depth of Moho discontinuity. Since 1980, the International Lithosphere Program emphasized the importance of investigating the Moho depth variation. Thereafter a lot of results have been published in the world (Braitenberg et al, 2000; Kaban et al,…  相似文献   

8.
Inversion of resistivity in Magnetic Resonance Sounding   总被引:3,自引:0,他引:3  
Magnetic Resonance Sounding (MRS, or Surface Nuclear Magnetic Resonance - SNMR) is used for groundwater exploration and aquifer characterization. Since this is an electromagnetic method, the excitation magnetic field depends on the resistivity of the subsurface. Therefore, the resistivity has to be taken into account in the inversion: either as a priori information or as an inversion parameter during the inversion process, as introduced in the presented paper. Studies with synthetic data show that water content and resistivity can be resolved for a low resistive aquifer even using only the amplitude of the MRS signal. However, the inversion result can be significantly improved using amplitude and phase of the MRS signal. The successful implementation of the inversion for field data shows that the resistivities derived from MRS are comparable to those from conventional geoelectric methods such as DC resistivity and transient electromagnetic. By having information about both the resistivity and the water content, MRS inversions give information about the quality of the water in the aquifer. This is of utmost interest in hydrogeological studies as this specific information cannot be determined solely by geoelectric measurements, due to the nonunique dependence of resistivity on water content and salinity.  相似文献   

9.
Karstic conduits play a crucial role for water supply in many parts of the world. However, the imaging of such targets is generally a difficult task for most geophysical methods. Magnetic Resonance Sounding (MRS) is a geophysical method designed for imaging of water bearing structures. Initially, MRS was developed for characterizing horizontally stratified aquifers. However, when applying a 1D MRS measuring setup to the imaging of 2D–3D targets, the size of which may be much smaller than the loop, the accuracy and the lateral resolution may not be sufficient. We have studied the possibility of simultaneously processing several MRS aligned along a profile to perform a Magnetic Resonance Tomography (MRT). This work emphasizes the gain of resolution for 2D–3D imagery of MRT versus the interpolation of 1D inversion results of MRS along the same profile. Numerical modelling results show that the MRT response is sensitive to the size and location of the 2D target in the subsurface. Sensitivity studies reveal that by using the coincident transmitting/receiving (TX/RX) setup and shifting the loop around the anomaly area, the depth, section and position of a single karstic conduit with a size smaller than the MRS loop size can be resolved. The accuracy of the results depends on the noise level and signal level, the latter parameter being linked to the depth and volume of the karstic conduit and the water content in the limestone matrix. It was shown that when applying MRT to the localization of 2D anomalies such as karstic conduits, the inclination of the geomagnetic field, the orientation of the MRT profile and the angle of crossover of the conduit by the MRT profile must be taken into account. Otherwise additional errors in interpretation should be expected. A 2D inversion scheme was developed and tested. Both numerical and experimental results confirm the efficiency of the developed approach.  相似文献   

10.
In this study, a regional scale gravity data set has been inverted to infer the structure (topography) of the top of the basement underlying sub‐horizontal strata. We apply our method to this real data set for further proof of concept, validation and benchmarking against results from an earlier forward modelling done elsewhere. Our aim is to carry out implicit structural inversion, i.e., to obtain a geologically reasonable model, without specifically solving for structure. The 2.5D volume of interest is parametrized with homogeneous horizontal prisms and a two‐lithology medium is assumed. A possible regional linear trend and a general floating reference are also inverted for. Using a gridded parametrization, linear programming is used to minimize the L1 ‐norm of the data misfit, relative to a floating reference. Given a known density contrast between the lithologies, an inversion using linear programming has the intrinsic advantage that a relatively sharp image of the sub‐surface is retrieved instead of a smooth one. The model recovered is almost bi‐modal and its general features seem to be robust with respect to several parametrization scenarios investigated. The floating reference and a linear trend in the data were also retrieved simultaneously. The inversion results, indicating two depressions in the basement, are robust and agree with those obtained earlier based upon detailed 2D forward modelling using many narrow, near‐vertical prisms.  相似文献   

11.
2D magnetic resonance tomography applied to karstic conduit imaging   总被引:1,自引:1,他引:1  
Karstic conduits play a crucial role for water supply in many parts of the world. However, the imaging of such targets is generally a difficult task for most geophysical methods. Magnetic Resonance Sounding (MRS) is a geophysical method designed for imaging of water bearing structures. Initially, MRS was developed for characterizing horizontally stratified aquifers. However, when applying a 1D MRS measuring setup to the imaging of 2D–3D targets, the size of which may be much smaller than the loop, the accuracy and the lateral resolution may not be sufficient. We have studied the possibility of simultaneously processing several MRS aligned along a profile to perform a Magnetic Resonance Tomography (MRT). This work emphasizes the gain of resolution for 2D–3D imagery of MRT versus the interpolation of 1D inversion results of MRS along the same profile. Numerical modelling results show that the MRT response is sensitive to the size and location of the 2D target in the subsurface. Sensitivity studies reveal that by using the coincident transmitting/receiving (TX/RX) setup and shifting the loop around the anomaly area, the depth, section and position of a single karstic conduit with a size smaller than the MRS loop size can be resolved. The accuracy of the results depends on the noise level and signal level, the latter parameter being linked to the depth and volume of the karstic conduit and the water content in the limestone matrix. It was shown that when applying MRT to the localization of 2D anomalies such as karstic conduits, the inclination of the geomagnetic field, the orientation of the MRT profile and the angle of crossover of the conduit by the MRT profile must be taken into account. Otherwise additional errors in interpretation should be expected. A 2D inversion scheme was developed and tested. Both numerical and experimental results confirm the efficiency of the developed approach.  相似文献   

12.
自适应非结构有限元MT二维起伏地形正反演研究   总被引:5,自引:1,他引:4       下载免费PDF全文
在山区进行MT勘探时,用规则网格有限元方法模拟起伏地形会受到限制.本文采用非结构三角网格可以有效地模拟任意二维地质结构,如起伏地形、倾斜岩层和多尺度构造等.正演引入自适应有限元方法,其在网格剖分过程中能根据单元误差自动细化网格,保证了正演结果的精度.将自适应有限元与Occam算法结合,且引用并行处理技术提高正反演计算速度.通过对比两个理论模型,讨论了地形对MT正演响应的影响;其次进行了不同地电模型带地形反演展示了本文算法的正确性和适用性;最后将该方法应用于实测MT数据处理,证明了自适应非结构有限元方法是复杂地形下处理MT数据的有力工具.  相似文献   

13.
基于有限差分正演的带地形三维大地电磁反演方法   总被引:4,自引:4,他引:0       下载免费PDF全文
本研究实现了一套基于有限差分(FD)方法的大地电磁测深数据带地形三维反演算法及代码.其中,在大地电磁场正演数值模拟方面,开发了起伏地形条件下基于交错网格剖分、有限差分方法的大地电磁测深三维正演代码;在满足平面波场假设的前提下,使用长方体网格剖分模拟三维起伏地形,实现了带地形三维正演计算;并设计理论模型进行试算,经试算结果与前人的有限元法计算结果对比,验证了所研发的带地形三维正演计算的正确性与可靠性.在反演方面,本研究基于非线性共轭梯度方法编写了大地电磁测深带地形三维反演代码,试验了不同的共轭梯度搜索因子β,避免了目标函数对海森矩阵(参数二次导数矩阵)的显式计算和存储,初步实现了大地电磁资料的带地形三维反演.最后,对一系列理论模型进行正演计算,利用其生成的合成数据模拟实测数据进行反演,并与现有的不带地形大地电磁测深三维反演结果比较,检验了所研发的带地形三维反演计算的可靠性与稳定性.  相似文献   

14.
To improve the knowledge of the regionally important Continental Terminal 3 (CT3) aquifer in south-western Niger, fifteen magnetic resonance soundings (MRS) were carried out in December 2005 in the vicinity of wells and boreholes. The output MRS geophysical parameters, i.e. water content and decay constants versus depth, were compared to hydrogeological characteristics, i.e. water table depth, total porosity, specific yield and transmissivity estimated from direct measurements, pumping tests and transient groundwater modelling. The MRS-determined parameters were then used to estimate the rates of groundwater recharge.Contained in poorly consolidated Tertiary sandstones, the CT3 aquifer's water table has continuously risen by 4 m in total over the past four decades. Additionally, a significant portion of this increase has occurred in the past decade alone, with an annual rise now ranging between 0.1 and 0.3 m depending on the monitored well. Increase in groundwater recharge due to land clearance and deforestation explains this situation. According to previous estimations, the pre-clearing recharge ranged from 1 to 5 mm per year in 1950–60 s, while more recent recharge rates (1990s–2000s) range from 20 to 50 mm per year. These recharge values are directly affected by estimated aquifer specific yield value, while the spatial variation of rates of water table rise can be attributed to large scale hydrodynamic heterogeneities in the aquifer. However, few field measurements were available to confirm these assumptions.The main results of this study are: (1) The water table depth and aquifer transmissivity are estimated from MRS output parameters with an average accuracy of ± 10% and ± 9% respectively. (2) The MRS-determined water content is linked to both the total porosity and the specific yield of the aquifer, but no quantitative formulation can be proposed as yet. (3) Using the average MRS-determined water content over the investigated area, i.e. 13%, the groundwater recharge rates can be estimated to be ~ 2 mm per year in the 1950–1960s (pre-clearing period), and ~ 23 mm per year for the last decade. (4) The variations in specific yield and transmissivity cannot explain by themselves the spatial variability of the rise of the water table. (5) The ranges in transmissivity and water content obtained from MRS are more realistic than the groundwater modelling outputs. Therefore, MRS could be used to better constrain the aquifer parameters in groundwater modelling with a dense site network.Finally, this work illustrates how MRS can successfully improve characterisation and transient multi-year groundwater balance of commonly found sedimentary aquifers, particularly when integrated with well observations and pumping tests.  相似文献   

15.
Over the last years, full-waveform inversion has become an important tool in the list of processing and imaging technologies available to the industry. For marine towed-streamer data, full-waveform inversion is typically applied using an acoustic approximation because S-waves do not propagate in water and elastic effects in recorded data are generally assumed to be small. We compare acoustic and elastic modelling and full-waveform inversion for a field data set acquired offshore Angola over sediments containing a salt body with significant topology. Forward modelling tests reveal that such geological structures lead to significant mode conversions at interfaces and, consequently, to significant relative amplitude differences when elastically and acoustically modelled traces are compared. Using an acoustic approach for modelling in full-waveform inversion therefore leads to problems matching the synthetic data with the field data, even for recorded pressure data and with trace normalization applied. Full-waveform inversion is unable to find consistent model updates. Applying elastic full-waveform inversion leads to more consistent and reliable model updates with less artefacts, at the expense of additional computation cost. Although two-dimensional marine towed-streamer data are least favourable for the application of full-waveform inversion compared to three-dimensional data or ocean-bottom data, it is recommended to check on the existence of elastic effects before deciding on the final processing and imaging approach.  相似文献   

16.
三维反演是磁测数据定量解释的重要方法,在金属矿勘探中扮演着重要的角色.但是在实际矿区的应用中,传统的磁总场异常反演方法依然存在两个问题:一是地面磁异常反演的深度分辨率较低,深部场源体的成像效果差;二是金属矿中可能包含强剩磁,反演结果可能是完全错误的.尽管前人对上述两个问题分别进行了广泛的研究,但尚未尝试同时解决这两个问题.本文在前人研究的基础上,提出了一种井地磁异常模量联合反演方法,该方法需要的控制参数少,无需加入额外的地质信息,且可用于多场源复杂磁异常的反演,具有较强的适用性.本文方法首先将地面和井中磁异常转化为模量数据,然后利用基于核函数或距离的加权函数将井地模量数据结合起来,使得该方法适用于联合反演.我们利用井地多种异常参量进行反演的模型试验表明,在强剩磁存在时,本文方法的效果优于其他方法,在减少剩磁影响的同时,也改善了深部成像效果,具有良好的应用前景.  相似文献   

17.
Transverse isotropy with a vertical axis of symmetry is a common form of anisotropy in sedimentary basins, and it has a significant influence on the seismic amplitude variation with offset. Although exact solutions and approximations of the PP-wave reflection coefficient for the transversely isotropic media with vertical axis of symmetry have been explicitly studied, it is difficult to apply these equations to amplitude inversion, because more than three parameters need to be estimated, and such an inverse problem is highly ill-posed. In this paper, we propose a seismic amplitude inversion method for the transversely isotropic media with a vertical axis of symmetry based on a modified approximation of the reflection coefficient. This new approximation consists of only three model parameters: attribute A, the impedance (vertical phase velocity multiplied by bulk density); attribute B, shear modulus proportional to an anellipticity parameter (Thomsen's parameter ε−δ); and attribute C, the approximate horizontal P-wave phase velocity, which can be well estimated by using a Bayesian-framework-based inversion method. Using numerical tests we show that the derived approximation has similar accuracy to the existing linear approximation and much higher accuracy than isotropic approximations, especially at large angles of incidence and for strong anisotropy. The new inversion method is validated by using both synthetic data and field seismic data. We show that the inverted attributes are robust for shale-gas reservoir characterization: the shale formation can be discriminated from surrounding formations by using the crossplot of the attributes A and C, and then the gas-bearing shale can be identified through the combination of the attributes A and B. We then propose a rock-physics-based method and a stepwise-inversion-based method to estimate the P-wave anisotropy parameter (Thomsen's parameter ε). The latter is more suitable when subsurface media are strongly heterogeneous. The stepwise inversion produces a stable and accurate Thomsen's parameter ε, which is proved by using both synthetic and field data.  相似文献   

18.
Magnetic resonance sounding (MRS) is an electromagnetic method designed for groundwater investigations. MRS can be applied not only for studying fresh-water aquifers, but also in areas where intrusion of saline water is rendering the subsurface electrically conductive. In the presence of rocks with a high electrical-conductivity attenuation and a phase shift of the MRS signal may influence the efficiency of the MRS method. We investigated the performance of MRS for allowing us to propose a procedure for interpreting MRS data under these conditions. For numerical modeling, we considered a subsurface with a resistivity between 0.5 and 10 Ω m. The results show that the depth of investigation with MRS depends upon the electrical conductivity of groundwater and surrounding rocks, on the depth of the saline water layer, and on the amount of fresh water above the saline water. For interpreting MRS measurements, the electrical conductivity of the subsurface is routinely measured with an electrical or electromagnetic method. However, due to the equivalence problem, the result obtained with these methods may be not unique. Hence, we investigated the influence of the uncertainty in conductivity distribution provided by transient electromagnetic measurements (TEM) on MRS results. It was found that the uncertainty in TEM results has an insignificant effect on MRS.  相似文献   

19.
We present a new approach to analyse the subsurface water content distribution obtained by inversion of MRS data in terms of resolution and penetration depth. It is based on a singular value decomposition (SVD) of the MRS forward operator to derive the model resolution matrix including regularisation parameters, i.e. including noise conditions. The approach takes loop size, subsurface resistivity distribution and noise conditions as input parameters affecting MRS into account and allows an assessment on resolution and penetration.The application of the new approach shows that the loop diameter must be carefully chosen depending on the investigation site to obtain optimal resolution and penetration. Using the introduced resolution measures the quality and reliability of the estimated model can be assessed.  相似文献   

20.
The T/P altimeter data 1993 – 1997 (cycles 11 – 194) has been analyzed with emphases on seasonal variations in sea surface topography (SST). The amplitude of the annual variations amounted to (5.9±0.3) mm when inverted barometer (IB) corrections were applied and (2.0±0.4) mm without any IB corrections. The amplitude of the semi-annual variations in SST was small with IB corrections applied: (0.6±0.3) mm. However, when no IB corrections were applied, it was (1.8±0.4) mm, i.e. the semiannual variations are at the same level as the annual variations with no IB corrections. The phase angle offset of the annual term has shifted by about 180° when IB correction was applied. The dynamics of the ocean-atmosphere system is discussed and it is concluded that it could, at least partly, be responsible for the above observed effects.  相似文献   

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

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