1 研究区地质概况
齐家—古龙凹陷位于松辽盆地中央坳陷西部，属于盆地的次级构造单元，面积约5 500 km2，主要地层包括下白垩统、上白垩统、第三系及第四系（见图1 a）。2021年大庆油田在古龙地区白垩系青山口组一段（简称青一段）发现页岩油储集层，落实含油面积1 413 km2，预测储量为12.68×108 t，标志着中国页岩油勘探再次取得重大突破。
前人研究表明，松辽盆地齐家—古龙凹陷下白垩统青山口组为最大湖泛期沉积[1⇓-3]，发育广泛分布的富有机质暗色泥页岩。图1 b为研究区青一段厚度等值线图，显示研究区青一段沉积时期位于深湖相沉积中心，整体呈脚掌形分布，沉积厚度大、TOC值高，为本区主力烃源岩发育层。青山口组发育互层、夹层和富黏土页岩3种类型页岩油储集层，其中青一段发育富黏土页岩油储集层，其主要岩石矿物组成包括长英质、黏土质和碳酸盐3大类，黏土矿物含量一般为30%~40%，纵向上又进一步细分为Q1—Q9共9个层组；页岩储集层的主要孔隙空间包括有机质孔缝、溶孔、黏土矿物晶间孔[4-5]，占总面孔率的80%以上。重点取心井岩心资料显示中等有机质含量的纹层状泥岩相具有较好生烃潜力，顺层发育的水平页理缝长从几微米到几十微米，宽几十纳米到几百纳米，总面孔率为22%~79%，既是重要的油气储集空间[3]，同时也为油气渗流提供了良好通道，使得这种富含黏土的页岩储集层覆压条件下水平渗透率平均值可达0.58×10-3 μm2。
针对青一段页岩储集层的上述特征，测井岩石物理研究和解释评价的主要难题包括各类孔隙空间的准确识别、以含油富集程度为核心的储集层参数定量计算等，如何确定适合此类储集层特征的饱和度计算模型无疑是其中的核心问题。富含黏土及发育页理缝使得常规的柱塞岩样制备、驱替岩电实验无法开展。围绕这一难题，前人开展了大量地质和地球物理评价方法研究。邹才能等[5-6]系统比较了致密油和页岩油的内涵与特征，总结了国内主要盆地页岩油的地质参数，但其仅从TOC、Ro和S1、氯仿沥青“A”等地球化学指标分析了含油性；还有学者[7⇓⇓⇓⇓⇓⇓⇓⇓-16]针对研究区青山口组泥页岩储集层开发了测井评价的新方法，但仅讨论了TOC、矿物组分、孔隙压力和脆性等参数计算模型；郑建东等[12]采用一维和二维核磁共振测井计算页岩含油饱和度，并简述了基于数字岩心技术确定含油饱和度模型的参数，但没有给出详细的方法原理。刘国强[13]认为现阶段针对非常规储集层的含油饱和度计算精度还很难满足国家储量规范要求，亟需发展新的岩石物理研究和计算方法。综上，笔者在前人研究基础上，以古龙凹陷青一段富含黏土和页理缝页岩油储集层为例，通过构建四维数字岩心和数值模拟，提出了1种新的页岩油含油饱和度计算新方法。
2 储集层特征
2.1 岩性特征
进一步采用QemScan（Quantitative Evaluation of Minerals by Scanning Electron Microscopy，简称QemScan）岩心分析系统对青一段储集层的矿物颗粒分布特征进行成像分析（见图2 b，仅列出含量大于1%的主要矿物类型），结果表明样品中伊利石含量最高（41.9%），含少量绿泥石，黏土矿物颗粒呈分散状分布，构成分散泥质分布形式；长英质矿物主要由石英和钠长石组成，与黏土矿物颗粒均匀混杂堆积。由此可见，富含以伊利石为主的黏土矿物且以分散泥质的形式发育是本区青一段页岩的主要特征。
2.2 物性与孔隙结构特征
根据研究区158块柱塞样孔隙度测试结果统计，青一段储集层孔隙较为发育，孔隙度主要为4%~12%，平均值为7.9%，有效孔隙度（不包括黏土矿物晶间孔和有机孔部分）主体为2%~4%，平均值为3.7%。
根据场发射扫描电镜实验结果，青一段页岩油储集层的孔隙类型可进一步分为残余粒间孔、颗粒溶孔、有机质孔、页理缝和黏土矿物晶间孔5种类型，其孔径依次减小[3]。总体上，青一段页岩孔隙具有3级结构特征，其中最大的孔隙直径为500~600 nm，中间一级的孔隙直径为200~250 nm，第3级孔隙直径为10~20 nm。进一步分析不同矿物纹层的孔隙占比和孔径大小（见图3 ）表明，长英质纹层主要以残余粒间孔、溶蚀微孔为主，等效孔径主体为100~600 nm（见图3 a），而有机质孔占比一般小于10%且以小于100 nm孔为主（见图3 b），不是主要孔隙类型。
电成像测井具有高达5 cm的纵向分辨率，可以揭示纹层特征，核磁共振测井则可以反映储集层的孔隙结构。图4 是GY2HC井青一段测井曲线综合图，第1道井径曲线显示本井段井壁规则，贴井壁偏心式核磁共振测井资料质量可靠，第6道核磁共振T2谱显示以小于10 ms的孔隙组分为主、峰值一般小于3 ms（取横向表面弛豫速率约为10 μm/s，换算为孔径约为50 nm[18]，与图3 b结果一致），也反映出以黏土晶间孔为主的纳米级微孔是本区青一段页岩储集层的主要孔隙类型。根据何文渊等对本区重点取心井的精细描述[19]，页理缝发育段在井径规则时核磁共振测井的T2组分也相对发育，高分辨率电成像测井也揭示纹层状结构发育。图4 最右一列是基于元素扫描测井计算的全岩矿物剖面，其矿物组合特征与岩心实验结果一致。
2.3 含油性特征
微纳米孔隙和页理缝是古龙页岩油的基本储集空间，纳米级孔隙是基础，微米级孔隙是“桥梁”，页理缝则成为页岩油气聚集的末端，三者构成了古龙页岩油储集层的运聚系统[3]。但是对比分析不同类型孔隙空间中原油富集程度来看，差异明显。曾维主等[16]通过对青山口组页岩岩心样品开展有机地球化学、低温氮气吸附和高压压汞实验分析发现，古龙页岩含油饱和度与汞孔隙度之间具有正相关关系，表明大孔径的页理缝、粒间孔隙中页岩油相对更加富集。图5 为GY2HC井青一段两块样品扫描电镜揭示的页理缝及粒间孔的含油特征，可以看出，在地表条件下由于应力释放，页理缝中的原油外溢（见图5 a），是本区主要的含油储集空间，图5 b为纹层状页岩扫描电镜照片，显示黏土晶间孔中含油特征明显。
3 岩电参数和饱和度计算方法
3.1 模型优选与岩电参数确定
从岩石物理特征看，页岩储集层与富含泥质的粉砂岩储集层类似，其电阻率高低不仅与地层水电阻率、含水饱和度、有效孔隙度有关，还与黏土的电阻率、相对含量及矿物颗粒的分布形式有关。根据前文对青一段页岩储集层的岩性特征分析，由于本区有机质含量一般不超过2%，且有机孔占比低，含油饱和度计算可以选择含分散状黏土的泥质砂岩模型。
$\frac{1}{{{R}_{\text{t}}}}=\frac{{{V}_{\text{cl}}}}{{{R}_{\text{cl}}}}{{S}_{\text{w}}}^{0.5n}+\frac{{{\phi }^{m}}}{a{{R}_{\text{w}}}\left( 1-{{V}_{\text{cl}}} \right)}{{S}_{\text{w}}}^{n}$
公式（1）中，Rt、Vcl、ϕ值可分别由电阻率测井、自然伽马测井、中子孔隙度测井等获取。由于研究区青一段沉积总体稳定，纯泥岩段Rcl可取常数2.0 Ω·m。根据地层水分析资料换算到地下，Rw取值为0.2 Ω·m。另外a值一般取1。因此，应用公式（1）计算饱和度的关键是确定岩电参数m、n值，这两个岩电参数的确定是实验研究的重点。其中m值的确定较为简单，n值确定的基本原理是：在真实岩心内部通过人工方式建立不同含水饱和度状态，测量对应的电阻率变化情况并拟合确定n值。为了实现不同饱和度状态，常用方法有油/气驱法、离心法、烘/风干法、半渗透隔板法和自吸增水法等。油/气驱法分别采用模拟油或真实原油、空气驱替饱含水样品，在不同驱替压差下测量电阻增大率、含水饱和度。离心法是利用离心机产生的离心力达到驱水或驱油的过程，该方法对离心机的最高转速、温度控制要求很严，而且存在离心后样品中水的分布状态不确定、影响电阻率测量精度等难题。隔板法在岩样一端用半渗透隔板夹持，采用两相驱替，在一定压力下计量通过的润湿相体积和电阻率。这3种方法都是采用两相流驱替的手段在岩心中建立不同饱和度状态，但影响测量结果可靠性的关键因素既包括计量装置的精度，也包括驱替过程中压力是否平衡、岩心内部的流体分布是否均匀等。自吸增水法是利用岩心的毛管压力自动吸水的原理，将岩心浸泡在配制的地层水溶液中，在不同时间节点重复测量岩心的重量和电阻率，计算含水饱和度变化情况。
根据前文对青一段储集层特征描述，富含黏土和页理缝的特征，使得常规柱塞岩样的加工和洗油等预处理过程均难以实现，本文采用线切割并包裹热缩膜的工艺加工完成了15块直径为25.4 mm的柱塞样。一方面，富含黏土矿物及大量晶间孔的存在使得有机溶剂洗油过程很难彻底完成，且容易形成次生裂纹从而破坏岩心完整性；另一方面，由于页理缝的大量存在，两相驱替过程中润湿相流体将会沿页理缝方向形成优势通道，其他孔隙中所含流体得不到有效驱替，因此常用的两相驱替岩电实验不适用于本文研究。
根据以上分析，本文研究过程中加工制作两组平行样品各15块，将自吸增水法测量结果与数值模拟结果进行对比，其中自吸增水岩电实验结果如图6 所示。图6 表明，根据自吸增水法测量结果确定的胶结指数ml=1.14，饱和度指数nl=2.17。但是，图6 中的两个图版数据点发散严重，相关系数偏低，主要是因为多数页岩样品的渗透率和自吸增水过程中吸水量极低、计量误差大（能够加工成型的样品仅含少量页理缝）。另外，对于页岩等复杂孔隙结构的储集层，自吸增水法建立含油饱和度的能力有限，最高含水饱和度一般不超过50%[23]，而且用于保护岩心的热缩膜会影响增水饱和过程中流体分布的均匀性，因此该实验结果并不能用于实际井资料处理。实际上近年来与国内页岩油地质勘探不断突破形成鲜明对比的是，针对此类储集层岩电参数的研究并没有取得实质性进展，多数工作仍然遵循上述常规流程，相关的企业、行业标准也尚未建立，亟需探索新的研究思路。
3.2 数字岩心模拟确定关键岩电参数
针对页岩等非常规储集层的岩石物理研究，近年来国内外很多学者都提出了“数字岩心+数值模拟”的新思路，基本流程就是借助CT扫描成像构建数字岩心格架进而开展电阻率等参数的数值模拟。CT图像分辨率与样品尺寸成反比，页岩储集层广泛发育纳米孔，纳米级成像要求样品尺寸不超过0.1 cm，另一方面，页岩等非常规储集层具有较强的微观非均质性，基于小尺寸样品的模拟结果与储集层宏观岩石物理性质相差较大[24]。针对这一问题本文提出构建全尺寸四维数字岩心的新方法，在不减小柱塞样品尺寸的前提下准确描述并定位各类纳米孔。
3.2.1 页岩储集层四维数字岩心建模
为兼顾页岩储集层非均质性和实验配套性等因素，选择直径为2.54 cm柱塞样作为研究对象，对应CT图像分辨率约为7.9 μm（见图7 a）。以C5#样品为例，采用阈值分割和二值化处理后的图像孔隙度为3.1%（见图7 b）。该样品气测孔隙度为9.02%，这说明由于分辨率限制，大约2/3左右的纳米孔没有被识别，这样的数字岩心不能用于电阻率模拟[25]。
对柱塞样端面进行抛光和镀碳层处理后，选择1块代表性方形区域并划分成若干个小的方形子区域，对每一子区域采用聚焦离子束扫描（FIBSEM）获得分辨率为10 nm的图像，再将这些图像拼接成1张超高分辨率、超大像素的图像，这一过程称为MAPS成像，可用于识别如图3 所示的各类孔隙。
为明确青一段不同矿物的纳米孔发育程度，对上述MAPS成像区域同步开展QemScan矿物测量，将两种测量结果进行对比（见图8 ），可以确定测量区域不同位置的组分类型及其孔隙率（面孔率），还可以界定MAPS图像上不同组分（包括有机质在内的各种矿物、孔、缝等）对应的灰度级别，结果如表1 所示。表1 中，阈值1、阈值2分别表示不同组分的灰度值下限、上限；体积比例是指基于某一组分的灰度阈值区间在2.54 cm（1 in）柱塞样微米CT三维图像中统计所占比例，孔隙率是指基于MAPS图像确定的该组分孔隙发育程度。数字岩心总孔隙度就是将不同组分所占的体积比例乘以其对应的孔隙率并累加，计算公式为：
${{\phi }_{\text{t}}}={{V}_{\text{p}}}+0.2{{V}_{\text{org}}}+0.15{{V}_{\text{il}}}+0.07{{V}_{\text{ca}}}$
表1 C5#样品不同组分的灰度级别映射及体积占比统计表 |
| 组分类型 | 阈值1 | 阈值2 | 体积 比例/% | 孔隙率 | 数字岩心总孔隙度/% | 柱塞样气测孔隙度/% |
|---|---|---|---|---|---|---|
| 孔隙 | 0 | 69 | 1.25 | 1.00 | 8.44 | 9.02 |
| 有机质 | 70 | 75 | 2.53 | 0.20 | ||
| 伊利石 | 76 | 90 | 38.80 | 0.15 | ||
| 石英/钠长石 | 91 | 115 | 44.40 | 0 | ||
| 方解石 | 116 | 186 | 12.40 | 0.07 | ||
| 黄铁矿等 | 187 | 255 | 0.69 | 0 |
从表1 可以看出，采用这种方法得到的四维数字岩心，其孔隙度约为8.44%，与气测孔隙度9.02%基本接近，说明已经将微米CT图像丢失的纳米级孔基本上找到并合理归位，据此可以开展数值模拟。需要强调的是，表1 中不同组分的分类及灰度阈值划分方案并非固定，对于不同地区、甚至同一地区不同层位都需要通过上述配套实验分析来具体确定。
3.2.2 基于四维数字岩心数值模拟确定饱和度模型参数
采用上述四维数字岩心构建方法完成6块数字岩心建模，获得类似图7 c的数字岩心格架，据此开展有限元法电阻率模拟，分析其导电特征。
有限元是1种高效能的数值计算方法，其基本原理是将偏微分方程转化为代数方程组求解，可应用于任何微分方程所描述的各类物理场模拟。针对全尺寸的四维数字岩心，假设在两端施加一定的电压并且电流保持恒定，按顺序依次扫描每一个像素单元，根据其对应的组分类型和孔隙率计算对应的等效电阻率，进一步拼装可得到整个岩心的等效电阻率[25]。限于篇幅，这里不再详细论述模拟计算过程。
首先假设所有孔隙完全饱含水（对应电阻率为Rw），模拟计算结果记为R0sim，其数值大小与孔隙结构指数m相关，并应与岩心完全饱含水测量的电阻率R0严格相等或误差满足精度要求。如果不相等则调整参数，通过反复循环直至达到要求后就相应地确定了该岩心的胶结指数和地层因素值。
本文研究的实验流程中，考虑到页岩油样品洗油效果差，采用烘干方法去除样品中的残余油，并结合来样和烘干后两种状态的二维核磁实验判断是否存在残余油，然后加压饱和盐水测量R0，其目的就是为了准确确定每块岩心的地层因素值。这一过程一方面是为了保证模拟的电阻率参数能够反映实际岩心电学特征，另一方面也说明采用基于小尺寸样品更高分辨率CT图像开展数模的做法行不通，因为其模拟结果缺乏实验刻度。6块岩心数值地层因素-数字岩心孔隙度关系模拟结果见图9 a。
进一步利用图像的腐蚀和膨胀运算模拟油驱水的过程，具体原理见文献[26]。就青一段储集层而言，页理缝的尺寸相对最大，也是重要的储集空间，模拟过程首先假设页理缝中的饱和盐水不同程度地被油替代，然后依次考虑残余粒间孔、晶间孔，按照孔径从大到小模拟差异化富油的过程，对每块岩心完成最多11个饱和度节点的电阻率模拟，6块岩心的电阻增大率-含水饱和度关系模拟结果见图9 b。
3.3 饱和度模型的适用性与精度分析
综合以上分析，针对研究区青一段页岩储集层的孔隙结构和岩石矿物特征，选择公式（1）计算含水饱和度，并取m=1.43，n=1.42，首先选择取心井的数据来验证该模型的精度。图10 为GY2HC井青一段饱和度处理结果，图中1—7道的曲线与图4 相同，第8、9两道中，Sws和Swl是基于Simandoux模型分别采用数模法和自吸增水实验确定的岩电参数计算的含水饱和度。分析二者的形态特征可以看出，在青一段纵向上沉积背景、岩石物性和电性变化不大的前提下，其含油性在该井段也不应该有显著变化，但Swl的数值变化区间在30%~90%，局部达到100%，这样大的变化显然不符合该段储集层的地质沉积背景，而基于数值模拟结果的Sws曲线数值分布在30%~60%，纵向变化相对较小，更符合研究区的地质背景。通过与本井段225块岩心含水饱和度数值结果（红色杆状线）对比，Sws、Swl曲线的平均绝对误差分别为4.6%和11.2%，说明采用Simandoux模型和数模方法确定的岩电参数可以准确计算青一段页岩含油饱和度，其精度满足现行国家储量规范关于含油饱和度计算精度的要求。
通过以上两口实际井资料的处理结果分析表明，Simandoux模型能够准确描述青一段富含泥质的页岩储集层岩电响应规律，采用数值模拟确定的岩电参数可以准确计算青一段页岩储集层的含油饱和度。
4 青一段含油饱和度分布规律
页岩储集层的含油性受沉积环境、TOC等多重因素控制，因此测井评价的结果应当与地质规律相符。关键井的取心实验测试结果已验证该含油饱和度模型的计算精度，进一步通过系统取心井的地质研究，可以判断本文方法模型计算结果的合理性。图12 是研究区全井段取心井GY8HC井青一段Q6—Q9层组的处理结果，图中第9道为本文方法计算的饱和度结果，第10道为层组划分结果。何文渊等[19]对该井取心段的储集层特征和实验测试数据统计分析揭示，Q6—Q9层组含油饱和度主要为30%~70%，平均为48.4%，并且Q8层组含油性差，Q9层组含油性较好。由图12 可见，第9道的含油饱和度数值整体也主要分布在30%~70%，并且在Q8层组数值最低，纵向含油分布规律与前人认识完全一致，说明采用上述方法模型计算的饱和度是合理的。
通过多井处理和对比，结合地质研究绘制了研究区青一段页岩储集层含油分布平面图。图13 是基于古页1井区16口井青一段Q1层组测井计算的含油饱和度平均值绘制的平面等值线图，对比图1 b和图13 可以发现，位于沉积中心的古页30和古页16井Q1层组含油饱和度较高，平均约为47%，而相对偏离沉积中心的古页1井Q1层组平均含油饱和度仅36%，这说明沉积中心的优质烃源岩更为发育、含油性更好。
5 结论
松辽盆地齐家—古龙凹陷青一段发育富含黏土和页理缝的页岩油储集层，黏土矿物含量高，并以分散状伊利石和长英质、碳酸盐岩颗粒混杂堆积为主。页理缝、残余粒间孔、粒内溶孔和黏土晶间孔是主要的储集空间和孔隙类型，有机孔占比小。为准确计算其页岩油储集层含油饱和度，基于Simandoux模型，通过构建四维数字岩心和数值模拟，提出了确定页岩岩石电性参数和含油饱和度计算的方法。首先通过聚焦离子束扫描成像及扫描电镜矿物定量评价分析配套测试和对比分析，对纳米孔准确描述和定位，构建四维数字岩心格架，采用有限元模拟确定页岩储集层岩石电性参数。然后通过与岩心数据和二维核磁共振测井饱和度结果对比表明，采用四维数字岩心+数值模拟得到的富含黏土型页岩油储集层含油饱和度公式误差符合现行国家储量规范要求。经多井对比分析表明，在单井纵向上和多井平面上，计算的含油饱和度与取心样品分析化验结果及区域沉积规律吻合，进一步验证了该方法模型较高的精度和良好的适用性。
符号注释：
a，m，n——Archie公式中的参数，无因次；ml，ms——采用自吸增水法实验和数字岩心数值模拟法确定的Simandoux模型参数，无因次；nl——自吸增水法岩电实验测量确定的饱和度指数，无因次；R0——岩石完全饱和水时测量的电阻率，Ω·m；R0sim——岩石完全饱和水时模拟计算的电阻率，Ω·m；Rcl——纯黏土的电阻率，Ω·m；Rt——地层真实电阻率，Ω·m；Rw——地层水电阻率，Ω·m；S1——游离烃含量，mg/g；Sw——地层含水饱和度，%；Sw_cpre——岩心实验分析的含水饱和度，%；Swl——基于Simandoux模型采用自吸增水法实验测量的岩电参数计算含水饱和度，%；Sws——基于Simandoux模型采用数值模拟的岩电参数计算含水饱和度，%；Sw_T1T2——二维核磁共振测井计算的含水饱和度，%；Tl，T2——核磁共振纵、横向弛豫时间，ms；TOC——总有机碳含量，%；Vcl，Vp，Vorg，Vil，Vca——黏土矿物、粒间孔隙、有机质、伊利石、方解石的体积含量，%；ϕ——地层岩石有效孔隙度，%；ϕt——总孔隙度，%。