李大潛蔡志杰陳娓譚永基王文娟王敬農(nóng)

（1.復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院，上海200433；2.上海立信會計學(xué)院數(shù)學(xué)與信息學(xué)院，上海201620；3.中國石油集團測井有限公司技術(shù)中心，陜西西安710021）

在不同子區(qū)域的交界面（記為γk，k＝1，…，5）上，由電流的連續(xù)性，應(yīng)滿足交界面條件

自然電位測井?dāng)?shù)學(xué)模型與求解方法

李大潛1蔡志杰1陳娓2譚永基1王文娟1王敬農(nóng)3

（1.復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院，上海200433；2.上海立信會計學(xué)院數(shù)學(xué)與信息學(xué)院，上海201620；3.中國石油集團測井有限公司技術(shù)中心，陜西西安710021）

編者按：李大潛先生，1937年生，復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院教授，中國科學(xué)院、發(fā)展中國家科學(xué)院及歐洲科學(xué)院院士，法國科學(xué)院及葡萄牙科學(xué)院外籍院士。曾為電阻率測井方法建立了統(tǒng)一的基本理論框架，據(jù)此制作的微球形聚焦測井儀器，20多年來一直成功地在大慶等多個油田使用。曾獲國家自然科學(xué)獎二等獎、何梁何利基金科學(xué)與技術(shù)進步獎、華羅庚數(shù)學(xué)獎、上海市科技功臣獎、蘇步青應(yīng)用數(shù)學(xué)獎等科技獎勵。謹(jǐn)以此文鼓勵測井界廣大同仁重視基礎(chǔ)理論和基礎(chǔ)方法研究。感謝李大潛院士對測井理論的突出貢獻！

對石油勘探開發(fā)中常用的自然電位測井方法，在以往散見于各處的研究成果及近期研究的基礎(chǔ)上，系統(tǒng)地總結(jié)了其相應(yīng)的數(shù)學(xué)模型，并補充了眾多算例，為自然電位測井在軸對稱情況下建立了完整的理論框架。有關(guān)內(nèi)容包括自然電位測井?dāng)?shù)學(xué)模型及其適定性，解的極限狀態(tài)及相應(yīng)的簡化模型，有效求解的數(shù)值方法，并揭示了有關(guān)自然電位測井方法的算例及應(yīng)用。

自然電位測井；軸對稱情形；數(shù)學(xué)模型；數(shù)值方法

0 引 言

在石油勘探開發(fā)中，自然電位測井是最常見也是最重要的測井方法之一。由于電化學(xué)等因素的作用，地層中正負(fù)離子具有不同的遷移速度，且泥巖顆粒經(jīng)常吸收正離子，因此在不同地層的交界面上會產(chǎn)生穩(wěn)定的電位差，稱為自然電位差，從而在地層中形成一個電場，為與人為供電產(chǎn)生的人工電場相區(qū)別，稱其為自然電場。自然電場的分布和巖性有密切關(guān)系，特別是在砂／泥巖剖面中能以明顯的曲線異常變化顯示出滲透性地層。因此，通過研究井眼內(nèi)自然電場中的電位變化即可反映井壁巖性，這種測井方法稱為自然電位測井。自然電位測井具有方法簡單、實用價值高等特點，是劃分巖性和研究儲集層性質(zhì)、求取測井參數(shù)以及其他地質(zhì)應(yīng)用中不可缺少的基本方法之一［1－3］。

在進行自然電位測井時，得到的自然電位曲線（SP曲線）隨著地層的改變而變化，可以指示滲透層的位置所在［4－8］。

1 自然電位測井的數(shù)學(xué)模型

式中，ρ為地層電阻率，在各個子區(qū)域Ωi中，ρ＝ρi均為常數(shù)（i＝1，2，3，4）。記σ＝1／ρ為地層電導(dǎo)率，在各子區(qū)域Ωi中，σ＝σi＝1／ρi（i＝1，2，3，4）。本文總假設(shè)σi＞0（i＝1，2，3，4）。

圖1 地層剖面示意圖

在軸對稱假設(shè)下，方程（1）可改寫為軸對稱形式

式中，n為Γ2的單位外法向量，而?u／?n為u的外法向?qū)?shù)。

注1 在實際應(yīng)用中，電極被套在一有相當(dāng)長度的絕緣套管上，此時Γ27應(yīng)表示此絕緣套管之表面，而以下整個的討論除特別說明外，均不受影響。

對于測量電極，它主要用來測量電極表面的電位值，而不發(fā)射電流。由于電極是良導(dǎo)體，其表面上的電位值應(yīng)處處相等，為一個（待定）常數(shù)，于是在電極表面（記為Γ0）上應(yīng)滿足等位面邊界條件［9］

