亚星游戏官网-www.yaxin868.com

一种水陆检波器数据海底反射系数反演的方法
【专利摘要】一种水陆检波器数据海底反射系数反演的方法，特点是包括以下步骤：用地震人工震源激发和采集水中检波器与陆地检波器地震数据并做预处理；计算水中检波器数据与陆地检波器数据傅立叶变换；计算上行波场以及上行波场一阶和二阶延迟数据；计算上行波场以及上行波场一阶和二阶延迟时间序列数据；计算海底反射系数特征方程系数；构造并求解海底反射系数特征方程；计算上行波数据能量；确定最佳海底反射系数；绘制海底反射系数剖面和存储海底反射系数于地震数据道头中；本发明实现水中检波器数据与陆地检波器数据的合并处理，消除地震数据中海水鸣震多次波干扰影响，有效提高地震数据信噪比和分辨率，具有一定的抗噪能力。
【专利说明】一种水陆检波器数据海底反射系数反演的方法

【技术领域】
[0001]本发明涉及油田的勘探、开发、开采技术，具体是为反映地下地层层位、油藏描述提供高分辨率的地震图形和数据的一种水陆检波器数据海底反射系数反演的方法。

【背景技术】
[0002]随着地震勘探技术的发展，海上油气勘探的难度和深度也越来越大，对地震资料的信噪比和分辨率要求也越来越高。三维地震可以精确地描述储层。OBC(Ocean BottomCable，海底电缆)，是一种联合海上和陆地地震数据采集技术，把检波器固定在海底，可以获得高分辨率三维地震数据。它使用一个固定的排列，把检波器固定在海底接收地震波，而一艘航船仅仅拖拽震源进行地震波激发。在OBC数据采集中，至少三艘船:一艘震源船，仅仅拖拽着气枪震源排列，进行地震波激发；一艘接收船，是固定不动的，连接着海底电缆，接收地震波；一艘船或者几艘船，铺设海底电缆和回收海底电缆。OBC目前仅仅局限于水深不超过150米的海域。但是随着电缆和电缆回收系统的改进，OBC技术正在向更深的海域推进。在一些特殊的地区，为了方便油藏监控试验，设计试验系统，把接收电缆保留在海底多年，以方便多次数据采集。为了得到地下介质高精度的三维图像，必须精确地知道所有炮点和接收点的位置。GPS (Global Posit1ning System,全球定位系统)数据或者基于海岸的无线电定位，提供炮船精确位置。同时也要精确定位电缆船，把每一个检波器放置到预先设计的位置。
[0003]在OBC数据采集时，由于海底和海面都是较强的反射界面。随着一个反射地震子波从地下达到海底，海底电缆中的检波器，感应并记录下这个反射地震子波。这个反射子波继续向上前进达到海面，受到海面的反射，然后改变方向向下传播，达到海底。海底电缆中的检波器，再一次感应并记录下这个地震子波。同时这个地震子波受到海底的反射，然后改变方向向上传播，达到海面，受到海面的反射，然后改变方向向下传播，达到海底。这个循环重复进行。而这些原始反射地震子波不希望的二次和后续达到，就是海水鸣震多次波(交混回响)。海水鸣震多次波是海上地震勘探数据中最大的噪声干扰。消除海水鸣震多次波噪声干扰，是海上地震数据处理中最为重要的步骤。
[0004]在浅水区，使用反褶积方法，可以有效消除海水鸣震多次波干扰，以恢复一次反射地震子波。但是对于海水深度超过10米的地区，原始反射地震子波与后续海水鸣震多次波干扰之间的时差很大，使得反褶积算法不能有效去除海水鸣震多次波干扰产生的同相轴。如果不能有效地从地震数据中去除海水鸣震多次波干扰，那么几个子波代表着一个反射界面，这样就会模糊地质断层界面。OBC采集的数据，提供了同一位置水中检波器数据和陆地检波器数据两种数据。这两种数据分别使用水中检波器和陆地检波器记录。水中检波器是一种压力检波器，记录的是地震波产生的压力变化；陆地检波器是一种质点速度检波器，记录的是质点速度变化。由于这两种检波器的记录机理不同，对于同一位置处海水鸣震多次波干扰，表现出不同特征。与水中检波器记录的海水鸣震多次波干扰相比，陆地检波器记录的海水鸣震多次波干扰表现出极性和振幅特征差异。两种检波器记录的海水鸣震多次波干扰，其极性是相反的，振幅是不同的，且相差一个与海底反射系数成比例的常数，这个常数值就是标定因子。因此利用这种振幅和极性特征差异，可以有效消除海水鸣震多次波干扰。这样消除海水鸣震多次波干扰的步骤包括:(I)在每一个接收点位置处，记录水中检波器数据和陆地检波器数据两种数据；(2)利用两种检波器的传感器灵敏度(传导常数)，调整陆地检波器数据的振幅，以匹配水中检波器数据的振幅；(3)计算确定标定因子、海水深度、海底反射系数参数，(4)利用标定因子，标定调整振幅后的陆地检波器数据；(5)把标定后的陆地检波器数据，与对应的水中检波器数据相加，得到水陆地检波器标定数据，(6)利用海水深度、海底反射系数参数，消除海水鸣震多次波干扰。这样，计算海底反射系数，构成了室内海上地震数据处理消除海水鸣震多次波干扰方法的关键步骤。
[0005]在包含海水鸣震多次波的OBC数据中，由于数据的能量是由地下反射波的能量和海水鸣震多次波的能量构成。消除海水鸣震多次波后，地震数据中，已经消除了多次波能量，剩下地下反射波的能量。因此为了有效估算海底反射系数，通常使用能量最小作为确定海底反射系数的准则。
[0006]目前常规处理方法是，采用扫描方法计算确定海底反射系数。采用预先设定一个海底反射系数范围值和扫描步长，采用扫描方法给出一系列的反射系数值，对一系列海底反射系数，计算一系列Backus海水鸣震多次波逆滤波器算子，对频率域水陆检上行波场数据(即，水陆检波器数据之和)乘以Backus海水鸣震多次波逆滤波器算子，得到滤波处理后的上行波场数据，然后使用傅立叶逆变换，将上行波场数据变换到时间域，在时间域计算一系列波场数据能量，其中最小能量所对应的海底反射系数，就是所求最佳海底反射系数。以上方法采用预先设定一个海底反射系数范围值和扫描步长，采用扫描方法，因此计算费时，且计算的海底反射系数数值必定存在误差。而且计算精度低，计算费时，不能满足实际数据处理需要。
[0007]发明目的
[0008]本发明提供一种目的在于快速、准确地进行水中检波器数据和陆地检波器数据合并处理，消除地震数据中海水鸣震多次波干扰，计算精确、效率高的水陆检波器数据海底反射系数反演的方法。


