魏 達 孫章慶

(①中國石油遼河油田分公司勘探開發(fā)研究院，遼寧盤錦 124010；②吉林大學地球探測科學與技術學院，吉林長春 130026)

0 引言

鈾礦是關系國家安全、國計民生的重要戰(zhàn)略資源，對其的需求與日俱增。近年來，隨著低成本、低污染的地浸采鈾技術的成熟，砂巖型鈾礦逐漸成為中國能源部門重要的找礦類型[1-4]。在勘探不斷深入的同時，面臨的地質問題也日趨復雜。為此，進入21世紀以來,有關部門也開展了針對砂巖型鈾礦的地震勘探工作，主要目的是查清含礦層的形態(tài)和構造特征。在地震反演方面，主要應用疊后波阻抗反演預測有效儲層[5-9]。但是，有關鈾礦體的預測技術的相關研究報道十分有限，其中一個主要原因是鈾礦體的厚度普遍較小(一般小于10m),常規(guī)確定性反演的縱向分辨率很難滿足精度需要。

地質統(tǒng)計學反演的概念最早由Haas等[10]在1994年提出，到了1998年Dubrule等[11]將該技術用于三維地震數據體反演。地質統(tǒng)計學反演技術由地質隨機建模與地震數據共同驅動，可以將各類地質信息和測井資料融入反演，突破了地震頻帶寬度的限制，可實現縱向高精度表征；同時利用地震資料橫向信息豐富的優(yōu)勢，反演結果充分展示了儲層等信息的橫向變化及非均質性[12-16]。該技術廣泛用于油氣勘探、開發(fā)領域，預測結果與井上地質信息吻合更好，大幅提高了薄儲層的表征精度[17-20]。在鈾礦勘探領域，人們也嘗試應用地質統(tǒng)計學反演預測鈾礦體的展布[21]，但是由于其預測結果的隨機性較強，實際應用效果局限性較大。地質統(tǒng)計學反演的可靠性主要受控于變差函數等關鍵參數的選取[22-24]，傳統(tǒng)的求取變差函數橫向變程方法的隨機性較強。

為了降低地質統(tǒng)計學反演的隨機性，提高有效性和反演精度，需要在求取橫向變程時融入更多的先驗信息并賦予變差函數更明確的地質含義。在油氣儲層反演領域，樊鵬軍等[25]認為對變差函數賦予明確的地質含義可以大幅提高反演可靠性和精度。鑒于此，本文將該思想引入QJD大型典型砂巖型鈾礦的礦體預測。該區(qū)在早期的石油勘探中采集了大量的三維地震資料，在多年的勘探過程中還積累了豐富的地質、鉆井、測井以及含礦性等資料。因此，完全可以充分利用該地區(qū)勘探程度高、井數多、井網分布均勻等先驗信息和三維地震資料豐富的橫向信息預測鈾礦體，尤其是在測井數據、礦體規(guī)模和含礦性等數據的充分控制下提出準確的變差函數求取方法，可降低地質統(tǒng)計學反演的隨機性，進而更精確、有效地反演礦體形態(tài)。為此，首先闡述地質統(tǒng)計學反演的基本原理；其次，闡述變差函數的定義及地質意義；然后，提出新的變差函數求取方法；最后，分析反演結果的可靠性和實際應用效果并給出相關結論。

1 地質統(tǒng)計學反演的基本原理與實現 流程

面向目標參數估計的地質統(tǒng)計學反演問題實際上可以歸結為一個貝葉斯參數估計問題，即：在一些觀測信息的基礎上，通過不斷更新先驗信息得到目標參數估計問題的條件最優(yōu)解。在貝葉斯推論框架下，地震數據反演問題實際為測井數據約束下基于地震數據的目標參數反演問題，其表達式為

P(R|L,S)=P(R|L)P(S|R)/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)方法進行概率評價。面向砂巖型鈾礦的地質統(tǒng)計學反演步驟如下(圖1)。

圖1 地質統(tǒng)計學反演方法的實現流程

