宋 滔，王緒本

（成都理工大學(xué) 地球物理學(xué)院，成都 610059）

0 前言

在使用有限單元法進(jìn)行地球物理的數(shù)值模擬時，最終需要求解一個大型稀疏線型方程組Ax＝B，而這個線性方程組往往是對稱正定的。在解點源問題時，如果采用齊次邊界條件，那么矩陣A只和地下電性參數(shù)有關(guān)，不包含源信息，這個時候只有列向量B包含有源信息。如果使用高斯消元法求解，那么對于每一個電源點都要進(jìn)行消元和迭代，這大大增加了計算量。阮百堯［1］使用Cholesky分解法，首先將矩陣A分解為兩個對稱矩陣，對于不同的列向量B只需要進(jìn)行一次順代和一次回代便可以求解方程組。但是通過對計算過程的分析發(fā)現(xiàn)，其算法的回代過程幾乎占用了整個求解過程時間的一半，沒有充分利用分解后矩陣的對稱性和稀疏性。在此基礎(chǔ)上，作者對該算法進(jìn)行了一定的改進(jìn)。

1 大型對稱稀疏矩陣的變帶寬存儲

對于變帶寬的對稱稀疏矩陣，采用二個一維數(shù)組來存儲［1］，用一個一維數(shù)組G存儲矩陣的下三角元素，并且以當(dāng)前行的第一個非零元素開始，以另外一個數(shù)組ID存儲矩陣中對角元素在矩陣中的位置，那么通過對對腳元素的索引，就可以實現(xiàn)對整個矩陣元素的訪問。以下面的5×5對稱正定稀疏矩陣為例：

這樣G＝（1，2，5，1，1，2，5，1，5），ID＝（1，3，4，7，9），以下的計算均以這個矩陣為例。

2 Cholesky分解法

矩陣A是正定對稱的，那么可以使用Cholesky分解法進(jìn)行計算，即可以得到

從式（3）與式（4）中可以看出，下三角矩陣L中的元素對應(yīng)矩陣A中下三角的元素，所以分解的結(jié)果矩陣L可以使用存儲A的二個一維數(shù)組來存儲索引。由于LT是L的轉(zhuǎn)置，所以LT中的元素也可以使用這二個一維數(shù)組來訪問。

由于A是稀疏的，含有大量的零元素，所以L和LT中也含有大量零元素，在順代和回代的過程中應(yīng)該充分利用這一特性。

上式（2）中矩陣A分解得到：

對于線型方程組Ax＝B，即LLTx＝B，進(jìn)行順代和回代便可求解。

3 順代和回代

順代的求解過程很簡單，對于Ly＝B，L是下三角矩陣，只需要依次求解y1、y2、…、yn即可。下一步是求解LTx＝y(tǒng)，由于LT是L的轉(zhuǎn)置，同時共用一個存儲空間：

LT是上三角矩陣，通過回代便可完成方程求解。方法（1）：按行依次求解xn、xn－1、…、x1，矩陣元素均可通過矩陣L的索引方式得到；方法（2）：回帶的同時，消去方程組中包含當(dāng)前求解的xi。這兩個方法在L非稀疏時，計算時間是一樣的，但是當(dāng)L是稀疏矩陣時，方法（2）的計算效率是非常高的，如下矩陣所示。

使用方法（2）回帶時，首先解出x5，然后消去方程組中的x5，即將第5列中的元素消去，因為LT中的列是L中的行，在數(shù)組G中是按行存儲的，所以這一步驟是很容易實現(xiàn)的。同時，G是壓縮存儲的，所以在此例中，消去最后一列，僅需要計算一步（消去LT中的l45）。在大型稀疏矩陣中，這種消去列的實現(xiàn)方法，對計算速度的提升是非?？捎^的。

4 Cholesky分解法的求解程序

！N 整變量，輸入?yún)?shù)，方程組的階數(shù)。

！N P 整變量，輸入?yún)?shù)，為系數(shù)矩陣A壓縮存儲的元素個數(shù)。

！A 輸入、輸出參數(shù)，N P個元素的一維實數(shù)組，輸入時存放系數(shù)矩陣的壓縮存儲元素；輸出時存放。

！Cho lesky 分解得到的下三角陣中變帶寬內(nèi)的元素。

！B 輸入、輸出參數(shù)，N 個元素的一維實數(shù)組。輸入時存放方程組右端的n維常向量；輸出時，存放解。

！向量。

！ID 輸入?yún)?shù)，N 個元素的一維整數(shù)組。存放系數(shù)矩陣A的各個對角線元素在壓縮的一維數(shù)組中的。

！位置坐標(biāo)。

5 算例分析