【发明内容】

[0009]本发明包括以下步骤:
[0010]I)用地震人工震源激发和采集水中检波器与陆地检波器地震数据并做预处理；
[0011]步骤I)所述的预处理是指对地震数据置标签、定义观测系统、水中检波器与陆地检波器地震数据分离、速度分析、叠加预处理工作；
[0012]2)计算水中检波器数据与陆地检波器数据傅立叶变换；
[0013]按照以下公式计算水中检波器数据傅立叶变换Η」(ω):
[0014]",(ω) = ￡>,.["]产(!)

H=I
[0015]式中，hj[n]为水检数据；i为虚数单位，且i2 = _1 ; ω为角频率,单位为弧度/秒(Ω/s) ;n表示共检波点道集数据时间样点顺序号，η = I, 2, 3，…，IN, IN表示共检波点道集数据时间样点总个数表示共检波点道集数据道顺序号= 1，2，3，…，IL，IL表示共检波点道集数据总道数。
[0016]按照以下公式计算陆地检波器数据傅立叶变换(^(ω):
[0017]^/(ω)-Χ-/[/?>/,"Η(2)

?-1
[0018]式中，gj [η]为陆检数据。
[0019]3)计算上行波场以及上行波场一阶和二阶延迟数据；
[0020]按照以下公式计算上行波场数据Cj ( ω ):
[0021]Cj ( ω ) = Hj ( ω ) +Gj ( ω )(3)
[0022]按照以下公式计算上行波场一阶延迟数据Dj(Co):
[0023]Dj ( ω ) = -2RS (Hj ( ω )+Gj ( ω )) θ?ω τ(4)
[0024]式中，Rs为海面反射系数，τ为水层双程旅行时间，单位为秒(S)。
[0025]按照以下公式计算上行波场二阶延迟数据Ε」(ω):
[0026]￡.(ω) = R; (//; (ω)+Gj (ω) )β?2ωτ(5 )
[0027]4)计算上行波场以及上行波场一阶和二阶延迟时间序列数据；
[0028]按照以下公式计算上行波场时间序列数据Cj[η]:
[0029]CM = -V ?~'(ω)β-?απ?ω(6)

2nJ
[0030]按照以下公式计算上行波场一阶延迟时间序列数据Clj[η]:
[0031]d [//] = —i D丨Uii)e "''"dω(7)

-
[0032]按照以下公式计算上行波场二阶延迟时间序列数据θ」[η]:
[0033]=.丄.1"''"cho(8)

2nJ-x
[0034]5)计算海底反射系数特征方程系数；
[0035]按照以下公式计算海底反射系数特征方程系数％、&1、a2、a3和a4:

IL IN
[0036]αο=ΣΣΓ/[/?]

J=I H=I

IL IN
[0037]α{ =[η]

j=\ η=\

IL IN
[0038]<h = ΣΣ("/2[/?]+ 2c:/[叱["])

7=1 h=1

IL IN
[0039]Cl、Clj [/?>,_ [η]

j=l H=I

IL IN
[0040]?4 = Σ Σe) [" ]( 9 )

7=1 ?=1
[0041]6)构造并求解海底反射系数特征方程；
[0042]按照以下公式构造海底反射系数特征方程:
[0043]a1+2a2R+3a3R2+4a4R3 = O(10)
[0044]式中，R为海底反射系数特征方程的变元。求解特征方程(10)，可以得到三个海底反射系数数值RpR2和R3。
[0045]7)计算上行波数据能量；
[0046]按照以下公式计算上行波数据能量Qp Q2和Q3:
[0047]Qy — ci{) + 6/| /?] + u-, /?| u'R' + 6/^/?,
[0048]Qi — + ct'R, + Rn + a i +
[0049]Qit — Cii^ + ?,χ R.' + u.0 R-, + U1R1 +(11)
[0050]8)确定最佳海底反射系数；
[0051]按照以下公式确定最佳海底反射系数Rbsst:
[0052]Rbest = 彳化 P J β , 02, ft }( 12)

