劉少華，呂瑞龍，伍東

1.長(zhǎng)江大學(xué)地球科學(xué)學(xué)院，湖北 武漢 430100 2.中國(guó)石油集團(tuán)工程技術(shù)研究院有限公司信息中心， 北京 102206

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

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

1 建模流程

基于整體法的三維斷層模型構(gòu)建需要先構(gòu)造不含斷層的三維地質(zhì)模型，再根據(jù)斷層數(shù)據(jù)對(duì)三維地質(zhì)模型進(jìn)行修改，恢復(fù)地質(zhì)體發(fā)生斷裂的真實(shí)狀態(tài)。生成斷層模型算法主要流程如下：

1)使用分層數(shù)據(jù)進(jìn)行地層劃分，對(duì)目標(biāo)地層生成六面體角點(diǎn)網(wǎng)格，根據(jù)層位點(diǎn)使用反距離加權(quán)插值算法計(jì)算六面體網(wǎng)格點(diǎn)的高程數(shù)據(jù)，構(gòu)建不含斷層的地質(zhì)模型；

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

3)將三維地質(zhì)模型沿?cái)鄬泳€劃分為上盤和下盤，對(duì)其進(jìn)行分區(qū)插值；

4)根據(jù)斷層線計(jì)算斷層在三維地質(zhì)模型中的位置，對(duì)斷層線經(jīng)過的六面體單元角點(diǎn)進(jìn)行移位，完成對(duì)包含斷層的地質(zhì)體三維模型的構(gòu)建。

2 地層模型構(gòu)建

2.1 模型介紹

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

2.2 模型構(gòu)建流程

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

2)確定建模的目標(biāo)地層；

3)對(duì)目標(biāo)地層構(gòu)建六面體網(wǎng)格，根據(jù)建模范圍及在X方向、Y方向的步長(zhǎng)確定每個(gè)六面體元角點(diǎn)的X、Y方向坐標(biāo)。

4)使用鉆孔頂層控制點(diǎn)對(duì)當(dāng)前層的所有六面體網(wǎng)格的頂面角點(diǎn)進(jìn)行插值計(jì)算，使用底層控制點(diǎn)對(duì)六面體網(wǎng)格底面角點(diǎn)進(jìn)行插值計(jì)算。

5)遍歷所有目標(biāo)地層，重復(fù)步驟3)和4)，即可完成模型構(gòu)建。

3 斷層數(shù)據(jù)及計(jì)算

3.1 斷層數(shù)據(jù)說明

1)控制點(diǎn)。所需控制點(diǎn)為斷層上盤頂面斷層線上的已知點(diǎn)的三維空間坐標(biāo)，且控制點(diǎn)個(gè)數(shù)應(yīng)不少于2個(gè)。使用前2個(gè)控制點(diǎn)計(jì)算斷層線在XOY平面投影線的直線方程，后續(xù)的點(diǎn)僅用來作為分區(qū)插值時(shí)的控制點(diǎn)。

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

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

4)地層斷距。斷層上下兩盤對(duì)應(yīng)位置的垂直距離，上盤上升為正，下盤上升為負(fù)。

3.2 計(jì)算公式

對(duì)于一個(gè)地層模型來說，不含斷層時(shí)只用考慮地層的頂、底2個(gè)面，包含斷層時(shí)斷層將地層斷裂為上盤和下盤2個(gè)部分，因此斷層模型就需要考慮上盤的頂面和底面、下盤的頂面和底面共4個(gè)面，導(dǎo)入的斷層線數(shù)據(jù)中的控制點(diǎn)為上盤頂面斷層線上的點(diǎn)，因此還需要根據(jù)斷層傾角、層厚、斷距來計(jì)算控制點(diǎn)在其余各面的對(duì)應(yīng)點(diǎn)坐標(biāo)。

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

(1)

(2)

z=z1-h

(3)

式中:h為各面對(duì)應(yīng)點(diǎn)與上盤頂面控制點(diǎn)在垂直方向的距離，當(dāng)計(jì)算上盤底面時(shí)，A點(diǎn)與上盤底面對(duì)應(yīng)點(diǎn)在垂直方向上距離為地層厚度，即h=h1；當(dāng)計(jì)算下盤頂面時(shí)，A點(diǎn)與下盤頂面對(duì)應(yīng)點(diǎn)在垂直方向上距離為斷距，即h=h2；當(dāng)計(jì)算下盤底面時(shí)，A點(diǎn)與下盤底面對(duì)應(yīng)點(diǎn)在垂直方向上距離為層厚加斷距，即h=h1+h2。

4 斷層入網(wǎng)

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

4.1 逐行剖切

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

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

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

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

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

3)重復(fù)步驟1)、2)，直至對(duì)模型最后一行進(jìn)行移點(diǎn)；

4)按斷層線將地層模型分區(qū)，對(duì)網(wǎng)格所有角點(diǎn)進(jìn)行分區(qū)插值，完成斷層建模。

計(jì)算六面體元網(wǎng)格中第i行被斷層線剖切的六面體元的列號(hào)Ir的計(jì)算公式如式(4)所示：

(4)

式中：a、b、c為斷層線方程ax+by+c=0的系數(shù)；Δx為六面體元在X軸方向的長(zhǎng)度；Δy為六面體元在Y軸方向的長(zhǎng)度；i為該行在六面體網(wǎng)格的行號(hào)。

4.2 逐列剖切

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

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

(5)

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

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

式中:j為該列在六面體網(wǎng)格的索引。

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

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

5 算法應(yīng)用

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

6 結(jié)論

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

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

1)保證數(shù)據(jù)結(jié)構(gòu)穩(wěn)定。傳統(tǒng)的方法[8,13,15]是六面體直接被斷層面切割，會(huì)導(dǎo)致一個(gè)六面體被切割為2個(gè)不同類型的體，不能保證一分為二后的體都為六面體(8個(gè)角點(diǎn))，其結(jié)果不僅會(huì)導(dǎo)致六面體元結(jié)構(gòu)不統(tǒng)一，造成體元表達(dá)困難，而且增加了體元個(gè)數(shù)，不便于模型的數(shù)據(jù)管理和運(yùn)算處理。該算法在不修改每個(gè)六面體元具有8個(gè)角點(diǎn)的數(shù)據(jù)結(jié)構(gòu)條件下，通過對(duì)角點(diǎn)坐標(biāo)進(jìn)行移位，能夠全面分析并展示平面剖切六面體元的情況，同時(shí)不修改體元數(shù)據(jù)量，為后續(xù)模型的分析應(yīng)用提供了方便。

2)降低時(shí)間復(fù)雜度，提高建模效率。該算法通過對(duì)已經(jīng)生成的三維模型六面體網(wǎng)格進(jìn)行處理，計(jì)算斷層線與角點(diǎn)網(wǎng)格的交點(diǎn)及交點(diǎn)高程，對(duì)相應(yīng)六面體元角點(diǎn)坐標(biāo)進(jìn)行移位，使之自適應(yīng)斷層面，從而能夠快速地完成斷層模型的生成。