0 引言

频率域可控源电磁法广泛应用于地下水勘探、环境调查、地质填图和油气勘探(汤井田和何继善, 2005; Constable, 2010; Streich, 2016).频率域电磁法可以分为两类：频率电磁测深法和频率电磁剖面法(李金铭, 2005).频率电磁测深法是通过改变发射频率达到测深目的的一类电磁法，典型的方法为可控源音频大地电磁法(CSAMT)，采用接地长导线或不接地回线作为场源，在波区测量不同频率的相互正交的电磁场分量，并计算卡尼亚电阻率，CSAMT的工作频率范围为0.125~8000 Hz(汤井田和何继善, 2005).频率电磁剖面法是保持发射频率不变，通过改变收发距得到某一深度范围内地下电导率分布，典型的方法为海洋可控源电磁法(CSEM)，采用长度为几百米的水平导线源随船拖拽移动发射低频信号，布设于海底表面的接收机阵列测量水平电场和三分量磁场响应，海洋CSEM的工作频率范围为0.1~10 Hz(Constable, 2010).

对于频率电磁测深法，在正演和反演解释过程中，为了得到一个三维模型的电磁响应，需要计算多个频率的三维正演；对于频率电磁剖面法，在勘探设计时需要计算多个频率的三维正演，用于确定一个特定三维模型的最优化发射频率.可见，频率域电磁法的勘探最优化设计和数据反演解释都需要计算多个频率的三维正演响应.目前模拟复杂三维介质的频率域电磁响应的正演算法主要有：有限差分法(Streich, 2009；殷长春等, 2014)，有限元法(徐志锋和吴小平, 2010；Um et al., 2013；李建慧等, 2016；Cai et al., 2017)和有限体积法(张烨等, 2012；杨波等, 2012；Haber and Ruthotto, 2014；彭荣华和胡祥云, 2016).不论采用哪种空间离散算法，最后的离散控制方程都可以表示为以下形式：

其中K表示离散形式的双旋度算子，M表示离散形式的电导率，u为未知的电场矢量，q为离散形式的源项.一般采用迭代法(殷长春, 2014)或直接法(Streich, 2009；彭荣华和胡祥云, 2016；Cai et al., 2017)求解上述方程.由于该方程的系数矩阵K+iωM是频率的函数，因此不论采用迭代法或直接法，不同频率的正演响应都只能独立求解，从而导致正演计算时间随发射频率个数呈线性倍数增加.

一种改进方法是针对迭代法的预处理重复利用技术(Um et al., 2013).在采用迭代法求解离散控制方程时，ILU预处理矩阵的求解需要耗费较长的计算时间.Um等(2013)提出只求解最低频率的ILU预处理矩阵，并将该预处理矩阵重复应用于其他频率的正演计算，从而减少总的正演计算时间.但该算法的改进幅度有限，而且该算法对于大的频率间隔(比如CSAMT)的正演模拟的有效性还有待验证.另一种改进方法是采用多频率并行计算(Plessix et al., 2007).由于各个频率的正演响应是相互独立的，因此可以采用多CPU的计算设备，同时计算多个频率的正演响应(彭荣华和胡祥云, 2016; Cai et al., 2017).当计算设备的CPU数量不少于发射频率数量时，其所有频率的总的正演计算时间与一个频率的正演计算时间相当.上述两种改进算法本质上都是基于不同发射频率正演响应的独立求解，当发射频率数量增加时，方程求解次数或CPU数量也需要相应增加.