Kbest -*M，-^2，九3 J
[0053]9)绘制海底反射系数剖面和存储海底反射系数于地震数据道头中。
[0054]发明效果
[0055]本发明实现水中检波器数据与陆地检波器数据的合并处理，消除地震数据中海水鸣震多次波干扰影响，有效提高地震数据信噪比和分辨率，具有一定的抗噪能力。
[0056]采用直接计算水陆数据上行波场时间序列以及其一阶和二阶延迟时间序列，实现直接计算确定出海底反射系数。计算量小、计算速度快、稳定性好和计算精度高。

【专利附图】

【附图说明】
[0057]图1水陆检合并共检波点道集数据对比；
[0058](a)水检数据
[0059](b)陆检数据
[0060](c)现有水陆检合并数据
[0061](d)本发明水陆检合并数据
[0062]图2水陆检合并共检波点道集数据频谱对比；
[0063](a)水检数据
[0064](b)陆检数据
[0065](C)现有水陆检合并数据
[0066](d)本发明水陆检合并数据
[0067]图3水陆检合并共炮点道集数据对比；
[0068](a)水检数据
[0069](b)陆检数据
[0070](C)现有水陆检合并数据
[0071](d)本发明水陆检合并数据
[0072]图4水陆检合并共炮点道集数据频谱对比；
[0073](a)水检数据
[0074](b)陆检数据
[0075](c)现有水陆检合并数据
[0076](d)本发明水陆检合并数据
[0077]图5水陆检合并叠加数据对比；
[0078](a)水检数据
[0079](b)陆检数据
[0080](c)现有水陆检合并数据
[0081](d)本发明水陆检合并数据
[0082]图6水陆检合并叠加数据频谱对比；
[0083](a)水检数据
[0084](b)陆检数据
[0085](C)现有水陆检合并数据
[0086](d)本发明水陆检合并数据
[0087]图7反演反射系数结果。

