面向復雜構造的油藏數值模擬網格剖分算法

孫鵬曹凱張慧賈志賓

（北京大學地球與空間科學學院北京100871）

面向復雜構造的油藏數值模擬網格剖分算法

孫鵬曹凱張慧賈志賓

（北京大學地球與空間科學學院北京100871）

隨著油氣勘探開發不斷進行，目標油氣田構造越來越復雜，然而當下商業建模軟件對復雜構造的表達均不夠理想，大多選擇簡化處理，誤差較大。基于此，通過研究對Cookie-cutter和Stair-Step算法進行了改善，引入地質界面的無限切平面，通過切割求交等處理，實現了對大部分復雜構造的表達，包括逆斷層、盲斷層、多級Y斷層、多級λ斷層、尖滅等。網格模型單元格保持六面體形態，地質界面以“近階梯”形式表達。一方面誤差較小，構造描述精細；另一方面網格正交性較好，有利于油藏數值模擬計算。此外，算法對斷層面、地層面、不整合面等地質界面采用無差別處理，邏輯簡單，易于實現。為驗證算法可行性，對局部復雜構造數據進行了剖分，網格質量較好。

網格；六面體；油藏數值模擬；復雜構造

Class NumberTP391

1 引言

網格剖分是油藏數值模擬的第一步，網格的質量對于數值模擬至關重要，直接影響模擬結果的可靠性［1～2］。隨著油氣資源勘探開發不斷進行，目標油氣田地質構造越來越復雜，對網格模型的復雜度、精確度要求越來越高。構建一種能保持地質三維特性、模擬精度高、生成速度快的網格，已成為油藏數值模擬最具挑戰性的難題［3～5］。

目前大多數商業油藏模擬軟件對結構化網格和非結構化網格都提供支持。但由于四面體網格具有單元形狀不理想、剛度過高和體積鎖死等缺點［6］，而PEBI、Radial等其他非結構化網格雖然克服了這些缺陷，但其構建、處理、可視化、附屬性都非常困難，故非結構化網格沒有得到廣泛應用［7］。結構化網格仍是市場主流，其中角點網格使用最為廣泛，尤以斯倫貝謝公司的“Pillar Grid”為代表，是油藏數值模擬和地質建模框架的行業標準［8］。

Pillar Gird基本滿足了目前大部分需求［8］，但也有致命缺陷：1）對于傾斜斷層，pillar會發生扭曲，使地質統計學算法產生誤差，進而可能導致流體模擬錯誤；2）對于復雜構造（如出現盲斷層、多級Y型斷層、多級λ斷層等），pillar方法無法剖分，只能進行大量簡化，如垂直化斷層、移除部分斷層等［10～11］。因此，提出一個能完整表達任何復雜構造，并且適用于油藏數值模擬的網格，是十分有必要的。

針對上述問題，本文構建了一種基于復雜構造模型的油藏數值模擬網格，并在北京大學信息地質研究實驗室開發的Creatar3.0三維可視化平臺上進行了實現。該網格可較精確表達逆斷層、盲斷層、多級Y斷層、多級λ斷層、尖滅等絕大部分復雜構造，且方便油藏數值模擬計算，在實際應用中效果較好。

2 網格概述

為了克服Pillar Grid的缺陷，學者提出了幾種新的網格形式，主要有Cookie-cutter和Stair-Step兩種（如圖1）［12］。這兩種網格與pillar方法最大的不同在于，允許格子在k方向上切割斷層，避免了大量扭曲格子的出現，從而保證了模型質量。Stair-Step網格以階梯形式逼近地質界面，各網格為準立方體，有利于油藏數值模擬，但斷層表達不夠精確。Cookie-cutter網格在斷層處采用cookie-cutter算法進行切割，可以精確表達斷層，但切割時會發生退化或者產生多面體，無法保證正交性。Cookie-cutter網格模型進行油藏數值模擬時，一般需要先將斷層面轉化為Stair-Step形式。本文在上述網格基礎上進行了改進，構建網格兼具二者優勢。地層面、斷層面、不整合面等地質界面采用“近階梯”形式表達，既可以較精確表達復雜構造，又維持了六面體特性，正交性較好。

圖1 “Stair-Step”（A）和“Cookie-cutter”（B）模式圖［12］

