0 引言

铀矿是关系国家安全、国计民生的重要战略资源，对其的需求与日俱增。近年来，随着低成本、低污染的地浸采铀技术的成熟，砂岩型铀矿逐渐成为中国能源部门重要的找矿类型^[1-4]。在勘探不断深入的同时，面临的地质问题也日趋复杂。为此，进入21世纪以来, 有关部门也开展了针对砂岩型铀矿的地震勘探工作，主要目的是查清含矿层的形态和构造特征。在地震反演方面，主要应用叠后波阻抗反演预测有效储层^[5-9]。但是，有关铀矿体的预测技术的相关研究报道十分有限，其中一个主要原因是铀矿体的厚度普遍较小(一般小于10m), 常规确定性反演的纵向分辨率很难满足精度需要。

地质统计学反演的概念最早由Haas等^[10]在1994年提出，到了1998年Dubrule等^[11]将该技术用于三维地震数据体反演。地质统计学反演技术由地质随机建模与地震数据共同驱动，可以将各类地质信息和测井资料融入反演，突破了地震频带宽度的限制，可实现纵向高精度表征；同时利用地震资料横向信息丰富的优势，反演结果充分展示了储层等信息的横向变化及非均质性^[12-16]。该技术广泛用于油气勘探、开发领域，预测结果与井上地质信息吻合更好，大幅提高了薄储层的表征精度^[17-20]。在铀矿勘探领域，人们也尝试应用地质统计学反演预测铀矿体的展布^[21]，但是由于其预测结果的随机性较强，实际应用效果局限性较大。地质统计学反演的可靠性主要受控于变差函数等关键参数的选取^[22-24]，传统的求取变差函数横向变程方法的随机性较强。

为了降低地质统计学反演的随机性，提高有效性和反演精度，需要在求取横向变程时融入更多的先验信息并赋予变差函数更明确的地质含义。在油气储层反演领域，樊鹏军等^[25]认为对变差函数赋予明确的地质含义可以大幅提高反演可靠性和精度。鉴于此，本文将该思想引入QJD大型典型砂岩型铀矿的矿体预测。该区在早期的石油勘探中采集了大量的三维地震资料，在多年的勘探过程中还积累了丰富的地质、钻井、测井以及含矿性等资料。因此，完全可以充分利用该地区勘探程度高、井数多、井网分布均匀等先验信息和三维地震资料丰富的横向信息预测铀矿体，尤其是在测井数据、矿体规模和含矿性等数据的充分控制下提出准确的变差函数求取方法，可降低地质统计学反演的随机性，进而更精确、有效地反演矿体形态。为此，首先阐述地质统计学反演的基本原理；其次，阐述变差函数的定义及地质意义；然后，提出新的变差函数求取方法；最后，分析反演结果的可靠性和实际应用效果并给出相关结论。

1 地质统计学反演的基本原理与实现流程

面向目标参数估计的地质统计学反演问题实际上可以归结为一个贝叶斯参数估计问题，即：在一些观测信息的基础上，通过不断更新先验信息得到目标参数估计问题的条件最优解。在贝叶斯推论框架下，地震数据反演问题实际为测井数据约束下基于地震数据的目标参数反演问题，其表达式为

$ P(R\left| {L, S) = } \right.P(R\left| {L)} \right.P(S\left| R \right.)/P(S) $

(1)

式中：R为待反演的目标参数分布；L为测井数据；S为地震数据。式(1)的含义为：在同时满足L和S条件下的R的后验概率P(R│L, S)为已知L的条件概率P(R│L)与S的似然函数P(S│R)/P(S)的乘积。由于概率分布的复杂性，式(1)无法求得解析解，但可以采用MCMC(Markov Chain & Monte Carlo)方法进行概率评价。面向砂岩型铀矿的地质统计学反演步骤如下(图 1)。

(1) 基础数据分析，包括对地震、测井、矿体等数据的分析。地震数据分析包括资料面积分布、资料规则化、资料主频和频宽等分析，确定地震资料是否满足目标预测的需要，同时可以对成果数据进行提频、去噪等预处理；测井资料分析包括井的数量、空间分布、曲线完整性、标准化等分析，核心是分析测井属性的岩石物理敏感性，进而确定目标反演参数，其中面向砂岩型铀矿矿体的目标参数是自然伽马(GR)数据。

