0 引言

大地电磁三维反演方法有众多实现，但相对成熟且应用更为广泛的有二类：基于数据空间的OCCAM方法(Siripunvaraporn et al., 2005)和非线性共轭梯度(NLCG)法(Newman and Alumbaugh, 2000；Egbert and Kelbert, 2012)，主要原因在于这两种方法的实现代码WSINV3DMT (Siripunvaraporn and Egbert, 2009)和ModEM(Kelbert et al., 2014)都曾经是或仍是开放源码，学者可以基于源码做更多的测试和二次开发工作.另外，有限内存拟牛顿法(LBFGS)等其他最优化方法也逐渐受到重视(Avdeev and Avdeeva, 2009).

近年来通过对大量实用数据集和理论模型的反演研究，我们发现目前的三维反演算法仍然存在一些缺点：1)作为最稳定方法的OCCAM实现，WSINV3DMT没有使用针对大规模问题的高斯牛顿隐式求解方案而采用了显式存储雅可比矩阵的方式，每次迭代中至少需要和数据个数相同的正演计算次数，需要巨量的内存和CPU时间，同时还必须多次搜索最佳正则化因子，因此在实际问题中的应用有较大限制；2)NLCG和LBFGS采用了隐式雅可比方案，内存消耗小、效率高，但目前的做法是对整个目标函数直接做线性化逼近，下降方向直接用上一次或前几次的总目标函数梯度和模型的改变量逼近.这一做法在正则化因子动态变化的场景下是不精确的，有可能导致反演发散；3)在NLCG和LBFGS反演中，目前的实现方法是直接忽略前文提到的下降方向不精确性，当拟合差下降比较慢时直接将正则化因子按固定比例缩小(如ModEM的实现)，并且对下降方向再次进行一个拟高斯牛顿的预条件化处理(Newman and Alumbaugh, 2000；Avdeeva, 2008)，但这种预条件技术并没有解决正则化因子变化时下降方向不准确的问题.上述后两点决定了自适应正则化技术很难在NLCG和LBFGS传统框架下有效应用，这也是为何NLCG和LBFGS反演相对更依赖初始模型，稳定性比OCCAM反演弱的根本原因.

为了达到减小内存需求和快速计算的实用化目标，本文在当前的NLCG和LBFGS框架下引入OCCAM思想，首先构建更精确的不依赖前几次正则化因子的下降方向，然后引入一种稳定且接近于OCCAM的，让每次迭代中模型尽量光滑的更新正则化因子方案，最终实现和当前NLCG或LBFGS反演计算效率相当的更稳定反演方法.

1 自适应正则化有限内存拟牛顿(AR-LBFGS)反演的基本理论

根据前文分析，若想基于目标函数梯度信息构建更精确的迭代下降方向，并实现自适应策略，必须对正则化约束下的反问题进行更深入的分析.

一般来说，对于三维大地电磁问题，正则化反问题的目标函数如下所示：

(1)

其中，Δd为当前模型响应和观测数据的残差向量，C_d为观测数据的协方差矩阵，一般为对角阵，且每个元素为相应观测数据的方差，m为反演网格电阻率组成的向量，C_m为模型协方差矩阵，一般可以取为梯度或拉普拉斯算子的内积得到的对称正定矩阵，Ψ_d和Ψ_m被称为数据目标函数和模型目标函数，正则化因子λ控制Ψ_m在总目标函数Ψ中的权重.一般来说, 在一次迭代中λ越大，下降方向就更倾向于让模型更光滑，反之则倾向于让数据拟合差更小.

为了让目标函数适应不同规模的网格剖分和数据集大小，我们将其归一化为以下形式：

(2)

其中，N_d和N_m分别为参与拟合的数据个数和参与反演的模型单元个数，从而保证不同规模的反问题在数值计算中都有较好的数值计算精度.上式也可用C_d^－1和C_m^－1的二范数进行归一化，但从计算复杂度考虑，取(2)式的形式已经能够满足AR-LBFGS反演的要求，求新的目标函数的极小与(1)式的极小等价.在正则化因子自适应场景中，反演迭代格式可以统一写为

