位场向下延拓系数矩阵性质及Barzilai-Borwein向下延拓法

张志厚 ，廖晓龙 ，姚 禹 ，范祥泰 ，路润琪

（西南交通大学地球科学与环境工程学院，四川 成都 611756）

重磁位场向下延拓技术可以将海平面磁测数据大深度向下延拓至海底，将航空重磁数据延拓至地形线，不仅突出了有用信息，还可以弥补条件恶劣地区重磁场资料的不足，并且重磁位场的向下延拓运算是许多重磁数据处理的核心步骤[1-5]. 此外，重磁位场的向下延拓技术在军事领域具有非常重要的作用[6-9].

重磁位场向下延拓的主要方法从类型上大体可以分为泰勒级数法[10-13]、等效源法[14-20]、积分迭代法[1-4,21-23]、Tikhonov正则化法[24-29]、Landweber正则化法[26,28,30]、奇异值分解法[26]等.

泰勒级数法的主要缺点是更高阶次的导数计算会导致数据的高频成分被无限放大，曾小牛等[31]对求导的频谱函数进行了改进，但阶次越多误差也被叠加得越多，最终导致下延误差反而很大；等效源法的主要缺点是计算时间长，且延拓效果受到“源”的类型、个数、深度及其组合的影响.

位场向下延拓的积分迭代法[21,32]计算效果好，诸多学者研究了其鲁棒性等性能[33-36]，该方法在数学上也被称为对称核第一类Fredholm积分方程的逐次逼近解法[37]. Tikhonov正则化法相比Landweber迭代法具有最好的综合下延性能[26,28]. 然而Tikhonov正则化法存在正则化参数的选取问题.

空间域的位场向下延拓实质上是求解一个第一类Fredholm型积分方程，该方程最突出的特性是其“不适定”性. 向下延拓方程离散化后的系数矩阵为实对称矩阵. 本文首先证明了该系数矩阵为对称的双重Toeplitz系统矩阵 （blak-Toeplitz-Toeplitzblock, BTTB）；其次假定在其正定的条件下，采用Barzilai-Borwein （BB）法[38-39]求解该线性方程组，并对迭代步长进行限制，保证了假定系数矩阵为正定情况下算法的收敛性，实现了位场的向下延拓计算，在迭代计算中使用BTTB与向量乘积的快速算法[40]对大数据进行计算，一定程度上节约了计算成本，并对该方法进行了理论模型无噪声数据与实际资料的检验，以及对比分析了该方法的收敛速度与位场向下延拓的积分迭代法的收敛速度.

1 空间域位场向下延拓

假设位场观测面为笛卡儿直角坐标系o-xyz的水平面 （z=0 ），该平面上的位场为f0(x,y)，计算z=h（h> 0）平面的位场fh(x,y)，且该平面为观测面以下至场源以上空间（z轴向下）.

根据向上延拓公式，可得

表示成褶积的形式为

2 位场延拓系数矩阵的性质

陈龙伟等[41]较为详细地证明了位场向下延拓系数矩阵为对称的分块Toeplitz矩阵，本文提供一种简洁的证明方法. 已知f0(x,y) 通过式（1）求解fh(x,y)即为向下延拓，由于已知场与待求场皆为离散化数据，且f0(x,y) 与fh(x,y)离散化网格规格一致，因此将

式（1）离散化[42]为

式中：Δx、Δy分别为方向x、y的离散化网格间距；M、N分别为方向x、y的离散化网格点总数，m=1,2,3,···,M，n=1,2,3,···,N.

由于离散化网格间距固定，为了方便表示，将f0(mΔx,nΔy) 与fh(iΔx,jΔy) 分 别 写 成f0(m,n)与fh(i,j) ，将式（4）中的f0和fh重排为列向量，分别记为

式（4）用矩阵表示为

式中：A为位场向下延拓的系数矩阵.

根据式（4），位场向下延拓的系数为

借鉴“格架函数”概念[43]，发现Q也具有平移等效性、互换对称性，等价关系为