更为有效的改进方法是将正演响应表示为发射频率的传递函数(蒋耀林, 2010)，从而实现一定频度范围内所有频率正演响应的同时求解.但由于该传递函数的分母中包含大维度的离散系数矩阵K，难以直接求解，一般采用模型降阶方法求解该传递函数(蒋耀林, 2010).所谓模型降阶，是指在某种情况下将一个较大系统转化为一个近似的较小系统的过程，以便降低大型复杂系统的理论分析难度和减少数据运算量，同时保证降阶系统与原始系统的近似误差足够小、降阶系统保持原始系统的某些性能(蒋耀林, 2010).目前模型降阶方法因其高效快速的计算性能被广泛应用于各类大型复杂问题的数值模拟与参数最优化计算中(蒋耀林, 2010).典型的模型降阶方法是基于m维多项式Krylov子空间投影的模型降阶方法(Druskin and Knizhnerman, 1994; Druskin et al., 1999)，但该算法的收敛速度依赖于计算频率范围和模型的电导率反差，对于较长的求解频率范围、存在电导率大反差模型的正演，多项式Krylov子空间的维度需要达到数千甚至上万，严重降低了该方法的计算效率.针对这一问题，Knizhnerman等(2009)提出采用基于m维有理函数Krylov子空间投影的模型降阶方法，有理函数相比多项式表现出更好的近似特性，不依赖于模型电导率反差，该算法的关键在于(1)选取最优化的m个有理函数极点，极点选取是否合适直接决定了有理函数近似特性的好坏(即m的大小)；(2)构建m维有理函数Krylov子空间需要逐次求解m个关于系数矩阵和极点的方程，如何有效求解这些方程直接决定了模型降阶方法的计算效率.Knizhnerman等(2009)通过求解第三类Zolotarev问题得到了准最优化的m个(m一般为几十)不同有理函数极点，采用迭代法逐次求解m次关于系数矩阵和不同极点的方程，该算法能有效计算一定频率范围内大反差电导率模型的正演响应，但该算法构建子空间的计算量较大.

本文研究快速实现一定频度范围内所有频率的可控源电磁法三维正演响应的同时求解.首先采用基于Yee氏交错网格的拟态有限体积算法(Haber and Ruthotto, 2014)实现空间离散，将离散控制方程的解表示为传递函数的形式.然后采用基于有理函数Krylov子空间投影的模型降阶方法求解该传递函数.针对有理函数Krylov子空间构建需要多次求解关于系数矩阵和极点的线性方程组的问题，本文提出采用构建m维单个重复极点的有理函数Krylov子空间，引入直接法求解器PARDISO (Schenk and Grtner, 2004)，只需要一次系数矩阵分解和m次矩阵回代即可实现有理函数Krylov子空间的构建，极大地减少了正演算法的计算量.针对最优化有理函数极点选取问题，由于本文将有理函数的极点限制为单个重复极点，极大地缩小了极点搜索范围，简化了最优化极点搜索过程.本文根据传递函数的有理函数Krylov子空间投影算法的误差分析理论(Güttel, 2010)，引入关于单个重复极点的收敛率函数，通过求解有理函数的最大收敛率直接给出最优化的单个重复极点公式.最终实现了传递函数的快速计算.最后通过多个理论模型的正演计算测试验证本文算法计算多频可控源电磁法三维正演响应的精确性和高效性，并通过一个三维海洋CSEM勘探设计最优化发射频率和接收区域选取的例子进一步说明本文算法的优点.

1 正演理论 1.1 拟态有限体积空间离散算法

设电磁场随时间的变化关系为e^iωt，忽略位移电流，则频率域可控源电磁法对应的控制方程为

(1a)

(1b)

其中E是电场强度矢量，B是磁感应强度矢量，s是外加源项，σ为电导率，ω=2πf为圆频率，μ是磁导率.采用自然边界条件(Haber and Ruthotto, 2014)

(2)

与有限元方法一样，MFV方法处理弱形式的控制方程：

(3a)

(3b)

其中W与E位于相同的Sobolev空间H(Curl; Ω)，F与B位于相同的Sobolev空间H(Div; Ω).内积定义为 ，其中A和G为Sobolev空间H(Curl; Ω)或H(Div; Ω)中的任意参数.

采用Yee氏交错网格进行空间离散.设定整个网格剖分区域的6个面均为规则矩形，采用正交规则矩形网格，在x、y和z三个方向上离散网格单元数为n_x、n_y和n_z.定义网格中心为(i, j, k)，网格单元如图 1所示.定义E在网格棱边中心，三个方向的场点分别为E_{i, j±1/2, k±1/2}^x、E_{i±1/2, j, k±1/2}^y和E_{i±1/2, j±1/2, k}^z，定义B在网格面中心，三个方向的场点分别为B_{i±1/2, j, k}^x、B_{i, j±1/2, k}^y和B_{i, j, k±1/2}^z.

