焦石壩地區(qū)頁巖Cole-Cole模型參數(shù)的快速計算方法

劉智穎，張柏橋，許 巍，羅 兵，李銘華，黎 偉

（1.中國石化江漢油田分公司勘探開發(fā)研究院，湖北武漢430223；2.長江大學地球物理與石油資源學院，湖北武漢430100）

Cole-Cole 模型是由PELTON 等人于1932年提出的泥頁巖孔隙電流傳導模型[1-2]，主要用于描述泥頁巖復電阻率隨交流電源頻率變化。標準的Cole-Cole 模型包含4 個參數(shù)，即ρ0、m、τ和c。這些參數(shù)可以反映頁巖的孔隙結構、礦物組分等特征[3-4]。如何利用巖心復電阻率頻譜實驗數(shù)據(jù)快速、準確地求出這4個參數(shù)，是重要的研究課題。

通常，最小二乘法程序的尋優(yōu)算法有模擬退火、遺傳算法及梯度下降算法等。雖然模擬退火、遺傳算法等全局尋優(yōu)算法可以有效避免陷入局部極小值的問題[5-7]，但其運算速度難以保證。相比之下，采用梯度下降算法更能保證運算速度和收斂性[8-9]。然而，過多的參數(shù)不僅造成巨大的運算量，而且增加了迭代算法陷入局部極小值的可能性。為了降低Cole-Cole 模型的復雜程度，PELTON 和MADDEN等[1,10]根據(jù)常見儲層的極化作用的類型及其特點，對標準Cole-Cole模型進行了簡化，將參數(shù)c給定為1/4、1/2 等值。這樣，最小二乘法的尋優(yōu)參數(shù)只剩下ρ0、m和τ這3個值，提高了運算速度。大量的實踐表明，前人提出的Cole-Cole 簡化模型在許多地區(qū)的應用效果良好。然而，實際巖心復電阻率資料表明，焦石壩地區(qū)c值的變化范圍較大，不能將參數(shù)c設為恒定值。同時，該地區(qū)參數(shù)τ的數(shù)量級相較ρ0、m和c很小，造成最小二乘法中的梯度下降算法難以保證一致收斂性，嚴重降低了收斂速度，甚至導致無法收斂。也就是說，此時梯度下降算法是病態(tài)的[11]。具體表現(xiàn)在：尋優(yōu)過程中，τ的迭代值會出現(xiàn)負值，超出其定義域，導致程序計算錯誤而無法收斂。如果簡單地將搜索步長s減小，雖然可以解決該問題，但程序計算速度嚴重減慢。

針對上述問題，通過從極值和漸近性質兩個方面對標準Cole-Cole 模型函數(shù)進行深入分析，給出了參數(shù)τ的一種簡便計算方法，以及ρ0、m和c這3 個參數(shù)的迭代初始值預估算法。不僅使得最小二乘法尋優(yōu)參數(shù)由原來的4個減少為3個，而且給出了ρ0、m和c的初始值，從而有效改善了程序收斂性。同時，還對梯度下降法的迭代步長進行了改進，給出了最優(yōu)步長。通過以上兩項優(yōu)化，程序的運算性能明顯提高。以焦石壩地區(qū)某井的巖心復電阻率頻譜數(shù)據(jù)為例，對改進前和改進后的程序運算速度進行了對比。對比結果表明，在相同的計算精度要求下，改進后的計算程序的迭代次數(shù)下降幅度較大，計算速度顯著提高。

1 參數(shù)的迭代初始值預估算法

選取合適的迭代初始值，對于改善程序的收斂性有很大的幫助。本文給出了ρ0、m、τ和c這4 個參數(shù)的初始值預算法，其中參數(shù)τ的初始值即為其尋優(yōu)結果，因而以上4個尋優(yōu)參數(shù)簡化為3個。

1.1 標準Cole-Cole模型簡介

標準Cole-Cole 模型函數(shù)描述了巖心復電阻率的實部和虛部隨電源頻率的變化情況[12]。形式如下：

