基于序貫高斯模擬法的水利工程壩址區速度建模

李偉, 周懷來*, 劉興業, 王元君, 馬慧蓮

(1.成都理工大學地球物理學院, 成都 610059; 2.成都理工大學地球勘探與信息技術教育部重點實驗室, 成都 610059)

全球對水利工程設施的建設需求的不斷增加,勘探力度逐步加大。中國擁有眾多河流,發展前景良好。在追求發展速度的同時,越來越注重水利等工程項目的質量。水庫作為水利工程中的基礎設施之一,為區域居民的生活和生產提供了重要的水資源保障[1]。在修建水庫之前,壩址的選擇至關重要,必須對影響水庫后期投入使用的各種因素進行必要的分析,也包括地質災害分析工作等[2],在全面符合條件的情況下才能進行規劃設計和施工,以確保水庫的正常投入使用。在壩址的選擇中,分析壩址區的速度和巖性是重要的兩個參考因素[3]。通常,壩址區巖層速度的橫向變化范圍大,僅依靠測井資料的垂向速度信息很難作為依據,需要利用其他條件加以約束[4]。通過建立速度模型或巖性模型可用于分析壩址區地下結構的穩定性,然而這對建模成像精度提出了更高的要求[5]。

水庫堤壩進行選址勘察時,會對速度及地下巖石的性質進行重點勘察,速度分析是必不可少的一步。巖石越堅硬致密,測井上的速度響應相對越高,而在巖石松散或孔隙較大部位,速度呈現異常低值,這對水利工程的修建以及后期的運行形成垮塌的風險隱患。分析不同類型巖石的速度對堤壩類型和選址的影響也是重要的一步[6]。對壩址區建立地下巖石的速度模型,是分析壩址穩定性的一種方法。然而這對建模精度有很高的要求,常規的建模方法很難實現,此問題待解決。為得到更高精度的近地表速度模型,學者們開展了許多的研究。趙玲芝等[7]運用多信息融合的方法建立近地表速度模型,將近地表調查資料與大炮初值層析反演結果進行協克里格插值技術融合,得到更高精度的速度模型并在靜校正及疊前深度偏移中得到了很好的效果。劉靈等[8]提出一種相控速度建模方法,綜合地震速度、測井速度、地震屬性及沉積相等多種信息,通過多屬性融合及神經網絡聚類分析來建立低頻相控速度模型,并用于約束反演結果,極大提高了速度建模的精度及分辨率。

水利工程設施的建設對地下速度建模精度要求嚴格[9],而常規的速度建模思路是利用測井資料得到數據進行插值,多為線性插值方法[10],對于水利壩址區此類復雜地層結構,常規插值方法難以解決問題,此外井數據只包含縱向速度信息,對于井之間的橫向大范圍空白區域缺乏約束條件。綜上所述常規插值方法比較簡單[11],建模效果不佳,需尋找高精度的速度建模方法。為此,有學者提出了利用地震層位約束來進行高精度建模的方法。在地震勘探領域,孟會杰等[12]為解決地震偏移成像精度不高的問題,應用層位約束的網格層析建模方法,通過不斷迭代、反演最終建立地下深度域高精度速度模型,并消除了實際工區中的構造成像假象。后來,有學者在對地質統計學建模的研究中發現,相比于常規線性插值建模,基于地質統計學理論的建模方法在效果上更好,精度更高,能夠充分利用變量的空間變化特征[13]。

為提高建模精度,準確分析水利壩址區內部構造特征,服務于壩址的選取及水庫的建設,鑒于前人專家學者的研究,現提出一種基于地質統計學理論的地震與測井相結合的速度建模方法,對常規的測井插值建模的方法以及地質統計學建模方法進行完善和補充,應用地質統計學序貫高斯模擬隨機建模,以地震層位信息作為橫向約束條件聯合測井資料建立速度模型,為壩址區巖石穩定性提供數據支撐,對選壩提供重要參考。同時,在前人對速度建模方法研究的基礎之上,創新性地將本文提出的基于地質統計學理論的建模方法應用于水利工程領域的壩址選取工作上。

1 方法原理

巖石速度是壩址選取的重要影響因素,因此以速度作為關鍵切入點對水利壩址區建立速度模型,進而評估待選壩址區的穩定性。常規速度建模方法較為簡單,模型精度無法滿足水利工程建設的要求,而使用基于地質統計學理論的地震與測井相結合的方法建立高精度速度模型,旨在解決壩址區速度建模精度低的問題,在高精度建模結果上圈定有利于修建壩址的區域,更加準確且具有說服力,為選壩修壩提供重要依據。