式（8）左右两端分别代入式（7）即可得证. 根据式（8），4位数组Q(m,n,i,j)的前两个参量都是常数1，故其可表示为二维数组，即将Q(m,n,i,j)简记为Q|m−i|+1,|n−j|+1.Q为四维数组，以行列式的形式排列即可得到：

可以看出：A为对称矩阵，大小为MN×MN，内部包含了N2个分块的对称矩阵，大小为M×M，即矩阵A为双重Toeplitz系统矩阵.

3 BB法

BB法是一种特殊的最速下降法[44]，假设A正定，式（7）采用BB法求解，位场向下延拓BB法的迭代格式为

式中：k为迭代次数为第k次向下延拓的计算值；t(k)为第k次的迭代步长；p(k)为第k次迭代待延拓场向上延拓值与已知场的差值向量

迭代停机准则：当 (p(k))Tp(k)<τ，即Cτ 时可停止计算，其中 τ为给定的适当小正数，C为常数，此时

若式（11）中的t(k)为常数时，上述迭代方法即为位场向下延拓的积分迭代法. 而向下延拓BB法相当于变步长的积分迭代法；王顺杰等[36]采用不动点原理证明了位场向下延拓积分迭代法的收敛性，但迭代步长的取值范围在0～2. 因此在求解式（4）时，将其迭代步长约束在0～2范围内可确保向下延拓BB法的鲁棒性. 考虑到迭代过程中必然会存在计算数据量大的问题，采用对称Toeplitz矩阵与向量相乘的快速算法[40]对式（11）中的和进行计算，有效地提高计算效率.

4 模型试验

4.1 二维模型

二维模型为无限延伸的水平圆柱体，正演计算两个不同高度观测线的理论重力异常，其中观测线与模型体垂直，对计算结果和收敛速度进行分析.

圆柱体模型的半径为500 m，线密度 λ=1000kg/m3，模型的中心埋藏深度为3 km. 计算点距为100 m，计算了z=0 和z=2km的理论重力异常. 将z=0理论异常分别采用BB法与积分迭代法向下延拓至z=2km. 理论重力异常与不同方法向下延拓后的结果如图1所示.

z=2km 理论重力异常的最大值与最小值分别为10.477 × 10−5、0.064 × 10−5m/s2；z=0理论重力异常的最大值与最小值分别为3.492 × 10−5、0.182 ×10−5m/s2. 结果表明：1） 迭代10次BB法计算结果的最大、最小值分别为8.645 × 10−5、−0.015 × 10−5m/s2，积分迭代法计算结果的最大、最小值分别为7.715 ×10−5、−0.164 × 10−6m/s2；2） 迭代50次BB法计算结果的最大、最小值分别为9.542 × 10−5、−0.173 × 10−6m/s2，积分迭代法计算结果的最大、最小值分别为9.132 ×10−5、−0.171 × 10−6m/s2. 由此可见，相同迭代次数下，BB法的计算结果与理论重力异常更为接近.

为了定量地评价本文所提方法的计算精度、稳定性，采用不同迭代次数计算结果的均方误差来评估其效果，均方误差为

式中：fcon,l和ftheo,l分别为第l个计算点的延拓值和理论值；L为计算的点数.

图2为位场向下延拓的BB法与积分迭代法在不同迭代次数下计算结果的均方误差. 可见两种方法都具有较好的收敛性，但向下延拓BB法的收敛速度更快.

图1 理论重力异常与不同方法向下延拓结果对比Fig. 1 Theoretical gravity anomaly compared with downward continuation results of different methods

图2 重力模型计算结果的均方误差Fig. 2 Mean square error calculated by gravity model

4.2 三维模型

三维模型为组合圆球型磁异常体，通常被用于位场向下延拓效果的检验[45-47]. 首先通过正演计算两个观测面的磁异常，然后采用不同计算方法将距离场源较远处的平面磁异常向下延拓至另一平面，通过均方误差来评估计算效果. 三维组合圆球形磁异常体参数见表1. 表中：(x0,y0,h0)为圆球型异常体中心坐标；r0为圆球异常体半径；M0异常体磁化强度；I0和A0分别为地磁场倾角和剖面磁偏角；I′和A′分别为异常体磁化强度的倾角和磁北夹角. 根据圆球型磁异常体的正演计算解析式[48]，分别计算两个不同观测面上的磁异常，如图3所示，其中虚拟观测区网格点个数为401× 401，Δx与 Δy都为50 m.

