直流電阻率三維正演的預(yù)條件共軛梯度法研究

楊金鳳，鄧居智，陳 輝，匡海陽(yáng)

（東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西撫州 344000）

直流電阻率三維正演的預(yù)條件共軛梯度法研究

楊金鳳，鄧居智，陳 輝，匡海陽(yáng)

（東華理工大學(xué) 放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西撫州 344000）

電阻率三維正演；SSOR－PCG；異常場(chǎng)法；有限差分法；E－SCAN法

0 前言

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，需要巨大的計(jì)算量及計(jì)算機(jī)內(nèi)存的電阻率三維正演在PC機(jī)上得以實(shí)現(xiàn)，并逐漸步入實(shí)用階段。但要實(shí)現(xiàn)快速高效的全三維電阻率反演，仍然具有一定的難度，而正演是反演的基礎(chǔ)，因而電阻率三維正演仍然是當(dāng)前研究熱點(diǎn)和難點(diǎn)［2～5］。如何快速、準(zhǔn)確地求解離散后的線性方程組，則是直流電阻率三維正演的核心問題。

在常用的線性方程組的求解方法中，預(yù)條件共軛梯度法是較好的迭代方法，特別是不完全Cholesky共軛梯度（ICCG）法［2，6～8］。但是對(duì)于大型復(fù)雜問題，必然要增加網(wǎng)格剖分的數(shù)量，那么計(jì)算機(jī)內(nèi)存的占用量勢(shì)必要增加，且計(jì)算速度也會(huì)受到影響。這就需要在此基礎(chǔ)上研究其它一些更適合的預(yù)處理方法，進(jìn)一步優(yōu)化內(nèi)存占用量和計(jì)算速度［9］?；谟邢薏罘值臄?shù)值計(jì)算方法，Spitzer［10］采用對(duì)稱超松弛預(yù)條件共軛梯（SSOR－CG）法求解離散的線性方程組，實(shí)現(xiàn)了快速、精確、實(shí)用性強(qiáng)的正演計(jì)算。吳小平［7］引入不完全Cholesky共軛梯度（ICCG）算法及一維稀疏存儲(chǔ)模式，實(shí)現(xiàn)了電阻率三維有限差分正演的快速計(jì)算。作者在SPITZER的研究基礎(chǔ)上，利用SSOR－CG法研究直流電阻率三維有限差分正演問題，采用第一類邊界條件和混合邊界條件，并引入異常場(chǎng)法來(lái)消除總場(chǎng)中由點(diǎn)電源導(dǎo)致的奇異性。

1 直流電阻率法三維正演

1.1 基本方程及其邊界條件

將點(diǎn)電源置于水平地面上A點(diǎn)，I為電流強(qiáng)度。產(chǎn)生電源三維地電場(chǎng)電位u滿足微分方程：

式中 σ（x，y，z）為電導(dǎo)率；f為電流源項(xiàng)，是源位置的函數(shù)。

當(dāng)電流強(qiáng)度為＋I(xiàn)的點(diǎn)電流源時(shí)，源函數(shù)為：

式中 （x，y，z）、（x0，y0，z0）分別為觀測(cè)點(diǎn)和點(diǎn)電流源位置；δ為狄拉克函數(shù)。

u滿足如下邊界條件的電位函數(shù)。

其中 Γs為區(qū)域Ω的地面邊界；?！逓閰^(qū)域Ω的地下無(wú)窮遠(yuǎn)邊界；n是邊界外法線方向的坐標(biāo)變量；r為源點(diǎn)到邊界上點(diǎn)的向徑；θ是邊界外法線方向的單位矢量n與r的夾角。

1.2 異常場(chǎng)法

從上面的方程可以看出，電位函數(shù)存在電流源項(xiàng)，即該函數(shù)在點(diǎn)電流源上存在奇異點(diǎn)。那么當(dāng)求解區(qū)域存在電流源時(shí)，如果采用總場(chǎng)電位法求解，往往造成電源附近解的不可信。一種常見的做法是將場(chǎng)源所在點(diǎn)從求解區(qū)內(nèi)“挖除”，這樣便在被“挖除”的場(chǎng)源點(diǎn)周圍形成新的求解邊界。這些新邊界上的節(jié)點(diǎn)離場(chǎng)源很近，而場(chǎng)源在那里產(chǎn)生的正常場(chǎng)，一般都大大超過不均勻體產(chǎn)生的異常場(chǎng)，因而這些邊界節(jié)點(diǎn)位置可以近似為正常場(chǎng)值。但是由于場(chǎng)源附近的電場(chǎng)變化梯度很大，而用差商代替微商的作法是以節(jié)點(diǎn)之間的電位變化為基礎(chǔ)的，所以為提高差分格式的近似程度，必須在場(chǎng)源附近減小網(wǎng)格步長(zhǎng)。但這樣會(huì)使網(wǎng)格剖分和構(gòu)組差分格式變得復(fù)雜、繁瑣，而減小步長(zhǎng)又勢(shì)必增加總結(jié)點(diǎn)數(shù)，從而增加計(jì)算量和計(jì)算費(fèi)用［11］。