地質統計學作為一門發展于20世紀60年代的數學地質邊緣學科,包含空間相關性分析、克里金插值以及隨機模擬三部分[14]。將研究區域在未知位置的函數值作為一組相關隨機變量,通過隨機變量理論分析來預測未知點的信息。相較于線性的插值方法有獨特的優勢特征,且超越了普通的插值問題[15],這也是利用地質統計學進行速度建模的理論依據。

地質統計學提出區域化變量[16]這一概念,能用空間分布來表達具有結構性的隨機變量。同時,為度量區域化變量的空間變異性,提出了變差函數的概念,它兼顧區域化變量的隨機變化和空間結構性規律[17],表示空間中數據方差隨著兩點之間相對距離的變化而在特定方向上的變異性[18]。假設隨機變量ξ滿足本征假設,則在任一方向α,相對距離為|h|的兩個隨機變量ξ(u)和ξ(u+h)的增量的方差稱為變差函數γ(h),表示為

(1)

式(1)中:u為空間中的任意點;h為兩點之間的距離,在地質統計學中稱為滯后距;E為數學期望。在平穩假設條件下變差函數γ(h)的變化只與滯后距h有關,而γ(h)隨h的變化稱為變差函數曲線,如圖1所示。圖1中c0表示滯后距h很小時,兩點之間具有一定的變異性。a表示變程,指區域化變量在空間中的相關范圍,變程越大說明觀測數據在更大范圍內相關。超出變程范圍后,兩點之間無關;在變程內,變量越近,變異性越小。c為基臺值,代表變量在空間中總變異性的大小。在速度建模過程中,對變差函數參數進行分析是至關重要的一步,其結果直接影響建模效果,所以需要對參數進行細致的分析。

圖1 變差函數曲線Fig.1 Variation function curve

通常實際觀測獲取的樣品數目有限,例如在速度建模過程中能提供的鉆井曲線信息較少,而通過有限觀測值經式(1)計算獲得的變差函數作為實驗變差函數γ*(h),表示為

(2)

式(2)中:N(h)為滯后距為h的點對個數;ξ為隨機變量;ui為空間中的任意點。

以h為橫坐標,γ*(h)為縱坐標即可得到實驗變差函數圖。實際應用中,可供利用的點對數量越多,獲取的實驗變差函數就具有越強的代表性。

實際的地質統計學研究用到的都是理論變差函數,因此需要在獲得實驗變差函數之后進行變差函數擬合以得到理論變差函數的參數[19]。其中球狀模型、指數模型、冪模型及高斯模型等是一些基本理論變差函數。根據實際原始數據特點選取相應的變差函數模型,不同模型的插值結果也會不同。本文研究使用球型變差函數模型,其靈活性大、數據適應性強,模型在原點處呈線性變化,在變程處達到基臺值。圖2為球型、高斯型、指數型這三類常用理論模型的變差函數曲線。

圖2 常用變差函數模型Fig.2 Common variation function model

使用球形模型進行變量空間估計時,需定義球體變差函數的長軸、短軸和中長軸的變程以及塊金值、拱高和基臺值等。同時定義球體的空間走向,作為變差函數的方位參數。本文的速度建模方法將測井得到的速度數據視作地質統計學中的區域化變量,對其進行變差函數分析,研究其數據特征選擇球型變差函數模型,求取最佳的變差函數參數。

確定變差函數重要參數后,即可進行地質統計學序貫高斯模擬。它是隨機模擬的一種方法,以已知信息為約束,使用隨機函數建立多個不同地質模型的方法。在模擬過程中考慮區域化變量隨機性的同時,還考慮了其空間相關性。且要求模擬結果忠實于觀測點的值[20]。假設隨機變量符合高斯分布,結合高斯概率理論與序貫模擬算法對連續分布的區域化變量進行模擬,以像元為單位按照選定的隨機路徑進行,并且將已經模擬過的節點的狀態值也直接加入條件數據中。模擬步驟如圖3所示[21]。序貫高斯模擬方法建立速度模型,相比于常規的建模方法在原理上有明顯優勢,它能充分利用實際觀測數據的空間信息,并對其進行變差函數分析以此建立更為精準的速度模型,為壩區提供高精度的地下速度信息,作為實際選壩的主要參考。

序貫高斯模擬可將地震層位作為約束條件加入模擬結果進行隨機建模。而無層位約束的模擬結果無法反映地下巖層的橫向變化特征,會得到一個相對光滑、連續的模型結果。水利壩址區地下構造因其速度在橫向變化范圍大,若無橫向方向的約束將缺失建模的精度,最終模型結果無說服性。因此文章將地震層位作為約束進行速度建模,提高建模的精度與準確性。

序貫高斯模擬法在多次模擬時會生成多個等概率的模型,其屬于條件模擬方法,因此模擬結果忠于觀測點處的數據,而對于非觀測點區域,模型之間差異較大。可對多次模擬得到的速度模型之間有明顯差異的區域做不確定性分析,選壩時則應避開此類位置。

