利用ICCG迭代技术加快电阻率三维正演计算

2005/9/12 17:48:28

　　摘要　一般而言，有限差分法求解点源三维地电场正问题所形成的大型稀疏线性方程组Ax＝b，直接解法的计算效率极低。本文从系数矩阵A的不完全Cholesky分解及矩阵特征值的特点等角度，说明了不完全Cholesky共轭梯度(ICCG)迭代技术可大大提高电阻率三维正演速度的内在原因。结合矩阵A的稀疏存储模式，使得内存需求也大大减少。
　　关键词　电阻率法勘探　三维正演　预处理　不完全的Cholesky共轭梯度法
　　中国图书资料分类法分类号　P631.34

3-D RESISTIVITY FORWARD CALCULATION ACCELERATEDBY ICCG ITERATION TECHNIQUE

Wu Xiaoping　Xu Guoming　(University of Science and Technology of China)

　Abstract　For the large sparse linear equations:Ax＝b，which are formed from the finite difference method used to solve the 3-D forward problem of geoelectrical field,in general,the computational efficiency with direct method is quite slow.In this paper,from the view point of the characteristics of the incomplete Cholesky decomposition of matrix A and its eigenvalue,the internal cause of greatly increased 3-D resistivity forward speed using the incomplete Cholesky conjugate gradient (ICCG) iteration technique is explained .Introducing the row-indexed sparse storage mode to store matrix A,the internal storage demand is greatly decreased.
　Keywords　resistivity prospecting;3-D forward;incomplete Cholesky conjugate gradient;sparse storage mode

1　引言
　　随着近年来高密度电测技术的应用，传统的直流电阻率法又展现了新的活力，成为浅层水文、工程、环境、考古等与人类社会生活密切相关的探测领域的重要手段^{［1，2，3］}。而浅层地质目标多表现为复杂的三维电性结构，对解释带来很大困难。典型的如大地电磁测深(MT)中的静位移效应^［4］，其根本的解决方法就是三维反演解释，这也是高精度电法发展的趋势^［5］。
　　众所周知，正演是反演的基础和前提。而大多数情况下，地电三维正演只能通过有限差分^［6］、有限单元^{［7，8，9］}等求得数值解，因而探寻高效快速的三维正演计算就显得尤为重要。基于有限差分的数值计算方法，笔者曾引入不完全Cholesky共轭梯度(ICCG)方法^{［10，11］}，同时结合系数矩阵的稀疏存储模式，求解二、三维地电场的正问题，在提高计算速度和减少内存需求上均获一定成功，取得初步效果。本文从讨论系数矩阵A的不完全Cholesky分解及矩阵特征值的特点出发，探讨了利用ICCG迭代技术加快电阻率三维正演计算的内在原因，希望能为ICCG方法在地球物理三维正演计算中的更好应用和发展起到抛砖引玉的作用。
　　
2　点源三维地电场的有限差分计算
　　设点电源的电流强度为I，位于坐标点(x₀，y₀，z₀)处，则其产生的点源三维地电场电位φ(x，y，z)满足微分方程：

其边界条件为：

 115-62.gif (2654 bytes)

式中　σ(x，y，z)——地下介质的电导率；
　　　　　n——边界外法线方向的坐标变量；
　　　　　θ——n和r的夹角；
　　　　　n——边界外法线方向的单位矢量；
　　　　　r——源点到边界上点的向径。
　　上述边值问题用有限差分的数值计算方法求解^［6，8］，可获得复杂三维结构上的地电场分布。只是需对整个研究区域进行N_x×N_y×N_z的三维网格剖分，未知节点个数太多，最后形成一大型稀疏线性方程组

Aφ＝S。(2)

其中A为容量矩阵，是一大型稀疏对称正定带状矩阵，有如下形式：

96-62.gif (2786 bytes)