2.1 網格橫向組織

網格模型橫向采用“XY regular”組織形式。如圖2所示，水平面P為網格俯視投影，X方向線X(i)與Y方向線Y(j)正交，構成二維規則網格。網格交點稱為原始點，記為點集V={v1，v2，…}；網格中心稱為中心點，記為點集C={c1，c2，…}。網格X、Y方向步長保持一定，分別為dx，dy。dx，dy大小由網格分辨率決定，可根據實際需求進行調整。

圖2 網格模型俯視圖

2.2 網格縱向組織

網格模型縱向步長可變，單元格頂點為原始點在地質界面相應切平面的豎直投影，其縱坐標取決于相應切平面的空間位置。具體形式如圖3所示，單元格G(cx|k)頂點P1、P2、P3、P4為直線VP(cx|1)，VP(cx|2)，VP(cx|3)，VP(cx|4)與平面F(cx|t-1)交點，頂點P5、P6、P7、P8為直線VP(cx|1)，VP(cx|2)，VP(cx|3)，VP(cx|4)與平面F(cx|t)交點。

圖3網格模型模式圖

圖3 所示符號定義如下，直線VP(vx)為通過原始點vx鉛垂線，稱為邊柱。直線CP(cx)為通過中心點cx的鉛垂線，稱為中心柱。可以看出，每個中心柱CP(cx)周圍環繞4個邊柱，分別記為VP(cx|1)，VP(cx|2)，VP(cx|3)，VP(cx|4)。本網格劃分方法對地層界面、不整合面、斷層面等地質界面不做區分，采用相同處理方法，各界面統一表示為面集H={H1，H2，…}。曲面Ht-1、Ht為兩個鄰近地質界面，平面F(cx|t-1)、F(cx|t)為曲面Ht-1、Ht在點Ct-1、Ct處的切平面（無限連續面），Ct-1、Ct為中心柱CP(cx)與Ht-1、Ht的交點。

2.3 地質界面表達

按照前述方法構建的網格模型，單元格頂面或底面為地質界面相應位置的切平面，地層面、斷層面、不整合面等地質界面則為多個格子的頂面或底面的集合，也就是說地質界面由一系列切平面以“近階梯”狀拼接而成。如圖4所示，當分辨率足夠大，切平面集合就能夠精確近似地質界面。曲面Ht既可以是地層界面，也可以是斷層面等其他地質界面，只是幾何形態有所差異。單元格側面為互相垂直的豎直平面，經過后期調整，可以充分保證網格正交性。網格模型各單元格均保持了六面體形態，避免了退化及多面體（多于6個面）的產生，降低了后期計算的復雜度。

圖4 不同分辨率下地質界面剖面圖

3 構造表達

剖分算法對地層面、斷層面、不整合面等地質界面不做區分，采用相同方法處理，構建任何復雜構造模型都可以理解為對地質界面的表達。而地質界面在幾何上就是曲面，因此構建復雜構造模型就簡化為對一系列曲面的表達。切平面拼接方式，可以處理絕大多數曲面組合，即可以表達大部分地質構造，包括逆斷層、盲斷層、多級Y斷層、多級λ斷層、尖滅等。

3.1 斷層

斷層是地質空間中最主要的構造，斷層的表達是網格剖分算法的核心環節，其精細程度直接影響網格的質量，進而影響油藏數值模擬的準確度。本文構建網格模型對各種斷層的表達如圖5所示。

相比于傳統剖分方法，對于正斷層的表達，本文剖分算法優勢不明顯。但對于逆斷層、盲斷層，尤其是復雜構造的表達，本方法則具有顯著改進：通過切平面求交的方法，有效解決了逆斷層地層重復、盲斷層不斷透地層等問題；地層、斷層的無差異處理，簡化了復雜構造的表達。從而使得逆斷層、盲斷層的表達與正斷層無異，復雜構造變為簡單構造的簡單疊加。

