0 引言

重磁资料反演是解释人员从重磁资料中获取地质信息的重要手段，而绝大多数反演方法都是以正演方法为基础构建的。因此重磁数据的正演一直是地球物理学家研究的重点和热点之一。

重磁正演问题分两大类：空间域正演和波数域正演。空间域正演的优点是直观、精度高；缺点是公式推导较难，导出结果繁琐，对于较复杂的形体，则难以导出正演解析表达式。空间域正演的另一缺点是计算速度慢。波数域正演的优点是公式推导简单，导出结果简洁。对于很多空间域无法导出的正演公式，在波数域中很容易导出。波数域正演的另一优点是计算速度快。波数域正演的主要缺点是精度低。为了提高波数域正演计算精度，1988年，Chai等^[1]提出了波数域最佳偏移抽样法(偏移量为0.22倍波数域采样间隔)，将重磁波数域正演精度提高了几十倍。2014年，Wu等^[2]将柴玉璞^[3-4]提出的移样离散傅里叶变换(SFT)与高斯节点积分相结合，提出了一种高精度快速位场波数域正演方法，对重磁勘探技术的发展做出了很大贡献，但尚有以下不足。

(1) 作者对所提方法的论证是从位函数谱的性质出发，未能充分展现高斯FFT算法广泛的应用领域。

(2) 文中缺少将高斯节点积分引入傅里叶变换数值计算的数学演绎过程，不免让读者认为文中公式(19)仅源于观察、公式(20)仅源于猜想。而此种感觉乃数学论证欠严谨所致。

(3) 文中未给出偏移量与高斯节点的换算关系、权系数与高斯求积系数之间的换算关系，也未直接给出区间[0, 1]内大于2的高斯抽样点(即偏移量)和权系数，无法实现编程。

本文从移样离散傅里叶变换理论出发，对高斯FFT算法展开论证，将其应用范围扩展到任意有界函数的正、反傅里叶变换。文中给出了将高斯节点积分引入傅里叶变换数值计算的严谨数学演绎过程，并在此过程中导出了偏移量与高斯节点的换算关系以及权系数与高斯求积系数之间的换算关系。这样借助数学手册^[5]中的高斯型求积节点和求积系数数据表，即可实现编程。

1 DFT变换的推广

目前采用的离散傅里叶变换(DFT)表达式为