在不同子區(qū)域的交界面（記為γk，k＝1，…，5）上，由電流的連續(xù)性，應(yīng)滿足交界面條件

由于不同介質(zhì)的交界面上存在自然電位差，因而

其中，Ek為交界面γk（k＝1，…，5）上的自然電位差，在應(yīng)用中常設(shè)為常數(shù)，“＋”和“－”分別表示函數(shù)在交界面兩側(cè)的值，而在交界面上的單位法向量n規(guī)定對兩側(cè)有同一取向。

由與交界面條件（8）的相容性，電位函數(shù)u在地表面上Γ11部分應(yīng)取常數(shù)

作變換

其中，u0為分塊常數(shù)函數(shù)

其中，

即在垂直交界面上不再存在自然電位差。因此，在Γ1＝Γ11∪Γ12上，邊界條件（3）和（9）此時可統(tǒng)一地寫為

定理1.1 對任意給定的幾何結(jié)構(gòu)及地層電阻率，電極上、井軸上及區(qū)域Ω1和Ω2中的自然電位函數(shù)僅依賴于常數(shù)E5、E1＋E3和E1＋E2＋E4，而不是依賴于獨立的5個常數(shù)Ek（k＝1，…，5）。

這一定理將會極大地減少制作相應(yīng)測井解釋圖版的工作量。

這樣，通過變換（10），就得到自然電位測井滿足的數(shù)學(xué)模型：擬調(diào)和方程（2）滿足相應(yīng)的邊界條件（14）、（4）～（6）及交界面條件（7）、（12）的邊值問題。為后面敘述方便，將這一模型記為（SP）。

我們可以考慮（而且在下面的討論中也要用到）更為一般的情形，即交界面上的自然電位差Ek（k＝1，…，5）可以不是常數(shù)，而是函數(shù)，記為Ek（s），其中，s是γk上分別以A（k＝1，3，5）和B（k＝2，4）為起點計算的長度，并假設(shè)Ek（s）是H?lder連續(xù)的。為保證交界面條件（8）與邊界條件（4）的相容性，總假設(shè)

而為保證交界面條件（8）與邊界條件（3）的相容性，邊界條件（9）應(yīng)取為

同樣可將邊界條件（3）和（16）統(tǒng)一地寫為（14），且變換后相應(yīng)得

在今后對自然電位測井的求解中，只會用到在A、B附近Ek為一般函數(shù)，而其余處Ek仍為常數(shù)的情形，此時式（15）自然滿足。

2 自然電位測井?dāng)?shù)學(xué)模型的解的性質(zhì)

2.1 解空間

為求解自然電位測井?dāng)?shù)學(xué)模型（SP），考慮等價的變分問題，為此先引進幾個有關(guān)的函數(shù)集合。

對任意給定的p（1＜p＜∞），引入集合

其范數(shù)定義為

這樣，與模型（SP）等價的變分問題為：求u∈Vp，使得其中，

如果交界面上的自然電位差Ek（s）（k＝1，…，5）在交匯點A和B處的代數(shù)和分別為0，即滿足

則稱交界面上的自然電位差滿足相容性條件［12］。此時，可在分塊空間中求解上述問題（SP），并可用有限元素法求其數(shù)值解。

然而，通常交界面上的自然電位差并不滿足相容性條件（24），此時整個電場中的電位函數(shù)具有一定的奇性，不再屬于分塊空間，從而不僅在理論分析上有本質(zhì)困難，而且在數(shù)值計算時，不能直接用有限元素法求解。

下面先考慮解應(yīng)隸屬的空間［8，12］。從現(xiàn)在起，有關(guān)的引理及定理均只列出結(jié)論，而略去其證明。

此處及以后，C（p）均表示一個僅依賴于p的正常數(shù)。

引理2.2 對任意給定的p（1＜p＜2），設(shè)Ek（s）∈Cα（γk）（k＝1，…，5），其中α＞1－，則Vp非空，且存在w∈Vp，使得

2.2 滿足相容性條件時的分塊解的存在性

設(shè)交界面上的自然電位差滿足相容性條件（24），并記u是變分問題（22）在V2中的解。由引理2.1，V2非空，并任意取定w∈V2滿足式（25）。記ν＝u－w∈，則對任意給定的φ∈，成立

即

其中，在Ωi中，σ＝σi，f＝（f1，f2）＝σi▽w（i＝1，2，3，4）。顯然，f∈［（Ω）］2。

由Lax－Milgram定理［13－14］，變分問題（28）在中存在唯一解ν∈，從而變分問題（22）存在唯一解u＝ν＋w∈V2。