由公式(3)可知，空间离散主要包含两部分内容：旋度算子离散和空间内积离散.旋度算子离散采用积分形式的斯托克斯定理处理.考虑矩形网格单元(i, j, k)，在x、y和z三个方向上的单元网格长度记为h_xi、h_yj和h_zk.单元网格6个表面分别记为S_{i±1/2, j, k}、S_{i, j±1/2, k}和S_{i, j, k±1/2}.以表面S_{i+1/2, j, k}为例，根据中点积分规则，计算电场的旋度x方向的分量，即关于表面S_{i+1/2, j, k}的投影为

(4)

同理可以得到电场旋度y方向和z方向分量的离散表示，整理将旋度算子表示为矩阵形式(Haber and Ruthotto, 2014)

(5)

其中P为包含剖分网格所有表面面积的对角矩阵，L为包含剖分网格所有棱边长度的对角矩阵，C为包含0和±1的矩阵，e是网格中电场E的矩阵表示形式.

空间内积离散包含位于面中心点的内积(B, F)和位于棱边中心的内积(σE, W).根据内积定义公式，采用中点平均即可实现空间内积离散为

(6)

(7)

其中b和f为B和F的矩阵表示形式，w为W的矩阵表示形式，M_μ和M_σ分别为单元体积的差分矩阵，具体形式见(Haber and Ruthotto, 2014).

利用旋度离散和内积离散的具体形式，消去f和w，得到控制方程空间离散的矩阵表示为

(8a)

(8b)

将公式(8a)代人(8b)，消去b，得到关于e的离散控制方程为

(9)

令A=M_σ^-1CURL^TM_μCURL，X=M_σ^-1s，τ=iω，则方程(9)可以整理为

(10)

方程(10)的解可以表示为

(11)

其中 为频率域的传递函数(蒋耀林, 2010)，发射源X为输入，接收电场e为输出.

对于频率域电磁响应的正演计算，一个很重要的问题是外加源项的离散.源附近的场变化剧烈，会显著地影响正演计算精度，目前有多种方法来处理这一问题，包括散射场方法(Newman and Alumbaugh, 1995)，直接离散方法(Weiss, 2013)，差异场方法(Kong, 2008)，虚拟场源校正方法(彭荣华和胡祥云, 2016).本文考虑采用基于有理函数Krylov子空间投影的模型降阶方法求解方程(11)，实现多个频率的可控源电磁响应的快速计算.该算法要求输入X对于不同频率是保持不变的.对于不同频率的正演计算，散射场方法、差异场方法和虚拟场源校正方法都需要计算不同频率的场的解析解，导致离散源项是频率的函数.直接离散方法只是在发射源区域的网格对源进行直接离散和赋值，不需要针对不同频率进行不同的计算，当发射源的长度位置和离散网格固定时，不同频率的发射源X是不变的，本文采用Weiss(2013)的方法对发射源进行直接离散.

1.2 有理函数Krylov子空间模型降阶算法

方程(11)中的传递函数g^τ(A)是频率的函数，将不同的频率代入到方程(11)，即可求得不同频率的正演响应，但由于传递函数的分母中含有大型稀疏矩阵A，直接求解方程(11)是困难的.本文采用有理函数Krylov子空间模型降阶算法，通过构建有理函数Krylov子空间，将大维度的传递函数投影到小维度的krylov子空间，从而实现方程(11)的快速求解.

有理函数Krylov子空间模型降阶算法本质上是一种近似算法.方程(11)的有理函数近似为

(12)

其中r_m(A)为有理函数，其一般形式为r_m=q_m^-1(A)p_m(A)，p_m是m阶非零多项式，(A-ξ_jI)，ξ_i为有理函数的m个极点.给定有理函数的阶数m，得到最小误差的有理函数近似等价为求解A在m维有理函数Krylov子空间Q_m+1=q_m(A)^-1κ_m+1(A, X)上的投影(Börner, 2010).本文采用Gram-Schmidt方法(Ruhe, 1994)构建有理函数Krylov子空间Q_m+1的正交基V_m+1，利用正交基得到矩阵A在有理函数Krylov子空间Q_m+1上的投影矩阵：