(1)基礎數據分析，包括對地震、測井、礦體等數據的分析。地震數據分析包括資料面積分布、資料規(guī)則化、資料主頻和頻寬等分析，確定地震資料是否滿足目標預測的需要，同時可以對成果數據進行提頻、去噪等預處理；測井資料分析包括井的數量、空間分布、曲線完整性、標準化等分析，核心是分析測井屬性的巖石物理敏感性，進而確定目標反演參數，其中面向砂巖型鈾礦礦體的目標參數是自然伽馬(GR)數據。

(2)結合井震數據對反演目標實施精細井震標定、精細構造解釋和構造地質建模。該地質模型構成整個反演的模型框架，并采用三維網格模型將其離散，構成模型參數的分布空間，進而獲取反演的初始模型。

(3)求取縱向、橫向變差函數。統(tǒng)計、分析測井數據求取縱向變差函數，統(tǒng)計、分析已開發(fā)區(qū)礦體規(guī)模和含礦性信息求取橫向變差函數，分別作為地質統(tǒng)計學的縱向和橫向約束條件，構成了反演的核心參數。

(4)利用MCMC算法獲得統(tǒng)計意義下的儲層參數隨機樣點分布。MCMC算法步驟為：①根據步驟(1)～步驟(3)確定初始模型和反演參數；②通過后驗概率密度函數采樣隨機生成一個新的模型；③對新模型實施正演計算得到合成地震數據，并利用實測地震數據計算上述新模型的似然函數；④求取新模型的先驗概率，再乘以步驟③得到的似然函數，得到新模型的后驗概率；⑤對比新模型的后驗概率與當前目標模型的后驗概率，若前者大于后者則接受新生成的模型，并將Markov鏈移動到新模型位置，否則隨機決定是否保留；⑥重復步驟②～步驟⑤，直到后驗概率不再改變，循環(huán)結束，并輸出計算結果。

2 變差函數的求取方法

2.1 變差函數的地質意義

從上述的地質統(tǒng)計學的基本原理與反演實現過程發(fā)現，選擇合適、有效且最能刻畫目標參數的先驗信息尤其重要，涉及求取縱向和橫向變差函數。同時，在地質統(tǒng)計學反演中，變差函數具有重要的地質意義。在鈾礦體預測方面，變差函數實際上反映了礦體在三維空間的變化特征，表征礦體的空間各向異性。其中，縱向變程反映礦體垂向厚度，其值決定反演的縱向分辨率；橫向變程則反映礦體的橫向發(fā)育規(guī)模，其不同方向取值反映儲層的空間各向異性，長軸方向代表礦體長度，短軸方向代表礦體寬度。鑒于變差函數具有實際地質指導意義，因此其求取方法非常重要。準確求取變差函數是降低地質統(tǒng)計學反演隨機性的重要途徑，是反演結果有效性和精度的重要保證。

2.2 縱向變程的求取

砂巖型鈾礦體厚度較小且變化快，在預測礦體時過大的縱向變程會使反演結果不能精確地反映礦體厚度，而過小的縱向變程則會成倍增加計算時間并增大反演結果的隨機性。為此，本文在求取縱向變程時，提出了一種結合測井資料與實際礦體形態(tài)的綜合統(tǒng)計方法，分析工區(qū)內均勻分布的20口井的GR曲線樣本點。目的層姚家組下段地層厚度普遍大于60m，測井數據采樣間隔為0.125m，因此樣本點個數遠大于(一般要求的)50(圖2)。QJD礦區(qū)31個主要礦體的最大厚度為7.70m，最小厚度為2.25m，平均厚度為4.48m，綜合測井數據和實際礦體厚度，確定研究區(qū)的縱向變程為4.5m。通過分析測井數值得到的變差函數(圖3)可以看出：數據點間相關性很明顯，且當縱向變程大于4.5m時，變差函數曲線趨于平緩。上述計算縱向變程的方法以大量實測數據統(tǒng)計、分析結果作為反演的縱向約束，盡可能在已有先驗井數據控制下得到符合地質意義的反演結果，進而最大限度地降低反演的隨機性。