如果相容性條件（24）不滿足，由引理2.1，V2為空集，因而變分問題（22）不存在分塊的解。而由引理2.2，對任意給定的p（1＜p＜2），Vp非空，從而可以考慮變分問題（22）的分塊解的存在性。

2.3 不滿足相容性條件時的分塊解的存在唯

一性及估計［6，8］

定理2.1 存在ε0＞0，對于任意給定的p（2－ε0＜p＜2），問題（SP）存在唯一的解u∈Vp，且成立

3 極限性態(tài)

圖2 橢圓型邊值問題

圖3 等位面邊值問題

及通常的邊值問題（見圖2）

其中，

相應(yīng)地，u∈W是下述變分問題的解：對任意給定的φ∈W0，成立

其中，

定理3.1 假設(shè)Ek（s）（k＝1，…，5）滿足相容性條件（24）及E5（C）＝0。lε（ν）及l(fā)（ν）分別是ν在空間（Ωε）及W0上的線性連續(xù)泛函，則問題（30）和（31）存在唯一的分塊解uε和u，且當(dāng)ε→0（其中ε為電極的高度）時，

在實際應(yīng)用中，交界面上的自然電位差Ek（k＝1，…，5）往往不滿足相容性條件（24），由定理2.1，此時對任意給定的p（2－ε0＜p＜2），問題（30）和（31）存在唯一的分塊解，而相應(yīng)的等位面邊界Γ0趨于原點時的極限性態(tài)則由如下定理給出。

定理3.2 當(dāng)相容性條件（24）不滿足時，存在ε0＞0，對任意給定的p（2－ε0＜p＜2），當(dāng)電極的高度ε→0時

下面用數(shù)值模擬驗證上述收斂性的結(jié)論。由于相容性條件（24）一般不滿足，解存在奇性，不能直接使用有限元素法求解，需要采用一些方法去掉這一奇性。這里采用過渡帶法（將在第4節(jié)中詳細(xì)介紹），這是最為簡便有效的方法。

取定以下基本幾何參數(shù)：目的層的半層厚H＝2m，井筒半徑R0＝0.125m，侵入帶半徑Rxo＝0.65 m，所考慮區(qū)域的徑向深度R＝1 000m，縱向高度Z＝600m。選取4組不同的物理參數(shù)分別進行計算：

（1）ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4，E1＝10 mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV；

（2）ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100，E1＝10 mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV；

（3）ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4，E1＝20mV，E2＝0mV，E3＝130mV，E4＝130mV，E5＝10mV；

（4）ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100，E1＝20 mV，E2＝0mV，E3＝130mV，E4＝130mV，E5＝10mV。

固定測量電極的直徑D0＝0.03m，將對應(yīng)于不同電極高度ε（當(dāng)ε→0時，原電極面退化為直徑為D0的圓周）的相應(yīng)問題在（D0／2，0）處的自然電位值（其中k＝1，2，3，4對應(yīng)于上述4組不同物理參數(shù)取值情況，下同），并列出ε＝0（即無電極時）在（D0／2，0）處的自然電位值以及在原點處相應(yīng)的自然電位值，計算結(jié)果見表1。

表1 電極高度的影響

現(xiàn)在同比例改變測量電極的高度ε和電極直徑D的大小，將相應(yīng)問題（在ε→0及D→0時，電極退化為原點）在（D／2，0）處的自然電位值記為（k＝1，2，3，4），并列出無電極時在原點處的自然電位值（k＝1，2，3，4），計算結(jié)果見表2。

4 求解方法

4.1 相容性條件（24）成立時的有限元素法

先考慮相容性條件（24）成立時，變分問題（22）的求解方法。由本文第2節(jié)，此時變分問題（22）存在唯一的分塊解u∈V2。

不失一般性，對區(qū)域Ω作三角形剖分（為明確起見，這里我們假定對區(qū)域作三角形剖分，事實上，以下方法對其他形狀的剖分也同樣適用），則變分形式（22）可改寫為

其中，e為任一三角形單元，且σe＝σi，若e?Ωi。

對任意給定的三角形單元e，記其3個頂點按逆時針順序排列分別為（rp，zp）（p＝1，2，3），函數(shù)u在各頂點上的值為up（p＝1，2，3）。對函數(shù)u作線性插值

其中，

為形函數(shù)，其系數(shù)ap，bp，cp（p＝1，2，3）分別為

而

為單元e的面積。由此可以得到元素剛度陣

將元素剛度陣Ae組裝成總剛度陣A，它是一個對稱矩陣。