為了提高計(jì)算精度和計(jì)算速度，避免由源奇異性引起的差分計(jì)算復(fù)雜性，作者采用將解析計(jì)算和數(shù)值計(jì)算相結(jié)合的方法，即異常場(chǎng)法［12、13］。假設(shè)電流源附近的電導(dǎo)率為σM，由供電電流源在以電導(dǎo)率σM為均勻半空間產(chǎn)生的正常場(chǎng)為up，與供電電流源無(wú)關(guān)的由電阻率分布異常引起的異常場(chǎng)為us?？傻萌缦玛P(guān)系式：

將式（5）、式（6）代入基本方程，得異常場(chǎng)方程

其中 邊界條件為：

1.3 微分方程離散化

利用有限差分法對(duì)式（7）進(jìn)行數(shù)值求解時(shí)，由于方程是二階微分方程，含有對(duì)電導(dǎo)率空間分布的導(dǎo)數(shù)，要求電導(dǎo)率具有一定的連續(xù)性才可以求解因此對(duì)式（7）進(jìn)行體積分，由格林定理將體積分化為面積分，消除對(duì)電導(dǎo)率的導(dǎo)數(shù)要求，再利用有限差分法構(gòu)成求解方程［2］，即：

對(duì)式（10）應(yīng)用格林定理，有：

式中 si，j，k為包圍區(qū)域Δvi，j，k的表面。

對(duì)求解區(qū)域進(jìn)行六面體剖分，然后在網(wǎng)格點(diǎn)上利用七點(diǎn)差分格式離散，如圖1所示。

圖1 直流電阻率三維正演網(wǎng)格剖分和離散示意圖Fig.1 3－D finite difference grid discretization

從式（11）可見，對(duì)研究區(qū)域及其邊界，邊界條件式（8）和式（9）可直接應(yīng)用到該方程中。利用中心差商近似得代替us/n，并沿體積Δvi，j，k的邊界積分，則第（i，j，k）節(jié)點(diǎn)及其周圍六個(gè)節(jié)點(diǎn)上的電位滿足式（12）。

對(duì)于全部有限網(wǎng)格的所有節(jié)點(diǎn)，式（12）可表示成矩陣形式：

式中 A為耦合系數(shù)；u為待求節(jié)點(diǎn)的異常場(chǎng)列向量；b為已知向量。

2 基于SSOR－PCG法的大型稀疏線性方程組求解

對(duì)于電阻率三維正演問題，大型線性方程組的求解是決定計(jì)算速度和效率的關(guān)鍵一步。而由前面的分析可知，電阻率三維正演所形成的矩陣是含有大量零元素的稀疏矩陣。如果使用一維變帶寬存儲(chǔ)方式，其中仍包含大量的零元素，很大程度上浪費(fèi)了內(nèi)存空間；而使用直接法求解該方程組，其速度又不是很理想。為此，作者采用一維非零元素的壓縮存儲(chǔ)的方式對(duì)系數(shù)矩陣的存儲(chǔ)，并利用預(yù)條件共軛梯度法求解大型線性方程組，極大地降低了內(nèi)存占用量和提高了計(jì)算速度。

2.1 一維非零元素壓縮存儲(chǔ)模式

由系數(shù)矩陣A的特點(diǎn)可知，其每行或每列至多含七個(gè)非零元，即含有七個(gè)非零對(duì)角元。而從式（12）可看出，節(jié)點(diǎn)（i，j，k）的電位值u只與其相鄰六個(gè)節(jié)點(diǎn)的電位值有關(guān)，因此在程序中定義七個(gè)一維整型數(shù)組ipos0（μ）、…、ipos6（μ）來(lái)指示系數(shù)C0（μ）、…、C6（μ）中的每一個(gè)非零元的列號(hào)，這里μ＝1、…、imjmkm，用以指示非零元的行號(hào)，也是節(jié)點(diǎn)編號(hào)。其中離散的節(jié)點(diǎn)按照從左到右，從前到后，從上到下的順序編號(hào)，按編號(hào)i、j、k，i＝1、…、im；j＝1、…、jm；k＝1、…、km，首先沿剖分的x方向讀取，然后再沿y方向，最后沿z方向依次循環(huán)讀取。這樣很容易得到任意位置編號(hào)（i，j，k）的非零元的行號(hào)，按式（15）計(jì)算：

