基于角點移位的三維地質斷層快速建模算法

劉少華，呂瑞龍，伍東

1.長江大學地球科學學院，湖北 武漢 430100 2.中國石油集團工程技術研究院有限公司信息中心， 北京 102206

計算機相關技術的發展，使得如今能夠通過三維建模的方式，為地下復雜的地質構造建立對應的實體模型，并提供可視化方法，這不僅能夠大大提高非專業用戶對于地質知識的認識與理解能力[1]，而且還為地質學家等的研究提供了一種新的方式。斷層是在地殼中分布極為廣泛的地質構造類型，斷層活動影響著地層的沉積與剝蝕，控制著凹陷的形成和演化[2]，不僅影響著油氣藏的形成與分布，還對油氣開發、剩余油的分布具有重要的影響，因此構造清晰準確的斷層模型對油氣田勘探開發人員更準確地進行油氣藏評價及開發管理有著重要意義。

對三維地質模型及地質斷層模型的問題，國內外已進行了較為廣泛的研究。在三維地質建模時，主要通過使用鉆孔數據對地層模型進行插值，在插值算法方面，CHEN等[3]對深度變化時剪切波速度的變化進行研究，結合鉆孔數據，使用克里金插值算法構建了蘇州城區的三維地質模型；ARFAOUI等[4]等研究了克里格算子在研究地區重力面預測方面的優勢; GONCALVES等[5]對地質結構隱式建模的機器學習算法進行了探討。目前，國內外對于三維斷層建模的方式主要分為整體法和局部法。整體法是先得到無斷層的三維地質模型，再在此基礎上根據斷層數據生成包含斷層的三維地質模型，適用于斷層斷距較小、斷層數量較少的地質體；局部法是基于分區插值的方法生成斷層模型，適用于斷層數據較多、斷距較大的地質體。三維地質模型的建模方式主要有表面數據模型、體元數據模型、混合結構數據模型，其中體元數據模型主要有三棱柱、四面體和六面體網格模型等[6]。表面數據模型建模所需數據量小，因此建模速度較快，但在表達地質體內部屬性方面不如體元模型[7]，針對體元模型的剖切是當前斷層模型構建的熱點研究內容[8]。SADEGH等[9]分析了地震測量地質構造的相關參數，刻畫了研究區的地層及地層的斷裂情況；LIU等[10]采用VC++結合OPENGL，構建了一個地質建模的平臺，通過研究地層的斷裂等對礦區的塌陷進行預測；馬鈞霆等[11]采用剖切廣義三棱柱生成多個四面體的方法，提出了對廣義三棱柱的剖切算法；王海起等[12]采用垂向剖切的方式對六面體網格模型進行剖切生成剖切面，效率較高，然而不能生成斜剖面；張文東等[13]采用了分層投影求交點的方法解決六面體網格求剖切面的問題，需要多次判斷是否相交并標記每一行的剖切六面體元的位置，實際剖切過程中會耗費較多資源；李定林等[14]利用三角網中三角網格的空間拓撲關系提高了獲取三角網格的效率，但是在大數據量網格的情況下，維護三角網格的拓撲關系需要花費較多的時間，影響剖切效率；李昌領等[15]使用切割三棱柱生成多面體，再將多面體轉化為多個四面體，在進行建模時需要進行多種數據結構的轉換，較為復雜；張適[16]基于面元對地層進行建模，對地層模型剖切后會產生空洞，還需對空洞進行縫合。成熟的商業軟件中，Petrel使用最優的基于階梯化網格封堵性分析算法，RMS使用非階梯化網格算法，控制條件繁雜，建模流程繁瑣，適用于復雜斷層建模。筆者提出一種基于角點移位的斷層快速建模算法，基于整體法[17]思路，先不考慮斷層對地層進行建模，計算斷層面剖切的角點網格單元，將這些角點網格的角點按照一定方法移至斷層面上，實現角點自適應斷層模型。算法的優點在于維護了數據結構的完整性和體元個數的不變性，從而提高建模效率及穩定性。

1 建模流程

基于整體法的三維斷層模型構建需要先構造不含斷層的三維地質模型，再根據斷層數據對三維地質模型進行修改，恢復地質體發生斷裂的真實狀態。生成斷層模型算法主要流程如下：

1)使用分層數據進行地層劃分，對目標地層生成六面體角點網格，根據層位點使用反距離加權插值算法計算六面體網格點的高程數據，構建不含斷層的地質模型；

2)融入斷層數據，對斷層數據進行計算處理，生成斷層上盤的頂面和底面、下盤的頂面和底面的斷層線；

3)將三維地質模型沿斷層線劃分為上盤和下盤，對其進行分區插值；

4)根據斷層線計算斷層在三維地質模型中的位置，對斷層線經過的六面體單元角點進行移位，完成對包含斷層的地質體三維模型的構建。