圖5（b）為低角度逆斷層網格模型剖面，其中H1、H2、H4為相鄰地層，被逆斷層H3錯斷。水平面P為模型底面。為表達方便，假設各地質界面為平面。點d(cx|k)為中心柱CP(cx)與各地質界面的交點，平面F(cx|k)為相應地質界面在點d(cx|k)處的無限切平面，點P(cx|1，k)為中心柱cx左側邊柱與切平面F(cx|k)的交點，點P(cx|2，k)為中心柱cx右側邊柱與F(cx|k)的交點（其中x=1，2，3，k=1，2，3…）。圖中除F(c1|1)、F(c3|2)外，其余切平面均與地質界面重合，沒有畫出。在邊柱與切平面求交過程中，若相鄰交點空間關系與對應切點空間關系不一致，則應根據切點空間關系，對交點進行調整。如點P(c1|1，2)處，交點P′(c1|1，2)在P(c1|1，1)之下，而對應切點d(c1|2)則在d(c1|1)之上，此時應將P′(c1|1，2)向上適當微調，至P(c1|1，1)之上P(c1|1，2)處。此外，由于一個邊柱被兩個中心柱共享，所以同一邊柱會求出兩套交點，若對應交點坐標不一致，如點P(c2|2，2)與點P(c3|1，3)，則需進行同一化處理，統一坐標值。得到單元格頂點P(vx|k)后（其中x=1，2，3，4，k=1，2，3…），連接相應頂點，即可完成網格模型構建。若考慮三維空間，一個中心柱則關聯四個邊柱，按上述方法依次進行處理即可。盲斷層、多級斷層構建方法同上。

3.2 尖滅

尖滅是指地層在橫向空間展布上變薄以至消失的現象，與地層堆積時的地理環境、地殼運動情況有關。油氣盆地先后經歷沉積及后期構造運動，尖滅廣泛存在。準確表達尖滅有助于判斷儲層位置，對于油氣開發具有重要意義。尖滅具體表達如圖6所示，構建方法與上述斷層一致。尖滅處六面體沒有發生退化，而是以一個“薄面”（如圖6加粗線段所示）進行模擬。既保持了六面體特性，又較精確地表達了尖滅。

圖5 斷層網格模型剖面

圖6 尖滅網格模型

4 網格剖分

剖分算法以地質界面一致性處理和切平面拼接為核心思想。實現過程以線面求交為主要技術手段。具體剖分流程如圖7所示，包括構造模型準備，初始參數設置，中心柱求交，求取無限切平面，邊柱求交，同一化處理及網格連接七個步驟。

4.1 構造模型

圖7 網格剖分流程

構造模型是本文剖分方法的基礎，構造模型復雜度及精確度決定了網格模型的質量。本文選用Creatar3.0可視化平臺進行構造建模，該軟件可處理任意復雜構造［12］，為本方法的實現提供了保障。如圖8所示，地層面、斷層面采用不規則三角網進行表達，以下步驟對地質界面進行處理即為對不規則三角網進行處理。

圖8 構造模型

4.2 參數設置

在進行網格剖分之前，需要對坐標系、分辨率等關鍵參數進行初始化。坐標系Z軸方向定為豎直向上，X軸、Y軸方向需要根據實際情況設定。選擇原則為盡量減少無效網格，盡可能體現油藏構造特性。實際應用中，可依據地質構造走向、地質體形態進行設定，如選取構造主應力方向或地質體俯視圖最小外包矩形的兩邊方向作為X軸、Y軸方向。網格分辨率通過步長dx、dy進行描述，dx、dy取值需綜合考慮原始數據分布情況和實際需求，如一般以井間有兩三個以上網格單元為佳。

4.3 中心柱求交

根據X軸、Y軸方向，dx、dy數值，可求出原始點及中心點坐標，進而完成水平方向二維規則網格構建。三維空間中邊柱、中心柱z坐標為0，x、y坐標分別與原始點、中心點坐標相同。獲得中心柱坐標后，可依次求出每個中心柱與各地質界面的交點。將計算結果存放在Column類型對象中，以便后續計算。求交過程中，可先行判斷邊柱與地質界面外包矩形空間關系，倘若二者相交，再詳細計算邊柱與不規則三角網的交點坐標，否則就是沒有交點，不必再具體求交，這樣可減少大量計算。

Column類用于記錄網格剖分過程中的重要信息，其主要屬性如表1所示。Column類幾何含義同圖3中Column，每個Column有一個中心柱及四個邊柱，邊柱坐標存放在Fiber中，中心柱坐標不需要直接記錄，可通過切點x、y坐標表達。每個Column中的切點按空間位置從下到上的順序進行記錄，此外還要記錄屬性信息，即與其相交的地質界面的屬性。Column按I、J進行索引，其中I方向平行于X軸方向，J方向平行于Y軸方向。