【具体实施方式】
[0088]本发明通过计算水中检波器数据和陆地检波器数据的傅立叶变换，进而计算上行波场数据以及其一阶延迟和二阶延迟数据，然后计算上行波场时间序列数据以及其一阶延迟和二阶延迟时间序列数据，再计算海底反射系数特征方程系数，构造并求解海底反射系数特征方程，由特征方程根计算上行波场时间序列能量，再由上行波场时间序列能量确定最佳海底反射系数，最后对水中检波器数据和陆地检波器数据进行合并处理，消除海水鸣震多次波干扰，以有效提高地震数据信噪比和分辨率。
[0089]以下结合附图和实例详细说明本发明。
[0090]本发明采用如下技术方案实现，包括以下步骤:
[0091](I)用地震人工震源激发和采集水中检波器与陆地检波器地震数据，并对地震数据置标签、定义观测系统、水中检波器与陆地检波器地震数据分离、速度分析、叠加等预处理工作；
[0092](2)计算水中检波器数据与陆地检波器数据傅立叶变换；
[0093]Backus海水鸣震多次波逆滤波器算子表达式为
[0094]B(Z) = (1-RsRZ)2 = 1_2RSRZ+RS2R2Z2(I)
[0095]式中，R为海底反射系数，是海底反射系数特征方程的变元，为本发明所求取的参数；RS为海面反射系数，由用户通过试验给定的参数，通常取值为-1 ;Z为延迟算子，其表达式为
[0096]Z = eiaT(2)
[0097]式中，i为虚数单位,且i2 =-1, ω为角频率,单位为弧度/秒(Ω/s) ； τ为水层双程旅行时间，单位为秒(S)，在本发明中作为已知参数，通过其它方法计算确定，与海水深度之间的关系为
2Ζ),、
[0098]T =——(3)
Y
[0099]式中，V为海水速度,单位为米/秒(m/s), D为海水深度,单位为米(m)。
[0100]根据震源鸣震和微屈反射鸣震多次波模型，为了有效消除震源鸣震和微屈反射鸣震多次波，使用算子BackUS海水鸣震多次波逆滤波器算子，乘以上行波场(水陆检波器数据之和)。即
[0101]Sj ( ω ) = (Hj ( ω ) +Gj ( ω )) B ( ω )(4)
[0102]式中，Sj(Co)为消除海水鸣震多次波干扰波场傅立叶变换数据，Η」(ω)为水检波场傅立叶变换数据，Gj(O)为陆检波场傅立叶变换数据，j表示共检波点道集数据道顺序号；j = I, 2，3，…，IL, IL表示共检波点道集数据总道数。且
[0103]A(CO) = tsW


W=I
[0104]/IjUo) ^hjInyim

n=\
[0105]0/(ω)-Χ-/[/7]^,',Η(5)

