CN107748834A

CN107748834A - 一种计算起伏观测面磁场的快速、高精度数值模拟方法

Info

Publication number: CN107748834A
Application number: CN201711170885.XA
Authority: CN
Inventors: 李昆; 戴世坤; 陈轻蕊; 陈龙伟; 张钱江; 赵东东
Original assignee: Central South University
Current assignee: Central South University
Priority date: 2017-11-22
Filing date: 2017-11-22
Publication date: 2018-03-02
Anticipated expiration: 2037-11-22
Also published as: CN107748834B

Abstract

本发明提供一种计算起伏观测面磁场的快速、高精度数值模拟方法，其通过复杂磁性体模型表示、高斯参数设计、离散偏移波数计算、磁化强度计算、波数域加权系数计算、二维离散傅里叶反变换等步骤，实现了起伏观测面磁场数值模拟在效率和精度上的统一。本发明解决了目前起伏观测面磁场数值模拟方法计算精度低、计算时间长，无法满足大规模航空磁测数据精细反演成像的问题，有助于开展大规模航空磁测数据三维磁化率精细反演成像、人机交互建模和解释的研究。

Description

一种计算起伏观测面磁场的快速、高精度数值模拟方法

技术领域

本发明涉及一种面向航空磁法勘探的数值模拟方法，特别是任意磁化率分布、起伏观测面磁场快速、高精度数值模拟方法。

背景技术

揭示地球内部物性分布和地质构造情况是地球物理学研究的基础科学问题。近年来，随着传感器技术和航空测量技术的发展，对地球磁场的测量范围越来越大，空间分辨率越来越高，获取的数据量随之剧增，这为揭示地球内部精细结构提供了数据基础。针对海量磁测数据的反演成像对传统反演成像方法提出了巨大挑战。高效、高精度、高分辨率的反演成像方法成为航空磁法勘探领域的研究重点和难点。正演是反演成像的基础，正演的计算精度与效率决定着反演精度与效率的高低。

针对复杂磁化率分布条件下磁场数值模拟，众多国内外学者进行了研究。文献(姚长利,郝天珧,管志宁,张聿文.重磁遗传算法三维反演中高速计算及有效存储方法技术.地球物理学报,2003.46(2):252-258.)提出了“格架分离”技术和“格架等效计算方案”，较好解决了计算效率和计算精度问题。但对于大规模剖分情形，该文献所给出的数值模拟方法的计算效率仍然比较低。文献(Tontini,F.C.,L.Cocchi,C.Carmisciano.Rapid 3-Dforward model of potential fields with application to the Palinuro Seamountmagnetic anomaly(southern Tyrrhenian Sea,Italy).Journal of GeophysicalResearch,2009.114.)采用三维傅里叶变换，给出了任意磁化率分布情形下磁场的波数域表达式，借助三维快速傅里叶变换算法，实现了快速数值模拟，该方法效率极高，但由于存在截断效应，该方法数值模拟精度低。文献(Wu,L.,Tian,G.High-precision Fourierforward modeling of potential fields.Geophysics,2014,79(5):G59-G68.)提出了一种重磁数值模拟的Gauss-FFT方法，该方法有效克服了传统傅里叶变换方法的截断效应问题，提高了数值模拟精度，但计算效率降低。

已有的磁场数值模拟方法，实现的是某一水平观测面磁场的计算。而实际的磁测，尤其是航空磁测，观测面都是起伏曲面，目前对于起伏观测面磁场快速、高精度数值模拟方法研究较少。因此寻找一种计算效率与计算精度高的起伏界面磁场的数值模拟方法，对于实现航空磁测数据三维反演成像具有重要的现实意义。

发明内容

本发明的目的在于提供一种计算起伏观测面磁场的快速、高精度数值模拟方法，该发明解决了目前起伏观测面磁场数值模拟方法计算精度低、计算时间长，无法满足大规模航空磁测数据精细反演成像的问题，有助于开展大规模航空磁测数据三维磁化率精细反演成像、人机交互建模和解释的研究。

本发明的技术方案是：

一种计算起伏观测面磁场的快速、高精度数值模拟方法，包括以下步骤：

步骤一：复杂磁性体模型表示：