这里C_t、C_b0、C_l、C_r、C_f、C_ba、C_p分别是节点(i，j，k)和上、下、左、右、前、后及自身节点的连接系数。可见矩阵A每行最多只有7个非零元素。φ＝(φ₁，φ₂，…，φ_{N_xN_yN_z})^T为节点组上的电位值；S是与供电电流源有关的右端项，只在供电节点p上有值，即S＝(0，……，0，S_p，0，…，0)^T，而S_p＝I。解以上差分方程求得电位φ(x，y，z)。
　　电阻率三维正演的速度基本取决于解此大型方程组的计算效率。直接解(2)式的计算效率非常低；首先系数矩阵A半带宽为N_x×N_z＋1(或N_x×N_y＋1或N_y×N_z＋1，与节点编号有关)，直接解法需用二维带状压缩存储其下三角矩阵带宽内的元素，存储量最少为(N_x×N_y×N_z)×(N_x×N_z＋1)。对于40×40×20的三维网格剖分，即使是单精度，其内存要求亦需103 MB，可见需要巨大的机器内存；再者直接解法对A做完全Cholesky分解；A＝LL^T，必然要对A带宽内的大量零元素进行操作，非常费时，因而计算速度也极慢。
　　近年来发展的预条件共轭梯度迭代算法，同时引入一维按行索引的稀疏存储模式存储系数矩阵A，则可有效避免直接解法存在的以上问题^{［10，11］}。
　　
3　共轭梯度(CG)方法
　　共轭梯度法是50年代初由Hestense and Stiefel^［12］提出的解对称正定线性方程组的迭代方法。现今预条件共轭梯度方法已成为解大型稀疏矩阵极为有效的方法。
　　线性方程组　　　　　Ax＝b　　　　　　　　(3)

的共轭梯度算法：
　　 97-63.gif (1998 bytes)

　　从上可见，整个计算流程只要求矩阵A与一个列向量的乘积，而矩阵A中各行的零元素对于A与一个列向量乘积是没有贡献的。根据按行索引稀疏存储模式，可将大型稀疏矩阵A压缩成两个一维数组sa和ija^［10］，sa是实型数组，存储矩阵A下三角非零元素，维长不超过4N_x×N_y×N_z，加上共轭梯度求解的4个一维辅助数组也不超过10N_x×N_y×N_z，大大减少了内存要求；ija是整形数组，存储的整型值是对数组sa的元素在矩阵A中位置的索引。因此由一维存储sa和ija，不难求得A与任意一个列向量的乘积。
　　对于(3)式的线性系统，当矩阵A接近单位矩阵时，CG迭代方法收敛快；而通常情况下，由于网络大小、物性参数可能相差几个数量级，所形成的大型方程组的系统矩阵A的特征值λ变化范围必然很大，A的条件数cond(A)＞10⁴，CG迭代就非常缓慢。而不完全Cholesky共轭梯度(ICCG)方法，正是运用矩阵A的不完全Cholesky分解进行预条件因子化，改善矩阵A的条件数，使得解大型差分方程的ICCG迭代大大加快。
　　
4　ICCG方法及其解释
4.1　不完全Cholesky因子化
　　不完全Cholesky分解^{［13，14］}的计算快速简单，其分解如下：

A≈CC^T。　　(5)

其中C是下三角矩阵，可从对角矩阵D求得。D的对角元素由下式定义：

98-63.gif (621 bytes)

显然求和中a_jk＝0的元素是无贡献的，D可简单求得。C由下式确定：

C＝UD^－1／2。

其中U是下三角矩阵，对角元素u_jj＝d_jj，非对角元素u_jk＝a_jk(k＜j)，即U的非对角元素与A一样。亦即不完全Cholesky分解因子C和矩阵U一样，其非对角元素与A相同并具同样的稀疏性，因此无需开辟另外空间来存储C，只需存储对角矩阵D即可，大大地节省了机器内存。
4.2　解大型稀疏线性方程组的ICCG方法
　　由(5)式的不完全Cholesky分解将(3)式重新写成：

［C^－1A(C^T)^－1］C^Tx＝C^－1b，　　(6)

如果(CC^T)^－1是矩阵A的逆的近似，则C^－1A(C^T)^－1将是一个近似的单位矩阵。由上述可知，CG方法应用于矩阵C^－1A(C^T)^－1将大大加快收敛速度。
　　将共轭梯度(CG)方法应用于改进后的(6)式，比较(3)式作以下替换：