(13)

则传递函数的有理函数Krylov子空间模型降阶解为(Güttel, 2010)

(14)

其中e₁=[1, 0, 0, …, 0]_(m+1)×1^T.由于m远小于系数矩阵A的维度，因此求解g^τ(T_m+1)的计算量比求解原传递函数g^τ(A)的计算量要小的多，小维度传递函数可以直接求解.采用方程(14)即可快速求得不同频率的电场响应.有理函数Krylov子空间模型降阶的具体算法如表 1所示.由表 1算法可知，该算法的关键在于：(1)选择最优化极点ξ_j, j=1, …, m；(2)算法的主要计算量在于逐次求解m个线性方程组(A-ξ_jI)x=v_j，如何有效求解这些线性方程组直接决定了该算法的计算效率.

根据有理函数Krylov子空间模型降阶算法可知，待求的m个线性方程组的系数矩阵主要由极点确定.如果选取m个极点都相同，即采用单个重复极点，则待求的m个线性方程组的系数矩阵完全相同，结合直接法求解器PARDISO(Schenk and Grtner, 2004)，只需要1次矩阵分解和m次矩阵回代即可实现m维有理函数Krylov子空间的正交基的构建.

对于最优化单个重复极点选取，根据传递函数的有理函数Krylov子空间投影算法的误差分析理论Güttel (2010)，可以引入关于单个重复极点的收敛率函数：

(15)

其中τ是正虚数，寻找最优化的实数极点ξ₀＜0，保证R在传递函数的参数变化范围T=[τ_min, τ_max]=2iπf_min, f_max内的最小值尽可能地大.对于给定的τ，在ξ＜0时是一个Cayley变换，当ξ沿负实轴(-∞, 0)移动时，ζ在复平面中为一个单位圆的上半圆.当传递函数的参数为单个参数(即求解单个频率时)T=[τ₀]，根据公式(15)可知当ζ₀=i，即ξ₀=iτ时，函数R取最大值：

(16)

当传递函数的参数变化范围为T=[τ_min, τ_max]时，当τ沿T范围变化时，ζ在复平面中为一个单位圆的一段圆弧.根据前面的推导已知ζ₀=i时R取最大值，因此为了保证整个参数变化范围内的R_min尽可能大，最优化单个重复极点选取为

(17)

对应的ζ在复平面中为关于ζ₀=i对称的单位圆的一段圆弧，在τ=ξ₀/i具有最大的收敛率R_max=1+ ，在τ_min和τ_max时具有最小的收敛率：

(18)

其中

对于有理函数Krylov子空间维度m的确定，注意到采用最优化极点的有理函数Krylov子空间模型降阶算法不依赖于模型电导率反差，只依赖于求解频率范围(Güttel, 2010).因此可以通过替代算法(Börner et al., 2015)快速确定有理函数Krylov子空间的维度m，即设置一个对角矩阵W，其本征值足够密集地等间隔分布在足够大的负半轴内，逐渐增大m，分别采用直接计算方法和有理函数Krylov子空间模型降阶算法计算各个频率点对应的传递函数，从而得到有理函数Krylov子空间模型降阶算法的误差，误差不再减小或误差小于给定的误差限时对应的m值，即为有理函数Krylov子空间的维度m.

2 模型计算 2.1 算法验证

首先采用典型的一维CSAMT模型和海洋CSEM模型，通过与解析解的正演结果进行比较，验证本文有理函数Krylov子空间模型降阶算法的计算精度和计算效率.本文计算所有模型三维正演的设备为32 G内存、4核主频3.6 GHz的Intel i7 CPU的台式电脑.

CSAMT方法一般采用接地长导线源作为发射源，在某一固定偏移距的接收点处，接收随频率变化的相互正交的电场和磁场，计算视电阻率和相位差(Routh and Oldenburg, 1999)：

(19)