那么對(duì)應(yīng)的列號(hào)也隨即可以確定。

從上面的存儲(chǔ)過程可以看出，其比按行壓縮存儲(chǔ)更為簡(jiǎn)單易行。由于每個(gè)節(jié)點(diǎn)都有一個(gè)位置（i，j，k），這樣就不需要保存每一個(gè)節(jié)點(diǎn)對(duì)應(yīng)的連接系數(shù)在矩陣中的行和列。而根據(jù)節(jié)點(diǎn)編號(hào)可以很容易地反向計(jì)算出節(jié)點(diǎn)的具體位置（i，j，k），同時(shí)其周圍六個(gè)節(jié)點(diǎn)的位置和它編號(hào)μ也可以求得。在計(jì)算矩陣與列向量的乘積時(shí)，只有七個(gè)節(jié)點(diǎn)相對(duì)應(yīng)的列向量中非零元與列向量進(jìn)行相乘，這樣就可以消除與零元素相乘的運(yùn)算，提高了解方程的速度，所需的存儲(chǔ)空間也得以大大減少。

2.2 SSOR－PCG迭代算法

共軛梯度迭代算法，對(duì)于系數(shù)矩陣接近單位陣時(shí)其對(duì)應(yīng)條件數(shù)很小，計(jì)算速度很快；而當(dāng)A的條件數(shù)cond（A）＞102時(shí)，則收斂速度就很慢。所以對(duì)于電阻率三維正演計(jì)算的類似問題，由于網(wǎng)格大小和物性參數(shù)可能相差幾個(gè)數(shù)量級(jí)，所形成的大型稀疏線性方程組的稀疏矩陣A的特征值的變化范圍很大，因而對(duì)應(yīng)的條件數(shù)必然很大，所以就要對(duì)式（14）的線性系統(tǒng)進(jìn)行預(yù)處理，改善矩陣A的條件數(shù)。然后再運(yùn)用CG法解該式，計(jì)算速度將會(huì)大大提高。

作者在本文采用SSOR預(yù)條件方法［14］，先將系數(shù)矩陣A分解為：

其中 E是下三角矩陣；I是單位矩陣。

其中 L＝I＋wE；w是松弛因子。

解式（14）的SSOR預(yù)條件共軛梯度迭代過程如下：

其中 r0、p0是共軛梯度迭代的初始化向量；u0是解的初值；α、β為常數(shù)。

從以上的計(jì)算流程可以看出，預(yù)條件矩陣很容易求得，且不增加計(jì)算內(nèi)存需求。這樣在整個(gè)計(jì)算過程中，主要就是系數(shù)矩陣A與任一列向量的乘積，以及預(yù)條件矩陣與任一列向量的乘積。而A中每一行的非零元素不參加計(jì)算，再加上用一維稀疏存儲(chǔ)，所以可以快速求得。關(guān)鍵是如何快速的求得預(yù)條件矩陣M的逆與任一向量的乘積，直接求解很困難。注意到L是下三角陣，而LΤ是上三角陣，可先將ti＝（LLΤ）－1ri轉(zhuǎn)化為求方程組LLΤti＝ri，然后用高斯消去法先順帶、后回帶解得ti。由于矩陣L中僅有少量的非零元素，并以稀疏存儲(chǔ)方式存儲(chǔ)，因此順帶和回帶過程僅涉及少量非零元素的操作，均可以非常簡(jiǎn)單、快速地完成。

松弛因子w的選擇是控制迭代收斂的重要參數(shù)，對(duì)模型給予不同的w值試算結(jié)果如下頁(yè)圖2所示。從圖2可看出，w在1.3～1.7之間迭代次數(shù)變化較小，對(duì)迭代過程影響不是很大；而w取“0”時(shí)，該迭代是發(fā)散的；最終取w值為“1.5”，此時(shí)的收斂速度最快。SSOR－PCG迭代的收斂標(biāo)準(zhǔn)是是第一次迭代后殘差的L2范數(shù)。

由此可見，單獨(dú)的SSOR－PCG迭代技術(shù)或系數(shù)矩陣稀疏存儲(chǔ)技術(shù)均不能有效的解決問題，但將兩者有機(jī)地結(jié)合，再加上合適的松弛因子的選擇就可形成快速高效的迭代算法。

3 模型算例