99-64.gif (790 bytes)

代入(CG)算式(4)稍作整理，便得到改进方法，即不完全Cholesky共轭梯度(ICCG)算法：
　　 100-64.gif (2300 bytes)

4.3　ICCG方法的解释
　　从上述的讨论中看到，ICCG迭代收敛速度的快慢完全取决于(5)式不完全Cholesky分解的近似程度。近似程度越好，C^－1A(C^T)^－1越接近单位矩阵，则迭代收敛速度越快。为了解释(5)式的近似性，我们以三维模型的有限差分计算为例。
　　设地下有一个40 m×40 m×9 m的低阻长方体，顶部埋深为4 m，围岩电阻率为100 Ω*m，异常体电阻率为10 Ω*m，三维网格为15×15×10。将其形成的系数矩阵A做完全Cholesky分解：A＝LL^T，可以观察到完全分解因子L^T中非零元素的大小(均为绝对值)在图1所示的方向上快速衰减至零。图2为某一半带宽内非主对角元素大小的变化情况，就是很好的证明。众所周知，完全Cholesky分解是稳定的，因此将其中较小的元素置为零，形成如式(5)的不完全Cholesky分解则和完全Cholesky分解应该是很接近的。

101-64.gif (2636 bytes)

图1　Cholesky分解矩阵的元素大小示意图

102-64.gif (1393 bytes)

图2　完全Cholesky分解因子中的元素值在远离主对角线
元素方向而衰减

　　从矩阵特征值的角度看，由于(5)式的近似，(CC^T)^－1A将接近单位矩阵，进一步而言就是，矩阵(CC^T)^－1A的所有特征值应近似为1.0。
　　图3引自参考文献［13］，用来说明ICCG方法中矩阵(CC^T)^－1A的特征值变化。这里的A是用五点差分格式解二阶椭圆方程形成的五对角系数矩阵，具体参数见文献［13］。

103-64.gif (2100 bytes)

图3　矩阵A、(L₀L₀^T)^－1A及(L₃L₃^T)^－1A的特征值

　　需要说明的是，图3中(L₀L₀^T)^－1A及(L₃L₃^T)^－1A分别是两种不同的不完全Cholesky共轭梯度法ICCG(0)和ICCG(3)的情况。本文的ICCG方法相当于ICCG(0)方法，即不完全Cholesky分解因子L₀较A的下三角矩阵多0条非零的次对角元素，也就是L₀与A具同样的稀疏性。同理，ICCG(3)是指L₃较A的下三角矩阵多3条非零的次对角元素，它是根据五对角系数矩阵的完全分解因子L^T性质(类似图1)，在L₃中多放置3条非零的次对角元素，以期使L₃和完全Cholesky分解更接近。从图3看到，对未知节点数不多的这样一个二维问题，形成系数矩阵A的条件数cond(A)＝7.503／0.058≈130，而cond(L₀L₀^T)^－1A)＝1.231／0.119≈10，条件数大大降低；并且(L₀L₀^T)^－1A的特征值绝大部分在1.0附近，Kershaw^［15］也得到类似的结果。由此不难理解运用ICCG方法解大型稀疏线性方程组的高效快速。当然文中涉及ICCG(3)方法，可参见文献［13］。
　　同理于二维问题，对于三维问题形成的七对角系数矩阵，亦可发展相应的ICCG(5)方法，即指L₅较A的下三角矩阵多5条非零的次对角元素。它也是根据上述完全分解因子L^T的性质(见图1)，在L₃中多放置5条非零的次对角元素，以期使L₅L₅^T和完全Cholesky分解更接近，而使三维正演的ICCG迭代收敛更快。
　　
5　ICCG方法与其它方法的比较
5.1　ICCG方法和直接方法
　　首先我们大致估计一下几种方法的计算量：
　　a.　直接方法：完全Cholesky分解约需n(m＋1)(m＋2)／2＋2n(m＋1)次乘法运算，其中n是系数矩阵的阶数，m是半带宽。
　　b.　ICCG方法：每次迭代约需16n次乘法运算^［13］。
　　c.　CG方法：每次迭代约需10n次乘法运算^［13］。
　　可以看到，直接方法的计算量随系数矩阵的阶数(节点数)的增大呈指数上升，而ICCG方法则是线性上升。因此随着网格数的增多，ICCG方法的计算量必然大大少于直接方法。
　　表1是不同网格情况下，在Pentium133微机上用NDP-Fortran编译，上述模型三维正演的ICCG方法和直接方法计算时间的比较。