圖2 工區(qū)測井資料分析(a)測井解釋圖(GR代表伽馬測井，RT代表深雙側向電阻率測井，RS代表淺雙側向電阻率測井)；(b)部分深度段測井數據

圖3 縱向變程分析圖塊金值約等于0，基臺值為0.3，因此空間相關度遠小于25%

2.3 橫向變程的求取

利用常規(guī)方法求取橫向變差函數時，由于井網密度很難達到橫向采樣要求，通常提取疊后稀疏脈沖反演平面屬性求取一個大致的變差范圍。本次研究首先對該方法進行測試，結果表明：在礦體橫向預測方面，反演結果橫向連續(xù)性較差，與實際礦體展布情況差別較大。說明該方法雖然彌補了井曲線橫向樣本不足的缺陷，但是其橫向分辨率有限且沒有充分結合變差函數的地質意義，所以很難解決地質統(tǒng)計學反演的橫向不確定性問題。

針對上述問題，基于QJD鈾礦床勘探程度較高且對礦區(qū)的沉積模式、礦體規(guī)模和成礦規(guī)律等的較成熟認識，并結合變差函數的地質意義，提出一種在已開發(fā)區(qū)礦體成礦模式和礦體展布規(guī)律先驗信息約束下的橫向變程求取方法。該方法的核心思想是：首先在有大量開發(fā)井的鄰近已開發(fā)區(qū)域建立礦體規(guī)模和單礦點最高平米鈾量之間的擬合關系，然后在反演的目標區(qū)域利用該擬合關系和反演約束井的最高平米鈾量，計算反演目標區(qū)礦體規(guī)模的橫向變程作為反演的核心參數(圖4)。為了建立上述擬合關系，首先需要確定礦體合并原則并圈定礦體，然后統(tǒng)計、分析礦體規(guī)模與礦體含礦性關系，最終在反演目標區(qū)域求取橫向變程，詳述如下。

圖4 橫向變程求取流程圖

(1)礦體合并原則

結合研究區(qū)儲層和區(qū)域隔水層發(fā)育特征，將目的層姚家組劃分為6套含礦層系，在各含礦層系內將可滲透砂巖中的異常層合并，將合并后平米鈾量大于1kg/m2且品位大于0.1‰的礦體定為工業(yè)層，將平米鈾量小于1kg/m2且品位大于0.1‰的礦體定為礦化層。本次研究以工業(yè)層為主要對象。

(2)圈定礦體

在充分考慮控礦因素的基礎上，參考相關行業(yè)規(guī)定，將各工業(yè)井的井點位置作為礦體中心位置的水平投影點，礦體的邊界線由各鄰井的外推點直接相連而成，連接時走向點與走向點相連、傾向點與傾向點相連。工業(yè)井與礦化井之間按勘查工程間距的1/2平推；工業(yè)井與無礦井之間則按勘查工程間距的1/4平推。

(3)礦體規(guī)模與礦體含礦性的關系

分析礦體規(guī)模與礦體含礦性的關系，首先圈定礦體的空間展布形態(tài)并確定礦體規(guī)模(長度和寬度)，然后標注各礦體的含礦性(平均平米鈾量)，最后通過擬合兩者的數值關系即可建立礦體規(guī)模與礦體含礦性的關系(圖5)。統(tǒng)計、分析QJD地區(qū)上百個礦體的規(guī)模與含礦性關系表明：在儲層非均質性變化不大的情況下，礦體的規(guī)模(長度及寬度)與礦體的平米鈾量呈明顯的正相關，即礦體規(guī)模越大則含礦性越高。

圖5 礦體的平面分布圖中顯示了一個局部區(qū)塊圈定的礦體①～④及其含礦性

(4)求取橫向變程

在步驟(1)～步驟(3)的基礎上，即可求取橫向變程，具體為：

(a)通過步驟(3)建立礦體規(guī)模與礦體最高平米鈾量之間的擬合關系

l=169.9lnx+102.7

(2)

b=7.2x+60