$ \left\{ {\begin{array}{*{20}{l}} {\tilde F(md) = \varDelta \sum\limits_{n = - \frac{N}{2}}^{\frac{N}{2} - 1} {\tilde f} (n\varDelta ){{\rm{e}}^{ - {\rm{i}}2\pi mn/N}}}\\ {\tilde f(n\varDelta ) = d\sum\limits_{m = - \frac{N}{2}}^{\frac{N}{2} - 1} {\tilde F} (md){{\rm{e}}^{{\rm{i}}2\pi mn/N}}} \end{array}} \right. $

(1)

式中：$ \tilde f\left( {n \mathit\Delta } \right)$表示函数f(x)以Δ为抽样间隔的抽样序列折叠成的长度为N的序列，n为序列序数；$\tilde F\left( {md} \right) $表示F(u)(f(x)的谱)以d为抽样间隔的抽样序列折叠成的长度为N的序列，m为序列序数。很明显，式(1)有较大的局限性，即只能由过函数原点的一组等间隔抽样值f(nΔ)，计算包含零频点的一组等间隔频点上的谱值F(md)。可以将式(1)推广为

$ \left\{ {\begin{array}{*{20}{l}} {{F_\xi }[(m + \eta )d] = \varDelta \sum\limits_{n = - \frac{N}{2}}^{\frac{N}{2}} {{f_\eta }} [(n + \xi )\varDelta ]{{\rm{e}}^{ - i2\pi (m + \eta )(n + \xi )/N}}}\\ {{f_\eta }[(n + \xi )\varDelta ] = d\sum\limits_{ - \frac{N}{2}}^{\frac{N}{2} - 1} {{F_\xi }} [(m + \eta )d]{{\rm{e}}^{i2\pi (m + \eta )(n + \xi )/N}}} \end{array}} \right. $

(2)

式(2)称为DFT_ξη变换，其中0≤ξ, η≤1。这一推广的正确性显而易见，把式(2)中的一个等式代入另一个等式，即可证明两式的互逆性。利用新确立的DFT_ξη变换可以表述一个重要定理——傅里叶变换离散化定理。

2 移样离散傅里叶变换理论和算法

柴玉璞在其发表在《中国科学》的论文^[3]中证明了一个重要定理——傅里叶变换离散化定理。该定理揭示了函数离散值与其谱离散值之间的关系，并据此给出了严谨的、突破性的傅里叶变换数值计算理论——移样离散傅里叶变换理论。

2.1 傅里叶变换离散化定理

设f(x)与F(u)为一傅里叶变换对，且均无无穷型间断点。依据下式

$ \left\{ {\begin{array}{*{20}{l}} {{f_\eta }[(n + \xi )\varDelta ] = \sum\limits_{k = - \infty }^\infty {{{\rm{e}}^{ - {\rm{i}}2\pi k\eta }}} f[(n + \xi + kN)\varDelta ]}\\ {{F_\xi }[(m + \eta )d] = \sum\limits_{l = - \infty }^\infty {{{\rm{e}}^{{\rm{i}}2\pi l\xi }}} F[(m + \eta + lN)d]} \end{array}} \right. $

(3)

将f(x)和F(u)离散化，并加权折叠成两个长度为N的复序列(dΔ=1/N)，则两序列恰好构成一DFT_ξη变换对

$ \begin{array}{*{20}{l}} {{\rm{DF}}{{\rm{T}}_{\xi \eta }}\left\{ {\sum\limits_{k = - \infty }^\infty {{{\rm{e}}^{ - {\rm{i}}2\pi k\eta }}} f[(n + \xi + kN)\varDelta ]} \right\}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{l = - \infty }^\infty {{{\rm{e}}^{{\rm{i}}2\pi l\xi }}} F[(m + \eta + lN)d]} \end{array} $

(4a)

或

$ \begin{array}{*{20}{l}} {{\rm{IDF}}{{\rm{T}}_{\xi \eta }}\left\{ {\sum\limits_{l = - \infty }^\infty {{{\rm{e}}^{{\rm{i}}2\pi l\xi }}} F[(m + \eta + lN)d]} \right\}}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = - \infty }^\infty {{{\rm{e}}^{ - {\rm{i}}2\pi k\eta }}} f[(n + \xi + kN)\varDelta ]} \end{array} $

(4b)

2.2 移样离散傅里叶变换理论和算法

将式(4a)中k=0和l=0的两项置于方程左边，其余项放在右边，可得

$ \begin{array}{*{20}{l}} {{{{\mathop{\rm DFT}\nolimits} }_{\xi \eta }}[f(n + \xi )\varDelta ] - F[(m + \eta )d]}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{l = - \infty ,l \ne 0}^\infty {{{\rm{e}}^{{\rm{i}}2\pi l\xi }}} F[(m + \eta + lN)d] - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{{\mathop{\rm DFT}\nolimits} }_{\xi \eta }}\left\{ {\sum\limits_{k = - \infty ,k \ne 0}^\infty {{{\rm{e}}^{ - {\rm{i}}2\pi k\eta }}} f[(n + \xi + kN)\varDelta ]} \right\}} \end{array} $

(5a)

上式即为DFT_ξη算法的误差方程。同理可由式(4b)导出IDFT_ξη算法的误差方程