2 建模步驟

(1)對測井速度曲線進行粗化處理,提高垂向采樣使數據更加連續,同時進行井曲線的標準化,使其更加符合地質規律,便于后續的速度插值與建模。

(2)分析變差函數參數。根據已有數據特征和精度,選取球形模型變差模型進行建模,確定球形模型主、次、垂向變程,以及最優的塊金值、拱高和基臺值等重要參數。

(3)對選取的變差函數參數進行序貫高斯模擬,首先進行無層位約束模擬,再以地震層位作為橫向約束進行模擬,完成不同層位的變差函數分析和分層速度建模。

(4)使用協克里金法將上述有、無層位約束情況下建立的速度模型進行加權[22],增強模型的縱向連續性以及準確性,形成最終的速度模型。并與常規的建模方法進行對比、評價。

(5)將最終得到的模型使用Garden公式計算出密度,得出反射系數,與雷克子波進行褶積得到合成地震記錄,與實際地震記錄相對比,評價速度建模的效果。

3 應用實例

在對以上方法研究的基礎上,對位于四川省蒙頂山某一處的近地表鉆井資料進行實際數據測試,該區域地層速度橫向變化范圍大,難以獲得精確的地下速度分布,因此需在高精度速度建模的結果上進行選址分析。工區范圍內分布11口可用井(W1～W11),井位信息如圖4所示,井位分布比較均勻。首先對工區內井曲線做批量標準化處理,得到的速度直方圖如圖5所示。圖6為井W6～W10 5口井標準化后的結果。此外對井速度數據做“粗化”處理,取速度值的算術平均對測井曲線進行重采樣,增強垂向上的連續性。

圖4 壩區井位信息Fig.4 Well location information of dam area

圖5 井曲線標準化前后直方圖Fig.5 Histogram of well curve before and after standardization

圖6 井曲線(W6～W10)標準化前后結果Fig.6 Well curve (W6～W10 ) results before and after standardization

完成上述兩步驟后進行序貫高斯模擬,若以地震層位作為橫向約束,則需要對每一層變差函數進行分析。表1為有層位無層位約束的變差函數參數的分析結果。本文變差函數模型選擇球形模型,由于水平方向變差函數參數對模擬結果影響很小,因此將其視為各向同性,只對垂向變差函數參數進行分析。

表1 有層位、無層位約束變差函數參數Table 1 Variation function parameters with horizon constraints and without horizon constraint

為分析不同變差函數對模擬結果的影響,下文使用無層位約束的模擬方式并通過設置不同的塊金值n(0、0.15、0.30、0.45)及垂向變程v(5、25、45、65)進行測試,如圖7和圖8所示。

圖7 不同塊金值的模擬結果Fig.7 The simulation results of different nugget values

圖8 不同垂直變程的模擬結果Fig.8 The simulation section results of different vertical ranges

通過分析,塊金值越大,變量的隨機性表現得越強烈,當塊金值為0時,模擬結果平滑,能較好體現出研究區域的速度變化趨勢,隨機性較小;對于垂向變差函數分析結果,不同垂直變程對模擬結果的差別較大。取較小的垂直變程能夠提高垂向分辨率,但對區域整體趨勢刻畫不足;若取較大的垂直變程,能體現出區域整體的變化趨勢,卻缺失垂向分辨率。因此選取最佳垂向變差參數才能達到較好的模擬結果,此處選取25作為最佳垂直變程。

完成變差函數分析后,使用序貫高斯模擬方法得到速度模型。首先圖9(a)展示無地震層位約束的速度模型過井W6～W10的剖面,圖中5條黑色橫線為地震層位,由上至下依次命名為S1～S5(S1、S5分別對應t=0 ms和t=90 ms時刻)。從圖中得出由于沒有層位的約束,在分層處速度橫向變化特征未能突出。圖9(b)為有地震層位約束后的模擬結果,在對每個層位的井數據都做了獨立的變差函數分析后,得到了相比無層位約束更高精度的結果,同時能反映出橫向速度變化特征,但其速度模擬結果在局部層位分界處之間變化劇烈,與常規地質規律認識不符。

圖9 無層位約束與有層位約束模擬結果過井W6～W10剖面Fig.9 Simulation section results without and with horizon constraint cross well W6～W10

將無層位約束的速度模擬結果與有層位約束的模擬結果進行協同克里金加權插值,加權步驟中將有層位約束的結果設置為主變量,無層位約束的模擬結果為次變量,通過分析確定權值后進行插值最終得到圖10所示的速度模型。對模擬結果進行分析,模型既能體現出速度的橫向變化的復雜性,也符合地質認識,且相較于上文結果精度進一步提升,為水庫的壩址選取分析以及安全性評價提供重要的數據支撐。

