2015, Vol. 58
2. 中国科学院地质与地球物理研究所油气资源研究重点实验室, 北京 100029
2. Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
Gravimetric density and seismic impedance are combined as a posteriori constraint to formulate the Tikhonov regularization model. Minimization of the Tikhonov regularization model can yield a reasonable inversion result of the density distribution of steam chambers in SAGD. First, a regression model of the seismic impedance and the density through well logging data is built, which can be used to calculate the density constraint data through the seismic P wave data and gives the density constraint. Then, an optimization model based on L2 norm is established. In solving the minimization model, Euler equations along with proper choice of the regularization parameter are simultaneously solved. Utilizing the micro-gravimetric Bouguer gravity anomaly data, the joint inversion scheme based on seismic impedance data and gravimetric data is performed, which yields the density distribution of the reservoir. By analyzing properties of the reservoir density, the density distribution of steam chambers in SAGD can be drawn.
Faced with the problem of development and monitoring during SAGD production, this paper investigates the joint inversion method based on the seismic and gravimetric data. The results obtained are as follows: (1) applying the P wave data as a constraint to invert the gravimetric density is proposed; (2) a posteriori regularization method is utilized to calculate the density, in particular, a posteriori choice of the regularization parameter is performed; (3) a quadratic regression model is built for the seismic impedance and the density, theoretical and field data applications reveal that the joint inversion algorithm is stable and can yield high precision of inversion results; (4) it is shown from the field data applications using our joint inversion method that the inverted low density distribution area coincides well with the known well gas injection position, and the connectivity of the low density area is similar to the forward model, hence proves the feasibility of the joint inversion algorithm in practice. We also remark that the joint inversion method relies on the proper initial model and high SNR data; this may be a requirement for practical usage of the joint inversion method.
For the development and monitoring problem of SAGD production, we develop a joint inversion method based on the seismic data and gravimetric data. Using a posteriori choice of the regularization parameter and a quadratic regression model for the seismic impedance and the density, an inversion procedure is performed by solving an Euler equation. The inversion results can well explain and draw the density distribution of steam chambers in SAGD.
地球物理数据的联合反演方法发源于20世纪70年代中期，Vozoff和Jupp(1977)开创了联合反演的先河，用迭代二阶马奎特阻尼最小二乘法实现了一维直流电测深(DC)和大地电磁测深(MT)资料的联合反演，解决了层状介质中的各向异性问题.Savino等(1980)利用地震P波走时和重力资料联合反演，研究华盛顿东部地区地壳上地幔密度、速度结构.Lines等(1988)使用地面地震数据、声波测井，地面重力及井中重力等资料研究了地震、重力同步反演及顺序反演方法，在反演过程中，充分利用了已有的地面地震、井下声波测井、VSP 数据、地面及井下重力数据等资料，从而大大缩小了模型的选择范围，减小了反问题的多解性，强化了解释过程，并分析了同步反演及顺序反演的各自的应用效果.Dobróka等(1991)对垂直地震剖面(VSP)走时数据、电法数据，采用基于最大频率值(MFV)的加权最小二乘算法进行联合反演，与阻尼最小二乘算法相比，该算法具有估计误差小以及初始模型选择对结果影响较小两个特点.Sun和Schuster(1992)提出了一个在层析成像反演中多个目标函数的极小化过程，该过程在层析成像反演中是十分有用的，特别是同时做几种类型数据模拟，该过程将分级的优化问题转成为等效约束优化的问题，从而使问题简单化.Zeyen和Pous(1993)在具有先验信息的基础上，如密度、磁化率、剩余磁化强度等，对重、磁场的联合反演问题进行了研究.而张贵宾等(1993)以BG理论为基础，在重磁异常线性反演中将该理论与吉洪诺夫正则化方法相结合求解地下密度源(或磁源)分布及质心(或磁质心)位置;在重、磁非线性反演中结合应用正则化方法和马奎特思想给出一种确定地下密度(或磁性)界面的稳定迭代算法——正则马奎特法.杨辉(1998)以地震资料解释的三维构造图作为先验信息，用重力三维正演剥离基底及基底以上界面所产生的重力效应，然后对分离后的基底岩性异常用稳健的SVD 算法来线性反演基底密度差.最后，利用重、磁、电、震、地面地质、钻井等资料综合解释了盆地的基底时代及岩性，取得了令人满意的地质效果.21世纪初，Parsons等(2001)应用重力数据及地震数据开展了顺序反演的工作，发展了重力数据及地震走时层析成像数据的顺序综合反演算法.Tondi等(2000，2003)综合利用折射波、宽角度反射波走时以及重力数据进行顺序综合反演来构造二维速度模型.随后又提出地震及重力数据顺序联合反演固体地球深-浅部分构造重建的方法，取得了很好的应用效果.Onizawa等(2002)给出了地震数据与重力数据同步反演较好的算法实例，Coutant等(2012)在其工作基础上进行改进，并在La Soufriere地区热液系统勘查中进行了应用.Hayashi等(2005)应用面波及微重力观测数据，通过建立横波速度与密度之间的关系进行联合反演，并应用阻尼最小二乘的方法进行求解，取得了很好的效果.Moorkamp等(2011)给出MT、重力及散射波数据联合反演的全局方法.
前人诸多的研究成果表明(Maceira and Ammon，2009;Gallardo and Meju，2004;Tikhotsky and Achauer，2008; Basuyau and Tiberi，2011)，联合反演在算法上可以大致分为两类：(1)顺序反演和(2)同步反演.在数据应用上可以分为两类：(1)基于相同物性地球物理观测数据之间的联合反演.如反射地震的走时和振幅联合反演，纵波和横波资料的联合反演;不同电法或电磁法所取得的观测数据的联合反演.这种联合反演有天然的合理性，因为它们均基于相同的岩石物性差异，物理基础相同，观测场之间必然存在着相关性；(2)基于不同岩石物性的地球物理观测数据之间的联合反演，如地震和重力，地震和MT，重力和MT等.这种联合的基础是不同物性之间存在着相关的内在联系，由相关的物性参数必然会诱发相关的物理异常，这是利用多种地球物理信息进行联合反演的地球物理基础(杨辉等，2002).
本文针对稠油开采中SAGD生产实际问题，分析单井与水平井组合情况下蒸汽腔体分布规律(Bulter，2004).基于密度反演后验约束正则化方法(Wang and Xiao，2001; 王彦飞，2007;王彦飞等，2011)，求解地震数据与重力数据联合反演的问题.通过测井数据，建立纵波阻抗与密度的直接关系，进行联合反演，得到储层含汽砂岩及含稠油砂岩的密度分布规律，最终分析SAGD生产中蒸汽腔体分布情况.
2 重力正演模型我们首先给出密度模型和观测系统：采用地表观测和接收数据，g(xi)为在第i个观测点的重力接收数据，dj为地下介质密度模型网格化后第j个网格的密度值(通常是未知的)，如图 1所示.
 | 图 1 密度模型及观测系统示意图 Fig. 1 Map of density model and acquisition geometry |