表1　ICCG方法和直接方法计算时间的比较　s

网格	直接法	ICCG法
9×9×10	1.5	0.5
15×15×10	10.0	1
19×19×10	46	1.5
25×25×10	205	2.5
39×39×20	—	15

　　实际计算也说明了随着网格节点数的增多，直接方法的计算量呈指数上升，而ICCG方法则是线性上升的。另外对于39×39×20网格，ICCG方法仅需1.2 MB的机器内存，而直接方法需近100 MB的机器内存，目前的机器配置还难于满足其内存要求。这些均表明对于三维问题，ICCG方法较直接方法在计算速度及内存要求上的优势是非常明显的。
5.2　ICCG方法和CG方法
　　虽然CG方法每次迭代所需的运算量较ICCG方法少，但由于ICCG方法是在对(3)式的线性系统预条件因子化后进行的，所以迭代次数较CG方法少得多。图4是网格15×15×10情况下，ICCG方法和CG方法的迭代收敛情况，ICCG迭代35次，CG迭代371次。

104-65.gif (5024 bytes)

图4　ICCG和CG迭代收敛速度比较

　　
6　利用ICCG方法计算的正演结果
　　由于从(3)式到(6)式的推导是精确成立的，因此(5)式不完全Cholesky分解的近似完全，不影响解的精确性。它只是用来改善(3)式的线性系统，使得C^－1A(C^T)^－1更接近单位矩阵，加快ICCG迭代的收敛速度。
　　一般来说，三维模型计算结果无解析解对比。不过可以根据地电场的互换性原理来判断三维地电场计算结果的正确与否。根据地电场的互换性原理，对于同一地下介质，无论其电性分布如何，在A点供电测得M点的电位值一定与在M点供电测得A点的电位值相等。同样以4.3节中三维模型的有限差分计算为例，只是网格加密为39×39×20。表2列出其20对互换性测量值的对比，互换性相对误差都非常小，大都仅0.01%左右。由此可见，文中利用ICCG迭代技术的点源三维地电场有限差分计算结果是完全可靠的。

表2　互换性测量值对比表

顺序	供电点	测量点	电位值	供电点	测量点	电位值	相对误差／%
1	10	11	1.605200	11	10	1.605199	0.0000668
2	10	12	0.7971908	12	10	0.7971889	0.0002393
3	10	13	0.5189223	13	10	0.5189380	0.0030209
4	10	14	0.3843988	14	10	0.3844308	0.0083267
5	10	15	0.3045164	15	10	0.3045579	0.0136232
6	10	16	0.2500160	16	10	0.2500578	0.0167121
7	10	17	0.2080722	17	10	0.2081102	0.0182619
8	10	18	0.1722966	18	10	0.1723047	0.0047048
9	10	19	0.1629066	19	10	0.1629062	0.0002470
10	10	20	0.1555632	20	10	0.1555598	0.0021936
11	10	21	0.1491902	21	10	0.1491798	0.0069716
12	10	22	0.1436742	22	10	0.1436583	0.0110664
13	10	23	0.1288378	23	10	0.1288225	0.0118665
14	10	24	0.1166769	24	10	0.1166637	0.0113090
15	10	25	0.1069725	25	10	0.1069626	0.0092564
16	10	26	0.0989768	26	10	0.0989832	0.0065114
17	10	27	0.0921943	27	10	0.0922173	0.0249149
18	10	28	0.0863158	28	10	0.0859638	0.4077905
19	10	29	0.0811427	29	10	0.0811449	0.0027271
20	10	30	0.0765403	30	10	0.0765413	0.0012946