表1 Column類主要屬性

4.4 求切平面

切平面記錄在TPlane數組中，通過切點和法向量進行表達，其中切點已由上步求出，接下來求出切平面法向量即可。根據切點與不規則三角網的空間關系，分為三種情況（如圖9所示）：1）切點位于三角形內部；2）切點位于三角形邊上；3）切點位于三角形頂點。情況1），三角形所在平面的法向量即為切平面的法向量，情況2）、3），可取關聯三角形法向量的平均值，作為切平面的法向量。本文采用面積加權平均法求取平均值，即

圖9 切平面法向量計算

4.5 邊柱求交

根據I、J索引遍歷Column，對于其中每一個邊柱，按照從下到上順序與切平面求交，并將交點坐標及其屬性存放在對應的InterPoint棧中。求交過程中，需要對每一個交點進行如下判斷：與前一交點縱坐標進行比較，如若該交點空間位置在前一交點之上，則繼續求下一交點；如若不然，則會發生交叉，需將該交點向上平移小段距離ds至前一交點之上。ds取值要綜合考慮地質體形態及模型分辨率。

相比于直接與地質界面求交，引入切平面具有兩方面優勢。一方面由于切平面無限大，每個單元格都可求出8個交點，確保了網格的六面體形態；另一方面，切平面具有單一屬性（由切點屬性決定），消除了同一單元格邊柱交點屬性不一致而引發的歧義。

4.6 同一化處理

對于每一個邊柱，由于為4個Column共享，交點坐標一共求出4套。然而地質含義及網格連接要求網格具有同一性，即同一邊柱同一屬性點只能有一個，故需要對4套交點進行同一化處理，本文采用均值法。具體而言，依次遍歷每一個邊柱，對每一個邊柱做如下處理：1）找出該邊柱關聯的4個Column，并取出相應的InterPoint棧；2）比較4個棧頂元素，選取屬性相同且空間位置相近的交點，求出其平均值；3）彈出計算過的棧頂元素，并用平均值對其交點坐標進行更新；4）重復1）～3），直至4個棧都為空。對于模型邊緣不足4套交點的，按實際數量處理。

4.7 連接網格

同一化處理后，Column中存放的交點就是網格模型單元格的頂點。按照Column中交點的組織方式依次進行連接，即可獲得單元格，進而完成網格模型構建。連接網格過程中，可同時對單元格I、J、K索引及地層屬性進行賦值。

5 模型舉例

為了驗證剖分算法正確可行，本文選擇局部復雜構造數據，在Creatar3.0可視化平臺上進行了剖分試驗。該區域由多級Y字形斷層限定，包括盲斷層、尖滅等構造現象。剖分結果如圖10所示，質量較好。

圖10 應用效果圖

6 結語

本文在Cookie-cutter和Stair-Step算法基礎上進行了改善，在一定程度上解決了復雜構造建模問題。

該方法通過切平面切割的方式，以近階梯形式表達地質界面，實現了對絕大多數復雜構造的精確表達，包括逆斷層、盲斷層、多級Y斷層、多級λ斷層、尖滅等。所建模型正交性較好，有利于進行油藏數值模擬。此外，算法對斷層、地層、不整合面采用一致性處理，邏輯簡單，易于實現。

但模型對傾角近垂直的地質界面表達不夠理想，變形較大。該模型不能直接用于油藏數值模擬，還需要進行細分、粗化等后期處理。網格質量的優劣尚需通過更多的實際應用，進行檢驗。

［1］Hocker C，Thom J.3-D Grid Types in Geomodeling and Simulation--How the Choice of the Model Container Determines Modeling Results［C］//2009 Annual Convention&Exhibition of the American Association of Petroleum Geologists，Denver，Colorado，2009.

［2］Bennis C，Sassi W，Faure J.One More Step in Gocad Stratigraphic Gris Generation：Taking into Account Faults and Pinchouts［C］//SPE European 3-D Reservoir Modeling Conference，Stavanger，Norway，1996.

［3］Aziz K.Reservoir Simulation Grids：Opportunities and Problems［J］.Journal of Petroleum Technology，1993，45（7）：658-663.