(2) 结合井震数据对反演目标实施精细井震标定、精细构造解释和构造地质建模。该地质模型构成整个反演的模型框架，并采用三维网格模型将其离散，构成模型参数的分布空间，进而获取反演的初始模型。

(3) 求取纵向、横向变差函数。统计、分析测井数据求取纵向变差函数，统计、分析已开发区矿体规模和含矿性信息求取横向变差函数，分别作为地质统计学的纵向和横向约束条件，构成了反演的核心参数。

(4) 利用MCMC算法获得统计意义下的储层参数随机样点分布。MCMC算法步骤为：①根据步骤(1)~步骤(3)确定初始模型和反演参数；②通过后验概率密度函数采样随机生成一个新的模型；③对新模型实施正演计算得到合成地震数据，并利用实测地震数据计算上述新模型的似然函数；④求取新模型的先验概率，再乘以步骤③得到的似然函数，得到新模型的后验概率；⑤对比新模型的后验概率与当前目标模型的后验概率，若前者大于后者则接受新生成的模型，并将Markov链移动到新模型位置，否则随机决定是否保留；⑥重复步骤②~步骤⑤，直到后验概率不再改变，循环结束，并输出计算结果。

2 变差函数的求取方法 2.1 变差函数的地质意义

从上述的地质统计学的基本原理与反演实现过程发现，选择合适、有效且最能刻画目标参数的先验信息尤其重要，涉及求取纵向和横向变差函数。同时，在地质统计学反演中，变差函数具有重要的地质意义。在铀矿体预测方面，变差函数实际上反映了矿体在三维空间的变化特征，表征矿体的空间各向异性。其中，纵向变程反映矿体垂向厚度，其值决定反演的纵向分辨率；横向变程则反映矿体的横向发育规模，其不同方向取值反映储层的空间各向异性，长轴方向代表矿体长度，短轴方向代表矿体宽度。鉴于变差函数具有实际地质指导意义，因此其求取方法非常重要。准确求取变差函数是降低地质统计学反演随机性的重要途径，是反演结果有效性和精度的重要保证。

2.2 纵向变程的求取

砂岩型铀矿体厚度较小且变化快，在预测矿体时过大的纵向变程会使反演结果不能精确地反映矿体厚度，而过小的纵向变程则会成倍增加计算时间并增大反演结果的随机性。为此，本文在求取纵向变程时，提出了一种结合测井资料与实际矿体形态的综合统计方法，分析工区内均匀分布的20口井的GR曲线样本点。目的层姚家组下段地层厚度普遍大于60m，测井数据采样间隔为0.125m，因此样本点个数远大于(一般要求的)50(图 2)。QJD矿区31个主要矿体的最大厚度为7.70m，最小厚度为2.25m，平均厚度为4.48m，综合测井数据和实际矿体厚度，确定研究区的纵向变程为4.5m。通过分析测井数值得到的变差函数(图 3)可以看出：数据点间相关性很明显，且当纵向变程大于4.5m时，变差函数曲线趋于平缓。上述计算纵向变程的方法以大量实测数据统计、分析结果作为反演的纵向约束，尽可能在已有先验井数据控制下得到符合地质意义的反演结果，进而最大限度地降低反演的随机性。

2.3 横向变程的求取

利用常规方法求取横向变差函数时，由于井网密度很难达到横向采样要求，通常提取叠后稀疏脉冲反演平面属性求取一个大致的变差范围。本次研究首先对该方法进行测试，结果表明：在矿体横向预测方面，反演结果横向连续性较差，与实际矿体展布情况差别较大。说明该方法虽然弥补了井曲线横向样本不足的缺陷，但是其横向分辨率有限且没有充分结合变差函数的地质意义，所以很难解决地质统计学反演的横向不确定性问题。

针对上述问题，基于QJD铀矿床勘探程度较高且对矿区的沉积模式、矿体规模和成矿规律等的较成熟认识，并结合变差函数的地质意义，提出一种在已开发区矿体成矿模式和矿体展布规律先验信息约束下的横向变程求取方法。该方法的核心思想是：首先在有大量开发井的邻近已开发区域建立矿体规模和单矿点最高平米铀量之间的拟合关系，然后在反演的目标区域利用该拟合关系和反演约束井的最高平米铀量，计算反演目标区矿体规模的横向变程作为反演的核心参数(图 4)。为了建立上述拟合关系，首先需要确定矿体合并原则并圈定矿体，然后统计、分析矿体规模与矿体含矿性关系，最终在反演目标区域求取横向变程，详述如下。