采用Routh和Oldenburg(1999)的1D层状模型，如图 2所示，为含有高阻和低阻层的5层模型，在背景电导率5×10^-4S·m^-1的半空间中，存在一个顶部埋深50 m、厚度150 m、电导率10^-4S·m^-1的高阻层，顶部埋深200 m、厚度300 m、电导率2×10^-2S·m^-1的低阻层，以及顶部埋深500 m、厚度700 m、电导率2×10^-3S·m^-1的地层.发射源(Tx)为长度1.5 km的线电流源，在偏移距2 km的位置放置接收机(Rx)，接收14个频率(2⁰~2¹³Hz)的电场和磁场，计算不同频率对应的视电阻率和相位差.

采用替代算法(Börner et al., 2015)确定有理函数Krylov子空间模型降阶算法的子空间维度m.引入一个本征值分布足够密集的对角矩阵，本文采用区间[10^-8, 10⁸]内对数等间隔分布的维度为500的对角矩阵A，定义一个维度为500的单位向量X，求解频率范围为[2⁰, 2¹³]Hz，则最优化单个重复极点选取为，根据公式(18)可知，在频率为 时具有最大收敛率，在频率为f_min=2⁰Hz和f_max=2¹³Hz时具有最小的收敛率，为R_min=1.1596.接下来分别采用直接计算方法和有理函数Krylov子空间模型降阶算法计算不同m值时各个频率点对应的传递函数，得到各个频率点的有理函数Krylov子空间模型降阶解的相对误差随子空间维度m的变化曲线如图 3所示.图 3中不同颜色的实线为不同频率对应的相对误差随m的变化曲线，红色虚线为最大收敛率R_max对应的相对误差随m的变化曲线，蓝色虚线为最小收敛率R_min对应的相对误差随m的变化曲线.由图 3可知，所有频率的的相对误差收敛率均位于R_max和R_min之间；在f=2⁶Hz和f=2⁷Hz处具有最大的收敛率，按照公式(18)可以求得其收敛率为2.3648，接近于最大收敛率R_max=2.4141；随着频率逐渐远离f₀=2⁶Hz，收敛率逐渐减小；在f_min=2⁰Hz和f_max=2¹³Hz处具有最小的收敛率R_min.各个频率对应的有理函数Krylov子空间模型降阶解的相对误差随m减小到一定数值后保持不变；在m足够大时，高频相比低频能够得到更小的相对误差，即在m足够大时，最低频率对应的相对误差是所有频率中最大的，因此我们选择以下策略来确定m：选择m为最低频率f_min对应的模型降阶解的相对误差不再减小时的最小m，根据图 3中的黑色实线可知，当m大于130之后，相对误差基本恒定，因此选取m=130.

接下来采用有理函数Krylov子空间模型降阶算法计算上述CSAMT模型，拟态有限体积法进行空间离散，空气层电导率设置为10^-6S·m^-1，采用非均匀网格剖分，最小网格长度为15 m，网格放大系数为1.3，总的网格单元数为31×51×58.同时采用一维正演程序Dipole 1D(Key, 2009)计算该模型.总的正演计算结果如图 4所示，图 4a为视电阻率，图 4b为相位差，其中浅色点线为1D正演计算结果，黑色点线为3D算法的计算结果.图 4c为三维正演计算的视电阻率的相对误差，图 3d为三维正演计算的相位误差.从图 4a和4b可知，三维有理函数近似算法计算的CSAMT响应与解析解吻合的很好.图 4c和4d表明，各个频率的振幅响应相对误差小于6%，相位误差小于1.5°.本文的有理函数近似算法计算所有14个频率正演响应的时间为113 s，即单个频率的平均正演响应时间仅为8 s，说明本文的三维有理函数Krylov子空间模型降阶算法能够快速精确地实现多频CSAMT三维正演响应.