2 地層模型構建

2.1 模型介紹

該研究基于六面體元結構生成三維地層模模型。每個體元由8個角點組成，根據建模范圍、六面體的長、寬設置六面體元8個角點的X、Y方向坐標，高程坐標需要根據鉆孔的層位數據，使用插值算法對每個六面體元的角點進行插值，適用的插值算法主要有克里金算法、反距離加權插值算法等。反距離加權插值算法具有很好的普適性，對斷裂起伏地層的刻畫更為明顯[18]。可以通過設置六面體元的長、寬、層數，來控制建模范圍內的六面體元個數，長、寬越小，層數越多，生成的三維地層模型中單位體積內體元數量越多，構建的模型就越精細。

2.2 模型構建流程

1)首先確定建模范圍、六面體網格中每個六面體元的長(X方向步長)和寬(Y方向步長)，以及建模的層數；

2)確定建模的目標地層；

3)對目標地層構建六面體網格，根據建模范圍及在X方向、Y方向的步長確定每個六面體元角點的X、Y方向坐標。

4)使用鉆孔頂層控制點對當前層的所有六面體網格的頂面角點進行插值計算，使用底層控制點對六面體網格底面角點進行插值計算。

5)遍歷所有目標地層，重復步驟3)和4)，即可完成模型構建。

3 斷層數據及計算

3.1 斷層數據說明

1)控制點。所需控制點為斷層上盤頂面斷層線上的已知點的三維空間坐標，且控制點個數應不少于2個。使用前2個控制點計算斷層線在XOY平面投影線的直線方程，后續的點僅用來作為分區插值時的控制點。

2)斷層傾角。斷層面與六面體網格頂面的二面角。

3)層厚。斷層處地層的厚度。

4)地層斷距。斷層上下兩盤對應位置的垂直距離，上盤上升為正，下盤上升為負。

3.2 計算公式

對于一個地層模型來說，不含斷層時只用考慮地層的頂、底2個面，包含斷層時斷層將地層斷裂為上盤和下盤2個部分，因此斷層模型就需要考慮上盤的頂面和底面、下盤的頂面和底面共4個面，導入的斷層線數據中的控制點為上盤頂面斷層線上的點，因此還需要根據斷層傾角、層厚、斷距來計算控制點在其余各面的對應點坐標。

若A(x1，y1，z1)、B(x2，y2，z2)為斷層上盤頂面斷層線上的2個點，層厚為h1，斷層傾角為θ，斷層斷距為h2，則A點在其他各面對應點(說明：如果A點位于上盤的頂面時，則要計算的是上盤底面、下盤的頂面和底面對應的點)的三維坐標P(x，y，z)計算公式如式(1)～(3)所示：

(1)

(2)

z=z1-h

(3)

式中:h為各面對應點與上盤頂面控制點在垂直方向的距離，當計算上盤底面時，A點與上盤底面對應點在垂直方向上距離為地層厚度，即h=h1；當計算下盤頂面時，A點與下盤頂面對應點在垂直方向上距離為斷距，即h=h2；當計算下盤底面時，A點與下盤底面對應點在垂直方向上距離為層厚加斷距，即h=h1+h2。

4 斷層入網

斷層入網即在三維地質模型的頂面和底面上，斷層線相交于網格線，把斷層線經過的六面體元角點修改為相應的斷層線與網格線的交點，把斷層線作為斷層上盤和下盤的邊界，最后再對網格角點進行插值計算，得出移位后角點的高程值，插值時考慮斷距及層厚數據，即可將斷層的上盤與下盤分裂開來，從而實現模擬地層的斷裂。在不考慮斷層線平行于坐標軸的情況下，將斷層入網問題分為逐行剖切和逐列剖切2種形式。

4.1 逐行剖切

將地層模型和斷層線垂直投影至XOY平面，如圖1(a)所示的斷層線位置，即在每一行上，斷層線切割六面體網格頂面不超過2個六面體元，此時斷層線在XOY平面的投影直線斜率k∈(-∞，-K)∪(K，+∞)(K為六面體元寬與長之比，下同)。

圖1 斷層線入網(逐行剖切)Fig.1 Fault line into the network (line by line division)

以六面體網格頂面為例，將六面體網格投影至XOY平面上得到如圖1(a)所示四邊形網格，逐行剖切方式算法步驟如下：

1)從四邊形網格的第1行(沿Y軸正方向)開始，計算斷層線m與該行網格的上下2個交點P1、P2；