由于在交界面γk（k＝1，…，5）上存在自然電位差，需對其作相應(yīng)的處理。為明確起見，規(guī)定僅保留上的電位函數(shù)值，而γ－k上的電位函數(shù)值需利用交界面條件（8）轉(zhuǎn)化為用γ＋k上的電位函數(shù)值來表示的形式。

式（50）右端的前一項是對元素剛度陣的貢獻，而后一項是對右端的貢獻，可寫為

其中，U＝（u1，…，uN）T為由電位函數(shù)在各節(jié)點的值構(gòu)成的列向量（對于交界面γ上的節(jié)點，表示γ＋上的電位函數(shù)值），而N為節(jié)點總數(shù)。

對Dirichlet邊界條件u｜Γ1＝0的一種處理方法是對系數(shù)矩陣A和右端項b作劃行劃列處理［9，16］：對Γ1上的任意一個節(jié)點k，將系數(shù)矩陣A的第k行和第k列劃去，同時將右端項b的第k行（即第k個分量）劃去。然而，這在計算機處理時需要對數(shù)據(jù)作大量的移動，會耗費較多的時間。

為了避免作這種數(shù)據(jù)移動，可采用另一種處理方法：對Γ1上的任意一個節(jié)點k，將系數(shù)矩陣A＝（apq）的第k行和第k列分別賦值為初始單位向量εk（即第k個分量為1，其余分量均為0的向量），即

同時，將右端項b的第k個分量賦值為0∶bk＝0。

對等位面邊界條件u｜Γ0＝常數(shù)作處理：固定Γ0上的一點，設(shè)其為節(jié)點k，對該等位面上的任意其他一個節(jié)點p，將系數(shù)矩陣A的第p行加到第k行，第p列加到第k列，然后再將第p行和第p列分別賦值為初始單位向量εp，即

這樣處理后仍保持矩陣A的對稱性。

然后，將右端項b的第p個分量加到第k個分量，再令第p個分量為0，即

經(jīng)過對邊界條件的處理，得到線性代數(shù)方程組

4.2 相容性條件（24）不成立時的求解方法

如果相容性條件（24）不滿足，由引理2.1，V2為空集，因而變分問題（22）不存在分塊解。而由定理2.1，存在ε0＞0，對任意滿足2－ε0＜p＜2的p值，存在唯一的分塊解u∈Vp。

由于奇性的存在，不能直接使用有限元素法求解（22），需要采用一些方法去掉這一奇性。有3種可供使用的方法：過渡帶法、移去奇性法和雙層位勢法。

4.2.1 過渡帶法

過渡帶法，又稱為交界面條件相容化方法［5］。這是最為簡便及有效的方法，特別值得推薦。對充分小的ε＞0，定義

其中，rA、rB分別為A、B的r坐標(biāo)，

容易看出，這里只是將原有的交界面上的自然電位差Ek（s）（k＝1，…，5）在交匯點A及B的附近適當(dāng)變化，使其滿足相容性條件（24），從而由本文第2節(jié)，對應(yīng)于（s）（k＝1，…，5）的問題存在唯一的分塊解uε，且可用有限元素法進行求解。當(dāng)ε充分小時，可以將uε作為原問題（22）的近似解。這一方法的合理性可由定理2.1保證，這就是定理4.1［8］。

定理4.1 存在ε0＞0，使對任意給定的p（2－ε0＜p＜2），過渡帶法的解uε當(dāng)ε→0時，在中強收斂于原問題的解u。

4.2.2 移去奇性法

在相容性條件（24）不滿足的情形下，解決奇性的另一種方法是移去奇性法［6］，其基本思想：構(gòu)造一個函數(shù)w，使得w在交界面上的跳躍值在A點和B點恰好滿足

然后記ν＝u－w，將變分問題（22）改寫為關(guān)于ν的形式，而ν滿足相容性條件（24）。

下面介紹w的構(gòu)造方法。

令

其中，

容易驗證

它滿足式（59）。

類似地，適當(dāng)選取

可使wB滿足式（60）。

記ν＝u－wA－wB，將變分問題（22）改為ν的方程，而ν滿足相容性條件（24），可用有限元素法求解。最后，由u＝ν＋wA＋wB得到原變分問題（22）的解。

4.2.3 雙層位勢法

第3種處理奇性的方法是雙層位勢法［17－18］。

在實際應(yīng)用中，常設(shè)交界面上的自然電位差Ek（k＝1，…，5）為常數(shù)。由前述，利用變換（10），可將交界面條件變?yōu)?/p>

回到三維空間的地層模型來考慮問題。若介質(zhì)在整個地層中是均勻的，由雙層位勢的定義［19－20］，具單位偶極矩面密度的偶極子雙層面Σ對空間中任一點M所產(chǎn)生的電位為