(3)

其中，k为迭代次数，Δm_k(λ_k)表示在正则化因子为λ_k的情况下求取的下降方向向量，α_k为线搜索步长.对于目标函数(2)式，其牛顿反演的下降方向满足：

(4)

其中，H为φ_d的海塞矩阵，和表示经过N_d和N_m归一化后的数据和模型协方差矩阵的逆，为避免海塞矩阵的求解并保证对称正定性，高斯牛顿法使用来近似难以求解且无法保证正定的海塞矩阵，得到：

(5)

其中，J_k为模型响应的雅可比矩阵，而OCCAM反演正是基于高斯牛顿法的一种变种，其最重要的思想是在每次迭代中反复搜索，使得φ_d下降到某个水平的最大λ_k值，从而得到当前迭代满足数据目标函数充分下降的最光滑模型.

分析公式(2)和(5)可知，OCCAM反演稳定的原因有二：1)总目标函数(2)式由一个非线性最小二乘和一个线性最小二乘组合而成，高斯牛顿方程仅对数据目标函数的海塞阵做了近似，线性最小二乘问题的模型目标函数是无需近似的，并且当前迭代的下降方向和前几次迭代的正则化因子无关，下降方向更精确的保证了光滑能力；2)通过一种搜索机制得到最佳正则化因子下的高斯牛顿步，从而让整个迭代过程的每一步既充分降低了数据拟合差，又保证了模型足够光滑.

反观传统的NLCG和LBFGS在解式(2)的极小问题时则无法保证以上两点，它们都是把总目标函数当成一个固定正则化因子的非线性最小二乘问题来直接求极小以达到最大的计算效率.对于NLCG和LBFGS，总目标函数梯度满足：

(6)

其中，g_d和g_m分别是数据和模型目标函数的梯度向量，总梯度向量g和当前迭代的正则化因子相关，对于第k次迭代，NLCG利用总梯度g_k－1和g_k以及下降方向Δm_k－1直接构建共轭梯度方向来逼近，避免了显式求取雅可比矩阵并求解方程组，可以认为是对海塞逆矩阵的某种求和式逼近，典型的逼近公式为

(7)

其中，被称为预条件阵，一般用一维或均匀空间的高斯牛顿总海塞逆矩阵来简化.由(7)式可知，这种近似虽然利用来做高斯牛顿预条件缩放，但上式仅在固定正则化因子的场景下能保证精确性，如果正则化因子发生了变化，g_k－1和g_k分别在不同的正则化因子下求得，那么下降方向的偏离是不可避免的，同时由于在下降方向的逼近中把模型目标函数这一线性最小二乘问题也当成非线性问题来处理，这本身也损失了迭代过程中正则化项的光滑能力.因此，尽管之前很多学者在二维共轭梯度框架下尝试过自适应正则化方案(Zhdanov, 2002；陈小斌等，2005)，但仍没有从根本上解决对上次正则化因子的依赖性问题，这些方案的正则化因子调整并不能保证反演每次迭代的下降方向为当前线性化的总目标函数下降方向.

对于传统LBFGS反演而言，情况和NLCG类似，不同的是用有限次(一般取3~20)之前迭代的总梯度和模型修正向量来构建对海塞逆的求积式更新，同NLCG一样，这种更新是基于总梯度的，而总梯度在不同迭代中若使用的是不同的正则化因子，那么下降方向将变得奇异，造成反演失败.

综上所述，要实现稳定的自适应正则化反演策略，首要条件是必须对目标函数(2)式的两项分开考虑，对于线性最小二乘的模型目标函数，其海塞矩阵是无需逼近的，若能够如高斯牛顿那样得到一个精度稍低的基于梯度信息的数据目标函数海塞阵逼近，则可以得到不依赖之前正则化因子的下降方向.

1.1 更精确的下降方向求取