确定目标区域以及目标区域内的异常体，建立包含所有目标区域的长方体模型，使得包含起伏地形的目标区域完全嵌入在该长方体模型中；将该长方体模型剖分成若干个小长方体，长方体模型的x、y、z方向小长方体的剖分个数分别为N_x、N_y、N_z；各小长方体的x方向边长均相同，各小长方体的y方向的边长均相同，在长方体模型的z方向上同一层的各小长方体z方向边长相同；

根据每个小长方体内嵌入的目标区域所对应的磁化率的分布，对每个小长方体磁化率进行相应的赋值，每个小长方体磁化率为常值，不同小长方体的磁化率取值不同，以此刻画任意磁化率分布复杂磁性体模型；将位于空气部分的小长方体的磁化率值设为零，以此刻画起伏地形。

步骤二：插值平面高度确定：

给定起伏观测面，在起伏观测面的最高点和最低点之间，确定N_obs个插值平面，其高度为Z_l，l＝1,2,…,N_obs。本发明所设的插值平面个数越多，最终计算得到的起伏观测面磁场的精度越大。在起伏观测面的最高点和最低点之间确定N_obs个处于不同高度的插值平面。进一步地，N_obs个处于不同高度的插值平面是均匀等间距设置的，均匀等间距设置的处理效果会让整个起伏观测面的计算精度更加平均。

步骤三：高斯参数设计：

给定x、y方向的高斯点个数L_x、L_y，区间[-1,1]上高斯点t_a、t_b及高斯系数A_a、A_b；其中，a＝1,2,…,L_x，b＝1,2,…,L_y。经测试，x、y方向的高斯点个数选取3或者4时即能够满足计算精度需求。t_a，t_b，A_a，A_b有确定的数值对应，可以查表得到，属于高斯数值积分常识性知识。

步骤四：根据空间域剖分参数和高斯参数，计算离散偏移波数，具体过程如下：

式中，

a＝1,2,…,L_x

b＝1,2,…,L_y

k_x，k_y分别表示x、y方向的偏移波数，Δk_x，Δk_y分别表示x、y方向基波数，N_x，N_y分别表示长方体模型其x、y方向小长方体的剖分个数，Δx，Δy分别表示小长方体其x、y方向的边长。

步骤五：磁化强度计算：

首先，根据公知的地球主磁场模型IGRF，计算各小长方体中心位置的地球主磁场三分量T_x(x_i,y_j,z_k),T_y(x_i,y_j,z_k),T_z(x_i,y_j,z_k)；其中：T_x(x_i,y_j,z_k)、T_y(x_i,y_j,z_k)、T_z(x_i,y_j,z_k)分别表示(x_i,y_j,z_k)处地球主磁场的x、y、z分量；(x_i,y_j,z_k)表示编号为(i,j,k)的小立方体几何中心坐标，i＝1,2,…,N_x(说明：此处i为斜体，斜体形式的i＝1,2,…,N_x)，j＝1,2,…,N_y，k＝1,2,…,N_z；N_z表示长方体模型其z方向小长方体的剖分个数。

其次，根据磁化率分布和地球主磁场，计算磁化强度

