我国煤炭资源储量丰富且开发前景广阔^[1,2]。煤层气含量是评价煤层含气储层参数之一^[3],其与煤层气井单井产量预测、资源开发前景及资源勘探息息相关^[4],在实际开采中为防止盲目开发带来的损失需进行区块储量预测^[5,6]。煤层气作为非常规天然气,与常规天然气在储层中的储集与渗流机理存在不同^[7],煤层气含量受控因素复杂繁多,与煤阶变质程度、温度、压力、有效埋深、厚度、构造特征,水文地质特征等地质因素相关联^{[8,9,10,11,12]}。对煤层气含量的预测自几十年前起就有学者进行了相关研究,煤层气含量的测定方法多样,相对直接且准确的是对煤心进行实验测量,但因煤心实测资料较少难以直接使用。继而可根据等温吸附理论、地质统计分析法等对气含量进行估算^[13]。国内外煤层气研究学者相继提出了多种煤层气含量评价预测方法,1997年由Kim A G.提出Kim法,将水分、灰分含量结合煤层压力、温度及平衡水状态校正量等进行煤层吸附气含量计算^[14],后基于这一方法将煤心工业组分含量引入并进行相关分析,提出改进Kim方程;Ahmed U^[15]等通过建立等温吸附模型对煤层气等温吸附线进行描述;Hawkins J M等^[16]根据Langmuir等温吸附理论,进一步提出用兰氏煤阶方程计算方法进行煤层气含量预测。

除地质方法外,利用测井手段求取煤层含气储层参数也逐步成为研究热点。由于煤层气含量受多方面因素影响,测井与煤层气含量的关系可能为线性或非线性。预测煤层气含量的方法多以回归分析、BP神经网络等算法为主。回归分析法多为通过分析常规测井曲线,建立测井与煤层气含量相关关系,预测煤层气含量,这类模型预测效果良好,有较好的应用前景^[17,18];邵先杰等^[19]基于煤岩工业组分间的关系及补偿密度、声波时差等的测井响应,对韩城矿区煤样应用多元回归分析和体积模型法建立了煤储集层含气量测井解释模型,计算结果匹配度好;黄兆辉等^[20]针对沁水盆地目标层,分析煤工业组分与相关测井响应特征,结合改进Langmuir煤阶方程对煤质与含气量关系进行建模;金泽亮等^[21]分析沁水盆地柿庄区块煤层气储层地质特征及测井资料,利用兰氏方程和多元线性回归导出了煤层气含量计算的经验公式,其结果具有可信度与有效性。在利用线性关系难以准确建模时,可利用神经网络等方法进行煤层气含量计算。将电阻率、体积密度、自然伽马等一系列测井曲线参数,通过神经网络训练预测值和期望值的吻合程度计算煤层气含量,该算法的非线性逼近能力强,对于已知存在某种潜在的联系但又无法用确切方程或算法表达的求解问题有其独到之处^[22]。侯俊胜等^[23]根据煤层气储层及其测井响应特征将神经网络方法应用于测井资料解释中,取得了满意的解释结果;潘和平等^[24]提出利用非线性方法预测煤层气含量,通过分析大量数据验证表明BP神经网络预测煤层气含量精度比较高;吴东平等^[25]应用BP神经网络对山西柳林杨家坪地区煤层进行煤层气含量预测,效果良好;连承波等^[26]将支持向量机引入煤层气含量的评价预测中。

在实际应用中,由于各测井参数对煤层气含量的灵敏度不同且组合方式存在差异,因而回归分析法得出的模型计算结果与实际数据存在一定偏差,误差分析数据不足导致应用效果受限,BP神经网络由于训练的复杂性并且对数据样本要求较高,使用存在局限性。针对上述问题,本文将灰色关联技术这一数学方法结合测井曲线参数进行煤层气含量预测,即通过计算各测井曲线参数与煤层气含量的灰色关联度并结合GM(0,N)模型进行煤层气含量预测,并用实际数据来验证本文方法的有效性和实用性。

对数据分析时使用灰色关联法所得的结果直观且对数据样本数量、数据形态及数据间隔要求较低,并且对于无典型分布规律的数据,定量计算与定性分析结果相符合,实用性强。最早且应用最为广泛的灰色关联法是由邓聚龙教授提出的一般关联度,基本思想为将各序列曲线的几何形状进行对比,根据其相似程度来比较判断各序列联系的紧密程度,曲线几何形状接近,则序列间关联度大,反之形状不接近则关联度小。

一般关联度的具体形式为规定一参考序列x₀与一对比序列x_i,其形式分别为

其中i∈[1,∞),将x₀与x_i之间的一般关联度定义为:

式中: maximaxk|x₀(k)-x_i(k)|为参考序列与对比序列间的两极最大差; minimink|x₀(k)-x_i(k)|为参考序列与对比序列间的两极最小差;ρ为分辨率系数,ρ∈[0,1],多取ρ=0.5。