(3)

式中：l為礦體長度；b為礦體寬度；x為礦點的最高平米鈾量。

(b)在反演的目標區(qū)域，利用20口井資料，基于式(2)和式(3)，在已知單礦點x的情況下，求取研究區(qū)各礦體的規(guī)模。

(c)對比步驟(b)的結果與初始GR 模型，反復修正橫向變程。

(d)判斷基于當前變程的GR模型與礦體展布規(guī)律的吻合程度。若吻合程度達到要求，則輸出橫向變程；否則重復步驟(a)～步驟(d)，直至吻合程度達到要求。

基于步驟(a)～步驟(d)，最終確定表征研究區(qū)礦體長度和寬度的橫向變程分別為700、300m。

3 鈾礦體預測精度分析及實際應用

為了評估基于本文提出的變差函數求取方法的鈾礦體反演精度及其實際應用效果，選取面積約為50km2的研究區(qū)展開試驗。研究區(qū)位于中國QJD礦區(qū)(圖6)，該區(qū)實際地震資料與井中合成地震記錄吻合良好(圖7)，其中礦區(qū)北部是近年QJD鈾礦床勘探的重點區(qū)域。隨著勘探程度的提高，礦體逐漸接近邊界，勘探成功率呈下降趨勢(近兩年的成功率不足40%)，迫切需要具有針對性的預測手段提高勘探效率，降低勘探成本。研究區(qū)構造北高南低，層間氧化帶由南向北發(fā)育，目的層姚家組主要為辮狀河三角洲平原相地層，砂體厚度大，儲層物性變化較小，目前共完鉆探井200余口，根據鉆井資料可以較準確地刻畫礦體的展布特征。

圖6 工區(qū)范圍示意圖

3.1 縱向變程對反演結果的影響

根據提出的縱向變程求取方法，求得該區(qū)的縱向變程為4.5m。為了驗證所取縱向變程的合理性，對比、分析了不同縱向變程的反演結果(圖8)，可見：①當縱向變程為15.0m時，反演結果分辨率很低，不能清晰地反映縱向的多個礦體，且礦體邊界較圓滑，與尖滅的礦體形態(tài)不符(圖8a)；②當縱向變程為10.0m時，反演結果縱向分辨率略有提高，礦體邊界逐漸清晰，但對尖滅的礦體形態(tài)的刻畫仍不夠清晰(圖8b)；③當縱向變程為4.5m時，反演結果的縱向分辨率較高，位于中間的主礦體形態(tài)清晰，且清晰地刻畫了其上、下的兩套小礦體，較真實地反映了礦體形態(tài)(圖8c)；④當縱向變程為2.0m時，反演結果的縱向分辨率提高不明顯。因此，縱向變程明顯影響反演結果，基于本文方法求取縱向變程的反演結果很好地刻畫了礦體(包括小礦體)形態(tài)。

圖8 不同縱向變程的反演結果(a)15.0m；(b)10.0m；(c)4.5m；(d)2.0m

3.2 橫向變程對反演結果的影響

為了測試本文求取橫向變程方法的效果，對比、分析了不同求取橫向變程方法的反演結果(圖9)，可見：①基于疊后稀疏脈沖反演平面屬性進行橫向約束方法(圖9b)、本文方法(圖9c)獲得的反演結果均刻畫了Y1旋回礦體，并且預測的礦體垂向厚度與礦體實際厚度(圖9a)基本吻合。②在橫向上，圖9b將旋回頂部的一套完整礦體描述為多個獨立小礦體，既不符合礦體合并原則也不符合礦體展布規(guī)律；圖9c將旋回頂部的礦體描述為同一套礦體，更好地反映了礦體的連續(xù)性，更符合礦體展布規(guī)律。圖10為實際礦體平面展布與本文求取橫向變程方法反演結果的GR均方根屬性平面圖。由圖可見：①兩者均呈現南北向的展布特征；②在紅色虛線圈定區(qū)域的反演結果中，不同規(guī)模礦體的主體空間位置、含礦性(圖10b)與已知井區(qū)的實際礦體空間位置、含礦性(圖10a)對應關系很好。綜上所述，基于本文方法求取橫向變程的反演結果具有很高的橫向分辨率，能準確地刻畫礦體形態(tài)。