$ \begin{array}{*{20}{l}} {{{{\mathop{\rm IDFT}\nolimits} }_{{\xi _\eta }}}\{ F[(m + \eta )d]\} - f[(n + \xi )\varDelta ]}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = - \infty ,k \ne 0}^\infty {{{\rm{e}}^{ - {\rm{i}}2\pi k\eta }}} f[(n + \xi + kN)\varDelta ] - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{{\mathop{\rm IDFT}\nolimits} }_{\xi \eta }}\left\{ {\sum\limits_{l = - \infty ,l \ne 0}^\infty {{{\rm{e}}^{{\rm{i}}2\pi l\xi }}} F[(m + \eta + lN)d]} \right\}} \end{array} $

(5b)

式(5a)、式(5b)中，左端第一项分别为DFT_ξη和IDFT_ξη算法表达式，第二项分别为离散真谱和函数真值表达式，右端为误差表达式。

实际上，式(5)以数学语言表述了一个严谨且具突破性的傅里叶变换数值计算理论——移样离散傅里叶变换。其严谨性在于它不仅给出了算法表达式，而且给出了误差表达式；其突破性在于它给出了更广义、更普遍的新算法——SFT算法(式(2))。利用这一新算法，既可由函数的任意一组等间隔抽样值计算任意一组等间隔频点上的谱值，也可以由任意一组等间隔频点上的谱值计算任意一组函数等间隔抽样值。这就为将高斯节点积分引入傅里叶变换数值计算创造了条件。因此，将高斯节点积分引入傅里叶变换数值计算，是移样离散傅里叶变换理论的一个重要应用。

3 将高斯节点积分引入傅里叶变换数值计算

高斯节点积分问题有严谨的数学证明^[6]，比较麻烦，也很深奥，涉及拉格朗日插值、多项式正交性、勒让德函数等问题。作为应用者，仅将其结论深入浅出地表述如下：对于区间[-1, 1]上有界函数f(x)的积分，可用其M个点的函数值加权求和高精度逼近，其条件是这M个点需是M阶勒让德多项式p_M(x)的M个根。高斯节点积分的数学表达为

$ \left\{ {\begin{array}{*{20}{l}} {\int_{ - 1}^1 f (x){\rm{d}}x \approx \sum\limits_{k = 1}^M {\lambda _k^{(M)}} f\left[ {x_k^{(M)}} \right]}\\ {\lambda _k^{(M)} = \frac{2}{{\left\{ {1 - {{\left[ {x_k^{(M)}} \right]}^2}} \right\}{{\left\{ {p_M^\prime \left[ {x_k^{(M)}} \right]} \right\}}^2}}}} \end{array}} \right. $

(6)

式中：x_k^(M)为M阶勒让德多项式的第k个根，也称高斯积分的第k个求积节点；λ_k^(M)为该节点对应的求积系数(x_k^(M)和λ_k^(M)的值无需计算，各种数学手册中均有x_k^(M)与λ_k^(M)数据表)；p′_M为M阶勒让德多项式的一阶导数。高斯节点积分的误差为

$ R = \frac{{{2^{2M + 1}}{{(M!)}^4}}}{{(2M + 1){{[(2M)!]}^3}}}{f^{(2M)}}(\theta ) $

(7)

式中：θ为区间[-1, 1]上的某一点；f^(2M)(θ)为函数f(x)的2M阶导数在该点的值。由式(7)可进一步分析高斯节点积分误差随节点M的增加而急剧变小的速率。例如：当M=2时，$ R = \frac{{{f^{(4)}}\left( \theta \right)}}{{135}}$；当M=4时，$ R = \frac{{{f^{(8)}}\left( \theta \right)}}{{3472875}}$。据此断定，如果将高斯节点积分引入傅里叶变换数值计算，随高斯抽样点数的增加，其计算精度将大幅度、递进式提高。

高斯节点积分是针对实函数的。由于复数谱的实部和虚部都是实函数，所以正、反傅里叶变换都可以引入高斯节点积分。所谓将高斯节点积分引入傅里叶变换数值计算，就是将傅里叶积分的N个等间隔子区间[nΔ, (n+1)Δ]上的积分，转换为区间[-1, 1]上的高斯积分。众所周知，对于变量t在区间[a，b]上的积分可通过变量替换

$ t = \frac{{b - a}}{2}x + \frac{{b + a}}{2} $

(8)

转换为区间[-1, 1]上的积分，即

$ \int_a^b f (t){\rm{d}}t = \frac{{b - a}}{2}\int_{ - 1}^1 f \left( {\frac{{b - a}}{2}x + \frac{{b + a}}{2}} \right){\rm{d}}x $

(9)

利用式(9)将子区间[nΔ, (n+1)Δ]上的积分转换为区间[-1, 1]上的高斯节点积分

$ \begin{array}{l} \int_{n\varDelta }^{(n + 1)\varDelta } {f(t){\rm{d}}t} = 0.5\varDelta \int_{ - 1}^1 f [(n + 0.5x + 0.5)\varDelta ]{\rm{d}}x\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \approx \varDelta \sum\limits_{k = 1}^M 0 .5\lambda _k^{(M)}f\left[ {\left( {n + 0.5x_k^{(M)} + 0.5} \right)\varDelta } \right]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \varDelta \sum\limits_{k = 1}^M {\lambda _k^{(M)}} f\left[ {\left( {n + \bar x_k^{(M)}} \right)\varDelta } \right] \end{array} $

(10)

进而可得

$ \begin{array}{l} \int_{ - \infty }^\infty f (t){\rm{d}}t \approx \sum\limits_{n = - \frac{N}{2}}^{\frac{N}{2} - 1} {\int_{n\varDelta }^{(n + 1)\varDelta } f } (t){\rm{d}}t\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{n = - \frac{N}{2}}^{\frac{N}{2} - 1} \varDelta \sum\limits_{k = 1}^M {\bar \lambda _k^{(M)}} f\left[ {\left( {n + \bar x_k^{(M)}} \right)\varDelta } \right]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = 1}^M {\bar \lambda _k^{(M)}} \varDelta \sum\limits_{n = - \frac{N}{2}}^{\frac{N}{2} - 1} f \left[ {\left( {n + \bar x_k^{(M)}} \right)\varDelta } \right] \end{array} $