式中：ρ為巖心復電阻率，Ω·m；ρ0為直流電阻率，Ω·m；m為極化率；i為虛數(shù)單位；ω為交流電源的角頻率，rad/s；τ為極化弛豫特征時間，s/rad；c是一個與極化類型有關的數(shù)[1]。

1.2 模型函數(shù)的極值及漸近性質

為了給出待求參數(shù)τ、ρ0、m、c初始值的預估方法，從極值和漸近性質兩個方面對式（1）進行了分析。

式中：ρR和ρI分別是巖心復電阻率ρ的實部和虛部。

現(xiàn)有研究表明，ρI（ω）曲線有一個極小值點[2]，首先對此極小值點的位置進行了分析，可以看出，用駐點法求式（3）極小值將導致求導運算煩瑣，為此，應用不等式法求出了ρI（ω）的極值及極值點的位置。

對于c∈（0，1）[1]，可以看出A和B均大于0。因此，應用重要不等式得：

當且僅當：x=1時取等號。

將式（5）代入式（4），并應用復合函數(shù)增減性原理可得：

當且僅當：x=1時取等號。

考慮到?c∈（0，1）有x是ω的單調增函數(shù)，則當ωτ=1 時ρI取極小值。即當時，式（1）的虛部ρI有且僅有一個極小值點。

顯然h（x）是單調遞增函數(shù)。應用復合函數(shù)增減性原理可得，ρR是x（或ω）的單調遞減函數(shù)，沒有極值點。此外，對于x=0的值和x→∞的漸近性質有：

當x→∞時，根據(jù)洛畢達法則不難得出：

最后，對式（2）和式（3）在x=0處的導數(shù)進行了計算。得出：

1.3 初值預估方法

根據(jù)前面的分析，結合本文使用的測量巖心復電阻率的阻抗分析儀TH2839A 的工作頻率，提出了τ、ρ0、m、c這4 個參數(shù)的初始值預估方法。文中將ρ0、m、c的初始值記作、m(0)、c(0)。

TH2839A 的工作頻率范圍是20 Hz ～4 MHz，可以輸出101個頻率下的巖心復電阻率數(shù)據(jù)。

1.3.1 參數(shù)τ的計算

先利用巖心復電阻率頻譜實驗數(shù)據(jù)繪制實部和虛部的頻譜曲線，找出虛部曲線的極小值對應的交流電角頻率值（記作：ω0），然后取ω0的倒數(shù)即可確定τ的初始值。根據(jù)式（6），此方法可以求出τ的嚴格解，因此，將這樣求得的初始值直接視為τ的尋優(yōu)結果。

求ω0的具體方法是：先在儀器輸出的101 個巖心復電阻率的虛部值中找出最小值，然后提取與該最小值相鄰的兩個工作頻率的復電阻率虛部值，再將這3 個虛部值及其各自對應的3 個工作頻率擬合為二次函數(shù)，最后利用二次函數(shù)的極值公式求出其極小值點對應的角頻率（圖1）。圖1中藍色點是測量的復電阻率虛部的最小值，該值不一定是嚴格的虛部最小值；綠色點是與之相鄰的2個數(shù)據(jù)點。虛線是二次擬合函數(shù)，其極小值對應的角頻率即為ω0。

實際算例表明，可以將該方法求出的τ值作為τ的尋優(yōu)結果。精度完全滿足要求，無需將τ代入最小二乘法尋優(yōu)程序中。

圖1 ω0 計算方法Fig.1 Computing method of ω0

1.3.2ρ0和m的迭代初始值

根據(jù)式（9）可以看出，可以選取20 Hz 時的復電阻率的實部值作為。算法為：

1.3.3c的迭代初始值

由式（11）可以看出，由于A和B都是c的函數(shù)，因此，只要知道和，即可確定c的值。和可以根據(jù)最小的2 個工作頻率（20 Hz和22.6 Hz）下測量的復電阻率值，應用差分法求出。這樣，式（11）即為一個關于c的方程組，方程組的解即可作為c的迭代初始值c(0)。

將式（14）和式（15）聯(lián)立為一個關于c的非線性超定方程組可以證明，在c∈（0，1）區(qū)間中，該方程組必存在一個最優(yōu)解。