根据图 1给出的重力观测模式，可以得到重力计算的物理模型如下：
上式可简写为如下形式：
在实际问题中，由于测量误差的存在，观测到的重力异常数据通常是带噪音的.假设g为重力异常理论值，实际数据记为gδ，并且假定
于是(3)式并不能精确求得，一般需要求解一个最小二乘问题：
然而直接求解由(5)式导出的正规化方程组是严重的病态问题，这是因为F*F的条件数要远远大于F的条件数.因此我们需要研究合适的正则化技巧.
3.1 Tikhonov正则化模型由于观测不足(n≤m)，重力反演问题通常是不适定的，因此直接求解(5)式相关的最小二乘问题通常会带来不稳定的计算效果(王彦飞等，2011).为克服问题的不适定性，本文建立如下基于2-范数的密度反演约束最优化模型：
利用偏差原理，公式(7)中α满足如下非线性方程(Wang and Xiao，2001)：
由以上推导，根据正则化反演原理，本文给出如下可以估计正则参数同时可以反演重力密度的正则化反演算法：
(1)取初始迭代值α0=0.01，δ=0.01，ε=0.0001，kmax=30，其中ε为算法终止迭代误差，kmax为终止迭代最大步数.给定位置距离算子 F，重力异常值gδ，初始先验密度模型d0，并令k∶ =0.
(2)利用高斯消去法解如下所示方程组:
(3)计算φ(αk)及φ′(αk):
(4)求取参数αk+1:
(5)终止准则：确定求解是否小于解的合理误差范围，若
可以采用SAGD稠油热采的简化模型进行试算(杨立强等，2007;孟巍等，2006;刘尚奇等，2007)，如图 2中的模型所示.正、反演模拟的数据参数为:
 | 图 2 SAGD生产二维密度模型:双直井与水平井组合二维模型(水平井垂直于纸面)(密度单位为g·cm-3) Fig. 2 2D forward density model in SAGD production: double straight wells with horizontal well joint 2D model (horizontal well perpendicular to the paper) (density unit: g·cm-3) |
