0 引言

20世纪60年代Rader^[1]提出了一种快速离散傅里叶变换算法，极大提高了离散傅里叶变换的计算效率，使得快速傅里叶变换(FFT)算法得到了更加广泛的应用，被评为20世纪十大最优秀算法之一。早期FFT算法要求数据维数具有2ⁿ形式，随后针对数据维数为3ⁿ、2ⁿK(n和K均为正整数)形式的FFT算法相继被提出^[2-5]。而基于质数分解的FFT算法的出现，解除了对数据维数的要求。Frigo等^[6]开发了FFTW软件包，集成各类FFT算法，实现了对任意维数数据的离散傅里叶变换快速计算。对于某些问题，空间或者频率域采样数据往往是不等间隔的，采用该采样方式的计算方法称为非规则采样数据离散傅里叶变换(NDFT，Non-Uniform Discrete Fourier Transform)。Dutt等^[7]、Anderson^[8]、Nguyen等^[9]、Potts等^[10]、Fessler等^[11]、Greengard等^[12]、Sørensen等^[13]通过改进FFT算法实现了NDFT。

由于傅里叶变换算子的振荡性，傅里叶变换可以归为振荡函数的积分。高振荡函数积分的计算是计算科学领域的研究热点。目前求解振荡函数积分的方法大致可以分为三类：Filon型方法，Levin型方法和渐进展开法。对于傅里叶变换而言，由于算子e^-ikx的矩已知，大多采用Filon型方法实现傅里叶变换^[14-18]。Clendenin^[14]采用的是线性函数；Piessens等^[19]采用了更为复杂的Chebyshev多项式；Arieh^[20]对基于Filon型方法的傅里叶变换算法中的积分点进行了分析，并在高震荡方程求解中取得良好的效果；Wu等^[21]将高斯积分应用于核函数的积分中，取得了较高的计算精度；Li等^[22]、Dai等^[23]和戴世坤等^[24]对高斯快速傅里叶变换与传统傅里叶变换在重、磁三维正演中的应用进行了对比研究。这些方法的区别在于剖分区间采用不同形式的多项式逼近核函数。

本文结合傅里叶变换积分与三次样条函数，提出一种快速傅里叶变换方法。该方法特点之一是能根据变换函数谱的变化趋势，任意地选取采样间距，提高了傅里叶变换剖分采样的灵活性，能兼顾计算精度和计算效率；特点之二是可根据离散后单元积分得到解析表达式，克服传统方法难以获得解析表达式的问题，进一步提高计算精度。

1 理论与方法 1.1 基于样条插值的一维傅里叶变换

函数f(m)的一维傅里叶正变换公式为

$ F(k) = \int_{ - \infty }^\infty {f(m)} {{\rm{e}}^{ - {\rm{i}}km}}{\rm{d}}m $

(1)

式中：k表示波数；F(k)为波数谱。对式(1)进行离散化

$ F(k) = \sum\limits_{i = 1}^M {\int_{{m_i}}^{{m_{i + 1}}} {f(m)} } {{\rm{e}}^{ - {\rm{i}}km}}{\rm{d}}m $

(2)

式中M表示总单元数。

对单元i的内核函数进行三次样条插值^[25]