为了让下降方向不依赖前几次迭代的正则化因子，我们需要将目标函数的非线性和线性最小二乘两项分开考虑，只利用前几次迭代的数据目标函数梯度及模型修正量信息去逼近数据目标函数海塞矩阵本身，然后如高斯牛顿法那样形成新的逼近方程.由于NLCG方法是求和式更新，且无法保证对牛顿步的尺度缩放(即保证步长为1)，因此难以在这种策略中应用.相对而言LBFGS方法可以采用一系列的对称正定矩阵求积式逼近海塞，因此，本文基于LBFGS的海塞更新公式进行下降方向的改造.

我们的基本思想是，只基于数据目标函数梯度和模型修正量，采用逼近得到的代替公式(5)中的，从而避免雅可比矩阵内积的求取，从而新的下降方向满足：

(8)

其中，g_{d, k}和g_{m, k}分别是公式(5)中右端表示数据目标函数和模型目标函数的两个梯度向量，我们这里不用雅可比矩阵表示是因为在反演中g_{d, k}并没有显式用雅可比矩阵计算而是通过解一次“拟正演”问题得到.对于g_{d, k}的求取则沿用NLCG反演的方法(Newman and Alumbaugh, 2000; Rodi and Mackie, 2001).

LBFGS是拟牛顿法解非线性最优化问题用得最广泛的一种更新公式，对于海塞矩阵和其逆，它有两种更新逼近形式，一种是传统LBFGS三维大地电磁反演应用的直接逼近海塞逆的方法(Nocedal, 1980)，

(9)

其中，m为需要利用的前有限次迭代梯度及模型修正量个数，一般取3~20之间，H_k表示海塞矩阵，我们这里没有使用最优化领域的用B_k表示海塞矩阵而用H_k表示的做法是为了和传统大地电磁反演的公式表示一致，V_k=I－ρ_ky_ks_k^T, ρ_k=1/y_k^Ts_k, s_k=m_k－m_k－1, y_k=g_k－g_k－1, 注意这里的g_k是逼近函数的梯度，并非对应公式(6)的总梯度.在实际反演中公式(9)还有一个非常快速的隐式双循环递归更新公式，非常利于编程实现，这也是为什么很多学者采用LBFGS直接逼近总目标函数梯度逆矩阵的原因.我们要实现的方法无法使用公式(9)而只能采用LBFGS的逼近海塞矩阵本身的另外一个等价形式(Byrd et al., 1994).为简化公式，首先定义两个矩阵存储m次向量对，

(10)

其中，每个列向量可以如下表示(g_k表示逼近函数的梯度)：

(11)

从而，海塞矩阵的LBFGS逼近可以用矩阵S_k和Y_k表示为

(12)

其中，, L_k和D_k分别是m×m矩阵S_k^TY_k的下三角(除主对角元素)和主对角矩阵，我们使用g_{d, k}和g_{d, k－1}代替公式(11)中的g_k与g_k－1即可得到数据目标函数φ_d的海塞矩阵逼近.

在三维大地电磁问题的实际编程实现中，显式存储也会消耗多次大规模向量的内积，因此显示存储不是最佳方案，同时要考虑到在光滑约束问题中正则化因子的影响有可能导致非正定，因此需要构建一个快速计算和一个任意x向量乘积的算法，用迭代法(如共轭梯度法)求解公式(8)，如果的求取返回了错误结果，则表示前几次迭代的梯度和向量已无法满足对称正定性逼近，此时需要进入“重启模式”丢弃之前存储的前几次迭代信息，让迭代从最速下降方向重新开始.

为了快速求取，定义F_k为矩阵的Cholesky分解，假设分解存在，定义一个新的2 m×2 m缓存矩阵，

(13)

引入2 m×1的缓存向量p来表示公式(12)中右边的矩阵和向量积，则p可以等价的通过解以下两次三角阵系统得到，

(14)

求出向量p后，可以很快速的通过下式计算：

(15)

从公式(13)和公式(14)可以看出，在迭代求解方程(8)之前，缓存矩阵W_k无需重复计算，并且可以很容易的通过对F_k矩阵的Cholesky分解计算是否存在来判断是否对称正定以决定反演迭代是否进入重启模式.如果进入重启模式，则重启后的下降方向满足：