由于無須嚴格求解c的值，因而采用枚舉法，從0.1，0.2，…，0.9 這9 個值中進行試探。這樣，僅需9次計算，即可找出最符合式（14）和式（15）的c值，作為迭代初始值c(0)。

2 用最小二乘法求ρ0、m和c

最小二乘法程序通過搜索誤差函數(shù)的最小值來確定待求參數(shù)的值。構造了如下誤差函數(shù)（式16），該式為不同電源頻率下式（1）算出的巖心電阻率的實部和虛部（ρRi、ρIi）與實際測量的巖心電阻率的實部和虛部（ηi、μi）之差的平方和。

然而，盡管尋優(yōu)的目標函數(shù)是式（16），但其不適合用于描述迭代的精度。這是由于對于完全相同的計算精度，如果N越大，則式（16）的值越大，這樣不利于檢驗算法的計算精度。因此，用相對誤差式（17）作為判斷迭代精度的標準。

式中：ρRi和ρIi分別表示式（1）計算出的第i個電源頻率下巖心復電阻率的實部和虛部，ηi和μi分別表示實際測量的第i個電源頻率下巖心復電阻率的實部和虛部，N為實驗儀器的頻率的個數(shù)。

2.1 改進前的梯度下降法

假設程序已經搜索至第k步，當前的參數(shù)值分別為ρ0(k)、m(k)、c(k)和τ(k)。此時，式（16）的值及其梯度分別為：

此時，式（19）的值為：

一般給定一個較小的搜索步長s（本文取7.0×10-5），即可使得迭代序列{:k=0,1,2,.....}收斂到4個待求參數(shù)ρ0、m、c和τ的精確結果。

2.2 對迭代步長的改進

為提高程序運算速度，直接用1.3.1 節(jié)中的方法求出τ，并對迭代步長s進行了改進。

仍然假設程序已經搜索至第k步，在的鄰域內將式（21）展開為泰勒級數(shù)，可用于表示的近似值。

其中：

精細化學品課程體系、課程內容的設置對化工專業(yè)學生專業(yè)能力和創(chuàng)新創(chuàng)業(yè)能力的培養(yǎng)起著至關重要的作用．隨著時代的進步，在全民創(chuàng)業(yè)、萬眾創(chuàng)新的時代背景下，高校的人才培養(yǎng)應該致力于滿足社會需求，培養(yǎng)學生的創(chuàng)新創(chuàng)業(yè)能力．精細化學品課程教學應該立足于專業(yè)的客觀實際，依據(jù)當前的發(fā)展現(xiàn)狀合理地安排、調整教學內容，及時了解精細化學品領域最新的發(fā)展動態(tài)、前沿領域，與時俱進，選擇合適的教學手段和方法，不斷更新、完善課程教學．精細化學品課程教學改革任重而道遠，作為高校教師，應該了解當下的實際情況，立足于創(chuàng)新、創(chuàng)業(yè)、應用型化工人才培養(yǎng)目標，積極探索精細化學品課程教學的新理念、新方法．

式中：x的定義與式（2）或式（3）相同。

由式（22）可以看出，第k+1 步的誤差函數(shù)值ε(r→(k+1))可以用一個關于搜索步長s的二次函數(shù)近似表示。這表明，只要s取值范圍足夠小，使得式（22）與ε(r→(k+1))的差別不大，式（22）將存在極小值。也就是說，選擇最優(yōu)的迭代步長，可使得誤差函數(shù)式（18）下降幅度達到最大。該最優(yōu)迭代步長（記作s0）可由式（22）的一階導數(shù)確定[15]。

顯然，如果矩陣H 是正定矩陣，則s0必然給出式（22）的極小值[16]，滿足最優(yōu)搜索步長的條件。

然而，如果直接將s0設定為搜索步長，可能存在兩個問題：①如果H 是非正定的，s0就有可能小于0，而完全根據(jù)數(shù)學定義檢驗H的正定性又是十分煩瑣的；②當s0過大時，式（22）的值將不再近似等于，造成s0不再是最優(yōu)的搜索步長。本文采用了一種較為簡單的確定搜索步長的方法。