作者將本算法分別用于點源二維電場和二維大地電磁的計算，與傳統(tǒng)的cholesky分解算法還有不帶平方根的cholesky算法進(jìn)行對比。在以下計算分析中，所使用的計算機(jī)為CPU：T6600，2.2GHz，內(nèi)存2G。

5.1 點源二維電場的計算

作者采用矩形網(wǎng)格剖分，電導(dǎo)率分塊均勻變化，雙線性插值有限元進(jìn)行模擬。x方向（水平沿測線方向）包括100個測區(qū)網(wǎng)格，兩邊分別有12個稀疏網(wǎng)格；y方向（垂直方向）一共設(shè)置20個網(wǎng)格，點源點移動100次，即分解后，進(jìn)行100次回代和順帶。地下空間設(shè)置為均勻半空間，分別使用本方法和cholesky分解法求解，如下頁表1所示。

從計算結(jié)果中可以看出，快速回代的cholesky方法對于點電源場求解速度的提升是非常大的，與傳統(tǒng)的順帶回代方法相對比，提升了五倍左右，而且單元數(shù)越多，其提升越明顯。

5.2 二維大地電磁正演模擬

地面測線長4km，測點間距100m，頻率范圍為1 000 Hz～0.01 Hz，對數(shù)等間隔采樣，一共41個頻點。

作者采用矩形網(wǎng)格剖分，電導(dǎo)率分塊均勻變化，雙二次插值有限元進(jìn)行模擬。x方向（水平沿測線方向）包括41個測區(qū)網(wǎng)格，兩邊分別有18個稀疏網(wǎng)格；y方向（垂直方向）一共設(shè)置56 個網(wǎng)格，TE模式下空氣網(wǎng)格設(shè)置14個。地下空間設(shè)置為均勻半空間，分別使用本方法和不帶平方根的cholesky分解法求解，如下頁表2、表3所示。

從表2、表3可以看出，使用快速回代的cholesky算法的速度，要比傳統(tǒng)的cholesky分解法快一倍。

6 結(jié)論

通過算例分析，快速回代的方法充分利用了Cholesky分解后矩陣的對稱性和稀疏性，大大加快了整體求解的速度。該方法對于使用有限元進(jìn)行點源場和大地電磁模擬，所形成的大型對稱稀疏線性方程組的求解是有效的。

表1 點源場求解時間對比Tab.1 The solution time of point source field

表2 TE模式求解時間對比Tab.2 The solution time of TE mode

表3 TM 模式求解時間對比Tab.3 The solution time of TM mode

［1］阮百堯，熊彬.大型對稱變帶寬方程組的Cholesky分解法［J］.物探化探計算技術(shù)，2000，22（4）：361.

［2］劉德貴，費景高，于江永，等.FORTRAN 算法匯編第一分冊［M］.北京：國防工業(yè)出版社，1980.

［3］彭國倫.Fortran 95程序設(shè)計（第一版）［M］.北京：中國電力出版社，2002.

［4］徐士良.計算機(jī)常用算法（第二版）［M］.北京：清華大學(xué)出版社，1995.

［5］阮百堯，徐世浙.電導(dǎo)率分塊線性變化二維地電斷面電阻率測深有限元數(shù)值模擬［J］.地球科學(xué)－中國地質(zhì)大學(xué)學(xué)報，1998，23（3）：303.

［6］熊彬，阮百堯.電位雙二次變化二維地電斷面電阻率測深有限元數(shù)值模擬［J］.地球物理學(xué)報，2002，45（2）：285.

［7］史明娟，徐世浙，劉斌.大地電磁二次函數(shù)插值的有限元法正演模擬［J］.地球物理學(xué)報，1997，40（3）：421.

［8］劉云，王緒本.大地電磁二維自適應(yīng)地形有限元正演模擬［J］.地震地質(zhì)，2010，32（3）：382.

［9］徐世浙.地球物理中的有限單元法［M］.北京：科學(xué)出版社，1994.

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

［11］樸化榮.電磁測深法原理［M］.北京：地質(zhì)出版社，1990.

［12］曾國.大地電磁二維有限元正演數(shù)值模擬［D］.湖南：中南大學(xué)，2008.

［13］王緒本，李永年，高永才.大地電磁測深二維地形影響及其校正方法研究［J］.物探化探計算技術(shù)，1999，21（4）：327.

［14］徐世浙.點電源二維電場問題的付氏變換的波數(shù)k的選擇［J］.物探化探計算技術(shù)，1988，10（3）：235.

［15］強(qiáng)建科，羅延鐘.三維地形直流電阻率有限元法模擬［J］.地球物理學(xué)報，2007，50（5）：1606.