(16)

其中，γ_k为与传统NLCG反演一样的预条件标量(Rodi and Mackie, 2001).公式(8)和公式(16)即为我们构建的不依赖前几次正则化因子的更精确下降方向，它是将目标函数中的非线性最小二乘和线性最小二乘分开考虑得到的结果，新的下降方向对于不同的正则化因子更新都是无偏的.这种更新也被作为预条件用于NLCG反演中(Newman and Alumbaugh, 2000)，但仍是基于总梯度更新的，没有解决正则化因子改变导致的下降方向偏差，在自适应正则化约束反演的场景，传统NLCG反演的下降方向可写为以下更简单形式：

(17)

其中f表示求取共轭方向的修正向量的函数，这里我们可以看到传统NLCG反演的两个问题，一是对于的逼近没有用(8)式而是基于总梯度，二是即使使用了(8)式作为预条件，但f的自变量λ_k－1 仍然会导致下降方向偏离真实共轭梯度方向.而传统LBFGS反演中，这一问题将更加严重，因为前m次迭代的正则化因子各不相同时，(9)式得到的方向会严重偏离真实的拟牛顿方向.总而言之，直接求解(8)式是自适应正则化反演的必要条件.在共轭梯度中应用本文的思想必须对(17)式进行改造，将f中的有关λ_k－1的项修正为λ_k下的值，但这样修改后是否仍然满足共轭梯度的理论仍需要更多的研究.

1.2 正则化因子自适应策略

自适应正则化反演的第二个重要问题是确定在每次迭代中最佳正则化因子的取值，目前主要有两类方法，一类将正则化因子的更新设定为随迭代次数衰减的光滑或分段函数的先验方法，如L曲线法(Hansen and O′Leary, 1993)、因子递减法(吴小平和徐果明, 1998；吴小平和汪彤彤, 2001)、目标函数比值法(陈小斌等，2005)、阶梯函数法(Zhdanov et al., 2007)等；第二类是在迭代中搜索一个能让数据拟合差减小到一定水平的最大正则化因子的后验方法，OCCAM反演(Constable et al., 1987)和阻尼最小二乘法都使用了这一思想.

对于三维大地电磁反演，正则化因子的先验更新策略在过去的研究中并未得到深入研究，以当前最流行的NLCG反演的ModEM程序为例，其更新策略是当数据拟合残差的减小速度太慢时就将正则化因子按固定比例缩小，这种自适应更新有两个缺点：第一是无法确定初值，如果初值设置得很小，反演一开始就处于几乎无正则化约束的状态，在数据目标函数降低到正则化项能起作用之前模型已经变得非常奇异了；第二是反演后期的正则化因子会不断的急剧变小，如果在反演中期模型不够光滑，那到了反演后期正则化项会完全不起作用.而后验方法(如OCCAM反演)需要在一次迭代中搜索多次最佳正则化因子，这对于海塞的高精度高斯牛顿逼近而言是可以接受的，但对于LBFGS逼近而言由于数据目标函数收敛效率低于高斯牛顿，并非每次迭代都是牛顿步(个别迭代的步长不为1，需线搜索确定)，从而需要内、外两层线搜索过程.

为了避免后验方法的低效缺点，借鉴目标函数比值法思想(陈小斌等，2005)，根据方程(8)的特点构建了一个基于梯度二范数的先验更新策略，保证每次反演迭代中方程(8)的右端两项的向量二范数处于同一数量级，即：

(18)

上式的含义在于让公式(8)的右端两项对于线性系统解的贡献大致相当，对均匀空间而言，||g_{m, k}||为0，因此在实现中正则化因子的取值可以用下式确定：

(19)