［4］Mallison B，Sword C，Viard T，et al.Unstructured Cut-Cell Grids for Modeling Complex Reservoirs［J］.Spe Journal，2013，19（2）：340-352.

［5］Lasseter T J，Jackson S A.Improving integrated interpretation accuracy and efficiency using a single consistent reservoir model from seismic to simulation［J］.The Leading Edge，2004，23（11）：1118-1121.

［6］徐能雄，武雄，汪小剛，等.基于三維地質建模的復雜構造巖體六面體網格剖分方法［J］.巖土工程學報，2006（8）：957-961.

XU Nengxiong，WU Xiong，WANG Xiaogang，et al.Approach to automatic hexahedron mesh generation for rock-mass with complex structure based on 3D geological modeling［J］.Chinese Journal of Geotechnical Engineering，2006（8）：957-961.

［7］Gringarten E J，Haouesse M A，Arpat G B，et al.Advantages of Using Vertical Stair Step Faults in Reservoir Grids for Flow Simulation［C］//2009 SPE Reservoir Simulation Symposium，Woodlands，Texas，USA，2009.

［8］King M J，Ballin P R，Bennis C，et al.Reservoir modeling：From RESCUE to RESQML［J］.Spe Reservoir Evaluation &Engineering，2012，15（2）：127-138.

［9］Filho Z M，Brazil E V，Sousa M C，et al.Cutaway Applied to Corner Point Models［C］//SIBGRAPI，Ouro Preto，Brazil，2012.

［10］Wu X H，Parashkevov R.Effect of Grid Deviation on Flow Solutions［J］.Spe Journal，2009，14（1）：67-77.

［11］Gringarten E J，Arpat G B，Haouesse M A，et al.New Grids for Robust Reservoir Modeling［C］//2008 SPE Annual Technical Conference and Exhibition，Denver，Colorado，USA，2008.

［12］Jayr S，Gringarten E，Tertois A L，et al.The need for a correct geological modelling support：The advent of the UVT-transform［J］.First Break，2008，26（10）：73-79.

［13］Mao P，Zhhaoliang L，Zhongbo C，et al.3-D Geological Modeling-Concept，Methods and Key Techniques［J］. Acta Geologica Sinica-English Edition，2012，86（4）：1031-1036.

［14］Wu Q，Xu H.An approach to computer modeling and visualization of geological faults in 3D［J］.Computers& Geosciences，2003，29（4）：503-509.

［15］Cherpeau N，Caumon G，Lévy B.Stochastic simulations of fault networks in 3D structural modeling［J］.Comptes Rendus Geoscience，2010，342（9）：687-694.

Partitioning Algorithm of Reservoir Numerical Simulation Grid for the Complex Structure

SUN PengCAO KaiZHANG HuiJIA Zhibin
（School of Earth and Space Sciences，Peking University，Beijing100871）

As the continuous exploration and development of oil and gas resources，the geological structure of target oil and gas field become more and more complicated.However，the commercial modeling softwares are not ideal，most of which simplify the complex structure and then model it，having a large error.Based on the Cookie-cutter and Stair-Step algorithm，a new hexahedron grid was proposed.This approach introduced the infinite tangent plane，and expressed the geological interface with“near Stair-Step”.On one hand，almost all the complex structure can be modeled，including reverse fault，blind fault，multi-Y faults，multi-λ faults，pinch-out and so on.On the other hand，cells preserve hexahedron and have a good orthogonality，which is favorable to calculate.Further more，the geological interfaces of fault，layer and unconformity are treated with the same method，lowing the algorithm complexity and easy to implement.The algorithm was verified with the local data of complex structure.The results are in line with expectation and meet the requirement.

grid，hexahedron，numerical simulation of reservoir，complex structure

TP391

10.3969/j.issn.1672-9722.2017.06.004

2016年12月2日，

2017年1月10日

國家自然科學基金項目“基于語義的多分辨率儲層數據組織與管理”（編號：41472113）資助。

孫鵬，男，博士研究生，研究方向：網格剖分算法。曹凱，男，博士研究生，研究方向：三維地質構造建模和相建模。張慧，女，碩士研究生，研究方向：三維地質建模。賈志賓，男，碩士研究生，研究方向：三維地質建模。