(1)网格分布：1)水平方向70个；2)垂直方向30个；
(2)网格大小为：10 m×10 m；
(3)观测间隔为10 m.
这样的网格分布与实际测量中的观测系统一致. 4.2 理论模型正演计算图 2中的模型为SAGD稠油生产中双直井与水平井组合方式下蒸汽腔体分布模型.模型给出了水平井垂直于纸面的情况.
模型中异常体(图 2蓝色部分)表示由于注入蒸汽、稠油被驱替的情况.产生的蒸汽代替稠油的空间位置，并造成储层区域重力值的变化.
模型的密度分别按照工区内含稠油砾岩岩石密度及含蒸汽砾岩岩石密度的平均值给出：(1)含稠油砾岩岩石密度为2.43 g·cm-3(图 2中褐色区域)；(2)含蒸汽砾岩岩石密度为1.98 g·cm-3(图 2中蓝色区域).
由于实际情况中稠油与蒸汽腔体有过渡带存在，为了使模型更为符合实际情况，文中对模型异常体边界进行线性平滑处理，以模拟油-水-汽的过渡变化关系.
对模型进行重力异常的正演计算，可以得到图 3a的重力异常正演结果.
 | 图 3 模型正演结果: (a) 模型重力正演结果； (b)蒸汽腔体密度分布模型(密度单位为：g·cm-3) Fig. 3 Forward simulation result: (a) gravity forward result; (b) density distribution model of steam chambers (density unit: g·cm-3) |
本文提出应用纵波阻抗数据作为约束进行重力联合反演的方法，模型的纵波阻抗数据如图 4所示.
 | 图 4 相对波阻抗模型(色标中数字无量纲) Fig. 4 Relative P-wave impedance model |
下面根据第2节给出的重力模型并利用我们提出的正则化方法进行反演计算.位置距离算子矩阵如图 5所示.计算表明，该位置距离矩阵是极度坏条件的，条件数达到3.1140×1021，因而必须应用正则化反演技巧.应用图 3a中理论正演重力数据进行反演，可得到图 6所示密度反演结果.为了验证反演算法的精 度，给出反演结果与真实模型的误差分布如图 7所示.文中在正演的观测数据中加入随机噪音，以模拟真实数据情况并说明提出算法的抗噪效果.信噪比计算公 式为r=10 log(S/N)，其中S和N分别是信号和噪声的平均功率，单位为dB.图 8给出信噪比为90时的重力正演数据.把初始密度模型(如图 3b)带入反演模型，应用上述正则化反演算法，经过迭代计算最终可以得到反演的模型密度如图 9所示.为了验证反演算法的精度，给出反演结果与真实模型的误差分布如图 10所示.
 | 图 5 模型的位置距离算子矩阵(纵横坐标分别为距离矩阵的维数大小，色标表示元素的大小) Fig. 5 Matrix of relative position operator |
 | 图 6 通过理论正演重力数据反演密度结果(密度单位为g·cm-3) Fig. 6 Density inversion result map via gravity data without noise (density unit: g·cm-3) |
 | 图 7 反演密度误差结果(密度单位为g·cm-3) Fig. 7 Error map of density inversion (density unit: g·cm-3) |
 | 图 8 信噪比为90的正演重力数据 Fig. 8 Forward gravity data curve with S/N=90 |
 | 图 9 反演密度结果(密度单位为g·cm-3) Fig. 9 Density inversion result (density unit: g·cm-3) |
 | 图 10 反演密度误差结果(密度单位为g·cm-3) Fig. 10 Error propagation of density inversion (density unit: g·cm-3) |
值得指出的是，正则参数α在密度反演中起着至关重要的作用.应用3.2节中的算法迭代得到的正则参数序列如图 11所示.正则参数分布曲线说明，随着迭代的进行，问题的条件数逐渐得到改善，正则解越来越逼近问题的真实解(王彦飞，2007).
 | 图 11 迭代得到的正则参数α Fig. 11 Regularization parameter α |
针对SAGD开发区，观测重力结果如图 12所示.利用阻抗Z和密度ρ对该区内两口观察井的测井曲线进行分析，建立二者之间的二次多项式函数关系表达式如下：
 | 图 12 观测线重力异常分布 Fig. 12 Distribution curve of observed abnormal gravity data |
 | 图 13 测井曲线中密度与波阻抗曲线二次回归关系曲线 Fig. 13 The square regression curve of density and P-wave impedance logs |
