利用偏移距矢量片地震数据预测储层裂缝的方法

文档序号:10723250阅读:424来源:国知局
利用偏移距矢量片地震数据预测储层裂缝的方法
【专利摘要】本发明提供了一种利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法,利用带有方位角信息偏移距矢量片地震数据直接反演出储层裂缝发育密度和裂缝方位。本发明克服了常规窄方位叠前地震裂缝预测技术存在的反演精度不高、空间分辨率不足等缺陷,能够保留更全面、完整的偏移距和方位角信息,为方位各向异性研究提供了更合适的地震数据,更有利于识别断层、裂缝和储层气水变化特征。
【专利说明】
利用偏移距矢量片地震数据预测储层裂缝的方法
技术领域
[0001] 本发明涉及地球物理勘探技术,具体地,涉及一种利用偏移距矢量片(Offset Vector Tile,简称0VT)宽方位地震叠前道集数据进行储层裂缝预测的方法。
【背景技术】
[0002] 在目前石油天然气勘探中,非常规、隐蔽性油气藏正在成为勘探的重点,裂缝性油 气藏成为当今全球油气增加储量、提高产量的重要领域之一。在储层预测和油气开发方面, 裂缝起到了至关重要的作用,裂缝增加了储集空间,改善了储层的基质渗透率和空隙连通 性,因此储层裂缝预测成为勘探开发中的一项关键技术。
[0003]地震勘探资料中包含了裂缝的信息,地震方法是识别裂缝型储层的重要手段,其 理论基础是各向异性理论。对于当前油气勘探中遇到的储层裂缝参数预测问题,目前的地 球物理预测方法主要分为两大类:叠后地震属性裂缝预测和叠前地震属性裂缝预测。
[0004]叠后地震属性裂缝预测常用的方法包括相干分析法、地震曲率分析法、方差分析 法、边缘检测分析法、分频数据分析法、吸收系数分析法等方法。叠后地震属性裂缝预测方 法未能充分挖掘叠前地震数据所包含的方位各向异性信息,预测精度和预测尺度不及叠前 预测方法。
[0005]叠前地震裂缝预测,目前普遍采用的技术是在常规窄方位叠前地震道集数据的基 础上,经过地震反演处理进行方位各向异性研究。由于常规窄方位地震数据在地震采集过 程中观测方位角信息受限、方位角分布不均匀,基于常规窄方位地震叠前道集数据裂缝预 测方法存在反演精度不高、空间分辨率不足等缺陷。

【发明内容】