3.1 垂直接觸帶模型

垂直接觸帶模型如下頁(yè)圖3（b）所示，左邊介質(zhì)的電阻率是1Ωm，右邊的電阻率為10Ωm。單

其預(yù)條件矩陣定義為：位點(diǎn)電源位于左邊介質(zhì)上，距界面6m，采用39× 39×20網(wǎng)格剖分，計(jì)算機(jī)配置為Intel Core 3.0GHz雙核處理器，內(nèi)存3.0GB，32為XP系統(tǒng)。SSOR－PCG算法迭代546次收斂，耗時(shí)3.69s。圖3（a）是單位點(diǎn)電源的總電場(chǎng)法，異常場(chǎng)法計(jì)算結(jié)果和理論電位曲線的比較。異常電場(chǎng)法計(jì)算的電位曲線和理論電位曲線幾乎完全吻合，平均誤差為0.69%；而總場(chǎng)法計(jì)算的電位曲線相差很大，其平均誤差為9.03%，尤其是在源點(diǎn)附近誤差很大，這是由源的奇異性引起的。

圖2 對(duì)不同的松弛因子w值SSOR－PCG算法所需的迭代次數(shù)Fig.2 The number of iteration required by SSOR－PCG for different relaxation factors

3.2 層狀介質(zhì)模型

模型2為二層介質(zhì)，第一層電阻率為10Ωm，第二層的電阻率為100Ωm，厚度為3m。采用39 ×39×20的網(wǎng)格剖分，SSOR－PCG迭代742次，耗時(shí)5.23s。

表1是單位點(diǎn)電源的總場(chǎng)法、異常場(chǎng)位法計(jì)算結(jié)果ut、us和理論電位值u的比較。異常場(chǎng)法計(jì)算的平均誤差△uS為0.16%，而總場(chǎng)法的平均誤差△ut為6.48%，精度提高一個(gè)多量級(jí)。特別是在電源附近，總場(chǎng)法計(jì)算誤差較大，異常場(chǎng)法計(jì)算精度較高。從表1還可看出，利用源點(diǎn)移除法，去除源點(diǎn)后總場(chǎng)法的平均誤差為3.45%，誤差減小了近一半，但還是比異常場(chǎng)法的平均誤差高，其平均誤差只有0.12%?？梢姰惓?chǎng)法對(duì)于提高計(jì)算精度，有顯著的效果。

3.3 低阻異常模型

為了說(shuō)明SSOR－PCG迭代算法在直流電阻率三維正演模擬中計(jì)算的可靠性，作者設(shè)計(jì)了一個(gè)典型的低阻異常模型，設(shè)地下有一個(gè)40m×40m× 20m，電阻率為1Ωm的低阻長(zhǎng)方體，如圖4所示其頂部埋深為20m，圍巖的電阻率為100Ωm；圖4（b）是電性分界線在地面（xoy平面）的投影。將計(jì)算區(qū)域剖分為37×37×21個(gè)單元。

作者采用的是高密度三維觀測(cè)方法ESCAN［15］，將所布設(shè)的電極網(wǎng)格化為21×21，其中x方向和y方向均布設(shè)21個(gè)電極，點(diǎn)距、線距均為5m。圖4（b）標(biāo)出地面上E－SCAN法測(cè)量的441個(gè)電極位置。

圖5為y＝0低阻異常體中心E－SCAN裝置視電阻率擬斷面圖，其中橫坐標(biāo)x為算術(shù)坐標(biāo)，縱坐標(biāo)AM/2為對(duì)數(shù)坐標(biāo)。由圖5可見，低阻模型會(huì)產(chǎn)生一個(gè)低阻異常區(qū)和高阻異常區(qū)，低阻異常區(qū)范圍為AM/2從2.5m至20m，x方向從－40m～40m變化。高阻異常區(qū)位于低阻異常體下方，橫向范圍與模型基本一致。

圖6為低阻異常體E－SCAN裝置在AM/2分別取3.54m、10.31m、15.2m、27.5m時(shí)的視電阻率等值線圖。圖6（a）中由于其探測(cè)深度還未到低阻異常體，視電阻率基本無(wú)變化；圖6（b）及圖6（c），由于受低阻異常體影響，在中心位置有一個(gè)明顯的低阻異常區(qū)，范圍比低阻異常體大；下頁(yè)圖6（d）其探測(cè)深度穿過異常體時(shí)，在中心位置出現(xiàn)一個(gè)明顯半徑大約為20m近似圓形高阻假異常區(qū)由此可見，E－SCAN裝置視電阻率受低阻異常體影響，隨著AM的增大由低阻異常轉(zhuǎn)向高阻異常異常范圍也隨之變小，并且E－SCAN裝置對(duì)異常體橫向分辨較好，縱向分辨率差。

