的)单元矩阵[K(E)]的各非零元素。而各矩形单元矩阵[K(E)] 之和便构成新的(或最终的)刚度矩阵[K]
?K???K(E) （2.72）
E此新的刚度矩阵[K]不再包括代表各矩形中心节点的行和列,它是一个阶数等于求解区内矩形网格节点数的方阵，其规模比先前包括矩形中心节点时大为减小。同时,按(2.61)式构成的线性方程组的节点函数值列矢量{V}和场源列矢量{I},也都不包括各矩形中心节点的对应量，其维数也等于求解区内矩形网格上的节点数。自然，解此线性方程组最终算出的结果,只包括这些矩形网格节点的函数值。
应该指出,上述消掉矩形中心节点函数值，使线性方程组降阶的方法，不同于直接对求解区进行矩形网格剖分的作法。它本质上仍然是作较细的三角网格剖分，因而对起伏地形及电性异常体的模拟和线性插值处理的近似程度都比直接作矩形网格剖分好。不过，它形成刚度矩阵[K]的计算量较后者大,并且存放各三角形单元电导率所占用的存贮量也较多。为了克服此缺点，可对内部电性均匀(?1??2??3??4)和不均匀（?1、?2、?3和?4不完全相同）的矩形单元分别进行处理，对导电性均匀的矩形单元（在整个求解区内占大多数），只存贮一个电导率值，并且可用简化的公式计算[K(E)]的元素。
2.2.2 线性方程组的解法
有限元法导出的线性方程组（2.61）式系数矩阵（刚度矩阵[K]）具有以下特点：
1）由于Krs=Ksr，即单元矩阵[K]是对称的，而总体的刚度矩阵[K]乃是单元矩阵[Ke]之和，所以刚度矩阵[K]也是对称的。此外，数学上还可以证明[K]是正定的； 2）刚度矩阵的阶数等于求解区的总节点数。通常，在求解高密度电阻率法二维正演问题时，节点数可达几千甚至更大，所以[K]又是大型的。不过，在刚度矩阵[K]中，只有对角元素及其附近的元素不为零，其余皆为零元素，故一般地说，[K]是稀疏矩阵。并且，如果对节点作恰当的编号，使中心节点与其相邻节点之
e
e
e
编号相差尽量小，则刚度矩阵[K]的非零元素将比较集中地分布在对角线附近的带状区域内，因此，[K]又是带状的。
总之，有限单元法在求解电法勘探正演问题时，所得到线性方程组的系数矩阵[K]是一个大型、稀疏和带状的对称、正定矩阵。有限单元法大部分时间要花在解线性方程组上，特别是在进行反演时，解方程组的耗时所占比例更大。因此，选用合适的方程组的解法，也是十分重要的。在有限单元法中，求解线性方程组最为常用的方法有迭代法（如SOR法和CG法）和直接法。实践表明，对这种大型、稀疏和带状的对称、正定系数矩阵，最好采用乔勒斯基（Cholesky）分解法求解。
乔勒斯基分解是将对称、正定的系数矩阵[K]分解为一个同阶的下三角矩阵[L]与其转置矩阵[LT]之积（故也称LLT分解），即
[K]?[L]?[L]T （2.73） 于是方程（2.61）变成
[L]?[L]T{U}?{I} （2.74） 它等价于以下两个方程
[L]T?{U}?{W} （2.75）
[L]?{W}?{I} （2.76）
先解第二个方程组，按如下“回代”公式计算中间变量{W}的各元素
Wi?(Ii??LipWp)/Liip?1i?1(i?1,2,???,N) （2.77）
然后，以算出的中间变量{W}作为已知右端项，求解第一个方程组，按如下“回代”公式计算{U}的各元素
U?(W??LU)/Liipipiip?i?1NN（2.78） (i?N,N?1,???,1)
以上两式中约定?(....)?0和?(....)?0
1N?10下三角矩阵[L]的非零元素Lij按下列递推公式算出
j?1??(Kij??LipLjp)p?1?1i?1?2Lij??(Kii??Lip)2p?1???0??j?i???j?i? （2.79）
??j?i??(i?1,2,??????,N)在按（2.79）式作乔勒斯基分解，计算下三角矩阵[L]的对角元素时，包含几次开方运算，它的计算量较大，计算精度稍差。为避开开方运算，可采用乔勒斯基分解的变形━LDLT分解，即将系数矩阵[K]分解为三个同阶矩阵之积 [K]?[L]?[D]?[L]T （2.80）
其中，[L]为单位下三角矩阵，[D]为一对角矩阵，它们的元素可按下列递推公式算出
j?1~~L?aij?Kij??aipjp?p?1? ~/d?Lij?aiji?~Ldi?Kii??aipipp?1i?1?j?i??? （2.81） ??j?i???(i?1,2,???,N)在对[K]作了LDLT分解之后，方程（2.61）变成
[L][D][L]T{U}?{I} （2.82） 它等价于以下两个方程
[L]?{W}?{I} （2.83） [L]T{U}?[D]?1?{W} （2.84） 它的解可按下列“回代”公式算出
Wi?Ii??LipWpp?1i?1(i?1,2,???,N) （2.85）
和
Ui?Wi/di??LpiUpp?i?1N(i?N,N?1,???,1) （2.86）
上述解线性方程组的直接法，相对于迭代法的缺点是要占用较多计算机储存单元。所以，在电子计算机上实现上述计算时，还应注意以下技巧，以尽量减少计算量和储存量：
（1）因为刚度矩阵[K]是对称的，所以只需计算和储存其一半元素，通常只需形成和储存下半个矩阵；
（2）尽量重复利用刚度矩阵[K]、场源列矢量{L}以及其中间矩阵来储存新生成的矩阵（如[L]等），而不是再开辟新矩阵；
（3）因[K]是稀疏的，因此，甚至并不需要储存其下半个矩阵的全部元素，而只储存其中的非零元素，因而可进一步大量节省计算机的储存量，这种储存方法称为“紧缩存贮方法”。数学上已证明，带状的对称、正定矩阵[K]作LLT或LDLT分解后，下三角矩阵[L]也是带状的，且其半带宽S不超过[K]的半带宽，故在采用紧缩存贮的条件下，[K]和[L]仍可先后占用相同的存贮单元。
如前所述，在采用图2.3所示矩形求解区和网格剖分时，实用的只是矩形网格上的节点，由于求解区通常沿水平（X）方向较长，节点数（N）较多，而沿垂直方向（Z）较短，节点数（M）较少，故将节点从左边开始，逐列由上至下进行编号。这样，刚度矩阵[K]的半带宽S最小（S=M+2）。
由于对称性，实际上只需要存贮[K]或[L]下三角中半带宽内的元素，故采用二维数组的形式存放这些元素。二维数组的列数等于半带宽S，而其行数与总节点数（N×M）相同。采用这种紧缩存贮方式，对用LLT或LDLT分解法解有限单元法的线性方程组是经济和方便的。并且，由于采用上述紧缩存贮，形成刚度矩阵[K]和分解、计算下三角矩阵[L]所涉及的元素数目都大为减少，所以不仅节省了存贮单元，而且也节省了计算量。
在计算高密度电阻率法正演问题时，需要计算不同供电电极的位置的电场分布。改变供电电极的位置不仅影响线性方程组的右端(场源)项{I}，而且使刚度矩阵[K]中与?2边界节点有关的元素发生变化，这是因为（2.69a?2.71c）各式中的????[IK?K1(?rK)?cos?k]k?2n?[Ik?1nK?K0(?rk)] 与供电点A位置有关。为了减
少计算量，在电极距不太大或求解区范围足够大时，可选用装置（主要是供电极）处于某一适中位置时的η值，固定不变。当装置移动或变化，供电电极位置改变时，只改变线性方程组的（2.61）的右端项{I}，而其系数（刚度）矩阵[K]保持不变。这样，在计算不同供电位置的节点函数值时，只需形成和分解一次刚度矩
阵[K]；而对每一供电位置，改变其右端项后由已分解好的[L]和?L?矩阵作回代，
T便可算出相应的节点函数值。因为与形成和分解矩阵[K]相比，回代的计算量甚小，所以采用上述处理方法计算多个供电位置的电场，增加的计算量并不大。这是提高电剖面法计算速度的一条重要措施。
2.2.3 视电阻率的计算
对于给定的供电电极位置经反傅氏变换（见附录B）计算出给定的各地面节点之电位后，便可进而计算给定装置给定位置上的视电阻率?s。
比如，对于温纳装置（在高密度电阻率法中有人称之为?装置，AM=MN=NB=na）。若给定的供电极A（I）和B（-I）之节点编号为nA和nB，则方程组（2.61）的右端列矢量{I}，除第nA和nB个元素分别为InA?I和InB??I之外，其它元素皆为零。由此算出A、B间各节点的电位值后，便可按下式计算测线上不同位置的视电阻率 ??sUnM?UnN （2.87)
?K?I?式中，UnM和UnN分别是测量电极M和N所在节点（编号为nM和nN）的电位值；装置系数K?按下式计算
K???(xA?xM)(xA?xN)?2?na (2.88)
xM?xN式中的xA,xB,xM,xN分别是节点nA,nB,nM,nN的横坐标值，a为最小电极距，n为电极隔离系数。
若使用偶极装置（在高密度电阻率法中有人称之为?装置，AB=MN=a，AM=BN=（n+1）?a）则
K???n(n?1)(n?2)a (2.89) 若采用三极装置（AM=na ，B->∞，MN=a），则
K??4?na (2.90)