(11a)

式中：$\overline x _k^{(M)} $表示第k个求积节点；$\overline \lambda _k^{(M)} $表示第k个求积节点对应的求积系数。高斯节点数取偶数，可充分利用高斯节点的对称性和对称高斯节点序列谱的共轭性，提高计算效率，故式(10)及之后的M规定为偶数。同理可得

$ \begin{array}{l} \int_{ - \infty }^\infty F (t){\rm{d}}t \approx \sum\limits_{m = - \frac{N}{2}}^{\frac{N}{2} - 1} {\int_{md}^{(m + 1)d} F } (t){\rm{d}}t\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = 1}^M {\bar \lambda _k^{(M)}} d\sum\limits_{m = - \frac{N}{2}}^{\frac{N}{2} - 1} F \left[ {\left( {m + \bar x_k^{(M)}} \right)d} \right] \end{array} $

(11b)

在设定Δd=1/N的条件下，若将式(11)中求和号下的f(t)和F(t)分别具体为傅里叶正、反变换的积分核，对照式(2)，即有如下重要结论：一个傅里叶积分可以M个移样离散傅里叶变换的加权求和高精度地逼近，其偏移量为$ \overline x _k^{(M)} = 0.5x_k^{(M)} + 0.5\left( {x_k^{(M)}为高斯节点} \right)$，权系数为$ \bar \lambda _k^{(M)} = 0.5\lambda _k^{(M)}\left( {\lambda _k^{(M)}{\rm{为该节点对应的求积系}}} \right)$。为了与移样离散傅里叶变换中的符号术语一致，在式(11)中f(t)和F(u)分别具体为正、反傅里叶变换的积分核后导出的结果中$\overline x _k^{(M)} $更名为偏移量，且在傅里叶正变换中以ξ_i^(M)替换，在傅里叶反变换中以η_i^(M)替换；$ \overline \lambda _k^{(M)} $更名为权系数，且在傅里叶正、反变换中，均以 $\overline \lambda _k^{(M)} $ 替换；下标由k换为i，其意义也由勒让德多项根的序号变为参与加权求和的移样离散傅里叶变换的序号。