圖9 不同求取橫向變程方法的反演結果(a)實際礦體剖面；(b)基于疊后稀疏脈沖反演平面屬性進行橫向約束；(c)本文方法U(鈾含量)、GR的量綱分別為0.001%、API，RT、RS的量綱均為Ω·m；Y1、Y2、Y3分別為姚一段、姚二段、姚三段。圖11同

3.3 反演結果驗證

前文分析了縱向、橫向變程對反演結果的影響，分析結果證明了本文提出的變差函數的有效性。但在將反演結果用于實際勘探之前，必須由驗證井驗證反演結果。為了驗證反演結果，本文采用去除已知井進行驗證的方法。反演時選取e井、f井、g井、h井所在區(qū)域進行測試，測試過程僅采用e井、g井、h井數據進行約束反演得到反演結果(圖11a)?？梢?，預測礦體與實際礦體的對應關系較好，反演結果準確預測了發(fā)育于f井目的層底部的礦體，預測礦體與e井連通，在g井附近尖滅，與實際礦體剖面(圖11b)吻合良好，充分證明了本文反演方法的可靠性。

圖11 反演結果的驗證(a)僅應用e井、g井、h井數據約束得到的反演結果；(b)實際礦體連井剖面

3.4 實際勘探效果分析

前文的精度分析表明，基于本文方法的反演結果，無論在縱向還是橫向都能準確刻畫礦體的空間展布特征，對于礦體的實際勘探具有重要意義。為了驗證本文方法的實際應用效果，通過綜合分析反演結果與控礦因素(氧化帶發(fā)育特征等)發(fā)現，在現有礦體北部的一個部署區(qū)域(圖10a藍色虛線框)內還存在礦體發(fā)育區(qū)，2019年底在該區(qū)部署并實施了12口探井，其中9口為工業(yè)礦井，見礦率為75%，遠高于QJD詳查勘探工業(yè)見礦率的平均水平。實際勘探成果充分證實了本文方法反演結果的可靠性。

圖10 實際礦體平面展布(a)與本文求取橫向變程方法反演結果的GR均方根屬性平面圖(b)

綜上所述，當賦予變差函數明確的地質含義后，反演結果與礦體的實際展布特征吻合更好，準確刻畫了鈾礦體的展布形態(tài)，同時也可以根據反演結果在未知區(qū)預測礦體分布，指導勘探部署。

4 結論

本文提出了在測井數據和地質意義控制下的變差函數求取方法，并用于QJD砂巖型鈾礦體預測的地質統(tǒng)計學反演，取得如下結論：

(1)地質統(tǒng)計學反演基于測井數據隨機建模，相對于常規(guī)確定性反演，繼承了地震資料橫向分辨率高的特點，同時具有更高的縱向分辨率。因此對于厚度小、橫向變化快的砂巖型鈾礦體的預測效果較好，其中變差函數的求取是影響反演結果的重要因素。

(2)基于本文提出的結合測井資料與實際礦體形態(tài)的綜合統(tǒng)計方法求取縱向變程的反演方法，在已有先驗井數據控制下得到了符合地質意義的反演結果，進而最大限度地降低了反演的隨機性，反演結果很好地刻畫了礦體(包括小礦體)的縱向形態(tài)。

(3)變差函數反映了礦體的三維空間展布規(guī)律，通過統(tǒng)計、分析礦體規(guī)模和礦點含礦性，擬合了兩者的關系，提出在已開發(fā)區(qū)礦體成礦模式和礦體展布規(guī)律等先驗信息約束下的橫向變程求取方法，獲得的反演結果與實際礦體分布吻合更好。

(4)利用本文方法得到的反演結果更準確地預測了礦體形態(tài)，參考反演結果部署的勘探井位顯著提高了礦體預測的可靠性和鉆探成功率。