其中，rPM是點M到Σ上某一點P的距離；dSP為雙層面Σ上的面積微元；n為負(fù)電荷所在面Σ—指向正電荷所在面Σ＋的單位法向量（見圖4）。

圖4 雙層位勢法

雙層位勢（71）具有如下的性質(zhì)［19－20］：

引理4.1 對任意給定的M?Σ，雙層位勢（71）在M點處是調(diào)和的，即

引理4.2 雙層位勢（71）在Σ上具有第1類間斷，且對任意給定的P0∈Σ，有

引理4.3 對任意給定的P0∈Σ，雙層位勢（71）的外法向?qū)?shù)在P0點是連續(xù)的，即

作極坐標(biāo)變換

易見有

定理4.2 設(shè)

則有

且有估計式

其中，R＞a2，Z＞z0。

現(xiàn)在利用前述雙層位勢的性質(zhì)，尋求問題（22）的數(shù)值解。在這里能夠得到簡便處理的，是應(yīng)用中在交界面上具常數(shù)自然電位差的情形，而且只針對測量電極尺寸較小、可以近似為對稱軸上一點的情形（見本文第3節(jié)）。此時的計算框架如下：

第1步 利用雙層位勢的性質(zhì)，找到滿足方程（2）及交界面跳躍條件w＋－w－＝Fk（k＝1，…，5）的函數(shù)w，其中F1＝F2＝F5＝0，而F3、F4由條件（70）式確定。

第2步 作平移ν＝?u－w，此時ν仍滿足原方程（2），但在交界面上連續(xù)，于是可利用相應(yīng)的變分問題，藉助于有限元素法求解。

第3步 令u＝u0＋w＋ν，即可求得原問題（22）的解。

根據(jù)上述計算框架，可給出具體的實施細(xì)節(jié)及有關(guān)公式［17，18］。

5 數(shù)值模擬

5.1 均勻地層

設(shè)整個地層中介質(zhì)是均勻的，即ρi＝ρj（i≠j）；又設(shè)自然電位差E4和E5一直延伸到無窮遠(yuǎn)處。此時，井軸上自然電位的精確解表達(dá)式［4，7］為

取H＝2m，R0＝0.125m，Rxo＝0.65m，E1＝10mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV，R＝1 000m，Z＝600m。將電極A0近似為井軸上的原點，分別利用過渡帶法和雙層位勢法在近20 000個節(jié)點的網(wǎng)格剖分下做數(shù)值計算，同時利用式（82）計算井軸上的自然電位精確值。比較發(fā)現(xiàn)它們在井軸上的值誤差很?。ㄒ妶D5）。再取E1＝20mV，E2＝50mV，E3＝130mV，E4＝140mV，E5＝30mV，其余參數(shù)不變，它們的結(jié)果亦非常接近（見圖6）。

圖5 均勻地層情形1

圖6 均勻地層情形2

5.2 非均勻地層

自然電位測井的主要目的之一是判定地層，也就是識別各地層的交界面。取H＝2m，R0＝0.125 m，Rxo＝0.65m，E1＝10mV，E2＝20mV，E3＝100 mV，E4＝120mV，E5＝30mV，R＝1 000m，Z＝600m?？紤]下面2種非均勻情形：

ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶2∶2

ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4

用雙層位勢法（設(shè)電極已按本文第3節(jié)中所述收縮為對稱軸上的1個點）繪出井軸上的SP測井曲線（見圖7及圖8中的點虛線）；同時用過渡帶法計算井軸上的SP測井曲線（見圖7及圖8中的實線）。從圖7及圖8中清楚地看到2種方法的計算結(jié)果很接近。

圖7 非均勻地層ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶2∶2

圖8 非均勻地層ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4

5.3 在交界面交匯點A附近自然電位及電場

固定物理參數(shù)ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30 mV，選取2組不同的E1值分別進行計算，繪制A點附近的等位線圖和電場圖。

（1）E1＝10mV；（2）E1＝200mV。

對此2種情形，均成立ΔB?E4－E3－E2＝0，即在交界面交匯點B處滿足相容性條件，而ΔA?E3－E1－E5≠0，故重點考慮在A點附近的情況。

圖9及圖11顯示2種情形下在A點附近的自然電位情況。由于交界面上電位有跳躍，我們在A點附近的3塊區(qū)域（Ω1、Ω2、Ω3）上分片繪制等位線，圖形顯示A點附近電位的變化比較劇烈。

圖9 情形1時的等位線圖

圖10 情形1時的電場分布圖