其中ε为一个很小的正数，可以认为是模型目标函数梯度二范数的门槛值，即目标函数梯度二范数小于ε的模型，我们认为等同于均匀半空间.在反演中如果ε取得太小会导致初始模型为均匀空间时第一次迭代不收敛，我们在均匀单位梯度算子正则化约束的情况下进行了大量的合成数据不同初始模型的试验，结果表明ε=0.001时即保证高效收敛又能保证算法稳定，这可大致表示我们认为相邻电阻率平均变化1.001倍时的模型与均匀半空间等价.必须指出的是，如果正则化约束项取其他形式时，ε的取值可以根据以上方法试验得出，即设定每个单元电阻率自然对数随空间变化的绝对值为0.001的模型，然后计算||g_{m, k}||即可以得到针对新的正则化约束形式最合理的ε取值.

考虑到AR-LBFGS的数据目标函数海塞逼近采用了相对高斯牛顿精度更低的拟牛顿形式，且反演网格不一定满足高精度的梯度计算，当数据目标函数计算精度较低时||g_{d, k}||随迭代的变化会不平滑，从而导致反演中||g_{m, k}||和λ_k剧烈的变化，使反演稳定性降低，我们可以对公式(19)进行多次迭代平均修正，以让正则化因子随迭代次数变化的曲线变得光滑，增强反演的稳定性，新的正则化因子更新公式可以定义为

(20)

这表示我们使用前三次迭代的数据目标函数和模型目标函数的梯度二范数的平均值之比作为正则化因子，试验表明这样的改动会稍微增加一些迭代次数，但稳定性加强了很多.我们也用高斯牛顿法进行了对比试验，结果表明高斯牛顿法中公式(19)是完全稳定的，在AR-LBFGS中更加的自适应正则化因子更新方案是(20).

1.3 AR-LBFGS反演的计算流程

根据公式(19)得到正则化因子后，可以用迭代法(如共轭梯度法)求解方程(8)，得到当前迭代的下降方向.与高斯牛顿法不同，LBFGS的迭代步长并不严格为1，所以需要经过一个插值线搜索(Moré and Thuente, 1994)过程来确定最佳步长，保证迭代稳定且收敛高效的步长，必须满足以下的Wolfe条件(Nocedal and Wright, 1999):

(21)

再考虑传统NLCG和LBFGS反演的做法，由于没有考虑正则化因子变化的情况，所以在反演实现中往往将条件错误的设定为

(22)

从而导致线搜索得到的满足充分下降和曲率条件(22)式的步长并非当前迭代的最佳步长，这也是造成反演不稳定的原因之一.在AR-LBFGS反演中，一般来说只有非常少数的迭代次数需要多次的线搜索尝试，大部分迭代初始步长1即可满足公式(21)而直接进入下一次迭代.

需要注意的是，线搜索条件能保证步长满足当前迭代的总目标函数和总梯度的充分下降和曲率条件，但并非意味着单独考虑数据目标函数和数据目标函数梯度条件也满足，因此正则化因子取得很大时AR-LBFGS对海塞的近似会失败，此时需进入重启模式并减小正则化因子(比如减小为原来的1/4).

综上所述，AR-LBFGS算法的下降方向是相对不依赖前几次迭代正则化因子的，为了让算法稳定收敛，使用了基于Wolfe条件的线搜索策略，在两种情况下算法会进入重启模式：一种是公式(8)求解得到的下降方向线搜索失败，这表示LBFGS逼近海塞的误差太大，另一种是迭代完成后更新拟牛顿缓存矩阵失败，无法求取下一次迭代的拟牛顿下降方向.当然也可以引入Levenberg-Marquardt阻尼最小二乘法的信赖域思想，取消线搜索而加入对正则化因子的后验更新，但鉴于AR-LBFGS的研究目前仍然在起步阶段，这一思想仍未实现.

AR-LBFGS的整个算法流程如图 1所示.

2 算例