圖10 協同克里金序貫高斯結果過井W6～W10剖面Fig.10 Co-kriging sequential gaussian section results cross well W6～W10

本文研究中使用的序貫高斯模擬屬于隨機模擬方法,多次模擬能得到多個等概率的模擬結果,但在鉆井位置均能忠實于觀測數據。多次模擬結果的差異體現了建模結果的不確定性。每次模擬后可對速度值差異明顯的區域進行分析和取舍,以此縮小范圍確定最佳的選壩位置。圖11展示使用協克里金序貫高斯再次模擬兩次的結果。

圖11 兩次協克里金序貫高斯模擬結果過井W6～W10剖面Fig.11 Twice co-kriging sequential gaussian simulation section results cross well W6～W10

結合圖10,在剖面上紅色區域圈定出了一處三次模擬結果未發生明顯變化且速度響應較高的區域,說明該位置的模擬結果響應更趨近于真實地下構造。在水庫修建中,基于安全考慮選壩位置在速度高值響應特征的巖層中最為合適,而此分析結果為壩址選取提供了一項重要的參考依據。

下文分析常規速度建模方法效果。圖12(a)為克里金插值法得到的速度模型,其結果得到的是一種相對較平滑的結果,局部的速度異常點未得到突出,只能用作反映全局特征。其次克里金插值是一種確定性方法,每次插值得到的模型相同,無法反映模型之間的差異性以及不確定性。

圖12 克里金插值與井速度散點插值結果過井W6～W10剖面Fig.12 Kriging interpolation and well velocity scatter interpolation section results cross well W6～W10

圖12(b)是井速度散點直接插值為速度體后的結果,圖中的黑色豎線為鉆井提供的速度觀測值。分析后得出插值效果并不理想,速度呈明顯條狀分布,且插值結果不完整,模型出現大片空白區域以及速度值為0的異常位置,不宜作為最終速度模型的參考。

通過以上分析,使用地震層位約束的地質統計學協克里金序貫高斯模擬得到的速度模型精度最佳,能較好地展示壩址區地下巖石的速度信息。圖13為t=63.60 ms時刻的速度模型平面圖,通過分析速度在平面上的分布特征,發現速度高值響應大多分布于壩區靠南區域,圖13中紅色實線為上文展示的過井剖面圖在此平面的投影,藍色框線部分圈定了一處有利修建水庫的壩址區,區域內速度值較高,同時在過井剖面上為高值特征響應,其結果為水利工程壩址區的選址提供了重要的參考依據。

圖13 平面上t=63.60 ms時刻協克里金序貫高斯法速度模型Fig.13 Co-kriging sequential gaussian velocity model on the plane at time t=63.60 ms

以圖10速度模型作為最終結果,利用Garden公式計算密度及反射系數,并與主頻為50 Hz的雷克子波進行褶積得到合成地震記錄,檢驗建模效果,如圖14(a)所示。

圖14 合成地震記錄與實際地震記錄過井W6～W10剖面Fig.14 Synthetic seismograms and the actual seismic record section cross well W6～W10

將計算得到的合成記錄與實際地震記錄圖14(b)進行對比,從中分析得出:合成記錄與實際資料的吻合程度較好,地震層位對應地震剖面上的同相軸部分連續性較好,通過本文的方法建立的速度模型對壩址區的穩定性分析及選址都能夠提供重要的依據,達到了速度建模的目的及實際效果。

本文提出的基于地質統計學理論的高精度速度建模方法圈定預測出的實際區域有利壩址區,已被當地相關水利施工單位認證,在實際現場勘測中,區域內地下巖性普遍致密,速度呈現高值響應。壩址位置選取正確,建設方案已被采取,現正進入施工階段。

4 結論

(1)地層的速度信息是影響壩址區選擇以及后期水利工程安全投入使用的重要因素。因此,對速度進行分析是必要的一步,有必要運用近地表速度模型的精細構建技術來確立精確近地表速度模型。

(2)針對水利工程壩址的選取,本文提出基于地質統計學隨機模擬的速度建模方法,并引入地震層位橫向約束建模過程。根據壩址區近地表速度分布的特點選取合理的變差函數模型是確保模型精度的重要前提,若參數選取不當則得到的速度模型無法用來評估壩址區的穩定性以及后期投入使用。利用測井數據作為硬約束,地震層位作為橫向約束的地質統計學序貫高斯模擬方法建立的速度模型符合地質規律。

(3)多次模擬可以得到多個等概率的速度模型,不同模型之間的差異可以反映建模結果的不確定性,為壩址選擇的風險決策提供依據。將提出的方法應用于實際數據資料上并取得了較好的成果,依據建立的速度模型分析并圈定出了有利于修建水利工程的壩址區,圈定方案被采用,實際正在施工當中,為水庫的選壩的合理性及后期的投入使用提供了重要的數據支撐與保障。