η=\
[0106]式中，Sj[η]为消除海水鸣震多次波干扰波场数据，hj[η]为水检数据，gj[n]为陆检数据，η表示共检波点道集数据时间样点顺序号；η = I, 2，3，…，IN, IN表示共检波点道集数据时间样点总个数。
[0107]按照公式(5)计算水中检波器数据傅立叶变换&(ω)和陆地检波器数据傅立叶变换 Gj(OJ)；
[0108](3)计算上行波场以及上行波场一阶和二阶延迟数据；
[0109]将方程(I)、(2)代入(4)，有
[0110]Sj ( ω ) = Cj ( ω ) +Dj ( ω ) R+Ej ( ω ) R2(6)
[0111]式中，
[0112]Cj ( ω ) = Hj ( ω ) +Gj ( ω )
[0113]Dj ( ω ) = -2RS (Hj ( ω )+Gj ( ω )) θ?ω τ
[0114]Ej(O)) = R；{/Ij(Oi) + Gj(M) y2"n( 7)
[0115]按照公式(7)计算上行波场数据Cj(Co)、上行波场一阶延迟数据Dj(Co)和上行波场二阶延迟数据& (ω)；
[0116](4)计算上行波场以及上行波场一阶和二阶延迟时间序列数据；
[0117]对方程(6)进行逆傅立叶变换，有
[0118]S」[η] = C」[n] +(Ij [n] R+e」[n] R2⑶
[0119]式中，
[0120]、L?J = —Si(O))^ ivwCho
[0121]c.[/7] = ^l Ci(Oy)C ",/ω
[0122]d: [//J =-1 D.{(?) β d (?.2KJ -CO
[0123]t [//J = -*—? Ε.(ω)β ,uwdO)(9)
J 2 nJ
[0124]按照公式(9)计算上行波场时间序列数据Cj [η]、上行波场一阶延迟时间序列数据Cij[η]和上行波场二阶延迟时间序列数据ej[n]；
[0125](5)计算海底反射系数特征方程系数；
[0126]消除海水鸣震多次波干扰波场数据Sj[n]的能量为

IL IN
[0127]Q ==+ + ChiR +U4Ra(10)

j=l η=\
[0128]式中，

IL IN
[0129]α0=ΣΣΓ/[/?]

J=I W=I

IL IN
[0130]Ci1 =[η]

j=l η=\
[0131]α2 = ΧΧ(?/；[/7] + 26^/7^^/7])

/=1 η=Λ

IL IN
[0132]α3 =2ΣΣdi [/?]e/ ["]./=1 η=1

IL IN
[0133]α4 =ΣΣβ/[/?](11)

7=1 η=1
[0134]按照公式(11)计算海底反射系数特征方程系数％、&1、a2、a3和a4 ；
[0135](6)构造并求解海底反射系数特征方程；
[0136]方程(10)两边对海底反射系数R求导数，有
[0137]-= i2j + R + 3ci^R + Aq^R(12)
SR"
[0138]令能量导数为0，得到海底反射系数R的特征方程
[0139]a1+2a2R+3a3R2+4a4R3 = 0(13)
[0140]按照公式(13)构造海底反射系数特征方程，求解特征方程(13)，可以得到三个海底反射系数数值RpR2和R3;
[0141](7)计算上行波数据能量；
[0142]将三个海底反射系数数值仏、R2和R3分别代入方程(10)，得到三个能量值Q1' Q2和Q3:
[01 43] Qi = “O + “I 及I + ^R' + “3 穴I + “4尺1
[0144]Q-ι — “() + (i\R^ + Rn + ci'R’ + u^Rt
[0145]Qnt — c1 + ο,χ + Ct0 + ci^R^ +( 14)
[0146](8)确定最佳海底反射系数；
[0147]三个能量值Q1A2和Q3中最小值所对应的海底反射系数，就是最佳海底反射系数。即最佳海底反射系数Rbest:
[0148]R = Min \Q^O,,Q,, }⑴)
/^[?,'馬'馬]-
[0149]按照公式(15)确定最佳海底反射系数Rbest ；
[0150](9)绘制海底反射系数剖面和存储海底反射系数于地震数据道头中。
[0151]本发明反演的海底反射系数参数，应用于随后水陆检数据合并处理。
[0152]本发明进行了实际数据处理验证。
[0153]图1是水陆检合并共检波点道集数据对比，横坐标是炮检距(Offset)，单位为米(m)，纵坐标是时间(Time)，单位为毫秒(ms)，(a)是水检数据，(b)是陆检数据，(c)现有水陆检合并数据，(d)是本发明水陆检合并数据；图2是水陆检合并共检波点道集数据频谱对t匕，横坐标是频率(Frequency),单位为赫兹(Hz),纵坐标是振幅(Amplitude), (a)是水检数据，(b)是陆检数据，(C)是现有水陆检合并数据，(d)是本发明水陆检合并数据；图3是水陆检合并共炮点道集数据对比，横坐标是道号(Trace)，纵坐标是时间(Time)，单位为毫秒(ms)，(a)是水检数据，(b)是陆检数据，(C)是现有水陆检合并数据，(d)是本发明水陆检合并数据；图4是水陆检合并共炮点道集数据频谱对比，(a)是水检数据，(b)是陆检数据，(C)是现有水陆检合并数据，(d)是本发明水陆检合并数据；图5是水陆检合并叠加数据对t匕，横坐标是共中心点号(CMP)，纵坐标是时间(Time)，单位为毫秒(ms)，(a)是水检数据，(b)是陆检数据，(c)是现有水陆检合并数据，(d)是本发明水陆检合并数据；图6是水陆检合并叠加数据频谱对比，(a)是水检数据，(b)是陆检数据，(C)是现有水陆检合并数据，
(d)是本发明水陆检合并数据；图7是反演反射系数结果，横坐标是共接收点号(ReceiverNumber),纵坐标是反射系数(Reflectivity)。。由图1、图2、图3、图4、图5和图6可以看出，利用本专利反演的海底反射系数，有效消除了水层产生的多次波干扰，提高了 OBC数据信噪比；由图2、图4和图6可以看出，利用本发明反演的海水深度，有效拓宽了海底电缆数据有效频带，提高了 OBC数据分辨率。因此本发明反演的海底反射系数满足了实际数据处理的需要。
【权利要求】
1.一种水陆检波器数据海底反射系数反演的方法，特点是包括以下步骤:1)用地震人工震源激发和采集水中检波器与陆地检波器地震数据并做预处理；2)计算水中检波器数据与陆地检波器数据傅立叶变换；3)计算上行波场以及上行波场一阶和二阶延迟数据；按照以下公式计算上行波场数据Cj(Co):Cj(GJ) = Η」(ω)+6」(ω)(3)按照以下公式计算上行波场一阶延迟数据Ddco):Dj(GJ) =-2Ιζ(Η」(ω)+6」(ω))θ?ωτ(4)式中，Rs为海面反射系数，τ为水层双程旅行时间，单位为秒(S)。按照以下公式计算上行波场二阶延迟数据1_(ω):Ej(M) = R； (// ,(ω) + G7(W) )ei2wt(5'、4)计算上行波场以及上行波场一阶和二阶延迟时间序列数据；按照以下公式计算上行波场时间序列数据cjn]:Cj[η] = 士J: Cy (O)e-—dω(6)按照以下公式计算上行波场一阶延迟时间序列数据+ [η]:".["] =丄厂 Di,:'"dco<?)按照以下公式计算上行波场二阶延迟时间序列数据ej[n]:e ；["] = j E丨((a) e "''"d ω(8 )5)计算海底反射系数特征方程系数；按照以下公式计算海底反射系数特征方程系数％、B1, a2、a3和a4:?ο-ΣΣ^;[?]
7=1 H=IO1=2丈 Y CvMiZy [η]
j=\ n=\?2 ^ΣΣ(β,/[/7]+2cJ ["知 / [/?])
J=I H=I
IL IN?3 =2ZZi7/
j=\ H=I



(9、
J=I n=\6)构造并求解海底反射系数特征方程；按照以下公式构造海底反射系数特征方程:a1+2a2R+3a3R2+4a4R3 = O(10) 式中，R为海底反射系数特征方程的变元，求解特征方程(10)，可以得到三个海底反射系数数值R!、M口 R3; 7)计算上行波数据能量； 按照以下公式计算上行波数据能量Qp Q2和Q3:
= + ?,χ + ci，R\ + a > /?| + ci^Ri
Q1 — Uii + ci'.R’ + ci0 + u'Rn + Q.、— u,^ + ^/| R' + U0 R.' + u.'R‘' + U^R ^(11) 8)确定最佳海底反射系数； 按照以下公式确定最佳海底反射系数Rbsst: Rbest = Min(12) 9)绘制海底反射系数剖面和存储海底反射系数于地震数据道头中。
2.根据权利要求1的方法，特点是步骤I)所述的预处理是指对地震数据置标签、定义观测系统、水中检波器与陆地检波器地震数据分离、速度分析、叠加预处理。
3.根据权利要求1的方法，特点是步骤2)所述的陆地检波器数据傅立叶变换是: 按照以下公式计算水中检波器数据傅立叶变换&(ω): = ^(I)

η=\ 式中，hj[n]为水检数据；i为虚数单位，且i2 = -1 ； ω为角频率,单位为弧度/秒(Ω/s) ；η表示共检波点道集数据时间样点顺序号，η = I, 2, 3，…，IN, IN表示共检波点道集数据时间样点总个数；j表示共检波点道集数据道顺序号；j = 1，2，3，…，IL, IL表示共检波点道集数据总道数； 按照以下公式计算陆地检波器数据傅立叶变换Gd ω):印⑴)=；^?高’^2)

?=1 式中，g彳η]为陆检数据。
【文档编号】G01V1/36GK104181586SQ201410379202
【公开日】2014年12月3日 申请日期:2014年8月4日 优先权日:2014年8月4日
【发明者】高少武, 赵波, 罗国安, 钱忠平, 颉冬莲, 张文栋 申请人:中国石油集团东方地球物理勘探有限责任公司, 中油油气勘探软件国家工程研究中心有限公司