事實上，通過觀察式（22）不難發(fā)現(xiàn)，即使H不是正定矩陣，只要s取得足夠小，序列{r→(k):k=0,1,2,.....}仍然是收斂的。

式（25）表明，只要s足夠小，就有ε(r→(k+1))＜ε(r→(k))。說明即使用式（24）求出的s0＜0，不能作為迭代步長，但此時只要給一個較小的搜索步長s即可，本文選擇1×10-7。

綜上所述，以1×10-7作為判斷s0大小的標準，若s0＞1×10-7或s0＜0則取s=1×10-7，若0＜s0＜1×10-7則取s=s0。

3 計算實例

將改進前的算法和改進后的算法用于焦石壩地區(qū)某井的巖心復電阻率頻譜數(shù)據(jù)處理，計算出了ρ0、m、τ、c4 個參數(shù)。下面以焦石壩頁巖氣田某區(qū)塊的8XX 井中2 號樣、18 號樣、23 號樣為例，說明不同算法的計算效果。巖心鉆切檔案資料見表1，實驗測量的復電阻率頻譜數(shù)據(jù)見圖2中的數(shù)據(jù)點。

編制了3個最小二乘法計算程序：第1個程序是將迭代步長設為定值（1×10-7），并隨機給定迭代初始值；第2個程序采用了2.2小節(jié)描述的最優(yōu)步長算法，但迭代初始值同第1個程序；第3個程序同時采用了文中的迭代初始值預估方法和最優(yōu)步長算法。3 個程序迭代計算過程中的相關參數(shù)見表2—表7。用第3 個程序求出的4 個參數(shù)繪制的復電阻率頻譜曲線見圖2中的實線。

表1 8XX井部分巖心鉆切檔案資料Table1 Records of several rock cores of well-8XX

可以看出，在相同的相對誤差條件下，雖然程序2（表4、表5）比程序1（表2、表3）在總迭代次數(shù)方面雖然有一定進步，但算法收斂性仍然不佳，無法滿足實際生產的要求。算法3（表6、表7）不僅收斂性很好，而且迭代次數(shù)大幅減少，程序運算速度顯著提高。

表2 程序1的參數(shù)及計算結果Table2 Parameters and computing results of program-1

表3 程序1的運行過程Table3 Processing of program-1

表4 程序2的參數(shù)及計算結果Table4 Parameters and computing results of program-2

圖2 巖心復電阻率頻譜數(shù)據(jù)（數(shù)據(jù)點）與用程序求出的參數(shù)繪制的頻譜曲線（實線）對比Fig.2 Comparison between real measured complex resistivity data and the calculations from program

表5 程序2的運行過程Table5 Processing of program-2

表6 程序3的參數(shù)及計算結果Table6 Parameters and computing results of program-3

表7 程序3的運行過程Table7 Processing of program-3

4 結論與討論

1）焦石壩地區(qū)頁巖儲層的參數(shù)c變化范圍較大，因而不能采用前人設計的簡化的Cole-Cole 模型，而必須使用標準Cole-Cole模型。同時，參數(shù)τ值過小，嚴重影響了最小二乘法尋優(yōu)程序的收斂性。

2）通過分析標準Cole-Cole 模型函數(shù)的極值點和漸近特性，發(fā)現(xiàn)參數(shù)τ值可以直接利用虛部曲線的極小值位置求出。同時，給出了ρ0、m和c這3個參數(shù)較好的迭代初始值計算方法，顯著改善了最小二乘法程序的收斂性。

3）對于復雜的非線性尋優(yōu)問題，如果將搜索步長設為固定值可能造成收斂速度過慢，而在迭代點的鄰域內，對尋優(yōu)目標函數(shù)進行泰勒級數(shù)展開，可獲得最優(yōu)的搜索步長，使得算法收斂速度加快。

4）本文并未考慮算法的魯棒性。事實上，由于TH2839A 阻抗分析儀的相敏部件的測量誤差，以及焦石壩地區(qū)的一些復雜的礦物，使得在某些頻率段巖心復電阻率與標準Cole-Cole 模型有較大差異。這些異常數(shù)據(jù)對算法魯棒性提出了一定要求，下一步工作將圍繞此問題展開。