圖10及圖12顯示2種情形下在A點附近的電場情況。由于各點電場強度之值不同，且在離交界面近的地方電場強度值顯著趨大，為了觀察電場的走向，圖10及圖12中各點的電場矢量均做了單位化處理。

我們觀察到在2種情形下A點附近的電場旋轉(zhuǎn)方向正好相反（情形1為逆時針，情形2為順時針）。注意到情形1下ΔAE3－E1－E5＝60＞0；而情形2下ΔA＝－130＜0。二者不同的符號決定了2種不同的電場走向：當(dāng)ΔA為正時，A點附近的電場按逆時針旋轉(zhuǎn)；而當(dāng)ΔA為負(fù)時，則按順時針旋轉(zhuǎn)。從其他類似的算例也可以看到這一點。

5.4 ρ1很大時套管壁上自然電位曲線

取定ρ1∶ρ2∶ρ3∶ρ4＝5000∶2∶3∶4，其中電阻率ρ1相對取值很大，可以近似地將井筒區(qū)域Ω1視為絕緣。分別取下面2組不同的自然電位差參數(shù)：

（1）E1＝10mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV；

圖11 情形2時的等位線圖

圖12 情形2時的電場分布圖

（2）E1＝20mV，E2＝0mV，E3＝130mV，E4＝130mV，E5＝10mV。

它們均滿足B點的相容性條件ΔB?E4－E3－E2＝0，而在A點ΔA?E3－E1－E5≠0。

圖13及圖14顯示對這2種情形套管壁（r＝D／2）上的自然電位曲線。從圖中可以觀察到在該曲線拐點處（它相應(yīng)于目的層與圍巖的分界處）自然電位值的跨度與相應(yīng)的ΔA值一致：在情形1時，ΔA＝60mV，而在曲線拐點附近自然電位的最大值為90mV，最小值為30mV；而在情形2時，ΔA＝100 mV，而曲線拐點附近自然電位的最大值為110 mV，最小值為10mV。這一事實與物理上的直觀結(jié)論吻合。

圖13 套管壁自然電位曲線（情形1）

圖14 套管壁自然電位曲線（情形2）

5.5 參數(shù)對自然電位曲線的影響

測井儀器工作時，其測井響應(yīng)不可能正好反映目的層的電參數(shù)信息。實際上，測井響應(yīng)將在不同程度上受到井筒、侵入帶和上下圍巖的影響。因此，測井的電阻率響應(yīng)稱為視電阻率，它與目的層的真電阻率是有差別的。為了讓測井視電阻率響應(yīng)盡可能地接近目的層真電阻率，要求對測井響應(yīng)作出校正和正確的反演、解釋，這就需要對大量的地層模型測井響應(yīng)數(shù)據(jù)進行比較研究。

影響自然電位測井響應(yīng)的因素多種多樣，可分別對不同的地層厚度、井筒半徑、侵入帶半徑、圍巖電阻率、侵入帶電阻率及目的層電阻率進行相應(yīng)的數(shù)值模擬，評估各因素對測井響應(yīng)的影響。

例如，可取自然電位差參數(shù)為E1＝20mV，E2＝0mV，E3＝130mV，E4＝130mV，E5＝10mV，所考察的地層尺寸為R＝1 000m，Z＝600m。此外，設(shè)定H＝2m，R0＝0.125m，Rxo＝0.65m，ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100為基本參考參數(shù)值，可通過數(shù)值試驗每次只考慮測井響應(yīng)隨其中一種參數(shù)變化的情況，其余的參數(shù)則取這組數(shù)據(jù)中相應(yīng)的值。具體計算結(jié)果從略。

6 將自然電位測井儀安裝在側(cè)向測井儀上對自然電位測井的影響

在實際測井中，常將自然電位測井與側(cè)向測井聯(lián)合使用，并將自然電位測井儀安裝在側(cè)向測井儀的上方（見圖15）。其中A為自然電位測量電極環(huán)，直徑為30mm，長度為250mm；B為一段絕緣表面的電纜，可視為絕緣體，直徑為30mm，長度為4.5 m；C為導(dǎo)體，直徑為90mm，長度有3種情況：5、10、15m；D為絕緣體，直徑為90mm，長度為800 mm；E為雙側(cè)向儀器＋聲波儀器，直徑為90mm，長度信息見圖16。

圖15 自然電位測井儀安裝在側(cè)向測井儀的上方

圖16 側(cè)向測井儀器數(shù)據(jù)