一般关联度于实际应用中在规范性与无量纲化保序性上存在缺陷,继而有学者提出斜率关联法。该关联法的基本思想为根据参考序列与对比序列在同一时段上各序列曲线斜率的接近程度来判定序列之间的关联度。若参考序列与对比序列在各时段上斜率方向相同且数值接近则两序列间关联程度大,反之则相对小。斜率关联度的计算公式为^[27]:

其中k∈[1,n-1],i∈[1,n]。

后对上述斜率关联法进行改进引入δ(k)函数,使其既可以反映序列间的正关联情况,也能反映负关联情况^{[28,29,30,31]}。本文在此基础上,结合实际情况,根据参考序列煤层气含量的分布情况进行分段加权平均,加权系数取决于各段上样本数量的多少。本文假定参考序列根据数值范围分选至m个区域,根据各区域样本数量的多少进行加权平均,计算公式可写为:

其中 Δ̅0w=1∑l=1mwl∑l=1m∑k=1n-1wl|x0(k+1)-x0(k)|;Δ̅iw=1∑l=1mwl∑l=1m∑1n-1wl|xi(k+1)-xi(k)|;δ(k)=±1,当(x0(k+1)-x0(k))(xi(k+1)-xi(k))≥0时值为1,当(x0(k+1)-x0(k))(xi(k+1)-xi(k))<0时值为0。

GM(0,N)模型为不含导数的静态模型,该模型与多元线性回归模型有着本质上的区别,一般的多元线性回归模型以原始数据为基础,GM(0,N)模型则是以原始数据一次累加的新数据序列为建模基础^[29]。GM(0,N)模型的建模形式弱化了原始数据间的随机性,将数据规律性提高^[32],且该模型对统计数据数量要求小,在实际生产中更能满足使用条件,这也是本文对工区数据使用GM(0,N)模型的原因之一。

规定 x1(0)={ x1(0)(1), x1(0)(2),…, x1(0)(n)}为系统特征序列,将相关因素序列规定为:

x2(0)={x2(0)(1),x2(0)(2),…,x2(0)(n)},x3(0)={x3(0)(1),x3(0)(2),…,x3(0)(n)},︙xN(0)={xN(0)(1),xN(0)(2),…,xN(0)(n)},

xi(1)为 xi(0)的一次累加生成序列,即1-AGO序列。可将GM(0,N)模型写成:

预测步骤为:

1) 求取各序列的1-AGO(一次累加)序列;

2) 构造数据矩阵:

B=1,x2(1)(2),…,xN(1)(2)1,x2(1)(3),…,xN(1)(3)︙1,x2(1)(n),…,xN(1)(n),Y=x1(1)(2)x1(1)(3)︙xN(1)(n);

3) 计算参数列: a^=[a,b₂,…,b_N]^T

参数列计算的最小二乘估计公式为:

4) 按照式(6)建立GM(0,N)模型;

5) 还原数据列得到预测值:

按照式(6)计算得到的数据无实际物理意义,通过序列内数据累减得到还原数据序列,公式为:

若要预测系统特征数据第n个值以后的值,仍可通过上述步骤建立GM(0,N)模型计算得到累加数据后进行累减逆生产即可^[33]。

本文使用的工区实例为山西省沁水煤田柿庄北区15号煤层气含量数据,结合本文上述方法对煤层数据进行处理并预测煤储层气含量并与邻近工区使用的多元回归法进行对比^[34]。

沁水煤田是中国目前产煤最多的大型石炭—二叠纪煤田^[35]。沁水煤田煤层气资源丰富,约占全国煤层气资源量的1/4,其勘探程度在国内最高、储量条件稳定、开发潜力巨大,新一轮资源评价结果显示,沁水盆地煤层气地质资源量大于1×10¹²>m^3[36]。本文所用数据其勘探区于沁水煤田柿庄北区15号煤层。勘探区共取得煤层气含量测试样品40个。将40个测试样品气含量及对应的标准化以后的自然电位、自然伽马、补偿密度、声波时差、补偿中子、深浅侧向电阻率及深度收集制表(表1)。

煤层气含量理论上受煤层有效埋深及煤层顶底板岩性等多个地质因素所控,而这些地质因素决定了煤层测井响应特征。因此,不同测井曲线上会因煤层气含量的变化出现不同响应。根据研究区柿庄北区15号层的40组煤储层气含量测试样品,能对煤层气测井响应与煤层气含量做相关性分析(图1),找到与煤层气含量相关性相对明显的测井曲线,根据拟合优度建立关联序(表2)。再应用改进的斜率关联法计算各测井参数与煤层气含量间的关联度与关联序(表3),为消除主控因素不同量纲的数据对预测结果的影响,进行关联度计算前需要对数据进行初值化处理^[18]。