$ \begin{array}{*{20}{l}} {f{{(m)}_i} = {d_i}{{(m - {m_i})}^3} + {c_i}{{(m - {m_i})}^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {b_i}(m - {m_i}) + {a_i}} \end{array} $

(3)

式中系数a_i、b_i、c_i、d_i可采用基于固定边界条件的样条插值计算得到。

将式(3)代入式(2)，得到

$ \begin{array}{l} F(k) = \mathop \sum \limits_{i = 1}^M \mathop \smallint \nolimits_{{m_i}}^{{m_{i + 1}}} [{d_i}{(m - {m_I})^3} + {c_i}{(m - {m_i})^2} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {b_i}(m - {m_i}) + {a_i}]{{\rm{e}}^{ - {\rm{i}}km}}{\rm{d}}m \end{array} $

(4)

令m′=m-m_i，则式(4)变为

$ F(k) = \mathop \sum \limits_{i = 1}^M \mathop \smallint \nolimits_0^L [{d_i}{m^{\prime 3}} + {c_i}{m^{\prime 2}} + {b_i}{m^\prime } + {a_i}]{{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}{\rm{d}}{m^\prime } $

式中L表示单元长度。令单元积分为

$ {I_{\rm{e}}} = (\int_0^L {{d_i}} {m^{\prime 3}} + {c_i}{m^{\prime 2}} + {b_i}{m^\prime } + {a_i}){{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}{\rm{d}}{m^\prime } $

(5)

据式(5)可得解析积分表达式为

$ {I_{\rm{e}}} = {d_i}{I_1} + {c_i}{I_2} + {b_i}{I_3} + {a_i}{I_4} $

(6)

其中

$ \begin{array}{*{20}{l}} {{I_1} = \int_0^L {{m^{\prime 3}}} {{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}{\rm{d}}{m^\prime }}\\ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{1}{{{k^4}}}\{ - 6 + k{m^\prime }[ - 6i + k{m^\prime }(3 + ik{m^\prime })]\} \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} {{\rm{e}}^{ - {\rm{i}}{k_x}({m^\prime } + {m_i})}}|_0^L \end{array} \end{array} $

$ \begin{array}{*{20}{l}} {{I_2} = \int_0^L {{m^{\prime 2}}} {{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}{\rm{d}}{m^\prime }}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{1}{{{k^3}}}[ - 2i + k{m^\prime }(2 + ik{m^\prime })]{{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}|_0^L} \end{array} $

$ {{I_3} = \int_0^L {{m^\prime }} {{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}{\rm{d}}{m^\prime } = \frac{1}{{{k^2}}}(1 + ik{m^\prime }){{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}|_0^L} $

$ {{I_4} = \int_0^L {{{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}} {\rm{d}}{m^\prime } = \frac{i}{k}{{\rm{e}}^{ - {\rm{i}}k({m^\prime } + {m_i})}}|_0^L} $

特别地，当波数k=0时，单元积分可简化为

$ {I_{\rm{e}}} = \int_0^L {({d_i}{m^{\prime 3}} + {c_i}{m^{\prime 2}} + {b_i}{m^\prime } + {a_i})} {\rm{d}}{m^\prime } $

(7)

对式(7)积分可得解析表达式

$ {I_{\rm{e}}} = {a_i}L + \frac{{{b_i}}}{2}{L^2} + \frac{{{c_i}}}{3}{L^3} + \frac{{{d_i}}}{4}{L^4} $

(8)

1.2 基于样条插值的二维傅里叶变换

对函数f(x, y)进行二维傅里叶正变换

$ F({k_x},{k_y}) = \int_{ - \infty }^\infty {\int_{ - \infty }^\infty {f(x,y)} } {{\rm{e}}^{ - {\rm{i}}({k_x}x + {k_y}y)}}{\rm{d}}x{\rm{d}}y $

(9)

式中：k_x为x方向的波数；k_y为y方向的波数；F(k_x, k_y)为波数谱。

式(9)的数值计算为二维积分，本文采取的策略是依次对x、y两个方向分别进行一维积分

$ {F({k_x},y) = \int_{ - \infty }^\infty {f(x,y)} {{\rm{e}}^{ - {\rm{i}}{k_x}x}}{\rm{d}}x} $

(10)

$ {F({k_x},{k_y}) = \int_{ - \infty }^\infty {F({k_x},y)} {{\rm{e}}^{ - {\rm{i}}{k_y}y}}{\rm{d}}y} $

(11)

2 正确性验证

为了分析误差，引入相对均方根误差^[26]，一维和二维的相对均方根误差公式分别为

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{rrms}} = \frac{{\sqrt {\sum\limits_{i = 1}^P {{{({{\hat B}_i} - {B_i})}^2}} } }}{{\sqrt {\sum\limits_{i = 1}^P {\hat B_i^2} } }} \times 100\% }\\ {{\rm{rrms }} = \frac{{\sqrt {\sum\limits_{i = 1}^P {\sum\limits_{i = j}^M {{{({{\hat B}_{ij}} - {B_{ij}})}^2}} } } }}{{\sqrt {\sum\limits_{i = 1}^P {\sum\limits_{i = j}^M {\hat B_{ij}^2} } } }} \times 100\% } \end{array}} \right. $

(12)

式中：P和M分别是x和y方向的节点数；B_i、B_ij分别表示表示一维和二维数值解；$ {\hat B_i}、{\hat B_{ij}}$分别表示一维和二维解析解。相对均方根误差能突出大异常值所占的比重，避免小异常和过零异常造成的误差失真。

为了验证该方法的正确性，对高斯函数进行样条插值傅里叶变换，并与解析公式进行对比。其中，一维和二维高斯函数及其傅里叶变换解析公式分别为

$ {f(x) = {{\rm{e}}^{ - \alpha {x^2}}} \Leftrightarrow F(k) = \sqrt {\frac{\pi }{\alpha }} {{\rm{e}}^{ - \frac{{{k^2}}}{{4a}}}}} $

(13)

$ {f(x) = {{\rm{e}}^{ - \alpha ({x^2} + {y^2})}} \Leftrightarrow F({k_x},{k_y}) = \frac{\pi }{\alpha }{{\rm{e}}^{ - \frac{{k_x^2 + k_y^2}}{{4a}}}}} $

(14)

式中α是任意非零常数，本文令α=0.001。

设定一维傅里叶变换的取值范围为[-100m, 100m]，采样点数为101，波数采样范围为[-0.3, 0.3]，采样点数为101；二维变换x、y方向点坐标范围均为[-100m, 100m], 采样点数均为101个，两个方向波数范围均为[-0.3, 0.3], 采样点数均为101。一维和二维傅里叶变换的波数域和空间域采样均采用等间隔剖分。

图 1为一维正、反傅里叶变换数值解与解析解曲线。可以看出，正、反傅里叶变换的数值解与解析解吻合度很高，正、反傅里叶变换的相对均方根误差分别约为4.5×10^-6和4.2×10^-7，验证了本文一维傅里叶正、反变换理论的正确性。

图 2为二维正、反傅里叶变换数值解与解析解等值线。可以看出，正、反傅里叶变换的数值解与解析解吻合度很高，正、反傅里叶变换的相对均方根误差分别约为6.0×10^-6和1.4×10^-5，验证了本文二维傅里叶正、反傅里叶变换理论的正确性。从计算结果可以看出，一维正反变换与二维正反变换的计算精度变化规律不同，这是由于不同的核函数谱的变化规律不同，因此其正、反傅里叶变换的计算精度不存在特定的规律性。

3 基于样条插值傅里叶变换的重磁场三维数值模拟

在频率域重磁场数值模拟中，傅里叶变换是关键步骤之一。设计连续介质模型，对基于样条插值傅里叶变换的重磁场三维数值模拟的计算效率与精度进行了研究，其中重磁场三维数值模拟采用空间波数混合域三维数值模拟方法^[22-24]。文献[25]表明基于高斯—FFT的重磁场数值模拟具有较高的计算精度，本文算例将高斯点数N_G=4的高斯FFT法重磁场数值模拟结果作为连续介质的解析解。本文算例中，令k_x=k_y=K。测试计算机配置为CPU-Inter Core i5-4590, 主频为3.3GHz, 内存为16GB。

3.1 重力场数值模拟

设计连续介质模型，其垂直方向密度不变，水平方向剩余密度ρ的空间变化可表示为

$ \rho = 2000 \times {{\rm{e}}^{ - 5 \times {{10}^{ - 8}}({x^2} + {y^2})}} $

(15)

模拟区域大小为20km(x)×20km(y)×10km(z)，坐标原点位于水平范围中心；异常体区域大小为20km(x)×20km(y)×3km(z)，顶面埋深为3km。模拟区域水平采样间隔均为200m，垂向采样间隔为100m，网格剖分数为101×101×101；异常体区域波数采样范围为-1.5×1.0^-2~1.5×1.0^-2。

图 3为不同波数下模拟的地面重力异常场及其梯度张量的相对均方根误差。可以看出，随着波数的增大，重力异常场及梯度张量的相对均方根误差逐渐减小，当波数K为71时，其计算精度与高斯FFT计算结果接近。当K继续增大时，梯度张量的计算误差有增大趋势，原因在于随着波数的增大，基于样条插值傅里叶变换的计算精度逐渐高于4点高斯FFT的计算精度，若把高斯FFT作为解析解，误差缓慢增大。

表 1所示为不同K时利用本文方法与高斯FFT方法计算地面重力场及其张量梯度耗时统计。可以看出，随着K的增大，样条插值FFT计算时间呈线性增长。在计算精度相近的情况下，即波数K=71时，样条插值FFT的计算时间为0.74s，高斯FFT的计算时间为4.85s，即前者约为后者的1/6。基于N_G=4的高斯FFT计算时间相当于传统FFT计算量的16倍，而本文算法仅需计算一次积分，因而在相同的计算精度下，基于样条插值的傅里叶变换计算速度高于高斯FFT算法。

图 4为K=71时地面重力异常场及梯度张量的数值解，图 5为图 4数据的绝对误差。从图 4和图 5可以看出，重力异常场的绝对误差最大约为0.05mGal，重力梯度张量的绝对误差最大约为0.06E，均比数值解(图 4)低3个数量级，表明其计算精度高，K选取71是合适的。

3.2 磁场数值模拟

设计一个连续介质模型，垂直方向磁化率不变，水平方向磁化率κ的空间变化为

$ \kappa = 0.02 \times {{\rm{e}}^{ - 5 \times {{10}^{ - 8}}x({x^2} + {y^2})}} $

(16)

模拟区域为20km(x)×20km(y)×10km(z)，坐标原点位于水平范围中心；异常区域大小为20km(x)×20km(y)×3km(z)，顶面埋深为3km。模拟区域网格剖分数为101×101×101；异常区域空间水平采样间隔均为200m，垂向采样间隔为100m。背景磁异常场为35A/m，磁倾角为45°，磁偏角为5°，波数K采样范围为-1.5×1.0^-2~1.5×1.0^-2。图 6所示为不同K时，采用样条插值FFT计算的地面磁异常场及其梯度张量的相对均方根误差。可以看出，随着K的增加，磁异常场及梯度张量的相对均方根误差逐渐减小；当K=101时，磁异常场及梯度张量的相对均方根误差小于0.1%，其计算精度与高斯FFT接近，当波数继续增大时，均方根误差降低的趋势变缓。

表 2所示为不同K时分别用本文方法及高斯FFT计算地面磁场及其梯度张量所耗时间对比。可以看出，随着K的增加，样条插值FFT计算时间呈线性增长。在计算精度相近的情况下，即K=101时本文方法耗时1.52s，高斯插值FFT算法耗时6.77s，即样条插值FFT耗时约为高斯插值FFT的1/4。

图 7为波数K=101时磁异常场及梯度张量的数值解，图 8为波数K=101时磁异常场及梯度张量的绝对误差。由图 7与图 8可以看出，磁异常场的绝对误差最大约为0.05nT，磁梯度张量的绝对误差最大约为2×10^-5nT/m，均比数值解(图 7)小2个数量级，表明本文方法计算精度高，且本实验中K选取101是合适的。

4 应用

为了进一步验证基于样条插值FFT方法对实测数据的应用效果，笔者利用“地理空间数据云”平台下载了安徽省数字高程数据(经度为117.6°~118.4°，纬度为30.5°~31.2°)，对此数据进行数值模拟。该区地形如图 9所示，境内中部丘陵、岗地起伏，呈北东向展布。假设该地区地下地层为水平层状介质^[26]，采用基于样条插值FFT的空间波数混合域三维数值模拟方法计算其重力异常场，结果如图 10所示。可以看出，重力高的区域与图 9所示起伏地形范围吻合很好，说明基于样条插值FFT方法应用于实际数据的处理是可靠的，具有较好的适应性。

5 结论

针对传统FFT存在截断效应的问题，本文结合FFT与三次样条插值算法，提出了一种基于样条插值的FFT算法。该方法充分利用样条插值拟合的高阶连续性及采样灵活的优点，兼顾计算精度和计算效率，实现了高效、高精度的FFT。得到如下结论：

(1) 利用高斯函数FFT，对比了数值解与解析解，验证了基于样条插值算法的一维、二维FFT的正确性。

(2) 将算法应用于频率域重、磁位场的计算，设计连续介质模型，对比基于样条插值FFT与高斯插值FFT的数值解，前者的计算精度更高。

(3) 对于连续介质模型，在计算精度相近的情况下，利用基于样条插值FFT进行重、磁场模拟计算，耗时统计数据表明，样条FFT的计算效率明显高于高斯FFT，前者的计算时间约为后者的1/6(重力数据)和1/4(磁力数据)。