由此可见，傅里叶积分

$ F(u) = \int_{ - \infty }^\infty f (x){{\rm{e}}^{ - {\rm{i}}2\pi ux}}{\rm{d}}x $

(12a)

可用下式高精度逼近

$ \bar F[(m + \eta )d] = \sum\limits_{i = 1}^M {\lambda _i^{(M)}} \mathrm{DFT} { _{\xi _i^{(M)}\eta }}\{ f[(n + \xi _i^{(M)})\varDelta ]\} $

(12b)

傅里叶积分

$ f(x) = \int_{ - \infty }^\infty F (u){{\rm{e}}^{{\rm{i}}2\pi ux}}{\rm{d}}u $

(13a)

可用下式高精度逼近

$ \bar f[(n + \xi )\varDelta ] = \sum\limits_{i = 1}^M {\lambda _i^{(M)}} \mathrm{DFT}{ _{\xi \eta _i^{(M)}}}\{ F[(m + \eta _i^{(M)})d]\} $

(13b)

由勒让德多项式的性质可知，高斯节点具有关于区间中心的对称性。式(4b)中的f为实函数，显然，利用该式(其左端只保留l=0项，即忽略截断误差)可证明，移样离散傅里叶反变换中，对称高斯节点序列的谱具有共轭性。因此，式(13b)可简化为

$ \begin{array}{*{20}{l}} {\bar f[(n + \xi )\varDelta ] = 2\sum\limits_{i = 1}^{M/2} {\lambda _i^{(M)}} Re \{ \mathrm{DFT}{ _{\xi \eta _i^{(M)}}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \{ F[(m + \eta _i^{(M)})d]\} \} } \end{array} $

(13c)

式(12b)和式(13c)即是一维傅里叶正、反变换数值计算的实用公式。

柴玉璞在其专著^[4]中将一维DFT_ξη和IDFT_ξη算法误差方程推广到二维，其形式与一维完全一致(式(5))。对于二维IDFT_ξη算法误差方程，略去有限效应项，令ξ=0，并将左端k₁=k₂=0项移到右端，得到

$ \begin{array}{*{20}{c}} { \mathrm{DFT} _{0\eta }^2\{ F[({m_1} + {\eta _1})d,({m_2} + {\eta _2})d]\} }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{{k_1} = - \infty }^\infty {\sum\limits_{{k_2} = - \infty }^\infty {{{\rm{e}}^{{\rm{i}}2\pi {k_1}{\eta _1}}}} } {{\rm{e}}^{{\rm{i}}2\pi {k_2}{\eta _2}}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f[({n_1} + {k_1}N)\varDelta ,({n_2} + {k_2}N)\varDelta ]} \end{array} $

(14)

式中：IDFT_0η²表示空间域不移样的二维移样离散傅里叶反变换；$f[({n_1} + {k_1}N) \mathit\Delta , ({n_2} + {k_2}N) \mathit\Delta ] $表示二维函数离散值；$ F[({m_1} + {\eta _1})d, ({m_2} + {\eta _2})d]$为前者的复谱；m₁、η₁分别表示u上的抽样序列号和偏移量；m₂、η₂分别表示v上的抽样序列号和偏移量；n₁、n₂分别表示x轴和y轴上的抽样序列号。

二维傅里叶变换相当于进行两次一维傅里叶变换，因此，若将高斯节点积分引入二维傅里叶变换(方法同一维，只是分别对两重积分实施)

$ F(u,v) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty f } (x,y){{\rm{e}}^{ - {\rm{i}}2\pi (ux + vy)}}{\rm{d}}x{\rm{d}}y $

(15a)

则其积分可用下式高精度逼近