表1 三维组合圆球型磁异常体参数Tab. 1 Magnetic anomaly parameters of 3D combination sphere

图3 不同高度面上的磁异常等值线（单位：nT）Fig. 3 Magnetic anomaly isolines at different planes (unit： nT)

将位于z= 0高度的理论磁异常分别采用BB法以及积分迭代法下延至z= 1 km的平面，计算结果见图4所示. 其结果除了边部的零等值线略有偏差，其余部分都与理论场（图3（b））高度吻合.

图4 理论磁异常向下延拓后的计算结果（单位：nT）Fig. 4 Calculation results of theoretical gravity anomaly after downward continuation (unit： nT)

两种方法计算结果的均方误差如图5所示，可以看出：随着迭代次数的增加，均方误差都在急剧衰减，BB法相比积分迭代法的收敛速度更快；在相同迭代次数的情况下，位场向下延拓BB法的计算结果的精度高于积分迭代法；同一收敛精度下，当均方误差在0.4 nT以下时，BB法需要迭代次数15次，而积分迭代法需要迭代35次以上，可见向下延拓BB法的计算速度提高了2倍.

图5 三维磁场模型计算结果的均方误差Fig. 5 Mean square error calculated by3D magnetic field model

此外，随着迭代次数增加，位场向下延拓BB法与积分迭代法的计算结果精度趋近，两种方法的本质都是迭代法，迭代法求解反问题时都会出现所谓的“半收敛”现象[44]，即在迭代的早期阶段，估计解可以得到很好的改善，但其长期行为会出现“饱和”效应，甚至会走向无序[49-50]，故而两种方法的计算精度极有可能会越来越相近.

5 计算速度

在主频2.30 GHz、内存1024 MB的计算机上，1024 × 1024个数据向下延拓迭代50次，BB法计算时间约为130 s. 另外，一般情况下依据经验取一定的迭代次数作为停机准则，BB法相比积分迭代法只需要较少的迭代次数便可达到相同的计算精度.

6 实际资料检验

通过某铜铁矿区实测网格化磁异常资料来评估下延计算效果，如图6所示. 图6（a）为已知实测的网格化磁异常，点距为50 m. 将该异常向上延拓10倍点距（向上延拓方法采用快速傅里叶变换方法），其结果如图6（b）所示. 随后再将该高度上的异常分别采用BB法与积分迭代法向下延拓500 m至原始高度，其结果如图6（c）、（d）所示（两种方法迭代次数相同），计算结果都与图6（a）基本一致.

图6 实例验证（单位：nT）Fig. 6 Validation by measured data (unit： nT)

原始场的最大、最小值分别为1518.3、−669.2 nT.对应的，向下延拓BB法计算结果的最大、最小值分别为1507.7、−640.7 nT，积分迭代法计算结果的最大、最小值分别为1497.3、−623.5 nT. 用式（11）计算BB法和积分迭代法下延结果的均方误差分别为24.4 nT和29.7 nT. 相比原始场的最大、最小值，两种方法的均方误差都较小，但BB法的精度更高.

为了进一步评价两种方法的计算精度，采用平均相对误差 ε进行评估，如式（13）所示.

用式（13）计算本文所提BB法和积分迭代法的平均相对误差分别为6.1%、7.7%，可见在实际应用中本文所提方法的效果优于积分迭代法.

7 结 论

1） 本文首先证明了位场向下延拓的系数矩阵为对称的双重Toeplitz系统矩阵，该性质可为今后的位场数值模拟节省成本.

2） 实现了一种位场向下延拓的计算方法，即位场向下延拓的BB法，该方法主要优点的计算精度高、速度快. 相比普遍采用的Tikhonov方法，文中方法不存在正则化因子的选择问题.