本文所有算例均在Inter XeonCPU@2.2 GHz，4处理器(64核)，内存为128 G的服务器上进行.文中设计的正演模型为三层模型，从地表往下，三层的电阻率分别为100 Ωm、10 Ωm和1000 Ωm，其中第一层厚度为14.618 km，中间层厚度为11.97 km；另外，在第一层中分布有A、B、C和D四个异常体(图 2)，图中各异常体用不同的灰度表示，电阻率分别为1 Ωm、1000 Ωm、10000 Ωm和10 Ωm，异常体详细数据如图 2所标识，图中白色部分表示第一层，电阻率为100 Ωm(图 2中蓝色字体所示)，需要说明的是图 2中四个异常体均按比例显示，G和F地层为示意图；NS和EW分别为北南和东西方向，图 2的Ⅱ中Z方向表示深度.在地表设计了“#”字型四条测线，每条线均匀分布有30个稀疏测点(图 3)，每个测点均包含48个频点，图 3表示正演模型不同切片图，切片具体位置如图中所述.网格模型剖分为42×42×61(南北向×东西向×垂向)，其中设置了7个空气层，模拟数据通过交错网格有限差分正演计算获得(Mackie et al., 1993; Siripunvaraporn et al., 2002; 谭捍东等, 2003).本节通过对不同的反演模型来测试算法的正确性和稳定性.

采用与正演相同的网格，反演初始模型分别给定为10 Ωm、100 Ωm、500 Ωm和1000 Ωm的四种均匀半空间，采用全张量反演，即Z_xx、Z_xy、Z_yx、Z_yy方向的阻抗以及倾子T_x和T_y的实部和虚部共12个分量.数据协方差矩阵C_d为主对角矩阵，其每个元素为该数据的方差，由于本文算例需要讨论正演合成的“无噪音”数据，在反演中将ZXY和ZYX的方差设定为其视电阻率相对误差为5%时的阻抗方差，ZXX和ZYY设定为ZXY和ZYX方差平均值的0.25%，这是根据假定ZXX和ZYY为0时，张量旋转任意角度时此两个元素能取到的最大幅值得到的，倾子的处理稍微复杂一些，首先设定其方差为幅值5%的平方和常数C(开始取1×10^-6)的最大值，在第一次正演完成后将C乘以一个系数重新调整方差，以保证初始的倾子RMS不大于阻抗RMS的一半.这种对无噪音数据的方差设定我们定义为门限方差，在反演实测的带噪音数据时，C_d的每个元素为门限方差与功率谱分析得到的方差的最大值.

使用自适应正则化策略，模型协方差矩阵为梯度算子内积，对每种模型进行100次反演迭代计算，其中迭代26次、60次和100次的结果见图 4、图 5和图 6，反演的部分参数及迭代计算时间如表 1所示，RMS、归一化的模型粗糙度和正则化因子随迭代次数的变化见图 7.

对于不同的初始模型，其反演迭代的时间是不一样的(表 1)，从图 4、图 5和图 6可以看出随着迭代次数的不断增加，测点覆盖区域及测点外区域的反演结果均越来越接近真实的模型.但也应该注意到在异常体的边缘，分辨率不高，这可能是由于测点比较稀疏的原因造成的，若想得到比较好的结果，最好测点能覆盖异常体的所有边缘，这样对异常体的约束会比稀疏测点要好，也更真实.在图 6中，经过100次迭代后，四种初始模型的反演结果在20 km以上的结构基本一致，但对于深部高阻层的反映是不同的，(B)、(C)、(D)三个结果其在深部有一部分高阻反映，(A)则没有，这可能由两个方面的原因造成：一方面，大地电磁响应本身对高阻相对不敏感，很难直接反演得到和真实模型电阻率一致的高阻背景；另一方面，AR-LBFGS虽然尽量压制了反演的初始模型依赖性，但对于数据约束不足的地区(我们的测点分布非常稀疏)，初始模型依赖性依然是存在的.