(1) 矿体合并原则

结合研究区储层和区域隔水层发育特征，将目的层姚家组划分为6套含矿层系，在各含矿层系内将可渗透砂岩中的异常层合并，将合并后平米铀量大于1kg/m²且品位大于0.1‰的矿体定为工业层，将平米铀量小于1kg/m²且品位大于0.1‰的矿体定为矿化层。本次研究以工业层为主要对象。

(2) 圈定矿体

在充分考虑控矿因素的基础上，参考相关行业规定，将各工业井的井点位置作为矿体中心位置的水平投影点，矿体的边界线由各邻井的外推点直接相连而成，连接时走向点与走向点相连、倾向点与倾向点相连。工业井与矿化井之间按勘查工程间距的1/2平推；工业井与无矿井之间则按勘查工程间距的1/4平推。

(3) 矿体规模与矿体含矿性的关系

分析矿体规模与矿体含矿性的关系，首先圈定矿体的空间展布形态并确定矿体规模(长度和宽度)，然后标注各矿体的含矿性(平均平米铀量)，最后通过拟合两者的数值关系即可建立矿体规模与矿体含矿性的关系(图 5)。统计、分析QJD地区上百个矿体的规模与含矿性关系表明：在储层非均质性变化不大的情况下，矿体的规模(长度及宽度)与矿体的平米铀量呈明显的正相关，即矿体规模越大则含矿性越高。

(4) 求取横向变程

在步骤(1)~步骤(3)的基础上，即可求取横向变程，具体为：

(a) 通过步骤(3)建立矿体规模与矿体最高平米铀量之间的拟合关系

$ l = 169.9\ln x + 102.7 $

(2)

$ b = 7.2x + 60 $

(3)

式中：l为矿体长度；b为矿体宽度；x为矿点的最高平米铀量。

(b) 在反演的目标区域，利用20口井资料，基于式(2)和式(3)，在已知单矿点x的情况下，求取研究区各矿体的规模。

(c) 对比步骤(b)的结果与初始GR模型，反复修正横向变程。

(d) 判断基于当前变程的GR模型与矿体展布规律的吻合程度。若吻合程度达到要求，则输出横向变程；否则重复步骤(a)~步骤(d)，直至吻合程度达到要求。

基于步骤(a)~步骤(d)，最终确定表征研究区矿体长度和宽度的横向变程分别为700、300m。

3 铀矿体预测精度分析及实际应用

为了评估基于本文提出的变差函数求取方法的铀矿体反演精度及其实际应用效果，选取面积约为50km²的研究区展开试验。研究区位于中国QJD矿区(图 6)，该区实际地震资料与井中合成地震记录吻合良好(图 7)，其中矿区北部是近年QJD铀矿床勘探的重点区域。随着勘探程度的提高，矿体逐渐接近边界，勘探成功率呈下降趋势(近两年的成功率不足40%)，迫切需要具有针对性的预测手段提高勘探效率，降低勘探成本。研究区构造北高南低，层间氧化带由南向北发育，目的层姚家组主要为辫状河三角洲平原相地层，砂体厚度大，储层物性变化较小，目前共完钻探井200余口，根据钻井资料可以较准确地刻画矿体的展布特征。

3.1 纵向变程对反演结果的影响

根据提出的纵向变程求取方法，求得该区的纵向变程为4.5m。为了验证所取纵向变程的合理性，对比、分析了不同纵向变程的反演结果(图 8)，可见：①当纵向变程为15.0m时，反演结果分辨率很低，不能清晰地反映纵向的多个矿体，且矿体边界较圆滑，与尖灭的矿体形态不符(图 8a)；②当纵向变程为10.0m时，反演结果纵向分辨率略有提高，矿体边界逐渐清晰，但对尖灭的矿体形态的刻画仍不够清晰(图 8b)；③当纵向变程为4.5m时，反演结果的纵向分辨率较高，位于中间的主矿体形态清晰，且清晰地刻画了其上、下的两套小矿体，较真实地反映了矿体形态(图 8c)；④当纵向变程为2.0m时，反演结果的纵向分辨率提高不明显。因此，纵向变程明显影响反演结果，基于本文方法求取纵向变程的反演结果很好地刻画了矿体(包括小矿体)形态。