[0006] 针对现有技术中存在的不足,本发明的目的之一在于解决上述现有技术中存在的 一个或多个问题。例如,本发明的目的之一在于克服常规窄方位叠前地震裂缝预测技术存 在的反演精度不高、空间分辨率不足等缺陷,提出了一种利用宽方位偏移距矢量片地震数 据进行叠前方位各向异性裂缝参数(例如,裂缝发育密度、裂缝方位)预测方法。
[0007] 为了实现上述目的,本发明的提供了一种利用偏移距矢量片地震叠前道集数据预 测储层裂缝的方法。所述方法包括以下步骤:
[0008] A、对于经过叠前时间偏移处理的偏移距矢量片地震叠前道集数据,根据目的储层 所在时间范围,选取裂缝预测时窗。
[0009] B、对于裂缝预测时窗内每一个成像点位置,计算偏移距矢量片地震叠前道集数据 的偏移距和方位角。
[0010] C、对于裂缝预测时窗内每一个成像点位置的偏移距矢量片地震叠前道集数据,构 建纵波反射系数R(0)。
[0012] 其中,Θ为实际地震炮检观测方位角,R(0)为实际地震炮检观测方位角Θ对应的纵 波反射系数,Ρ为裂缝介质对称轴的方位角表示炮检方向和裂缝走向之间的夹角,Ρ 为各向同性介质下的反射系数,Q为与裂缝发育密度相关的各向异性反射系数。
[0013] D、对于裂缝预测时窗内每一个成像点位置的偏移距矢量片地震叠前道集数据,构 建目标函数:
[0015] 其中,Θ,为第i个偏移距矢量片地震叠前道(简称道)的实际地震炮检观测方位角, R^i(0i)为第i个偏移距矢量片地震叠前道的实际地震炮检观测方位角01对应的纵波反射 系数,N为该成像点位置偏移距矢量片地震叠前道集数据的最大反射道数,Θ Ν为实际最大炮 检观测方位角。
[0016] 在Ε达到最小的情况下,反演求解得到储层的弹性参数值ρ、ο、?)·??(2州及 Q,sin(2<p)。
[0017] E、利用所述步骤D的反演结果,计算得到裂缝密度e和裂缝方位角#。
[0020] 在本发明的一个示例性实施例中,所述方法还可以包括步骤F:对所有成像点偏移 距矢量片地震叠前道集数据,针对目标储层时窗范围,重复以上步骤A至步骤E,以得到整个 三维工区地震数据体目标储层的裂缝密度e和裂缝方位角#反演结果。
[0021] 在本发明的一个示例性实施例中,其中,采用基于最小二乘解的计算方法来使步 骤D中的E最小,求解方程组得到储层的弹性参数值P、Q、0_ms'(2的及心〇),
[0022] 所述方程组为:
[0027]在本发明的一个示例性实施例中,在所述步骤B中,根据第i个偏移距矢量片地震 叠前道的横向偏移距信息〇FFx、纵向偏移距信息0FFy,计算得到所述第i个偏移距矢量片地 震叠前道的偏移距〇FFi和观测方位角0i,
[0029]在本发明的一个示例性实施例中,优选地,所述方法还可以包括在步骤B和C之间 进行地步骤G:对裂缝预测时窗内给定成像点位置的偏移距矢量片道集数据进行部分叠加 处理。所述步骤G可以包括:
[0030] G1、将给定成像点位置的偏移距矢量片地震叠前道集数据的方位角划分为X个扇 区,母个扇区内的道数为Mm,m = 1,2,…,X。
[0031] G2、对任一扇区m内的所有道进行叠加,得到所述任一扇区m内新生成地震道的偏 移距0FFm_new和方位角彳目息9m_new。
[0033]其中,OFFj为第m个扇区内第j道的偏移距,Θ」为第m个扇区内第j道的方位角,j = 1, 2,…,Mm; G3、重复所述步骤G2,直到所有扇区都被处理。
[0034]在本发明的一个示例性实施例中,所述X大于3且小于道集最大覆盖次数。例如,X 取值为18。
[0035]在本发明的一个示例性实施例中,将采集的三维宽方位地震数据在偏移距矢量片 域进行处理,得到偏移距矢量片地震叠前道集数据。
[0036] 在本发明的一个示例性实施例中,所述方法还包括将根据反演的储层裂缝发育密 度和裂缝方位参数进行储层裂缝发育带预测。
[0037] 在本发明的一个示例性实施例中,所述步骤(A)至(D)中,成像点位置可以是预测 时窗内任意一个成像点的位置,也可以是预测时窗内的一段成像点范围,也可以是全部三 维工区数据体的所有成像点范围,即该技术方案适用于三维工区中任意时窗范围、任意平 面范围内偏移距矢量片叠前道集数据的反演。具体反演时窗的选取是根据具体地震工区目 标储层的地质层位决定的。
[0038] 与现有技术相比,本发明的有益技术效果包括:
[0039] (1)本申请采用宽方位偏移距矢量片叠前道集数据进行反演计算,该道集数据保 留了更全面、完整的偏移距和方位角信息,通常其道集覆盖次数比常规窄方位地震道集数 据的覆盖次数更高,克服了常规窄方位叠前地震裂缝预测技术存在的反演精度不高、空间 分辨率不足等缺陷,为方位各向异性研究提供了更合适的地震数据,能够更精确、更高效地 研究振幅、速度等地震响应随方位角变化的特征,更有利于识别断层、裂缝和储层气水变化 特征。
[0040] (2)本申请所构建的目标函数表达式清晰简洁,易于推广应用,方程组求解采用了 基于最小二乘解的计算方法,保证了数值计算的稳定性和收敛性,提高了计算效率和反演 精度。
【具体实施方式】
[0041 ]在下文中,将结合示例性实施例详细地描述根据本发明的利用偏移距矢量片地震 叠前道集数据预测储层裂缝的方法。
[0042]本发明所阐述的偏移距矢量片地震叠前道集数据这一概念,是目前石油、天然气 勘探领域宽方位地震勘探技术中广泛采用的通用概念。所谓偏移距矢量片,是指在一个共 接收线炮集中按炮线距和检波线距等距离划分成许多小矩形,每一个矩形就是一个偏移距 矢量片。偏移距矢量片处理技术的优势在于:偏移距矢量片道集数据保留了更完整的偏移 距和方位角信息,为方位各向异性研究提供更精确的地震数据,有利于研究速度、振幅等地 震响应特征随方位角的变化特征,增强了断层、裂缝和地层岩性变化的识别能力。
[0043] 本发明利用同时带有方位角和偏移距信息的偏移距矢量片地震叠前道集数据直 接反演出储层裂缝发育密度和裂缝方位。在一个示例性实施例中,本发明利用偏移距矢量 片地震叠前道集数据预测储层裂缝的方法包括如下步骤:
[0044] (1)以本领域公知的方法进行三维宽方位地震勘探采集,以得到适合进行宽方位 处理的地震数据,然后,对采集得到的三维宽方位地震数据在偏移距矢量片(0VT)域进行处 理,得到能够进行储层裂缝密度和裂缝方位反演的偏移距矢量片地震叠前道集数据,根据 目的储层所在时间范围,选取裂缝预测时窗。
[0045] (2)导入经过叠前时间偏移处理的偏移距矢量片地震叠前道集数据,对于裂缝预 测时窗内每一个成像点位置,其最大道数为N,根据第i个偏移距矢量片地震叠前道的横向 偏移距信息〇FFx、纵向偏移距信息0FFy,计算该偏移距矢量片道的偏移距OFFjP观测方位角 9i,i = l,2,3,...,N,其计算公式如下:
[0048] ( 3)基于Ruger纵波HTI介质(横向各向同性介质,Hor i zontal Trans verse Isotropic)反射系数公式,对其进行简化,建立如下的计算公式:
[0049] R(0) = P + Q cos:(<7 - φ)
[0050] 对该公式进行三角函数展开,得到如下形式:
[0052] 其中Θ为实际地震炮检观测方位角,R(0)为实际地震炮检观测方位角Θ对应的纵波 反射系数,炉为裂缝介质对称轴的方位角,-识表示炮检方向和裂缝走向之间的夹角,P为 各向同性介质下的反射系数,而Q为与裂缝发育密度相关的各向异性反射系数。在上述公式 中,纵波反射系数随观测方位角的变化表现为椭圆特征。
[0053] (4)对于步骤(2)输入的裂缝预测时窗内每一个成像点位置的偏移距矢量片地震 叠前道集数据,构建目标函数:
[0055] 其中,Θ,为第i个偏移距矢量片地震叠前道的实际地震炮检观测方位角, 为实际地震炮检观测方位角应的纵波反射系数,Ν为该成像点位置偏移距矢量片地震 叠前道集数据的最大反射道数,ΘΝ为实际最大炮检观测方位角。
[0056] 在Ε达到最小的情况下,反演求解得到储层的弹性参数值Ρ、Q、&· c?(2#及 Q · ?η(2φ) 〇
[0057]为使E达到最小值,采用基于最小二乘解的计算方法,求解如下方程组:
[0062 ] 计算上述方程组的解,即可获得P、Q、β · ??(2的及0 · 的值。
[0063] (5)计算裂缝密度及裂缝方位角。在方位各向异性介质中,上述步骤(3)中纵波反 射系数随观测方位角的变化表现为椭圆特征,P+Q为椭圆特征的长轴,P-Q为椭圆特征的短 轴,椭圆长轴方向即为裂缝走向方位角,长轴与短轴之比指示了裂缝发育强度。
[0064]利用步骤(4)的反演结果,分别计算椭圆长轴P+Q,椭圆短轴P-Q,则可计算得到裂 缝密度e和裂缝方位角夢如下:
[0067] (6)对于提供的所有成像点位置的偏移距矢量片地震叠前道集数据,针对目标储 层时窗范围,重复以上步骤(1)至步骤(5),则可得到整个三维工区地震数据体目标储层的 裂缝密度e和裂缝方位角#反演结果。然后,即可根据整个地震工区目标储层的裂缝密度和 裂缝方位角反演结果进行裂缝性储层预测和识别。
[0068] 本发明步骤(3)至步骤(5)中所阐述的纵波反射系数随观测方位角的变化规律及 裂缝密度、裂缝方位计算方法,其理论基础为各向异性理论,即地震波在各向异性介质中平 行或垂直裂缝方向传播时具有不同的传播特征。当射线平面与裂缝发育带方向平行时,纵 波反射振幅最大;当射线平面与裂缝发育带方向垂直时,纵波反射振幅最小;其它方向,纵 波反射振幅介于二者之间,其振幅近似为周期180度的正余弦曲线。
[0069] 在上述步骤(2)至步骤(5)中,成像点位置可以是预测时窗内任意一个成像点的位 置,也可以是预测时窗内的一段成像点范围,也可以是全部三维工区数据体的所有成像点 范围,即该技术方案既适用于三维工区中任意时窗范围、任意平面范围内偏移距矢量片叠 前道集数据的反演。具体反演时窗的选取是根据具体地震工区目标储层的地质层位决定 的。
[0070] 以下将结合具体示例来详细说明本发明的示例性实施例。
[0071] 在本发明的一个具体示例中,利用偏移距矢量片地震叠前道集数据预测储层裂缝 的方法可以由以下步骤实现:
[0072] (1)宽方位纵波地震数据偏移距矢量片(0VT)域叠前时间偏移处理。
[0073]以本领域公知的方法进行三维宽方位地震勘探采集,以得到适合进行宽方位处理 的地震数据,然后,对采集得到的三维宽方位地震数据在偏移距矢量片(0VT)域进行处理, 得到能够进行储层裂缝密度和裂缝方位反演的偏移距矢量片地震叠前道集数据。
[0074] (2)预测时窗选取。
[0075]导入经过叠前时间偏移处理的偏移距矢量片地震叠前道集数据,根据目标储层所 在时间范围,选取裂缝预测时窗[Ta,Tb]。针对不同沉积特征的目标储层,预测时窗[Ta,Tb] 具有不同的取值范围。如针对浅层碎肩岩沉积储层须家河组的储层裂缝预测,[T a,Tb]分别 代表须家河组顶、底反射界面的反射时间;又如针对深层碳酸盐岩沉积储层龙王庙组的储 层裂缝预测,[T a,Tb]分别代表龙王庙组顶、底反射界面的反射时间。
[0076] (3)对于预测时窗内给定成像点位置的偏移距矢量片地震叠前道集数据,计算偏 移距和方位角信息。
[0077]导入经过叠前时间偏移处理的偏移距矢量片地震叠前道集数据,其最大道数为N, 根据第i个偏移距矢量片地震叠前道的横向偏移距信息〇FFx、纵向偏移距信息0FFy,计算该 偏移距矢量片的偏移距OFFi和观测方位角(实际地震炮检观测方位角)0i,i = l,2,3,. . .,N, 其计算公式如下:
[0080] (4)对给定成像点位置的偏移距矢量片地震叠前道集数据进行部分叠加处理。
[0081] 由于原始偏移距矢量片(0VT)叠前道集数据通常信噪比较低,一般需要对道集数 据进行部分叠加处理(所谓部分叠加处理,是指将道集分为若干个扇区,分别对每个扇区内 的地震道进行叠加,即分扇区叠加。相对于传统意义上的道集全部叠加,该处理过程被称为 道集分部分叠加),以提高道集的信噪比,进而提高后续反演算法的稳定性和精确度。具体 实施过程如下:
[0082]假设给定成像点位置的偏移距矢量片(0VT)道集数据方位角范围为0度至180度, 可以将方位角范围按照10度的间隔划分为18等份:0-10度为扇区1,10-20度为扇区2,…, 170-180度为扇区18,则对应每一个扇区内的数据道数为M m,m= 1,2,…,18。
[0083]对于第m个扇区,将该扇区内的地震道进行叠加处理,则该扇区内新生成地震道的 偏移距〇FFm_ne3w和方位角信息0mne3W由如下公式计算:
[0086] m=l ,2,---,18
[0087] 其中,OFFj为第m扇区内第j道的偏移距,Θ」为第m个扇区内第j道的方位角,j = 1, 2,…,Mm〇
[0088] 上述公式计算得到的偏移距0FFm_nejP方位角0mnew,即可作为新生成的部分叠加 偏移距矢量片地震叠前道集数据的更新后实际观测偏移距信息和实际观测方位角信息,并 对外输出,然后将已完成部分叠加的偏移距矢量片地震叠前道集数据作为后续反演计算的 输入数据。
[0089] (5)基于Ruger纵波HTI介质反射系数公式,对其进行简化,建立如下的计算公式:
[0090] R{0) = P + Q- co^-φ-φ)
[0091] 对该公式进行三角函数展开,得到如下形式:
[0093] 其中,Θ为实际地震炮检观测方位角,R(0)为实际地震炮检观测方位角Θ对应的纵 波反射系数,供为裂缝介质对称轴的方位角表示炮检方向和裂缝走向之间的夹角,P 为各向同性介质下的反射系数,而Q为与裂缝发育密度相关的各向异性反射系数。在上述公 式中,纵波反射系数随观测方位角的变化表现为椭圆特征。
[0094] (6)对于步骤(4)输入的裂缝预测时窗内每一个成像点位置的偏移距矢量片地震 叠前道集数据,构建目标函数:
[0096] 其中,0,为实际地震炮检观测方位角,ΙΓ-ΗΘΟ为实际地震炮检观测方位角Θ,对应 的纵波反射系数,i = l,2,3, ...,Ν,Ν为该成像点位置偏移距矢量片地震叠前道集数据的最 大反射道数,ΘΝ为实际最大炮检观测方位角。在Ε达到最小的情况下,反演求解得到储层的 弹性参数值p、Q、Ρ · mv(2p)及β1 · 〇
[0097] 为使E达到最小值,采用基于最小二乘解的计算方法,求解如下方程组:
[0102] 计算上述方程组的解,即可获得P、Q、β · ??'(2妁及C) _.咖叫V的值。
[0103] (7)利用步骤(4)的反演结果,分别计算椭圆长轴P+Q,椭圆短轴P-Q,则可计算得到 裂缝密度e和裂缝方位角梦如下:
[0106] (8)对于提供的所有成像点偏移距矢量片地震叠前道集数据,针对目标储层时窗 范围,重复以上步骤(1)至步骤(7),则可得到整个三维工区地震数据体目标储层的裂缝密 度e和裂缝方位角P反演结果。
[0107] 具体地讲,要达到目标储层段裂缝预测目的,需要针对目标储层段反演时窗[Ta, Tb]内的所有成像点进行反演计算。如需针对浅层碎肩岩沉积储层须家河组进行储层裂缝 预测,则需要针对须家河组顶、底反射界面的反射时间[Ta,Tb]内所有成像点进行反演计算; 又如需针对深层碳酸盐岩沉积储层龙王庙组进行储层裂缝预测,则需要针对龙王庙组顶、 底反射界面的反射时间[T a,Tb]内所有成像点进行反演计算。然后,即可根据整个地震工区 目标储层的裂缝密度和裂缝方位角反演结果进行裂缝性储层预测和识别。
[0108] 综上所述,本发明的优点包括:采用了宽方位偏移距矢量片地震叠前道集数据进 行反演计算,保留了更全面、完整的偏移距和方位角信息,为方位各向异性研究提供了更合 适的地震数据,克服了常规窄方位叠前地震裂缝预测技术存在的反演精度不高、空间分辨 率不足等缺陷,反演计算出的储层裂缝发育密度和裂缝方位信息更加精确、分辨率更高。另 外,构建的目标函数表达式清晰简洁,易于推广应用,基于最小二乘解的反演算法,保证了 数值计算的稳定性和收敛性,提高了计算效率和反演精度,反演结果更有利于断层、裂缝和 储层气水变化特征的有效识别。
[0109] 也就是说,本发明通过利用宽方位偏移距矢量片地震叠前道集数据能够直接反演 出裂缝密度和裂缝方位,进而可实现更加准确的储层裂缝发育带预测。本发明针对复杂岩 性油气藏勘探中裂缝性储层识别问题具有广泛的应用前景。
[0110]尽管上面已经通过结合示例性实施例描述了本发明,但是本领域技术人员应该清 楚,在不脱离权利要求所限定的精神和范围的情况下,可对本发明的示例性实施例进行各 种修改和改变。
【主权项】
1. 一种利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法,其特征在于,所述 方法包括W下步骤: A、 对于经过叠前时间偏移处理的偏移距矢量片地震叠前道集数据,根据目的储层所在 时间范围,选取裂缝预测时窗; B、 对于裂缝预测时窗内每一个成像点位置,计算偏移距矢量片地震叠前道集数据的偏 移距和方位角; C、 对于裂缝预测时窗内每一个成像点位置的偏移距矢量片地震叠前道集数据,构建纵 波反射系数R(9),其中,Θ为实际地震炮检观测方位角,R(0)为实际地震炮检观测方位角Θ对应的纵波反 射系数为裂缝介质对称轴的方位角,心表示炮检方向和裂缝走向之间的夹角,P为各 向同性介质下的反射系数,Q为与裂缝发育密度相关的各向异性反射系数; D、 对于裂缝预测时窗内每一个成像点位置的偏移距矢量片地震叠前道集数据,构建目 标函数:其中,θι为第i个偏移距矢量片地震叠前道的实际地震炮检观测方位角,real(0i)为第i 个偏移距矢量片地震叠前道的实际地震炮检观测方位角θι对应的纵波反射系数,N为该成 像点位置偏移距矢量片地震叠前道集数据的最大反射道数,Θν为实际最大炮检观测方位 角; 在Ε达到最小的情况下,反演求解得到储层的弹性参数值P、Q、^ms'U的及口抑; Ε、利用所述步骤D的反演结果,计算得到裂缝密度e和裂缝方位角Ρ,2. 根据权利要求1所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 其特征在于,所述方法还包括步骤F:针对目标储层时窗范围,对所有成像点的偏移距矢量 片地震叠前道集数据,重复W上步骤A至步骤E,W得到整个Ξ维工区地震数据体目标储层 的裂缝密度6和裂缝方位角|^反演结果。3. 根据权利要求1所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 其特征在于,其中,采用基于最小二乘解的计算方法来使步骤D中的E最小,求解方程组得到 储层的弹性参数值P、Q、資·娜(2約及0·.如口 W,所述方程组为:4. 根据权利要求1所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 其特征在于,在所述步骤B中,根据第i个偏移距矢量片地震叠前道的横向偏移距信息OFFx、 纵向偏移距信息OFFy,计算得到所述第i个偏移距矢量片地震叠前道的偏移距OFFi和实际地 震炮检观测方位角曰1,5. 根据权利要求1所述的偏移距矢量片地震叠前道集数据预测储层裂缝的方法,其特 征在于,所述方法还包括在步骤B和C之间进行地步骤G:对给定成像点位置的偏移距矢量片 道集数据进行部分叠加处理。6. 根据权利要求5所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 其特征在于,所述步骤G包括: G1、将给定成像点位置的偏移距矢量片地震叠前道集数据的方位角划分为X个扇区,每 个扇区内的道数为Mm,m=l,2,…,X; G2、对任一扇区m内的所有道进行叠加,得到所述任一扇区m内新生成地震道的偏移距 0FFm_new和方位角f目息0m_new,其中,OF門为第m个扇区内第j道的偏移距,为第m个扇区内第j道的方位角,j = l, 2 , · ··, Mm ; G3、重复所述步骤G2,直到X个扇区都被处理。7. 根据权利要求1所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 所述X大于3且小于所述偏移距矢量片地震叠前道集的最大覆盖次数。8. 根据权利要求1所述的利用偏移距矢量片地震叠前道集数据预测储层裂缝的方法, 其特征在于,将采集的Ξ维宽方位地震数据在偏移距矢量片域进行处理,得到偏移距矢量 片地震叠前道集数据。
【文档编号】G01V1/30GK106094029SQ201610715839
【公开日】2016年11月9日
【申请日】2016年8月24日
【发明人】谢万学, 巫芙蓉, 李德珍, 何光明, 雷虹, 陈丹, 王珑
【申请人】中国石油集团川庆钻探工程有限公司地球物理勘探公司
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1