$ \left\{ {\begin{array}{*{20}{l}} {\bar F\left[ {\left( {{m_1} + \eta } \right)d,\left( {{m_2} + \eta } \right)d} \right]}\\ {\quad = \sum\limits_{i = 1}^M {\lambda _i^{(M)}} {{{\mathop{\rm DFT}\nolimits} }_{\xi _{2,i}^{(M)}\eta }}\left[ {{{\bar F}^\prime }_{\xi _{2,i}^{(M)}\eta }} \right]}\\ {{{\bar F}_{\xi _{2,i}^{(M)}\eta }} = \sum\limits_{i = 1}^M {\lambda _i^M} {{{\mathop{\rm DFT}\nolimits} }_{\xi _{1,i}^{(M)}\eta }}\left\{ {f\left[ {\left( {{n_1} + \xi _{1,i}^{(M)}} \right)\varDelta ,} \right.} \right.}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\left. {{\kern 1pt} \left( {{n_2} + \xi _{2,i}^{(M)}} \right)\varDelta } \right]} \right\}} \end{array}} \right. $

(15b)

式中ξ_{1, i}、ξ_{2, i}分别表示第i个移样离散傅里叶变换中x轴和y轴上的偏移量。同理，将高斯节点积分引入二维傅里叶反变换

$ f(x,y) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty F } (u,v){{\rm{e}}^{{\rm{i}}2\pi (ux + vy)}}{\rm{d}}u{\rm{d}}v $

(16a)

则上式可高精度逼近为

$ \left\{ \begin{array}{l} \bar f\left[ {\left( {{n_1} + \xi } \right)\varDelta ,\left( {{n_2} + \xi } \right)\varDelta } \right]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{i = 1}^M {\lambda _i^{(M)}} {{\mathop{\rm IDFT}\nolimits} _{\xi \eta _{2,i}^{(M)}}}\left[ {{{\bar f}^\prime }_{\xi \eta _{2,i}^{(M)}}} \right]\\ {{\bar f}^\prime }_{\xi \eta _{2,i}^{(M)}} = \sum\limits_{i = 1}^M {\lambda _i^{(M)}} {{\mathop{\rm IDFT}\nolimits} _{\xi \eta _{1,i}^{(M)}}}\left\{ {F\left[ {{m_1} + \eta _{1,i}^{(M)}} \right)d} \right.,\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\left. {\left( {{m_2} + \eta _{2,i}^{(M)}} \right)d} \right]} \right\} \end{array} \right. $

(16b)

式中η_{1, i}、η_{2, i}分别表示第i个移样离散傅里叶变换中u轴和v轴上的偏移量。由式(14)可判断，二维傅里叶反变换的第一重变换中，其对称高斯节点序列的谱无共轭性，但第二重变换其对称高斯节点序列的谱具有共轭性。因此式(16b)可简化为