3.2 横向变程对反演结果的影响

为了测试本文求取横向变程方法的效果，对比、分析了不同求取横向变程方法的反演结果(图 9)，可见：①基于叠后稀疏脉冲反演平面属性进行横向约束方法(图 9b)、本文方法(图 9c)获得的反演结果均刻画了Y1旋回矿体，并且预测的矿体垂向厚度与矿体实际厚度(图 9a)基本吻合。②在横向上，图 9b将旋回顶部的一套完整矿体描述为多个独立小矿体，既不符合矿体合并原则也不符合矿体展布规律；图 9c将旋回顶部的矿体描述为同一套矿体，更好地反映了矿体的连续性，更符合矿体展布规律。图 10为实际矿体平面展布与本文求取横向变程方法反演结果的GR均方根属性平面图。由图可见：①两者均呈现南北向的展布特征；②在红色虚线圈定区域的反演结果中，不同规模矿体的主体空间位置、含矿性(图 10b)与已知井区的实际矿体空间位置、含矿性(图 10a)对应关系很好。综上所述，基于本文方法求取横向变程的反演结果具有很高的横向分辨率，能准确地刻画矿体形态。

3.3 反演结果验证

前文分析了纵向、横向变程对反演结果的影响，分析结果证明了本文提出的变差函数的有效性。但在将反演结果用于实际勘探之前，必须由验证井验证反演结果。为了验证反演结果，本文采用去除已知井进行验证的方法。反演时选取e井、f井、g井、h井所在区域进行测试，测试过程仅采用e井、g井、h井数据进行约束反演得到反演结果(图 11a)。可见，预测矿体与实际矿体的对应关系较好，反演结果准确预测了发育于f井目的层底部的矿体，预测矿体与e井连通，在g井附近尖灭，与实际矿体剖面(图 11b)吻合良好，充分证明了本文反演方法的可靠性。

3.4 实际勘探效果分析

前文的精度分析表明，基于本文方法的反演结果，无论在纵向还是横向都能准确刻画矿体的空间展布特征，对于矿体的实际勘探具有重要意义。为了验证本文方法的实际应用效果，通过综合分析反演结果与控矿因素(氧化带发育特征等)发现，在现有矿体北部的一个部署区域(图 10a蓝色虚线框)内还存在矿体发育区，2019年底在该区部署并实施了12口探井，其中9口为工业矿井，见矿率为75%，远高于QJD详查勘探工业见矿率的平均水平。实际勘探成果充分证实了本文方法反演结果的可靠性。

综上所述，当赋予变差函数明确的地质含义后，反演结果与矿体的实际展布特征吻合更好，准确刻画了铀矿体的展布形态，同时也可以根据反演结果在未知区预测矿体分布，指导勘探部署。

4 结论

本文提出了在测井数据和地质意义控制下的变差函数求取方法，并用于QJD砂岩型铀矿体预测的地质统计学反演，取得如下结论：

(1) 地质统计学反演基于测井数据随机建模，相对于常规确定性反演，继承了地震资料横向分辨率高的特点，同时具有更高的纵向分辨率。因此对于厚度小、横向变化快的砂岩型铀矿体的预测效果较好，其中变差函数的求取是影响反演结果的重要因素。

(2) 基于本文提出的结合测井资料与实际矿体形态的综合统计方法求取纵向变程的反演方法，在已有先验井数据控制下得到了符合地质意义的反演结果，进而最大限度地降低了反演的随机性，反演结果很好地刻画了矿体(包括小矿体)的纵向形态。

(3) 变差函数反映了矿体的三维空间展布规律，通过统计、分析矿体规模和矿点含矿性，拟合了两者的关系，提出在已开发区矿体成矿模式和矿体展布规律等先验信息约束下的横向变程求取方法，获得的反演结果与实际矿体分布吻合更好。

(4) 利用本文方法得到的反演结果更准确地预测了矿体形态，参考反演结果部署的勘探井位显著提高了矿体预测的可靠性和钻探成功率。