4 結(jié)論

作者在本文詳細(xì)推導(dǎo)了異常場(chǎng)法所形成的方程，并對(duì)方程進(jìn)行差分離散得到大型稀疏線性方程組，將求解該大型稀疏線性方程組的SSOR－PCG迭代算法、系數(shù)矩陣的稀疏存儲(chǔ)技術(shù)，以及消除點(diǎn)電源奇異性的異常場(chǎng)法三者有機(jī)的結(jié)合在一起。通過對(duì)二個(gè)簡(jiǎn)單模型的模擬結(jié)果與解析解對(duì)比，結(jié)果表明該方案的數(shù)值解與理論解基本吻合，說(shuō)明將解析計(jì)算和數(shù)值計(jì)算相結(jié)合，并利用該算法求該正演計(jì)算中的異常電位，能消除點(diǎn)電源奇異性影響，在保持快速計(jì)算的基礎(chǔ)上，顯著提高了計(jì)算精度，而且內(nèi)存需求大大減少。另外對(duì)低阻地電模型的正演模擬結(jié)果表明，SSOR－PCG是電阻率三維有限差分穩(wěn)定高效的計(jì)算方法，實(shí)現(xiàn)了快速精確的電阻率三維正演計(jì)算，為其三維反演奠定了基礎(chǔ)。

［1］ 傅良魁.電法勘探教程［M］.北京：地質(zhì)出版社，1983.

［2］ DEY A，MORRISON H F.Resistivity modeling for arbitrarily shaped three－dimensional structures［J］.Geophysics，1979，44（4）：615.

［3］ MA Q Z.The boundary elemental method for 3－D DC resistivity modeling in layered earth［J］.Geophysics，1998，63（2）：399.

［4］ 阮百堯，熊彬，徐世浙.三維地電斷面電阻率測(cè)深有限元數(shù)值模擬［J］.地球科學(xué)，2001，26（1）：73.

［5］ 魯晶津，吳小平，KLAUS SPITZER.直流電阻率三維正演的代數(shù)多重網(wǎng)格方法［J］.地球物理學(xué)報(bào)，2010，53（3）：700.

［6］ LI Y G，SPITZER K.Three－dimensional DC resistivity forward modeling using finite elemental in comparision with difference solutions［J］.Geophys.J.Int.，2002，151（3）：924.

［7］ 吳小平，徐果明，李時(shí)燦.利用不完全Cholesky共軛梯度法求解電源三維地電場(chǎng)［J］.地球物理學(xué)報(bào)1998，41（6）：838.

［8］ WU X P.A 3－D finite－element algorithm for DC resistivity modeling using the shifted incomplet Cholesky conjugate gradient method［J］.Geophys.J Int.，2003，154（3）：947.

［9］ 劉斌，李術(shù)才，李樹忱，等.基于預(yù)條件共軛梯度法的直流電阻率三維有限元正演研究［J］.巖土工程學(xué)報(bào)2010，32（12）：1846.

［10］SPITZER K.A 3－D finite－difference algorithm fo DC resistivity modeling using conjugate gradien methods［J］.Geophys.J.Int.，1995（123）：903.

［11］羅延鐘，張桂青.電子計(jì)算機(jī)在電法勘探中的應(yīng)用［M］.武漢：武漢地質(zhì)學(xué)院出版社，1987.

［12］LOWRY T，ALLEN M B，SHIVER P N.Singularity removal：a refinement of resistivity modeling techniques［J］.Geophysics，1989，54（6）：766.

［13］ZHAO S，YEDLIN M J.Some refinements on the finite－difference method for 3－d dc resistivity modeling［J］.Geophysics，1996，61（5）：1301.

［14］吳小平，汪彤彤.利用共軛梯度算法的電阻率三維有限元正演［J］.地球物理學(xué)報(bào)，2003，46（3）：428.

［15］黃俊革，王家林，阮百堯.三維高密度電阻率ESCAN法有限元模擬異常特征研究［J］.地球物理學(xué)報(bào)，2006，49（4）：1206.

book=147,ebook=147

1001—1749（2012）03—0303—07

P 631.3＋22

A

10.3969/j.issn.1001－1749.2012.03.11

楊金鳳（1984－），女，碩士，主要從事電法勘探正反演研究。

國(guó)家科技支撐計(jì)劃（2009BAB43B03；2011BAB04B03）

2012－01－05改回日期：2012－01－09