在進行側(cè)向測井時，自然電位測井儀不工作；而在進行自然電位測井時，側(cè)向測井儀不工作。由于側(cè)向測井儀由眾多大小不同的電極組成（包括1個發(fā)射電極，若干屏蔽電極及測量電極），而且為隔開這2類測井儀器，在其間還插有圓柱形的導(dǎo)體C（相當(dāng)于一個不發(fā)射電流的電極）及圓柱形絕緣電纜B及D。在進行自然電位測井時，盡管這些電極處于不工作的狀態(tài)，即它們不發(fā)射電流，但由于它們及絕緣體的存在，原先自然電位場所在區(qū)域被它們在井軸附近占據(jù)了一部分，原先的一部分本質(zhì)上滿足絕緣條件的井軸，要代之以這些電極Γk（k＝1，…，l）表面上的齊次等位面邊界條件

這樣，原先求解自然電位的等位面邊值問題（1）～（9）要代之以等位面邊值問題（1）～（9）及（83）～（85），它們具有不同的求解區(qū)域及一些不同的邊界條件，整個自然電位場要相應(yīng)地發(fā)生變化。特別地，在測量電極Γ0上的自然電位值U（這是自然電位測井主要關(guān)心的對象）要相應(yīng)地發(fā)生變化。一個迫切需要考慮的問題是：這樣做對測量電極上的自然電位值U會不會帶來重大的影響。

為了檢驗側(cè)向測井儀器的電極串對自然電位測井的影響，針對圖16所示的數(shù)據(jù)，我們?nèi)“雽雍馠＝2m，井筒半徑R0＝0.125m，侵入帶半徑Rx0＝0.65m，并對所考慮區(qū)域取徑向深度R＝1 000m及縱向高度Z＝600m，通過過渡帶法進行數(shù)值計算。記未加入側(cè)向儀器電極串時測量電極Γ0上的自然電位值為U，加入側(cè)向測井儀器電極串時相應(yīng)的電位值為?U，取下面的幾組數(shù)據(jù)進行計算，計算結(jié)果見表3。

（1）ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4，E1＝10mV，E2＝20mV，E3＝100mV，E4＝120mV， E5＝30mV；

（2）ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100，E1＝10mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV；

（3）ρ1∶ρ2∶ρ3∶ρ4＝1∶2∶3∶4，E1＝20mV，E2＝50mV，E3＝130mV，E4＝140mV，E5＝10mV；

（4）ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100，E1＝20 mV，E2＝50mV，E3＝130mV，E4＝140mV，E5＝10mV。

對LC的長度分別為5、10m及15m時分別進行計算，結(jié)果表明加入側(cè)向測井儀器電極串后對自然電位測井的影響很小。

下面我們再重點考察隔開2組電極系、具絕緣表面的電纜B的長度LB大小對測量自然電位值U的影響。取E1＝10mV，E2＝20mV，E3＝100mV，E4＝120mV，E5＝30mV，ρ1∶ρ2∶ρ3∶ρ4＝1∶10∶50∶100進行考慮，見表4。此種情況下，未加入側(cè)向測井儀器電極串時測量自然電位值U為126.5869mV。

表3 計算結(jié)果

表4 電纜B的長度LB大小對測量自然電位值U的影響

由此可見，絕緣電纜B的長度LB不能太長，特別是不能太短，否則將（可能嚴(yán)重）影響自然電位測井的結(jié)果。現(xiàn)有儀器設(shè)計中取LB的長度為4.5m，值的相對誤差（－U）／U約為0.9%，符合工程需要。而LB的長度最好選擇為6m或者說6.5m，其相對誤差分別為0.18%及0.27%。

7 結(jié) 論

（1）建立了自然電位測井的數(shù)學(xué)模型，并在交界面上的自然電位差滿足相容性條件的特殊情況及不滿足相容性條件的一般情況，分別在不同的函數(shù)空間中，揭示了該模型在數(shù)學(xué)上的適定性。

（2）在電極的尺度很小時，討論了自然電位測井問題的解的極限性態(tài)，并據(jù)此對該數(shù)學(xué)模型作了適當(dāng)?shù)暮喕?/p>

（3）對交界面上的自然電位差滿足相容性條件及不滿足相容性條件的2種情況，分別提出了如何有效地利用有限元素法求解自然電位測井問題的數(shù)值方法。

（4）根據(jù)實際的需要，提供了眾多的算例，驗證了前述的理論，并對有關(guān)的應(yīng)用提供了指導(dǎo)。

［1］ 張庚驥.電法測井［M］.東營：中國石油大學(xué)出版社，1996.

［2］ 楚澤涵，黃隆基，高杰，肖立志.地球物理測井方法與原理：上冊［M］.北京：石油工業(yè)出版社，2007.

［3］ 潘和平，馬火林，蔡柏林，牛一雄，等.地球物理測井與井中物探［M］.北京：科學(xué)出版社，2009.

