自然電位測井數學模型與求解方法

李大潛蔡志杰陳娓譚永基王文娟王敬農

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

在不同子區域的交界面（記為γk，k＝1，…，5）上，由電流的連續性，應滿足交界面條件

自然電位測井數學模型與求解方法

李大潛1蔡志杰1陳娓2譚永基1王文娟1王敬農3

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

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

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

自然電位測井；軸對稱情形；數學模型；數值方法

0 引 言

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

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

1 自然電位測井的數學模型

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

圖1 地層剖面示意圖

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

式中，n為Γ2的單位外法向量，而?u／?n為u的外法向導數。

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

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

在不同子區域的交界面（記為γk，k＝1，…，5）上，由電流的連續性，應滿足交界面條件

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

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

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

作變換

其中，u0為分塊常數函數

其中，

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

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

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

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

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

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

同樣可將邊界條件（3）和（16）統一地寫為（14），且變換后相應得

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

2 自然電位測井數學模型的解的性質

2.1 解空間

為求解自然電位測井數學模型（SP），考慮等價的變分問題，為此先引進幾個有關的函數集合。

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

其范數定義為

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

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

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

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

下面先考慮解應隸屬的空間［8，12］。從現在起，有關的引理及定理均只列出結論，而略去其證明。

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

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

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

設交界面上的自然電位差滿足相容性條件（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 極限性態

圖2 橢圓型邊值問題

圖3 等位面邊值問題

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

其中，

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

其中，

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

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

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

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

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

（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，將對應于不同電極高度ε（當ε→0時，原電極面退化為直徑為D0的圓周）的相應問題在（D0／2，0）處的自然電位值（其中k＝1，2，3，4對應于上述4組不同物理參數取值情況，下同），并列出ε＝0（即無電極時）在（D0／2，0）處的自然電位值以及在原點處相應的自然電位值，計算結果見表1。

表1 電極高度的影響

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

4 求解方法

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

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

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

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

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

其中，

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

而

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

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

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

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

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

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

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

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

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

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

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

經過對邊界條件的處理，得到線性代數方程組

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坐標，

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

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

4.2.2 移去奇性法

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

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

下面介紹w的構造方法。

令

其中，

容易驗證

它滿足式（59）。

類似地，適當選取

可使wB滿足式（60）。

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

4.2.3 雙層位勢法

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

在實際應用中，常設交界面上的自然電位差Ek（k＝1，…，5）為常數。由前述，利用變換（10），可將交界面條件變為

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

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

圖4 雙層位勢法

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

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

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

引理4.3 對任意給定的P0∈Σ，雙層位勢（71）的外法向導數在P0點是連續的，即

作極坐標變換

易見有

定理4.2 設

則有

且有估計式

其中，R＞a2，Z＞z0。

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

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

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

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

根據上述計算框架，可給出具體的實施細節及有關公式［17，18］。

5 數值模擬

5.1 均勻地層

設整個地層中介質是均勻的，即ρi＝ρj（i≠j）；又設自然電位差E4和E5一直延伸到無窮遠處。此時，井軸上自然電位的精確解表達式［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個節點的網格剖分下做數值計算，同時利用式（82）計算井軸上的自然電位精確值。比較發現它們在井軸上的值誤差很?。ㄒ妶D5）。再取E1＝20mV，E2＝50mV，E3＝130mV，E4＝140mV，E5＝30mV，其余參數不變，它們的結果亦非常接近（見圖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

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

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

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

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

固定物理參數ρ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塊區域（Ω1、Ω2、Ω3）上分片繪制等位線，圖形顯示A點附近電位的變化比較劇烈。

圖9 情形1時的等位線圖

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

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

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

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

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

（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）上的自然電位曲線。從圖中可以觀察到在該曲線拐點處（它相應于目的層與圍巖的分界處）自然電位值的跨度與相應的ΔA值一致：在情形1時，ΔA＝60mV，而在曲線拐點附近自然電位的最大值為90mV，最小值為30mV；而在情形2時，ΔA＝100 mV，而曲線拐點附近自然電位的最大值為110 mV，最小值為10mV。這一事實與物理上的直觀結論吻合。

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

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

5.5 參數對自然電位曲線的影響

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

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

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

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

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

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

圖16 側向測井儀器數據

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

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

為了檢驗側向測井儀器的電極串對自然電位測井的影響，針對圖16所示的數據，我們取半層厚H＝2m，井筒半徑R0＝0.125m，侵入帶半徑Rx0＝0.65m，并對所考慮區域取徑向深度R＝1 000m及縱向高度Z＝600m，通過過渡帶法進行數值計算。記未加入側向儀器電極串時測量電極Γ0上的自然電位值為U，加入側向測井儀器電極串時相應的電位值為?U，取下面的幾組數據進行計算，計算結果見表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時分別進行計算，結果表明加入側向測井儀器電極串后對自然電位測井的影響很小。

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

表3 計算結果

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

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

7 結 論

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

（2）在電極的尺度很小時，討論了自然電位測井問題的解的極限性態，并據此對該數學模型作了適當的簡化。

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

（4）根據實際的需要，提供了眾多的算例，驗證了前述的理論，并對有關的應用提供了指導。

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

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

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

［4］ 彭躍軍.自然電位測井中的數學方法［J］.應用數學與計算數學學報，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］.數學年刊，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］ 李大潛.有限元素法在電法測井中的應用［M］.北京：石油工業出版社，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］.葉其孝，王耀東，應隆安，等譯.北京：人民教育出版社，1981.

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

［13］Gilbarg D，Trudinger N S.二階橢圓型偏微分方程［M］.葉其孝，等譯.上海：上?？茖W技術出版社，1981.

［14］陳亞浙，吳蘭成.二階橢圓型方程與橢圓型方程組［M］.北京：科學出版社，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］李大潛.有限元素法續講［M］.北京：科學出版社，1979.

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

［18］陳娓，左立華.自然電位測井數學模型的一類求解方法［J］.應用數學學報，2009，32（4）：732－745.

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

［20］復旦大學數學系.數學物理方程：2版［M］.上海：上?？茖W技術出版社，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，Fudan 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 本文編輯 王環）

1004－1338（2012）03－0211－14

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