M_x(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_x(x_i,y_j,z_k) (3)

M_y(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_y(x_i,y_j,z_k) (4)

M_z(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_z(x_i,y_j,z_k) (5)

式中，χ(x_i,y_j,z_k)表示编号为(i,j,k)的小立方体的磁化率值；M_x(x_i,y_j,z_k)、M_y(x_i,y_j,z_k)、M_z(x_i,y_j,z_k)分别表示(x_i,y_j,z_k)处磁化强度的x、y、z分量。

步骤六：采用二维快速傅里叶变换算法，快速计算波数域磁化强度三分量：

其中，为z方向第k层波数域磁化强度的x、y、z分量，(k_x,k_y,z_k)表示磁化强度在波数域的坐标，k＝1,2,…,N_z，i表示虚数单位(说明：式中e的上标正体i表示虚数单位)。

步骤七：计算N_obs个插值平面各自对应的波数域加权系数

式中，μ表示真空磁导率，Δz_k表示长方体模型其第k层小长方体的z方向边长，表示x,y,z三个方向波数域加权系数，ζ表示异常体z方向的长度范围；sgn()为符号函数：

步骤八：根据波数域磁化强度和波数域加权系数，计算N_obs个插值平面上各自对应的波数域磁场

式中，表示高度为Z_l的插值平面上的波数域磁场。

步骤九：采用二维快速傅里叶反变换算法，快速计算N_obs个插值平面上各自对应的磁场

其中，A_a，A_b为高斯加权系数。

步骤十：计算起伏观测面上的磁场

采用三次样条插值算法，根据步骤九计算得到的N_obs个插值平面上各自对应的磁场，插值得到起伏观测面磁场。

本发明的步骤七中计算N_obs个插值平面各自对应的波数域加权系数，其过程为

(1)计算

将计算结果存储于计算机。

(2)计算

将计算结果存储于计算机。

(3)计算

将计算结果存储于计算机。

(4)采用查询的方式，根据存储的E_k、V_l和V，计算N_obs个插值平面对应的波数域加权系数

本发明是一个有机整体，包含三个关键环节：第一个是在特殊的模型剖分方式下，给出了精确的棱柱体组合模型磁场波数域计算公式；第二个是多个插值平面情形下波数域加权系数快速算法；第三个是单个插值平面磁场快速、高精度正演计算。基于上述基础，本发明实现了任意磁化率分布、起伏观测面磁场数值模拟在计算效率和计算精度上的统一。

与现有技术相比，本发明具有以下优点：

(1)模型剖分方法简单、灵活，很容易刻画任意磁化率分布复杂磁性体以及起伏观测面；

(2)能够实现起伏观测面磁场的快速、高精度数值模拟，可以满足大规模航空磁测数据三维精细反演成像、人机交互建模和解释的需求；

(3)大规模数值模拟时，占用计算机内存少，计算效率和计算精度高。

附图说明

图1为本发明的流程图；

图2为复杂磁性体模型表示；

图3为长方体磁性体模型平面图；

图4为起伏观测面磁场理论值

图5为起伏观测面磁场计算值

图6为计算值与理论值对比；

图7为随插值平面个数变化起伏观测面相对均方根误差图；

图8为随插值平面个数变化起伏观测面计算时间趋势图

图中符号说明如下：

N_z：表示长方体模型其z方向小长方体的剖分个数；

N_x：表示长方体模型其x方向小长方体的剖分个数；

N_y：表示长方体模型其y方向小长方体的剖分个数；

rrms：表示相对均方根误差；

表示N_obs个插值平面分别对应的插值平面高度。

本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

如图1所示，图1为本发明一种计算起伏观测面磁场的快速、高精度数值模拟方法的流程图。本发明包括以下步骤：

步骤一：复杂磁性体模型表示：

如图2所示，确定目标区域，建立包含所有目标区域的长方体模型，使得包含起伏地形的目标区域完全嵌入在该长方体模型中；将该长方体模型剖分成若干个小长方体，长方体模型的x、y、z方向小长方体的剖分个数分别为N_x、N_y、N_z；各小长方体的x方向边长均相同，各小长方体的y方向的边长均相同，在长方体模型的z方向上同一层的各小长方体z方向边长相同。本发明对长方体模型的z方向上不同层的各小长方体z方向边长的长度是否一致没有要求，可以相同也可以不同。

步骤二：插值平面高度确定：

如图3所示，给定起伏观测面，在起伏观测面的最高点和最低点之间，均匀等间距确定N_obs个处于不同高度的插值平面，其高度为Z_l，l＝1,2,…,N_obs。

步骤三：高斯参数设计：

式中，

a＝1,2,…,L_x

b＝1,2,…,L_y

步骤五：磁化强度计算：

首先，根据公知的地球主磁场模型IGRF，计算各小长方体中心位置的地球主磁场三分量T_x(x_i,y_j,z_k)，T_y(x_i,y_j,z_k)，T_z(x_i,y_j,z_k)；其中：T_x(x_i,y_j,z_k)、T_y(x_i,y_j,z_k)、T_z(x_i,y_j,z_k)分别表示(x_i,y_j,z_k)处地球主磁场的x、y、z分量；(x_i,y_j,z_k)表示编号为(i,j,k)的小立方体几何中心坐标，i＝1,2,…,N_x，j＝1,2,…,N_y，k＝1,2,…,N_z；N_z表示长方体模型其z方向小长方体的剖分个数。

其次，根据磁化率分布和地球主磁场，计算磁化强度

M_x(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_x(x_i,y_j,z_k) (24)

M_y(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_y(x_i,y_j,z_k) (25)

M_z(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_z(x_i,y_j,z_k) (26)

其中，为z方向第k层波数域磁化强度的x、y、z分量，(k_x,k_y,z_k)表示磁化强度在波数域的坐标，k＝1,2,…,N_z，i表示虚数单位。

步骤七：计算N_obs个插值平面各自对应的波数域加权系数

式中，μ表示真空磁导率，Δz_k表示长方体模型其第k层小长方体的z方向边长，表示x,y,z三个方向波数域加权系数，ζ表示异常体z方向的长度范围，sgn()为符号函数

式中，表示高度为Z_l的插值平面上的波数域磁场。

其中，A_a，A_b为高斯加权系数。

步骤十：计算起伏观测面上的磁场

(1)计算

将计算结果存储于计算机。

(2)计算

将计算结果存储于计算机。

(3)计算

将计算结果存储于计算机。

下面给出一具体实施例，验证本发明所提出的方法用于计算任意磁化率分布起伏观测面磁场的效率和精度。

目标区域有一个异常体，该异常体为一个棱柱磁性体，目标区域范围为：x方向从-1000m到1000m，y方向从-1000m到1000m，z方向从0m到640m(z轴向下为正)；异常体展布范围为：x方向从-500m到500m，y方向从-500m到500m，z方向从100m到600m；磁化率为0.01(SI)；目标区域地球主磁场为45000nT，磁偏角为5度，磁倾角为45度。将目标区域剖分成200×200×64个大小相同的小棱柱体。计算起伏地形的高度范围为(-500m,-300m)。

本发明方法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU-Inter Core i5-4590，主频为3.3GHz，内存为16GB。当插值平面个数为15时，磁场模拟算法理论值如图4所示，计算值如图5所示，从形态上看两者是一致的。图6所示为y＝0m时的计算值与理论值，两者吻合很好。为了进一步研究插值平面个数对计算精度的影响，图7所示为随着插值平面个数的变化整个起伏观测面的相对均方根误差，图8所示为随着插值平面个数的增大计算时间趋势图，从图7、图8可以看出随着插值平面个数的增大，计算精度逐渐增大，计算时间增加缓慢，说明该算法在模拟起伏地形磁场时具有较高的计算精度与效率。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims

1.一种计算起伏观测面磁场的快速、高精度数值模拟方法，其特征在于，包括以下步骤：

步骤一：复杂磁性体模型表示：

根据每个小长方体内嵌入的目标区域所对应的磁化率的分布，对每个小长方体磁化率进行相应的赋值，每个小长方体磁化率为常值，不同小长方体的磁化率取值不同，以此刻画任意磁化率分布复杂磁性体模型；将位于空气部分的小长方体的磁化率值设为零，以此刻画起伏地形；

步骤二：插值平面高度确定：

给定起伏观测面，在起伏观测面的最高点和最低点之间，确定N_obs个插值平面，其高度为Z_l，l＝1,2,…,N_obs；

步骤三：高斯参数设计：

给定x、y方向的高斯点个数L_x、L_y，区间[-1,1]上高斯点t_a、t_b及高斯系数A_a、A_b；其中，a＝1,2,…,L_x，b＝1,2,…,L_y；

<mrow> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>=</mo> <mo>(</mo> <mrow> <mi>m</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <mrow> <msub> <mi>t</mi> <mi>a</mi> </msub> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <msub> <mi>&Delta;k</mi> <mi>x</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>=</mo> <mo>(</mo> <mrow> <mi>n</mi> <mo>+</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mrow> <mo>(</mo> <mrow> <msub> <mi>t</mi> <mi>b</mi> </msub> <mo>+</mo> <mn>1</mn> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> <msub> <mi>&Delta;k</mi> <mi>y</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

式中，

<mrow> <msub> <mi>&Delta;k</mi> <mi>x</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> <mrow> <msub> <mi>N</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mi>&Delta;</mi> <mi>x</mi> </mrow> </mfrac> </mrow>

<mrow> <msub> <mi>&Delta;k</mi> <mi>y</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mn>2</mn> <mi>&pi;</mi> </mrow> <mrow> <msub> <mi>N</mi> <mi>y</mi> </msub> <mo>&CenterDot;</mo> <mi>&Delta;</mi> <mi>y</mi> </mrow> </mfrac> </mrow>

a＝1,2,…,L_x

b＝1,2,…,L_y

k_x，k_y分别表示x、y方向的偏移波数，Δk_x，Δk_y分别表示x、y方向基波数，N_x，N_y分别表示长方体模型其x、y方向小长方体的剖分个数，Δx，Δy分别表示小长方体其x、y方向的边长；

步骤五：磁化强度计算：

首先，根据公知的地球主磁场模型IGRF，计算各小长方体中心位置的地球主磁场三分量T_x(x_i,y_j,z_k),T_y(x_i,y_j,z_k),T_z(x_i,y_j,z_k)；其中：T_x(x_i,y_j,z_k)、T_y(x_i,y_j,z_k)、T_z(x_i,y_j,z_k)分别表示(x_i,y_j,z_k)处地球主磁场的x、y、z分量；(x_i,y_j,z_k)表示编号为(i,j,k)的小立方体几何中心坐标，i＝1,2,…,N_x，j＝1,2,…,N_y，k＝1,2,…,N_z；N_z表示长方体模型其z方向小长方体的剖分个数；

其次，根据磁化率分布和地球主磁场，计算磁化强度

M_x(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_x(x_i,y_j,z_k) (3)

M_y(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_y(x_i,y_j,z_k) (4)

M_z(x_i,y_j,z_k)＝χ(x_i,y_j,z_k)T_z(x_i,y_j,z_k) (5)

式中，χ(x_i,y_j,z_k)表示编号为(i,j,k)的小立方体的磁化率值；M_x(x_i,y_j,z_k)、M_y(x_i,y_j,z_k)、M_z(x_i,y_j,z_k)分别表示(x_i,y_j,z_k)处磁化强度的x、y、z分量；

<mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>y</mi> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>x</mi> </msub> </munderover> <msub> <mi>M</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>y</mi> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>x</mi> </msub> </munderover> <msub> <mi>M</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>y</mi> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>x</mi> </msub> </munderover> <msub> <mi>M</mi> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>

其中，为z方向第k层波数域磁化强度的x、y、z分量，(k_x,k_y,z_k)表示磁化强度在波数域的坐标，k＝1,2,…,N_z，i表示虚数单位；

步骤七：计算N_obs个插值平面各自对应的波数域加权系数

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>&mu;</mi> <mrow> <mi>i</mi> <mn>2</mn> <msub> <mi>k</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>&lsqb;</mo> <msubsup> <mo>&Integral;</mo> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>+</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> </msubsup> <mrow> <mo>(</mo> <mo>-</mo> <mi>s</mi> <mi>g</mi> <mi>n</mi> <mo>(</mo> <mrow> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>-</mo> <mi>&zeta;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>K</mi> <mo>|</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>-</mo> <mi>&zeta;</mi> <mo>|</mo> </mrow> </msup> <mi>d</mi> <mi>&zeta;</mi> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>&mu;</mi> <mrow> <mi>i</mi> <mn>2</mn> <msub> <mi>k</mi> <mi>x</mi> </msub> </mrow> </mfrac> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>&lsqb;</mo> <msubsup> <mo>&Integral;</mo> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>+</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> </msubsup> <mrow> <mo>(</mo> <mo>-</mo> <mi>s</mi> <mi>g</mi> <mi>n</mi> <mo>(</mo> <mrow> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>-</mo> <mi>&zeta;</mi> </mrow> <mo>)</mo> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>K</mi> <mo>|</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>-</mo> <mi>&zeta;</mi> <mo>|</mo> </mrow> </msup> <mi>d</mi> <mi>&zeta;</mi> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>&mu;</mi> <mi>K</mi> </mrow> <mrow> <mn>2</mn> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>k</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>(</mo> <mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> </mrow> <mo>)</mo> <mo>&lsqb;</mo> <mrow> <msubsup> <mo>&Integral;</mo> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> <mrow> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>+</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> </mrow> </msubsup> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>K</mi> <mo>|</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>-</mo> <mi>&zeta;</mi> <mo>|</mo> </mrow> </msup> <mi>d</mi> <mi>&zeta;</mi> </mrow> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>

步骤八：根据波数域磁化强度和波数域加权系数，计算N_obs个插值平面上各自对应的波数域磁场；

<mrow> <msub> <mover> <mi>B</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>z</mi> </msub> </munderover> <mo>&lsqb;</mo> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msub> <mover> <mi>M</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>

式中，表示高度为Z_l的插值平面上的波数域磁场；

<mrow> <msub> <mi>B</mi> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <mn>4</mn> <msup> <mi>&pi;</mi> <mn>2</mn> </msup> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>a</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>L</mi> <mi>x</mi> </msub> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>b</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>L</mi> <mi>y</mi> </msub> </munderover> <msub> <mi>A</mi> <mi>a</mi> </msub> <msub> <mi>A</mi> <mi>b</mi> </msub> <mo>&lsqb;</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>N</mi> <mi>x</mi> </msub> <mo>/</mo> <mn>2</mn> </mrow> <mrow> <msub> <mi>N</mi> <mi>x</mi> </msub> <mo>/</mo> <mn>2</mn> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>N</mi> <mi>y</mi> </msub> <mo>/</mo> <mn>2</mn> </mrow> <mrow> <msub> <mi>N</mi> <mi>y</mi> </msub> <mo>/</mo> <mn>2</mn> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mover> <mi>B</mi> <mo>~</mo> </mover> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>+</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

其中，A_a，A_b为高斯加权系数；

步骤十：计算起伏观测面上的磁场

2.根据权利要求1所述的计算起伏观测面磁场的快速、高精度数值模拟方法，其特征在于，步骤二中，在起伏观测面的最高点和最低点之间确定N_obs个处于不同高度的插值平面，所设插值平面个数越多，最终计算得到的起伏观测面磁场的精度越大。

3.根据权利要求1或2所述的计算起伏观测面磁场的快速、高精度数值模拟方法，其特征在于，步骤二中，在起伏观测面的最高点和最低点之间均匀等距确定N_obs个处于不同高度的插值平面。

4.根据权利要求3所述的计算起伏观测面磁场的快速、高精度数值模拟方法，其特征在于，步骤七中计算N_obs个插值平面各自对应的波数域加权系数，其过程为

(1)计算

<mrow> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>,</mo> <mi>k</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>N</mi> <mi>z</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mi>E</mi> <mrow> <msub> <mi>N</mi> <mi>z</mi> </msub> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>K</mi> <mrow> <mo>(</mo> <msub> <mi>z</mi> <msub> <mi>N</mi> <mi>z</mi> </msub> </msub> <mo>+</mo> <mfrac> <mrow> <msub> <mi>&Delta;z</mi> <msub> <mi>N</mi> <mi>z</mi> </msub> </msub> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>

将计算结果存储于计算机；

(2)计算

将计算结果存储于计算机；

(3)计算

<mrow> <mi>V</mi> <mo>=</mo> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>x</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>x</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>ik</mi> <mi>y</mi> </msub> <mrow> <mo>(</mo> <mo>-</mo> <mfrac> <mrow> <mi>&Delta;</mi> <mi>y</mi> </mrow> <mn>2</mn> </mfrac> <mo>)</mo> </mrow> </mrow> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow>

将计算结果存储于计算机；

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>x</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>&mu;</mi> <mrow> <mi>i</mi> <mn>2</mn> <msub> <mi>k</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mi>V</mi> <mo>&CenterDot;</mo> <msub> <mi>V</mi> <mi>l</mi> </msub> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>y</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mi>&mu;</mi> <mrow> <mi>i</mi> <mn>2</mn> <msub> <mi>k</mi> <mi>x</mi> </msub> </mrow> </mfrac> <mi>V</mi> <mo>&CenterDot;</mo> <msub> <mi>V</mi> <mi>l</mi> </msub> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>

<mrow> <msub> <mover> <mi>W</mi> <mo>~</mo> </mover> <mi>z</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>k</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>k</mi> <mi>y</mi> </msub> <mo>,</mo> <msub> <mi>z</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>Z</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mi>&mu;</mi> <mi>K</mi> </mrow> <mrow> <mn>2</mn> <msub> <mi>k</mi> <mi>x</mi> </msub> <msub> <mi>k</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mi>V</mi> <mo>&CenterDot;</mo> <msub> <mi>V</mi> <mi>l</mi> </msub> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <msub> <mi>E</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>E</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>21</mn> <mo>)</mo> </mrow> <mo>.</mo> </mrow>