［4］ 彭躍軍.自然電位測井中的數(shù)學(xué)方法［J］.應(yīng)用數(shù)學(xué)與計算數(shù)學(xué)學(xué)報，1988，2（2）：35－43.

［5］ Li Tatsien，Tan Yongji，Peng Yuejun，et al.Mathematical Methods for the SP Well－Logging［M］.Applied and Industrial Mathematics，Venice－1（edited by R.Spigler），Kluwer Academic Publishers，1991，343－349.

［6］ Li Tatsien，Tan Yongji，Peng Yuejun.Mathematical Model and Method for the Spontaneous Potential Well－Logging［J］.European Journal of Applied Mathematics，1994，5：123－139.

［7］ 李海龍.均勻地層中自然電位測井方程的“精確解”［J］.數(shù)學(xué)年刊，1996，17A（1）：87－96.

［8］ Zhou Yi，Cai Zhijie.Convergence of a Numerical Method in Mathematical Spontaneous Potential Well－Logging［J］.European Journal of Applied Mathematics，1996，7（1）：31－41.

［9］ 李大潛.有限元素法在電法測井中的應(yīng)用［M］.北京：石油工業(yè)出版社，1980.

［10］Li Tatsien，et al.Boundary Value Probelms with Equivalued Surface and Resistivity Well－Logging，Pitman Research Notes in Mathematics Series 382，Chapman ＆Hall／CRC，1998.

［11］Adams R A.索伯列夫空間［M］.葉其孝，王耀東，應(yīng)隆安，等譯.北京：人民教育出版社，1981.

［12］彭躍軍，一類偏微分方程的邊值問題適定的充分必要條件［J］.同濟大學(xué)學(xué)報，1988，16（1）：91－100.

［13］Gilbarg D，Trudinger N S.二階橢圓型偏微分方程［M］.葉其孝，等譯.上海：上海科學(xué)技術(shù)出版社，1981.

［14］陳亞浙，吳蘭成.二階橢圓型方程與橢圓型方程組［M］.北京：科學(xué)出版社，1991.

［15］Cai Zhijie.Asymptotic Behavior for a Class of Elliptic Equivalued Surface Boundary Value Problem with Discontinuous Interface Conditions［J］.Applied Mathematics，A Journal of Chinese Universities，Ser.B，1995，10（3）：237－250.

［16］李大潛.有限元素法續(xù)講［M］.北京：科學(xué)出版社，1979.

［17］左立華，陳娓.自然電位測井的雙層位勢法［J］.工程數(shù)學(xué)學(xué)報，2008，25（6）：1013－1022.

［18］陳娓，左立華.自然電位測井?dāng)?shù)學(xué)模型的一類求解方法［J］.應(yīng)用數(shù)學(xué)學(xué)報，2009，32（4）：732－745.

［19］Kellogg O D.Foundations of Potential Theory［M］.Berlin：Springer－Verlag，1929.

［20］復(fù)旦大學(xué)數(shù)學(xué)系.數(shù)學(xué)物理方程：2版［M］.上海：上?？茖W(xué)技術(shù)出版社，1961.

Mathematical Modeling and Numerical Method for the Spontaneous Potential Well－logging

LI Tatsien1，CAI Zhijie1，CHEN Wei2，TAN Yongji1，WANG Wenjuan1，WANG Jingnong3
（1.School of Mathematical Sciences，F(xiàn)udan University，Shanghai 200433，China；2.School of Mathematics and Information，Shanghai Lixin University of Commerce，Shanghai 201620，China；3.Technical Center，China Petroleum Logging CO.LTD.，Xi’an，Shaanxi 710021，China）

Spontaneous potential well－logging is a useful method for the petroleum exploration and exploitation.Based on a series of previous results published separately in diverse journals and a recent research，the authors provide a complete theoretical and numerical framework of the mathematical modeling and numerical method together with examples and applications for the spontaneous potential well－logging in the axi－symmetric case.The contents consist of mathematical model and its well－posedness，limit behavior of solutions and a reduced mathematical model，efficient numerical methods，examples and applications.

spontaneous potential well－logging，axi－symmetric case，mathematical modeling，numerical method

P631.321；O29

A

2012－04－05 本文編輯 王環(huán)）

1004－1338（2012）03－0211－14

1復(fù)旦大學(xué)數(shù)學(xué)科學(xué)學(xué)院，E－mail：dqli＠fudan.edu.cn；2上海立信會計學(xué)院數(shù)學(xué)與信息學(xué)院，E－mail：weichen418＠gmail.com；3中國石油集團測井有限公司技術(shù)中心，E－mail：wangjingnong123＠163.com