7　结束语
　　运用预条件共轭梯度迭代技术加快三维正演速度的关键在于其因子化方法。不完全Cholesky预条件因子化正是充分利用了差分方程系数矩阵A简单的带状、稀疏性质，大大改善了线性系统的条件数，形成的ICCG方法应用于电阻率三维有限差分正演取得良好效果。值得指出的是，对于有限元计算形成的线性系统则较为复杂，系数矩阵中的非零元素位置与剖分有关，用上述预条件方法将使方程组的解极不稳定。Manteuffel^［16］、Papadrakakis & Dracopoulos^［17］也分别提出了一些针对有限元计算形成的线性系统的预条件化方法，但应用到三维问题的有限元计算似乎还有待进一步研究。

作者简介　吴小平　男　32岁　博士　电磁法勘探

作者单位：吴小平　徐果明　(中国科学技术大学地球与空间科学系　合肥　230026)

参考文献
　1　White P A.Electrode arrays for measuring groundwater flow direction and velocity.Geophysics,1994;59(2):192～201
　2　底青云等.高密度电阻率法在珠海某防波坝工程中的应用.地球物理学进展，1997；12(2)：79～88
　3　Noel M,Xu B.Archaeological investigation by electrical resistivity tomography: a preliminary study.Geophys.J.Int,1991;107:95～102
　4　王家映.关于大地电磁的静校正问题.地质科技情报，1992；11(1)：69～76
　5　何继善.电法勘探的发展和展望.地球物理学报，1997：40(增刊)：308～316
　6　Dey A,Morrison H F.Resistivity modeling for arbitrarily shaped three-dimensional structure.Geophysics.1979;44(4):753～780
　7　Holcomble H T,Jiracek G R.Three-dimensional terrain correction in resistivity surveys.Geophysics,1984;49:33～59
　8　周熙襄.电法勘探数值模拟技术.成都：四川科学技术出版社，1986
　9　徐世浙.地球物理中的有限单元法.北京：科学出版社，1994：131～297
　10　吴小平，徐果明.不完全Cholesky共轭梯度法及其在地电场计算中的应用.石油地球物理勘探，1998；33(1)：89～94
　11　吴小平，徐果明，李时灿.利用不完全Cholesky共轭梯度法求解点源三维地电场.地球物理学报，1998；33(6)
　12　Hestense M,Stiefel E.Methods of conjugate gradients for solving linear systems.J.Res.Nat.Bur.Standards,1952;49:409～436
　13　Meijerink J A,Van Der Vorst H A.An iterative solution method for linear system of which the coefficient matrix is a symmetric matrix.Math.Comp.,1977;31:148～162
　14　Smith J T,Conservative modeling of 3-D electromagnetic fields,Part Ⅱ:biconjugate gradient solution as an accelerator,Geophysics,1996;61:1319～1324
　15　Kershaw D S.The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations.J.Comp.Phys.,1978;26:43～65
　16　Manteuffel T A.An incomplete factorization technique for definite linear systems.Math.Comp.,1980;34:473～497
　17　Papadrakakis M,Dracopoulos M C.Improving the efficiency of incomplete Cholesky preconditionings.App.Numer.Methods,1991;7:603～612

煤炭网版权与免责声明：

凡本网注明"来源：煤炭网www.coal.com.cn "的所有文字、图片和音视频稿件，版权均为"煤炭网www.coal.com.cn "独家所有，任何媒体、网站或个人在转载使用时必须注明"来源：煤炭网www.coal.com.cn "，违反者本网将依法追究责任。

本网转载并注明其他来源的稿件，是本着为读者传递更多信息的目的，并不意味着本网赞同其观点或证实其内容的真实性。其他媒体、网站或个人从本网转载使用时，必须保留本网注明的稿件来源，禁止擅自篡改稿件来源，并自负版权等法律责任。违反者本网也将依法追究责任。如本网转载稿件涉及版权等问题，请作者在两周内尽快来电或来函联系。