利用阻抗Z和密度ρ的关系式，得到纵波阻抗数据如图 14所示.利用纵波阻抗数据作为约束，把观测数据带入反演算法进行实际重力数据的反演计算，密度反演结果如图 15所示.我们从反演结果中可以看到，水平井注、采汽位置及垂直井注汽位置附近有低密度分布区，文中初步判断为SAGD生产中蒸汽替换部分.同时在低密度区周围有相对高密度区分布，文中推测这些相对 高密度区是由含水砂岩及含稠油砂岩混合相态介质组成.
 | 图 14 纵波阻抗剖面(图中虚线矩形为直井注汽位置； 空心圆圈为注、采汽水平井) Fig. 14 P-wave impedance section (dotted rectangle: steam injection range of straight well; hollow circles: produced horizontal well) |
 | 图 15 密度反演结果(图中虚线矩形为直井注汽位 置；空心圆圈为注、采汽水平井)(密度单位为g·cm-3) Fig. 15 Density inversion section (dotted rectangle: steam injection range of straight well; hollow circles: produced horizontal well) (density unit: g·cm-3) |
文中的密度反演结果直接反映了储层密度的变化，该结果可以用来分析含汽砂岩、含稠油砂岩、油汽水混合砂岩的空间分布区域.这说明我们提出的联合反演方法可以对调整SAGD生产中的注采方案起到辅助作用.
5 结论本文针对SAGD生产中开发监测问题，发展了综合应用地震及重力数据反演储层密度的联合反演算法.本文取得了如下研究成果：
(1)提出应用纵波阻抗数据作为约束进行重力联合反演的方法；
(2)提出密度反演后验约束正则化方法，特别是给出了后验选取正则参数的计算方法；
(3)建立了阻抗和密度的二次回归模型，针对理论模型和实际数据，通过波阻抗数据作为约束进行联合反演，在算法上提高了稳定性，同时得到较高的反演精度.
(4)通过联合反演在实际数据中的应用可以看到，反演得到的低密度分布区域与已知井注汽位置吻合很好，且反演得到的低密度区连通性特征与正演模型形态相似，这些特征均验证了联合反演算法在实际生产中的有效性.
同时，文中的联合反演方法对信噪比及初始模型的要求较高，在实际应用中需要较高的数据信噪比及较为精确的初始模型作为约束.
致谢 十分感谢审稿人提出的宝贵意见，使得论文内容更加充实.感谢辽河油田提供相关数据并允许此文的发表，同时感谢对于此工作给予重大帮助的同事们.
| [1] | Basuyau C, Tiberi C. 2011. Imaging lithospheric interfaces and 3D structures using receiver functions, gravity, and tomography in a common inversion scheme. Computers & Geosciences, 37(9): 1381-1390. |
| [2] | Bulter R M. 2004. The behavior of non-condensable gas in SAGD: A rationalization. Journal of Canadian Petroleum Technology, 43(1): 28-34. |
| [3] | Coutant O, Bernard M L, Beauducel F, et al. 2012. Joint inversion of P-wave velocity and density, application to La Soufrière of Guadeloupe hydrothermal system. Geophysical Journal International, 191(2): 723-742. |
| [4] | Dobróka M, Gyulai á, Ormos T, et al. 1991. Joint inversion of seismic and geoelectric data recorded in an underground coal mine. Geophys. Prosp., 39(5): 643-665. |
| [5] | Gallardo L A, Meju M A. 2004. Joint two-dimensional DC resistivity and seismic travel time inversion with cross-gradients constraints. J. Geophys. Res., 109(B3), doi: 10.1029/2003JB002716. |
| [6] | Hayashi K, Matsuoka T, Hatakeyama H. 2005. Joint analysis of a surface-wave method and micro-gravity survey. Journal of Environmental and Engineering Geophysics, 10(2): 175-184. |
| [7] | Lines L R, Schultz A K, Treitel S. 1988. Cooperative inversion of geophysical data. Geophysics, 53(1): 8-20. |
| [8] | Liu S Q, Wang X C, Gao Y R, et al. 2007. SAGD process with the combination of vertical and horizontal wells in super-heavy oil reservoir. Petroleum Exploration and Development (in Chinese), 34(2): 234-238. |
| [9] | Maceira M, Ammon C J. 2009. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure. J. Geophys. Res., 114(B2), doi: 10.1029/2007JB005157. |
| [10] | Meng W, Jia D, Xie J N, et al. 2006. Optimization of geological design for SAGD process in super heavy oil reservoirs using a combination of vertical and horizontal wells. Journal of Daqing Petroleum Institute (in Chinese), 30(2): 44-47. |
| [11] | Moorkamp M, Heincke B, Jegen M, et al. 2011. A framework for 3-D joint inversion of MT, gravity and seismic refraction data. Geophys. J. Int., 184(1): 477-493. |
| [12] | Onizawa S, Mikada H, Watanabe H, et al. 2002. A method for simultaneous velocity and density inversion and its application to exploration of subsurface structure beneath Izu-Oshima volcano, Japan. Earth, Planets and Space, 54(8): 803-817. |
| [13] | Parsons T, Blakely R, Brocher T. 2001. A simple algorithm for sequentially incorporating gravity observations in seismic travel time tomography. Int. Geol. Rev., 43(12): 1073-1086. |
| [14] | Savino J M, Rodi W L, Masso J F. 1980. Simultaneous inversion of multiple geophysical data sets for earth structure. 45th Annual International Meeting, SEG. |
| [15] | Sun Y H, Schuster G T. 1992. Hierarchic optimizations for smoothing and cooperative inversion. Annual Meeting Abstracts, Society of Exploration Geophysicists, 745-748. |
| [16] | Tikhotsky S, Achauer U. 2008. Inversion of controlled-source seismic tomography and gravity data with the self-adaptive wavelet parametrization of velocities and interfaces.Geophys. J. Int., 172(2): 619-630. |
| [17] | Tondi R, De Franco R, Barzaghi R. 2000. Sequential integrated inversion of refraction and wide-angle reflection travel times and gravity data for two-dimensional velocity structures. Geophys. J. Int., 141(3): 679-698. |
| [18] | Tondi R, De Franco R, Biella G. 2003. Deep and shallow solid-earth structures reconstructed with sequential integrated inversion (SII) of seismic and gravity data. Space Science Reviews, 108(1-2): 115-130. |
| [19] | Vozoff K, Jupp D L B. 1977. Effective search for a buried layer: An approach to experimental design in geophysics. Society of Exploration Geophysicists, 8(1): 6-15. |
| [20] | Wang Y F. 2007. Computational Methods for Inverse Problems and Their Applications (in Chinese). Beijing: Higher Education Press. |
| [21] | Wang Y F, Stepanova I E, Strakhov V N, et al. 2011. Inverse Problems in Geophysics and Solution Methods (in Chinese). Beijing: Higher Education Press. |
| [22] | Wang Y F, Xiao T Y. 2001. Fast realization algorithms for determining regularization parameters in linear inverse problems. Inverse Problems, 17(2): 281-291. |
| [23] | Yang H. 1998. Basement density inversion using gravimetric and seismic data and the integrative interpretation. OGP (in Chinese), 33(4): 496-502. |
| [24] | Yang H, Dai S K, Song H B, et al. 2002. Overview of joint inversion of integrated geophysics. Progress in Geophysics (in Chinese), 17(2): 262-271. |
| [25] | Yang L Q, Chen Y M, Wang H Y, et al. 2007. Physical and numerical simulation of steam assisted gravity drainage with vertical and horizontal well combination in extra heavy oil reservoir. Journal of China University of Petroleum (in Chinese), 31(4): 64-69. |
| [26] | Zeyen H, Pous J. 1993. 3-D joint inversion of magnetic and gravimetric data with a priori information. Geophysical Journal International, 112(2): 244-256. |
| [27] | Zhang G B, Shen N H, Wang X C, et al. 1993. The program system on generalized linear cooperative inversion of the potential anomalies. Journal of Changchun University of Earth Sciences (in Chinese), 23(2): 197-204. |
| [28] | 刘尚奇, 王晓春, 高永荣等. 2007. 超稠油油藏直井与水平井组合SAGD技术研究. 石油勘探与开发, 34(2): 234-238. |
| [29] | 孟巍, 贾东, 谢锦男等. 2006. 超稠油油藏中直井与水平井组合SAGD技术优化地质设计. 大庆石油学院学报, 30(2): 44-47. |
| [30] | 王彦飞. 2007. 反演问题的计算方法及其应用. 北京: 高等教育出版社. |
| [31] | 王彦飞, 斯捷潘诺娃I E, 提塔连科V N等. 2011. 地球物理数值反演问题. 北京: 高等教育出版社. |
| [32] | 杨辉. 1998. 重力、地震联合反演基岩密度及综合解释. 石油地球物理勘探, 33(4): 496-502. |
| [33] | 杨辉, 戴世坤, 宋海滨等. 2002. 综合地球物理联合反演综述. 地球物理学进展, 17(2): 262-271. |
| [34] | 杨立强, 陈月明, 王宏远等. 2007. 超稠油直井-水平井组合蒸汽辅助重力泄油物理和数值模拟. 中国石油大学学报(自然科学版), 31(4): 64-69. |
| [35] | 张贵宾, 申宁华, 王喜臣等. 1993. 位场广义线性综合反演系统的建立. 长春地质学院学报, 23(2): 197-204. |