$ \left\{ \begin{array}{l} \bar f\left[ {\left( {{n_1} + \xi } \right)\varDelta ,\left( {{n_2} + \xi } \right)\varDelta } \right]\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = 2\sum\limits_{i = 1}^{M/2} {\lambda _i^{(M)}} {\mathop{\rm Re}\nolimits} \left\{ {{{{\mathop{\rm IDFT}\nolimits} }_{\xi \eta _{2,i}^{(M)}}}\left[ {{{\bar f}^\prime }_{\xi \eta _{2,i}^{(M)}}} \right]} \right\}\\ {{\bar f}^\prime }_{\xi \eta _{2,i}^{(M)}} = \sum\limits_{i = 1}^{M/2} {\lambda _i^{(M)}} {{\mathop{\rm IDFT}\nolimits} _{\xi \eta _{1,i}^{(M)}}}\left\{ {F\left[ {\left( {{m_1} + \eta _{1,i}^{(M)}} \right)d,} \right.} \right.\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {\left. {{\kern 1pt} \left( {{m_2} + \eta _{2,i}^{(M)}} \right)d} \right]} \right\} \end{array} \right. $

(16c)

式(15b)和式(16c)即是二维傅里叶正、反变换数值计算的实用公式。式中两重一维移样离散傅里叶变换的高斯节点数理论上可以不同，但一般情况下无此必要，故本文未作区别。

柴玉璞在其专著^[4]中详细论述了SFT与DFT的关系，并指出SFT仍可借用FFT实现快速计算。其实，只需将DFT_ξη的定义式(式(2))展开，上述关系和结论便一目了然。

4 位场波数域正演中的应用

位场是有界函数，而且其谱呈负指数规律快速衰减，所以高斯FFT算法可用于位场反傅里叶变换，且效果特别好。位场反傅里叶变换，即位场波数域正演，也是重磁勘探的重要课题之一，也特别需要借助这一数学方法提高计算精度。所以本文选择位场波数域正演作为这一数学方法的应用领域，并沿用Wu等^[2]对这一快速、高精度位场波数域正演方法的命名——高斯傅里叶正演。

设计一个由5个同尺寸的直立方柱体组成的模型。5个直立方柱体的高度均为2km，顶面埋深为1km，顶面面积为10km×10km，剩余密度均取2.0g/cm³。这5个方柱体的平面位置：其中一个的中心与图 1的中心点重合，其余4个的中心分别与图 1的四条边的中点重合。正演计算面积为64km× 64km，网格距为0.5km，网格数据点为128×128。图 1是采用空间域方法计算的模型重力异常。图 2为高斯抽样点数分别为2、4、6时的高斯傅里叶正演结果及其对应的正演误差。高斯抽样点数为2、4、6时的均方根误差分别为13.300、0.060、0.001mGal，说明随着高斯抽样点数的增加，误差急剧减小，当高斯抽样点增加到4时，精度已可满足常规解释要求。

高斯傅里叶正演的一般规律是，场源离计算窗中心越远，其谱的波动频率越高，则需要更多的高斯抽样点才能达到期望精度，特别是在源跨越计算窗边、延伸到计算窗外的情况。当高斯抽样点为2时，中心方柱体的正演误差已经很小，以至于可以忽略，但周围四个方柱体的正演误差仍相当大(图 2a右)；当高斯抽样点为4时，周围四个方柱体的正演误差已微乎其微(图 2b右)；当高斯抽样点为6时，整个计算区域内的误差已经降至非常低(约10^-3mGal)(图 2c右)。实际上此时的误差是谱的截断误差，它不再随高斯抽样点的增加而减小。

理论和实践都证明，无论场源分布如何，高斯傅里叶正演都能通过增加高斯抽样点达到预期的精度，且仍能保持位场波数域正演计算速度快的优点。

5 结论

(1) 快速傅里叶变换(FFT)解决了傅里叶变换数值计算的效率问题，但其抽样点需是一个傅里叶积分被划分成的数个子区间的端点；高斯节点积分具有很高精度，但其抽样点须是子区间内的某些点。移样离散傅里叶变换(SFT)理论的确立，客观上为解决这一矛盾创造了条件。因此可以说，在地球物理学家将FFT和高斯节点积分相结合、提出位场高斯傅里叶正演算法的创新中，SFT功不可没。

(2) 本文以严谨的数学演绎，不仅导出了傅里叶积分可以数个移样离散傅里叶变换(SFT)的加权求和高精度地逼近的结论，而且导出了偏移量与高斯节点的换算关系、权系数与高斯求积系数之间的换算关系。这为位场高斯傅里叶正演方法的推广提供了严谨的理论支撑和编程实现上的可行性和方便条件。

(3) 高斯FFT算法已经在位场波数域正演中得到成功应用。因为高斯节点积分理论和移样离散傅里叶变换理论的充分条件，都是任意有界函数，因此有理由相信，该方法未来有可能在更广阔的领域获得更多应用。

感谢中国石油集团东方地球物理公司贾继军、邓玉友和孙喜明三位高级工程师在计算机操作和绘图方面提供的帮助。