通过详细分析各模型100次迭代后的RMS，模型粗糙度和正则化因子λ曲线(图 7)，可得到如下的结论：(1)对于不同初始的模型，起始的正则化因子是不一样的，但是最终都会收敛到一个比较稳定的值，不过最终收敛的值越小，反演得到的结果越光滑，越接近真实的模型；(2)模型的粗糙度在迭代到一定的次数后，反演的结果也越来越接近真实的模型，模型粗糙度值接近收敛于0.1，不同的是各模型的收敛速度不一致；(3)当半空间模型与真实模型完全不同时，AR-LBFGS会寻找一个有较小RMS的最平滑模型；(4)通过综合比较四种模型以及三个参数可以看出，对于反演而言，并不是迭代次数越多越好，这些参数在迭代到一定的次数后，其值都逼近于某一个常数，因此选择合适的迭代次数既能减少反演的时间，也能得到比较满意的结果.

3 不同大地电磁参数作用下的讨论

上一节中反演用到了所有的数据，也就是全张量反演，即xx、xy、yx、yy方向的阻抗以及T_x和T_y倾子的实部和虚部共12个分量.在实际测量的数据中，有时候数据某一个方向的阻抗向量或者倾子向量会受到较大程度的干扰，此时如果将受干扰的数据也进行反演，则很难得到理想的结果.另外，当全张量数据都参与反演时，内存需求极大，反演所需要的时间也会长于只用部分数据参数的反演(表 2).当然，在数据质量良好的情况下，参与反演的参数越多，对模型的约束就越大，得到的结果也会相对而言更精确.不过一个好的算法，如果仅对数据中的部分参数进行反演，也能得到相对一致的结果，这就说明算法的适应性较强，算法也比较稳定.为此，本节将对100 Ωm的均匀半空间模型，分别对三种不同参数组合下的数据进行反演：(1)全张量反演，包括Z_xx, Z_xy, Z_yx, Z_yy, Z_zx, Z_zy的实部和虚部共12个分量；(2)不包括倾子的阻抗张量数据反演，仅包括Z_xx, Z_xy, Z_yx, Z_yy的实部和虚部共8个分量；(3)仅有xy、yx方向(副对角方向)的阻抗张量的数据反演，即Z_xy, Z_yx的实部和虚部共4个分量.得到的结果见图 8—10.

通过图 8—10可以看出，即使仅有TE和TM数据，反演得到的结果也能很好的反映实际模型，因此在数据无干扰的情况下，使用较少的参数就能得到理想的结果.但是不同参数下得到的结果也有差别，比如对1 Ωm低阻异常体而言，全张量数据反演得到结果更准确，边界也更清晰；分析其原因为参与反演的数据越多越能相对降低反问题本身的不确定性，并且倾子数据对电阻率突变边界的分辨能力比阻抗数据更强.但带倾子反演的时候，数据误差和反演的门限误差的设定是非常值得讨论的问题，由于篇幅所限，本论文暂未讨论人工增加随机噪音的反演，我们初步的实验表明全张量数据反演中，对倾子门限误差的控制需保证初始阻抗RMS大致为初始倾子RMS的两倍，这样能让反演在开始阶段尽量拟合阻抗而非倾子，以避免得到一些看起来分辨率高的假象，这些研究将随后尽快发表.

4 结论

本文所提出的自适应正则化有限内存拟牛顿反演算法，在每次迭代中，只基于数据目标函数梯度本身(而非总目标函数梯度)来逼近数据拟合项的海塞矩阵，然后通过求解类高斯牛顿下降方向方程，得到了不依赖前几次迭代正则化因子的更精确下降方向.选取一个能保证求解下降方向方程右端项的二范数模中关于数据梯度和模型梯度两项的二范数处于同一个数量级的二范数因子，从而得到一个能使反演迭代更加稳定的自适应正则化因子.

AR-LBFGS反演算法可以选择不同的初始模型，且算法有很好的稳定性，反演中数据偏移不会下降太快，而正则化因子在反演迭代的后期需要缓慢下降，这样才会使粗糙度缓慢增长.利用不同的数据参数进行反演比较，可以看出，使用的参数越多，反演得到的结果越精确；而在数据无干扰的情况下，较少的参数就能得到可以接受的结果.因此数据采集应尽可能的精确，减少干扰，这样利用较少的数据参数即可得到相对准确的地下电性结构，节省了计算机的内存和时间需求，不需要大型工作站，利于三维反演的全面应用.