不同于陆地CSAMT方法，海洋CSEM方法接收不同频率沿多个偏移距的电场和磁场响应，一般采用沿测线(inline)方向的电场的振幅和相位作为响应数据.本文采用Constable和Weiss(2006)的1D海洋CSEM层状模型，如图 5所示，空气层电导率为10⁶ S·m^-1，海水层深度为1000 m，电导率为0.3 S·m^-1，地层电导率为1 S·m^-1，在海底下方1000 m处存在一个电导率为0.01 S·m^-1的高阻层.发射源位于海底上方100 m处，发射源长度为100 m，接收机位于海底表面，接收机间距500 m，计算发射频率分别为0.1 Hz、0.25 Hz、0.5 Hz、0.75 Hz、1.0 Hz时沿inline方向的电场振幅和相位响应.

分别采用1D程序Dipole 1D(Key, 2009)和3D有理函数Krylov子空间模型降阶算法计算上述模型.采用3D算法进行正演算法时，空气层电导率设置为10^-6 S·m^-1，采用非均匀网格剖分，最小网格长度为20 m，网格放大系数为1.3，总的网格单元数为102×42×58.最优化极点为-1.9869，对于极点数m选取，采用替代算法(Güttel, 2010)计算最低频率f_min=0.1 Hz对应的模型降阶解的相对误差与m的关系如图 6所示，由图 6可知，当m大于30之后，模型降阶解的相对误差基本恒定，因此选取m=30.

正演计算结果如图 7所示，其中黑色、红色、绿色、蓝色和青色分别为0.1 Hz、0.25 Hz、0.5 Hz、0.75 Hz、1.0 Hz的正演计算结果.图 7a为沿inline方向的电场的振幅，图 7b为沿inline方向的电场的相位，其中实线为1D正演计算结果，圆圈为3D有理函数Krylov子空间模型降阶算法的计算结果.图 7c为3D正演计算的电场振幅的相对误差，图 7d为3D正演计算的电场相位的误差.采用3D有理函数Krylov子空间模型降阶算法，只需要一次正演，即可得到所有频率的正演响应.从图 7可知，3D模型降阶算法计算的海洋CSEM响应与解析解吻合得很好，各个频率的振幅响应相对误差小于5%，相位误差小于3°.3D算法计算上述5个频率正演响应的时间为310 s，即单个频率的平均正演响应时间仅为62 s，说明本文的三维有理函数Krylov子空间模型降阶算法能够快速精确地实现多频海洋CSEM正演响应.

2.2 勘探最优化设计

接下来分析一个典型三维海洋CSEM模型，如图 8所示，海底1 km处存在一个厚度100 m，水平尺寸4 km×4 km的矩形高阻目标体，发射源Tx为发射电流强度为1A的线电流源，位于海底上方100 m处，发射源距离海底目标体水平方向边界3 km，接收机Rx位于海底表面，接收inline方向的电场.

采用三维有理函数Krylov子空间模型降阶算法计算典型的发射频率为0.1 Hz、1 Hz的三维正演响应.采用非均匀网格剖分，最小网格长度为20 m，网格放大系数为1.3，总的网格单元数为112×48×59.最优化极点为-1.9869，对于极点数m选取，由于计算频率范围与计算图 5的一维CSEM模型的频率范围是一致的，因此其极点数仍然选取为m=30.采用模型降阶算法同时计算2个频率的正演响应，其中不含目标体的正演时间为471 s，含有目标体的正演时间为478 s.

正演计算结果如图 9所示，图 9a中虚线为不含目标体时沿inline方向的电场振幅，实线为含有目标体时沿inline方向的电场振幅，绿色虚线为深水域接收器的噪音水平5×10^-16V/(Am²)(Mittet and Morten, 2013)，图 9b为归一化振幅.由图 9可知，发射频率为0.1 Hz时，接收信号幅值大，在0~10 km偏移距范围内的接收信号幅值都大于噪音水平，但异常响应小，最大归一化异常为2；发射频率为1 Hz时，异常响应大，最大归一化异常大于40，但接收信号幅值小，异常响应明显的主要区域，即在偏移距大于5 km区域的接收信号幅值小于噪音水平，而偏移距小于5 km的区域，归一化异常均小于2.因此采用0.1 Hz或1 Hz的发射频率均不能得到最优化的异常响应(即大于噪音水平的信号幅值存在最大的归一化场)，显然存在一个最优化的发射频率，保证接收信号大于噪音水平，同时信号的归一化幅值最大.

