辛普森公式的推廣形式及應用

劉智穎， 王向公， 任威龍, 龍耀萍

(1.長江大學 油氣資源與勘探技術教育部重點實驗室， 湖北 武漢 430100；2.長江大學 地球物理與石油資源學院， 湖北 武漢430100；3.北京吉奧特能源項目責任有限公司， 北京 100195)

當前, 簡單而常見的數值積分方法是無限細分法(亦稱“黎曼和”).在應用過程中，需將區間劃分成盡可能多的份數以保證計算精度.然而，較細的劃分會導致計算時間增加，而且當被積函數變化幅度較大時，其誤差增大,故辛普森公式應運而生.常見的辛普森公式大都只適用于直角坐標系.事實上很多積分需要在柱坐標和球坐標下進行，如果做坐標變換后再求解則十分繁瑣而且會造成誤差累積.因此，本文分析、研究并推導了柱坐標系和球坐標系下的辛普森公式，且成功應用于地球物理測井中自然伽馬射線強度函數的數值積分運算.

1 三維柱坐標系和球坐標系下辛普森公式的推廣

1.1 辛普森公式的常見形式及其適用范圍

常見的辛普森公式(公式(1))只能應用于直角坐標系下一元函數的數值積分

(1)

1.2 柱坐標系和球坐標系下的辛普森公式

為書寫方便，將柱坐標系下的坐標(ri,θj,zk)及球坐標系下的坐標 (ri,θj,φk)皆記作i,j,k，被積函數在節點處的值記作fi,j,k.腳標i,j,k可以是整數或半整數(其意義與式(1)同).

柱坐標和球坐標剖分體元內插值函數的體積分值有相似的矩陣[2]表達形式(公式2).

(2)

將式(2)中的腳標x替換為i即可得到柱坐標系下的公式.式中矩陣的意義如下：

其中A和B兩個矩陣只需將Z矩陣中符號zk分別換成θj和ri即可.

將式(2)中的腳標x替換為k即可得到球坐標系下的公式.式中矩陣的意義如下：

其中Φ、Θ和R三個矩陣只要將柱坐標系公式中Z矩陣中的zk分別換成φk、θj和ri即可.

1.3 復化辛普森公式

在實際應用中，由于求解區域的邊界不規則，有的節點可能會落到積分區域以外.嚴格地講，這時需要對求解區域重新剖分以適應不規則的邊界，但一般只要已經滿足了精度就可以省去這些繁瑣的步驟，只要將落到積分區域以外的節點的函數值設定為零[3]即可.

2 辛普森公式的推廣形式在地球物理測井中的應用

自然伽馬測井[4]是常見的地球物理測井方法之一.自然伽馬射線強度函數值的計算在分析、應用自然伽馬測井資料時十分重要，該函數描述了自然伽馬測井值(GR值)與地層放射性物質分布情況的關系，需用柱坐標系下的體積分表示.應用辛普森公式的推廣形式便可實現對該函數值的計算[5].

通常，為了便于解釋人員對相關數據的分析，儀器實測的GR值需要進行“井眼校正”[4].而校正系數是通過計算自然伽馬射線強度函數得到的.將不同的井筒直徑做為變量輸入自然伽馬射線強度函數就可以計算出一系列離散的GR井眼校正系數.應用這些離散的數據可以繪制出以井筒直徑為橫坐標，以GR井眼校正系數為縱坐標的圖版(圖1).

圖1 GR井眼校正圖版

3 結論和認識

通過分析和總結當前常用數值積分算法的不足，研究數值積分的基本原理，從而推導辛普森公式在三維柱坐標系和球坐標系下的推廣形式.對此做出如下總結.

(1)提出并推導了辛普森公式的推廣形式，拓展了辛普森公式的適用范圍，提高了三維柱坐標系和球坐標系下數值積分的運算速度和精度，彌補了常見數值積分方法的不足.

(2)將辛普森公式的推廣形式應用到地球物理測井領域，完成了自然伽馬射線強度函數值的計算，從而研制了GR井眼校正圖版.

[1] 關治,陸金甫.數值分析基礎[M].北京: 高等教育出版社, 1998.

[2] 錢椿林.線性代數[M].3版.北京: 高等教育出版社, 2000.

[3]錢能. C++程序設計教程[M].2版.北京: 清華大學出版社, 2003.

[4]黃隆基.放射性測井原理[M].北京: 石油工業出版社,2000.

[5] 商慶龍.新型自然伽馬測井響應的三維數值模擬以及數據合成與高分辨率處理技術[D].吉林：吉林大學，2013：43-56.