2)如果相交單元格是該行最后一個，如圖1(a)中單元格B，則將單元格B的左角點M1、M2分別移至交點P1、P2，由于相鄰單元格角點不共用，因此需將單元格A的右角點M1、M2分別移至P1、P2；如果相交單元格不是該行最后一個，如圖1(a)中單元格C，則將單元格C的右角點M2、M3分別移至P2、P3，將單元格D的左角點M2、M3分別移至P2、P3，移點效果如圖1(b)所示；

3)重復步驟1)、2)，直至對模型最后一行進行移點；

4)按斷層線將地層模型分區，對網格所有角點進行分區插值，完成斷層建模。

計算六面體元網格中第i行被斷層線剖切的六面體元的列號Ir的計算公式如式(4)所示：

(4)

式中：a、b、c為斷層線方程ax+by+c=0的系數；Δx為六面體元在X軸方向的長度；Δy為六面體元在Y軸方向的長度；i為該行在六面體網格的行號。

4.2 逐列剖切

將地層模型和斷層線投影至XOY平面，如圖2(a)所示斷層線位置，即在每一列上，斷層線切割六面體元不超過2個，此時斷層線在XOY平面的投影直線斜率k∈[-K，0)∪(0，K]。逐列剖切算法思路與逐行剖切相同，區別只是對六面體網格逐列進行分析，具體流程在此不進行贅述。

計算六面體元網格中第j列被斷層線剖切的單元格索引Ic的計算公式如式(5)所示：

(5)

圖2 斷層線入網(逐列剖切)Fig.2 Fault line into the network (column by column division)

圖3 斷層位置的特殊情形Fig.3 Special case of fault location

式中:j為該列在六面體網格的索引。

式(4)和式(5)組成了快速判定斷層線與六面體元網格相交的計算模型，加快了定位相交六面體元的速度，大大提高了斷層入網的速度。

特殊地，針對斷層線位于地質模型一角的情況，為了保證六面體元結構的統一性及六面體元個數的不變性，防止模型超出建模范圍，將如圖3(a)所示頂面為ABCD的六面體元N頂面的左角點A、B進行修改，將A點坐標修改為斷層線與線段AB的交點坐標，由于斷層線與線段DC所在直線交點位于建模范圍之外，因此將D點坐標修改為C點坐標，然后修改六面體元N同行的上一個六面體元M頂面的右側兩點坐標為修改后的A、D兩點坐標。如圖3(b)所示為斷層線位于一角時移點完成后效果圖，此時整個六面體元網格中，六面體元個數保持不變，且六面體元結構仍為8個頂點(六面體元N頂面的C、D點重合)。

5 算法應用

使用所提出的算法對包含斷層的地層進行三維建模，算法應用的計算機環境為：①操作系統：Windows 7 旗艦版 64位；②CPU規格： CPU @2.50GHz；③內存：金士頓 4G。算法應用選擇了A工區12口井的鉆井數據，建模地層有一條斷層穿過。經測試，當網格大小設置為71m×90m×1層時，網格數量有7396個(網格角點數量59168)，建模時間為1650ms。而文獻[13]對于六面體元的剖切，當數據點個數58275時，模型剖切時間為5684ms，在網格數量相等情況下，筆者所提出的算法效率比其高3倍。如圖4所示為對2個地層進行三維斷層建模的前后對比效果圖，其中地質模型建模設置為橫向范圍為6040m，縱向范圍為7680m，網格大小設置為120m×120m×1層。

6 結論

筆者提出了一種基于整體法的三維地質模型可視化表達方法，對含斷層的地層進行快速建模，該算法通過對斷層數據進行處理，判斷斷層面剖切六面體網格的剖切類型，設計了斷層表達數學模型以計算斷層線在斷層上、下盤的頂面和底面上的位置，提出了快速判定斷層線與六面體元網格相交計算模型，從而快速解算出逐行或逐列方向上，斷層線相交六面體元的索引，并對六面體元的角點進行移位，從而使得模型表達更符合地層斷裂的真實狀態。對于單斷層建模，該算法具有以下優勢：

圖4 三維地質斷層模型Fig.4 3D geological fault model

1)保證數據結構穩定。傳統的方法[8,13,15]是六面體直接被斷層面切割，會導致一個六面體被切割為2個不同類型的體，不能保證一分為二后的體都為六面體(8個角點)，其結果不僅會導致六面體元結構不統一，造成體元表達困難，而且增加了體元個數，不便于模型的數據管理和運算處理。該算法在不修改每個六面體元具有8個角點的數據結構條件下，通過對角點坐標進行移位，能夠全面分析并展示平面剖切六面體元的情況，同時不修改體元數據量，為后續模型的分析應用提供了方便。

2)降低時間復雜度，提高建模效率。該算法通過對已經生成的三維模型六面體網格進行處理，計算斷層線與角點網格的交點及交點高程，對相應六面體元角點坐標進行移位，使之自適應斷層面，從而能夠快速地完成斷層模型的生成。