为了得到该模型最优化的异常响应，需要模拟不同频率的正演.采用常规正演算法需要多次正演响应，而采用本文有理函数Krylov子空间模拟降阶算法，只需要一次正演，即可得到一个频率范围内所有频率的正演响应.对于该3D模型，本文选取频率范围为0.1~2 Hz，间隔为0.1 Hz的20个频率.分别计算3D模型含有目标体和不含目标体时的响应，网格剖分仍然采用非均匀网格剖分，最小网格长度为20 m，网格放大系数为1.3，总的网格单元数为112×48×59.由于发射频率范围发生变化，最优化极点变化为-2.8099，采用替代算法(Güttel, 2010)计算的最小频率的模型降阶解的相对误差与m的关系如图 10所示，当m大于35之后，模型降阶解的相对误差基本恒定，因此选取m=35.

采用有理函数Krylov子空间模型降阶算法计算20个频率的3D正演响应，不含目标体的正演时间为495 s，含有目标体的正演时间为514 s.对比可知，计算20个频率的正演时间相比计算2个频率(0.1 Hz、1 Hz)的正演时间仅有少量增加(不含目标体的正演时间增加24 s，含有目标体的正演时间增加36 s)，进一步说明了模拟降阶算法适合于计算多个频率的正演响应.

计算所有20个频率沿inline方向电场的的正演响应结果如图 11所示.图 11a为含有目标体相比不含目标体时的电场的归一化振幅，横坐标为偏移距，纵坐标为发射频率，剖面值为归一化电场振幅，图中剔除了含有目标体时测量电场振幅小于噪音水平的数据，即图中的白化区域.由图 11a可以清晰地看出，归一化电场振幅较大(大于4)的区域主要分布在发射频率为0.4~0.8 Hz、偏移距为5.5~8.5 km范围.图 11b为含有目标体相比不含目标体时的电场的相位差，横坐标为偏移距，纵坐标为发射频率，剖面值为相位差，图中同样剔除了含有目标体时测量电场振幅小于噪音水平的数据，即图中的白化区域.由图 11b可知，相位差较大(大于80)的区域主要分布在发射频率为0.1~1.0 Hz、偏移距为5.5~9.5 km范围.

进一步绘制了含有目标体时测量电场振幅的等值线(图 11a和11b中的实线)，偏移距越大，发射频率越高，测量电场的振幅越小.在实际采集数据时，需要保证测量电场的振幅尽可能大，从而保证测量数据质量，因此在确定最优化发射频率和偏移距时，归一化电场振幅、相位差和测量电场的振幅需要综合考虑.

3 结论

本文采用有理函数Krylov子空间模型降阶算法实现了多频可控源电磁法三维正演响应的快速计算.采用基于Yee氏交错网格的拟态有限体积法实现控制方程的空间离散，从而将任意频率的电场响应表示为关于频率参数的传递函数形式.本文提出采用单个重复极点的有理函数Krylov子空间模型降阶算法求解该传递函数，并根据传递函数的有理函数Krylov子空间投影算法的误差分析理论得到了最优化的单个重复极点的计算公式，进一步结合直接法求解器PARDISO，利用Gram-Schmidt方法，将有理函数Krylov子空间模型降阶算法的计算量减少为1次系数矩阵分解和m次矩阵回代，从而实现了多个发射频率的可控源电磁法三维正演响应的快速计算.

通过与多发射频率的CSAMT和海洋CSEM层状模型的一维解析解正演结果对比，说明了本文的三维有理函数近似算法能够快速精确地计算多发射频率的可控源电磁法正演响应.对比CSAMT模型正演和海洋CSEM模型正演可知，当发射频率越多时，本文算法的计算速度优势更显著.分析了一个海洋CSEM三维模型勘探最优化设计的例子，采用本文有理函数Krylov子空间模型降阶算法，只需要2次正演(含有目标体和不含目标体)，即可实现最优化频率选取.进一步说明了有理函数Krylov子空间模型降阶算法能够快速计算一定频率范围的三维正演响应，因此特别适合于勘探前最优化频率选取.