1948年，世界上首套固定翼航空电磁系统在加拿大北湾试飞成功，标志着航空电磁测量应用的开始。1954年，频率域航空电磁法(frequency-domain airborne electromagnetic, FAEM)应用于实际的矿产勘探中，成功地在加拿大New Braunswick省勘探到Health Steele矿床^[1]. 经过几十年的不断优化与改进，该方法愈发成熟，广泛应用于地质调查^[2-4]、矿产资源勘探^[5]、水资源环境调查^[6-7]及危险废弃物场地环境勘察^[8-9]等诸多领域.

目前，航空电磁二、三维正演较成熟，三维反演研究取得了很好的成效，但在实际运用中，由于航空电磁具有数据量大、反演对计算机硬件要求高等特点，暂未完全实用化^[10]. 航空电磁数据的一维反演是目前最实用的方法，常用的频率域航空电磁一维反演和快速成像方法主要有质心深度成像法^[11-12]、电导率-深度成像法^[13]等. 该类方法的优点是能够较快地获得地电参数，缺点是成像分辨率较低，准确性较差^[14].

21世纪迎来了机器学习(machine learning, ML)发展的热潮，机器学习在图像识别、语音识别、人脸识别等诸多领域均取得巨大的成功^[15-18]. 在航空电磁数据的反演中，以深度神经网络为代表的ML方法具有大数据非线性反问题的显著逼近优势，受到科研人员的密切关注. 廖晓龙等^[19]利用改进粒子群算法并结合深度神经网络模型对FAEM数据进行反演，避免了传统反演方法引起的病态性问题. Ahl^[20]利用神经网络实现不需要先验信息的FAEM数据的自动反演. Puzyrev等^[21]利用卷积神经网络，对一维频率域和时间域航空电磁数据进行联合反演，克服了传统线性方法过度依赖于初始模型的弊端，实现了完全由数据驱动的、端到端的反演方法. 以上方法虽然都取得了较好的应用效果，但仍存在不足：1）神经网络以经验风险最小化为优化目标，致使模型易陷入局部极小值，无法统观全局；2）在处理小规模样本数据集时，通过复杂的训练过程拟合数量有限的样本，致使训练得到的模型泛化能力不强；3）构建庞大的正演数据集需要大量的计算成本，耗时较长^[22].

支持向量机(support vector machine, SVM)是较传统的具有扎实理论基础及较强泛化性能的ML方法^[23]. SVM能够克服地球物理传统线性反演方法计算时间过长、求解稳定性较差的弊端，避免一些传统非线性反演方法的超参数调参问题，提高了计算速度，增强了模型的鲁棒性. SVM较深度神经网络的优势是所需的训练样本较少^[24]，节省了计算成本. 诸多科研工作者将SVM应用于地球物理数据处理，取得了突破性的进展. 钟仪华等^[25]结合主成分分析与最小二乘支持向量机，提出岩性识别的方法，该方法发挥了SVM擅长解决小样本数据集、非线性分类预测问题的优势，大大提高了岩性识别的准确率. Li等^[26]利用基于改进粒子群的支持向量机对沙漠地震随机噪声进行降噪，该方法的去噪效果明显优于小波去噪等传统方法，实际案例表明了该方法的有效性. 在处理回归问题方面，Li^[27]将支持向量回归(support vector regression, SVR)应用于海洋地震数据的多层波去除，提高了成像精度. 唐长江^[28]利用多输出支持向量回归对时间域航空电磁数据进行反演，取得了较好的反演效果. 以上研究为基于SVR的FAEM反演带来了新的应用前景.

本文提出基于多输出最小二乘支持向量回归(multiple-output least square support vector regression, MLS-SVR)的频率域航空电磁数据一维反演方法. 利用频率域航空电磁一维正演构建样本数据集；将航空电磁响应数据作为模型的输入端，将地电参数作为模型的输出端，搭建MLS-SVR框架，利用网格搜索法及K-折交叉验证法进行参数寻优；将训练得到的最优参数应用于模型中，对理论模型进行反演预测，与单输出支持向量回归(single support vector regression, S-SVR)和多输出支持向量回归(multiple-output support vector regression, M-SVR)的反演方法进行对比. 为了验证本文方法的有效性，针对实测航空电磁数据进行反演，与Occam反演方法及奇异值分解反演方法进行对比.

频率域航空电磁一维反演是根据地电模型的电磁场数据，求取相应的一维地电模型参数，因此必须确定电磁场数据与一维地电模型间的映射关系，使之能够根据给定的地电模型参数计算相应的电磁场数据，即正演模拟计算. 利用正演模拟构建样本数据集，正演计算是实现MLS-SVR反演至关重要的基础.

目前，频率域航空电磁测量系统按照飞行器类型，可以分为固定翼系统与直升机系统. 按照发射线圈与接收线圈的相对位置，分为水平共面(horizontal co-plane, HCP)与垂直同轴(vertical coaxial, VCX) 2种类型^[29]. 正演模拟所用飞行器类型为固定翼类型，装置类型为HCP(见图1). 由Nabighian^[30]的推导可知，垂直磁偶源接收线圈的垂直磁场强度表达式为

(1) $ \begin{split} {H_{\rm{T}}}{\rm{ = }}&\frac{m}{{4{\text{π}}}}\int_0^\infty {\left[{\rm{exp}} \;{\left( { - \lambda \left( {z - h} \right)} \right) + } \right.} \\ & \left. {{\gamma _{{\rm{TE}}}}{\rm{exp}}\;\left( {{\mu _0}\left( {z - h} \right)} \right)} \right]\frac{{{\lambda ^3}}}{{{\mu _0}}}{J_0}\left( {\lambda \gamma } \right){\rm{d}}\lambda . \end{split}$

式中：

(2) $ \begin{split} &{\gamma _{{\rm{TE}}}} = \frac{{{Y_0}\, {- }\, {{\hat Y}_1}}}{{{Y_0} + {{\hat Y}_1}}},\\ &{Y_n} = \frac{{{u_n}}}{{{z_n}}} = \frac{{{u_n}}}{{{\rm{i}}\omega {\mu _n}}}, \end{split}$

(3) $ {\hat Y_n} = {Y_n}\frac{{{{\hat Y}_{n + 1}} + {Y_n}\tanh \; \,({\mu _n}{h_n})}}{{Y{}_n + {{\hat Y}_{n + 1}}\tanh \; \,({\mu _n}{h_n})}},\,\,{\hat Y_N} = {Y_N}, $

(4) $ u_{n}=\sqrt{\lambda^{2}-\left(\omega^{2} \mu_{n} \varepsilon_{n}-{\rm{i}} \omega \mu_{n} \sigma _{n}\right)}\,. $

m为发射磁矩，是线圈内流经电流、线圈匝数与线圈面积三者之积；γ_TE为TE模式下的反射系数；Y_n为本征面波导纳； $\hat Y_n $为本征阻抗；r为接收线圈与发射线圈2个线圈之间的距离；h为线圈平面与地表之间的距离；ω为发射线圈的发射角频率；μ_n与σ_n分别为第n层地层的磁导率和电导率，通常认为大地磁导率等于真空中的磁导率，即μ_n=μ₀；ε_n为第n层地层的介电常数；J₀为零阶贝塞尔函数。

发射线圈激发的电磁波在空气中传播后由接收线圈接收到的电磁波信号称为一次场，经地面下介质感应后接收到的电磁波信号称为二次场，如图1所示。图中，H₁为一次场磁场强度，H₂为二次场磁场强度。在采集到的频率域航空电磁数据中，H₁远远大于H₂^[10]。在进行反演处理之前，将H₁与H₂进行分离，将分离得到的H₂与H₁相比得到两者的比值，之后乘以10⁶，作为反演所用数据，即^[31]

(5) $ d=\frac{{H}_{{\rm{T}}} - {H}_{1}}{{H}_{1}}\times {10}^{6} . $

式中： $ {H_{\rm{T}}} $为总场磁场强度， $ d $为归一化垂直磁场分量. 在层状介质电磁响应推导公式中出现了零阶贝塞尔函数的广义积分，为了提高运算速度与计算精度，采用快速汉克尔变换数值积分进行计算^[19].

对于样本集 $ ({{\boldsymbol{x}}}_{{m}},{{\boldsymbol{y}}}_{{m}})\in {{\bf{R}}}^{d}\times {\bf{R}}，{m}=1,2, \cdots ,N $，其中N为样本中点的个数，d为维度. 经过非线性映射 $ \varphi ({{\boldsymbol{x}}_{{m}}}) $将样本数据集从低维（d维）平面映射到高维平面（D维）（D>>d）. 在D维空间中，构造线性函数 $ f({\boldsymbol{x}})={\boldsymbol{\omega}}^{\mathrm{T}}\varphi ({\boldsymbol{x}})+b $对数据集进行拟合. 采用最小二乘思想解决以下规划问题^[32]：

(6) $ \left.\begin{array}{l} \min \,\, {\boldsymbol{\omega}}^{\mathrm{T}} {\boldsymbol{\omega}}/2+ c \displaystyle \sum_{m=1}^{N} e_{m}^{2}/2; \\ {\rm{s.\; t.}}\,\, y_{i}={\boldsymbol{\omega}}^{\mathrm{T}}\varphi\left({\boldsymbol{x}}_{m}\right)+b+e_{m}. \end{array}\right\} $

式中：c为惩罚系数， $ {e_m} $为误差变量， $ {\boldsymbol{\omega}} $为映射到的特征空间的权重向量， $ \varphi ({{\boldsymbol{x}}}_{{m}}) $为非线性映射函数， $ b $为阈值. 引入拉格朗日乘子 $ {\alpha _{{m}}} $，得到

(7) $\begin{split} L({\boldsymbol{\omega}}, b, e, \alpha)=\;&\frac{1}{2} {\boldsymbol{\omega}}^{\rm{T}} {\boldsymbol{\omega}}+\frac{1}{2} c \sum_{m=1}^{N}{e}_{m}^{2}-\\ &\sum_{m=1}^{N} \alpha_{m}\left({\boldsymbol{\omega}} \varphi\left({\boldsymbol{x}}_{m}\right)+b+e_{m}-y_{m}\right). \end{split} $

根据KKT条件(Karush-Kuhn-Tucker conditions），解决最优化问题的一种方法，用于求某函数在指定定义域上的全局最小值)求解约束优化问题，可以得到

(8) $ \frac{\partial L}{\partial {\boldsymbol{\omega}}}, \frac{\partial L}{\partial b}, \frac{\partial L}{\partial e_{m}}, \frac{\partial L}{\partial \alpha}=0 \rightarrow\left\{\begin{array}{l} {\boldsymbol{\omega}}=\displaystyle \sum_{m=1}^{N} \alpha_{m} \varphi\left({\boldsymbol{x}}_{m}\right), \\ \displaystyle \sum_{m=1}^{N} \alpha_{m}=0, \\ \alpha_{m}=c e_{m} ,\\ {\boldsymbol{\omega}} \varphi\left({\boldsymbol{x}}_{m}\right)+e_{m}-y_{m}. \end{array}\right. $

消去 $ \boldsymbol{\omega} $及e_m，可得

(9) $ \left[\begin{array}{cc} 0 & \boldsymbol{Q}^{\mathrm{T}} \\ \boldsymbol{Q} & \boldsymbol{P P}^{\mathrm{T}}+c^{-1} \boldsymbol{I} \end{array}\right]\left[\begin{array}{l} b \\ \boldsymbol{A} \end{array}\right]=\left[\begin{array}{l} 0 \\ \boldsymbol{Y} \end{array}\right] . $

式中： $ {\boldsymbol{Q}}=[1,1, \cdots, 1]^{\mathrm{T}} $，维度为N； $ \boldsymbol{P}=\left[\varphi\left({\boldsymbol{x}}_{1}\right), \varphi\left({\boldsymbol{x}}_{2}\right), $ $ \cdots, \varphi\left({\boldsymbol{x}}_{N}\right)\right] $；I为单位矩阵； $\boldsymbol{A}= \left[\alpha_{1}, \alpha_{2}, \cdots, \alpha_{N}\right]^{\mathrm{T}}$； $ \boldsymbol{Y}= $ $ \left[{y}_{1}, {y}_{2}, \cdots, {y}_{N}\right]^{\mathrm{T}}. $可得回归函数：

(10) $ f({\boldsymbol{x}})=\sum_{m=1}^{N} \alpha_{m} K\left({\boldsymbol{x}}, {\boldsymbol{x}}_{m}\right)+b . $

式中： $ K\left({\boldsymbol{x}}, {\boldsymbol{x}}_{m}\right) $为核函数.

最小二乘支持向量回归(least square support vector regression, LS-SVR)区别于传统SVR的特点如下：1）调整约束条件，将不等式约束更改为等式约束，从而将问题转化为线性方程组的求解；2）标准的SVR中Lagrange乘子为非零数值，在LS-SVR中乘子序列 $ {\alpha _m} $与误差序列 $ {e_{{m}}} $成正比，将 $ {\alpha _{{m}}} $称为支持数值谱^[33].

MLS-SVR是利用不敏感代价函数将LS-SVR推广到多输出^[34]，传统的LS-SVR无法同时对多个参量进行输出，在需要对多个参量进行输出时采用的处理方式是训练多个支持向量机模型针对每一个参量进行单独输出，在一定程度上忽略了参数之间可能存在的非线性交叉相关性^[35]，因此对于MLS-SVR的研究十分重要. 给定一组航空电磁数据样本 $M=\{({{\boldsymbol{x}}}_{{m}},{{\boldsymbol{y}}}_{{m}})\}, $ $ {m}=1,2,\cdots ,N,{x}_{i}\in $ $ {{\bf{R}}}^{{d}_{1}},\,{y}_{i}\in {{\bf{R}}}^{{l}}$，其中 $ {x_{{l}}} $和 $ {y_{{l}}} $分别为航空电磁正演数据和地电模型参数. MLS-SVR将正演数据映射到高维特征空间，基于结构风险最小化原理，须解决以下规划问题：

(11) $ \left.\begin{array}{l} \min\,\, \displaystyle \sum_{n=1}^{l}\left\|{\boldsymbol{\omega_{n}}}\right\|^{2}/2+ c \displaystyle \sum_{m=1}^{N} \displaystyle \sum_{n=1}^{l} e_{m, n}^{2}/2+c_{0} \displaystyle \sum_{m=1}^{N} \displaystyle \sum_{n=1}^{l}\left|{e_{m, n}}\right|, \\ \text { s. t. } y_{m, n}={\boldsymbol{\omega}}_{{{n}}}^{\mathrm{T}}\varphi\left({{{\boldsymbol{x}}}}_{m}\right)+b_{n}+e_{m, n}. \end{array}\right\} $

式中： $ {c_0} $为样本的惩罚系数， $ \eta $为样本的总体拟合误差， $ c $为单维输出拟合误差惩罚系数， $ e $为样本的输出误差. 引入拉格朗日乘子：

(12) $ \begin{split} & L({\boldsymbol{\omega}}, b, e, \alpha) = \\ & \sum_{n=1}^{l}\left\|{\boldsymbol{\omega_{n}}}\right\|^{2}/2 + {c} \sum_{m=1}^{N} \sum_{n=1}^{l} e_{m, n}^{2}/2 + c_{0} \sum_{m=1}^{N} \sum_{n=1}^{l} \left|e_{m, n} \right|- \\ &\sum_{m=1}^{N} \sum_{n=1}^{l} \alpha_{m, a}\left(y_{m, n}-{\boldsymbol{\omega_{n}}}^{\mathrm{T}}\varphi ({{\boldsymbol{x}}}_{{m}}) -b_{n}-e_{m, n}\right).\\[-18pt] \end{split} $

根据2.1节的KKT条件求解约束优化问题，消除变量后得到线性方程组：

(13) $ \left. {\begin{array}{l} \left[ {\begin{array}{*{20}{c}} {\rm{0}}&{{{\boldsymbol{Q}}^{\rm{T}}}}\\ {\boldsymbol{Q}}&{{\boldsymbol{T}} + {c^{ - 1}}{\boldsymbol{I}}} \end{array}} \right]\left[ \begin{array}{l} b_n\\ {\boldsymbol{A}}_n \end{array} \right] = \left[ \begin{array}{l} 0\\ {\boldsymbol{Y}}_n + {{{c_0}}}{\boldsymbol{Z}}_n/{c} \end{array} \right],\\ {\boldsymbol{T}} = \left[ {\begin{array}{*{20}{c}} {k\left( {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_1}} \right)}&{k\left( {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_2}} \right)}& \cdots &{k\left( {{{\boldsymbol{x}}_1},{{\boldsymbol{x}}_N}} \right)}\\ {k\left( {{{\boldsymbol{x}}_{2,}},{{\boldsymbol{x}}_1}} \right)}&{k\left( {{{\boldsymbol{x}}_2},{{\boldsymbol{x}}_2}} \right)}& \cdots &{k\left( {{{\boldsymbol{x}}_2},{{\boldsymbol{x}}_N}} \right)}\\ \vdots & \vdots &{}& \vdots \\ {k\left( {{{\boldsymbol{x}}_N},{{\boldsymbol{x}}_1}} \right)}&{k\left( {{{\boldsymbol{x}}_N},{{\boldsymbol{x}}_2}} \right)}& \cdots &{k\left( {{{\boldsymbol{x}}_N},{{\boldsymbol{x}}_N}} \right)} \end{array}} \right]. \end{array} } \right\}$

式中： $ {\boldsymbol{A}}_{n}=\left[\alpha_{1, n}, \alpha_{2, n}, \cdots, \alpha_{N, n}\right]^{\mathrm{T}} $；Y_n、Z_n为N维的列向量. 根据上述方程，可得 $ {\alpha _{{n}}} $、b_n，得到频率域航空电磁数据的多输出最小二乘回归反演算子：

(14) $ f^{(n)}({\boldsymbol{x}})=\sum_{m=1}^{N} \alpha_{m,n} K\left({\boldsymbol{x}}, {\boldsymbol{x}}_{m}\right)+b_{n}. $

式中： $ {\alpha _{{{m,n}}}} $为多输出参量对应的乘子序列， $ {b_{{n}}} $为多输出参量对应的阈值.

在利用SVR反演的过程中，核函数的选择是至关重要的. 核函数的应用能够在不改变原始数据的条件下，将其从低维空间变换至高维空间，实现线性可分. 核函数的使用可以将数据从低维空间转化至无限空间(见图2)，保证了在低维空间计算高维点积^[36]. 选择适应相关问题的核函数有利于优化反演过程，提高反演精度，加快反演速度.

常用的核函数有以下4种：线性核函数(linear kernel)、多项式核函数(polynomial kernel)、Sigmoid核函数(Sigmoid kernel)及高斯径向基核函数(radial basis function, RBF). 与多项式核函数相比，高斯径向基核函数所需要确定的参数数量少，在一定程度上降低了函数的复杂程度；当多项式阶数较高时，核矩阵的元素值会趋近于无穷，高斯径向基核函数会降低计算的困难程度. 对于来自神经网络的Sigmoid核函数，当参数满足特定的条件时，Sigmoid的核是半正定的. 对于频率域航空电磁一维反演的非线性问题，采用高斯径向基(RBF)核函数来进行反演. 高斯径向基核函数的表达公式为

(15) $ K\left({\boldsymbol{x}}_{i}, {\boldsymbol{x}}_{j}\right)=\exp \,\left(-\frac{\left\|{\boldsymbol{x}}_{i}-{\boldsymbol{x}}_{j}\right\|^{2}}{2 \sigma^{2}}\right) . $

式中：σ为高斯核带宽.

采用MLS-SVR求解频率域航空电磁的反问题，需要确定以下2个参数：惩罚系数c及核函数参数 $ \sigma $. 其中，惩罚系数表示训练过程中对错误的容忍度，决定拟合得到f(x)的复杂程度，若c过大，则模型的泛化能力会下降；核函数参数 $ \sigma $影响拟合曲线的光滑程度. 使用网格搜索的方法寻找最佳参数，使用K-折交叉验证确定最优参数的选取准则^[37]，本文选取的K为4. 训练的精度以均方根误差（root mean square error, RMSE）为评价标准，根据均方根误差淘汰不合理的参数组合，产生一组最佳参数，保证均方根误差最小. 将训练样本集随机地分为4个大小相同的子集，任意选出其中1个子集作为预测集，剩余的3个子集作为训练集. 选取一组参数( $ c_{1}^{*}, \sigma_{1}^{*} $)，对训练集进行训练得到拟合函数，用该拟合函数对预测集进行预测，得到RMSE值. 每一个子集都作为测试集，对于选取的这组参数( $ c_{1}^{*}, \sigma_{1}^{*} $)会得到4个RMSE，求平均数 $ \Delta {R_1} $，以反映该组参数对模型的训练效果. 选取另一组( $ c_{2}^{*}, \sigma_{2}^{*} $)参数，重复上述操作得到 $ \Delta {R_2} $. 以此往复，将每一组参数均进行以上操作，选取最小均方误差平均值对应的那组参数作为最优模型参数. RMSE的公式如下：

$ \text { RMSE }=\sqrt{{M^{ - 1}} \sum_{i=1}^{M}\left[{{{{y}}}}_{i}-f\left({\boldsymbol{x}}_{i}\right)\right]^{2}} . $

式中：M为样本数据的个数，y_i为第i个数据的理论值，f(x_i)为第i个数据的预测值.

利用网格寻优法主要有以下2个优点：1）避免启发式或穷举式长时间寻找最优参数；2）可以平行计算，节省运算成本. K-折交叉验证方法可以尽可能地利用样本数据集，除此之外，可以防止过拟合现象的产生.

分别采用S-SVR、M-SVR及MLS-SVR，对层状地电模型进行反演测试. 层状地电模型包含2大类参数：一种是物性参数，即每层地层的电阻率；另一种是位置参数，即每层地层的厚度.

对于2层地电断面模型，依据电阻率可以分为以下2种模式：低-高型(G型)及高-低型(D型). 考虑到目前频率域航空电磁法的有效勘探范围不超过180 m，且线圈发射电磁波频率有限^[38]，正演模拟参数设置如下. 频率参数共设置5个，地电断面模型参数共设置3个，分别为第1层地层电阻率 $ {\;\rho _{\text{1}}} $、第2层地层电阻率 $ {\rho _{\text{2}}} $及第1层地层厚度 $ \mathop h\nolimits_1 $. $ {\;\rho _{\text{1}}} $、 $ {\;\rho _{\text{2}}} $的取值为100~1 000 $ \Omega \cdot \mathrm{m} $，间隔为50，即100, 150 $,\cdots , $ 1 000 $ \Omega \cdot \mathrm{m} $(第1层与第2层电阻率不相等). 第1层地层厚度为15~150 m，间隔为15，即15，30 $,\cdots , $ 150 m. 通过正演计算，共得到3 420个样本数据.

将正演得到的航空电磁数据先取自然对数，再进行归一化处理之后作为输入端. 对于地电断面模型参数的归一化处理不同，由于这些参数均为等间隔取值，且量纲不同，对其进行独立尺度变换，作为SVM的输出端^[39].

将样本数据集分为训练集与测试集，随机取出总样本数据的95%作为训练样本集(3 420×95%=3 249)，剩下的5%作为反演验证集(3 420×5%=171). 超参数的设置是利用网格搜寻法及K-折交叉验证的方法，寻找最优参数. c与 $ \sigma $的取值范围的选取依赖于先验性经验，若步长选取太小，则计算成本过高，训练速度缓慢；若步长选取太大，则易遗漏模型的最佳参数^[40]. 令m = [−5, 9]，n = [−5, 5]，步长均为1；取模型参数c = 2^m， $ \sigma $ = 2ⁿ，由指数函数的性质可知，c与 $ 0 $是非负的. 参数设置完成后进行训练，训练得到各地电断面模型参数的RMSE，如表1所示.

对比表1中3个地电断面模型参数的RMSE可知，其中第1层地层电阻率 $ \mathop \rho \nolimits_{\text{1}} $的RMSE最小， $ \mathop \rho \nolimits_2 $及 $ \mathop h\nolimits_1 $的RMSE较大，说明利用SVR对第1层地层电阻率的反演精确度最高，第2层地层电阻率及第1层地层厚度的RMSE相较于第1层电阻率虽然较大，但反演结果准确. 利用MLS-SVR反演方法，可以降低各参数的RMSE，说明MLS-SVR提高了2层地电断面模型的反演精度.

利用训练完成的模型对反演验证集数据进行反演，由于样本数据集无法一一展示，随机抽取列出5组反演结果，如表2所示. 相对误差的对比结果如表3所示. 分析表中3种方法的相对误差可知，与S-SVR及M-SVR相比，MLS-SVR降低了反演的平均误差. 如图3所示为3种方法对抽取的G型及D型地电模型反演结果的对比图. 图中，ρ为电阻率，D为深度. 从图3可以看出，MLS-SVR的反演结果相较于S-SVR及M-SVR更接近真实模型.

对于3层地电断面模型，共有5个地电模型参数，分别为 $ \;{\rho _{\text{1}}} $、 $ \;{\rho _{\text{2}}} $、 $\; {\rho _{\text{3}}} $、 $ \mathop h\nolimits_1 $、 $ \mathop h\nolimits_2 $. 对这些参数进行如下的设置：3层地层的电阻率均为100~1 000 $ \Omega \cdot {\rm{m}} $，采样间隔均为50，且相邻2个地层的电阻率不等，如 $ {\rho _{\text{1}}} $为100, 150 $,\cdots ,$ 1 000 $ \Omega \cdot {\rm{m}} $. 第1层与第2层地层的厚度为15~150 m，采样间隔均为15，第1层与第2层的厚度相加不得超过150 m. 按照上述方法进行取样，共计29 160组数据. 取其中95%的数据(29 160×95%=27 702)作为训练样本集，以训练得到反演模型，剩余数据(29 160×5%=1 458)作为反演验证集.

数据处理步骤与上面2层模型部分类似，将模型的输出端设置为5个参数，分别是：第1层地层电阻率 $ {\rho _{\text{1}}} $、第2层地层电阻率 $ {\rho _{\text{2}}} $、第3层地层电阻率 $ {\rho _{\text{3}}} $以及第1层地层厚度 $ \mathop h\nolimits_1 $和第2层地层厚度 $ \mathop h\nolimits_2 $. 将各类参数设置完成后开始进行训练. 训练得到最佳参数，对应的RMSE如表4所示.

利用训练完成的模型对反演验证集数据进行反演预测，抽取列出5组反演结果，如表5所示. 相对误差的对比如表6所示.

分析表4可知，浅部地层地电参数的RMSE比深部地层小，说明MLS-SVR的反演精度随着深度降低，但能够较准确地反演出各地电模型参数. 分析表5、6可知，MLS-SVR对3层地电模型的反演精度高于M-SVR和S-SVR. 如图4所示为3种方法对于3层地电模型(A型、H型、K型、Q型)反演结果的对比图. 分析图4可知，MLS-SVR反演结果更接近真实模型.

无论是2层地电断面模型还是3层地电断面模型，SVR对于第1层地层电阻率的反演精度较高，反演结果较精确. 对于物性参数及位置参数，RMSE随着深度的增加而增大，说明反演精度随着深度的增加而降低，但反演结果较准确. 对于2层、3层地电模型，MLS-SVR的反演效果优于M-SVR及S-SVR. 2层地电模型的反演效果优于3层模型.

如图5所示为5层理论模型MLS-SVR反演结果. 可以看出，随着深度的增加，反演精度降低，但能够反映地下空间电阻率的变化趋势. 2层及3层地电模型的反演效果总体上优于5层地电模型.

为了研究噪声对该方法的影响，在正演得到的反演验证集数据中加入3%及5%的高斯随机噪声，利用MLS-SVR反演方法对含噪声数据进行反演，与理论数据的反演结果进行对比. 随机抽取5组2层模型反演结果和5组3层模型反演结果，无噪声及含噪声数据反演结果的相对误差如表7所示（序号1~5对应2层地电模型，序号6~10对应3层地电模型）. 含噪数据反演相对误差较无噪数据反演相对误差略有增大，相对于理论模型层参数的增加幅度较小，表明MLS-SVR反演方法具有较强的鲁棒性.

为了验证该方法的有效性，利用MLS-SVR对某地区近地表的实测航空电磁数据进行反演. 实测数据是由HCP装置测量得到的垂直磁场分量，采集过程中分别设置5个频点，分别为386、1 538、6 257、25 790、100 264 Hz. 考虑到实测数据的频点、电阻率以及深度与理论试验中各地电模型参数不同，且在实测数据的采集过程中夹杂噪声，根据该地区实地踏勘及已有的矿物岩性资料对正演参数进行重新设置，在样本数据集中添加0~5%的随机噪声. 确定好相关参数后进行模型训练，利用训练好的模型进行反演预测. 选取其中4个有代表性的测点，与Occam反演方法及奇异值分解反演方法(singular value decomposition, SVD)进行对比.

如图6所示为殷长春等^[31]利用Occam反演方法对实测数据进行反演. 图中，M为测线起点与终点间的距离. 从反演结果来看，该地区地层大致可以分为3层：第1层电阻率为8~18 $ \Omega \cdot {\rm{m}} $，厚度大致为10 m；第2层电阻率为3~8 $ \Omega \cdot {\rm{m}} $，厚度大致为10 m；第3层的电阻率为8~18 $ \Omega \cdot {\rm{m}} $.

利用MLS-SVR方法来对实测数据进行反演预测，预测输出的地电模型设置为3层地电模型，共有5个输出参数. 如图6所示为A、B、C、D 4个测点的MLS-SVR与Occam反演及SVD反演结果对比图.

图7中的虚线为A、B、C、D（4个点的坐标在图5中依次为(25.7，0)、(157.7,0)、(288.2,0)、(466.1,0)） 4个测点的MLS-SVR反演结果，做标记的实线为Occam反演结果，不做标记的实线为SVD反演结果. 可以看出，MLS-SVR反演出地下电性结果的变化趋势为：“高-低-高”，与Occam反演方法及SVD反演方法得到的电性变化趋势相符，与SVD反演方法相比更加接近真实地质模型，说明MLS-SVR在FAEM数据反演中的可行性. 在实际工程中，对于实测数据的反演往往要求得到更多的地电参数. 随着地电参数的增多，数据集规模呈指数性增大，这给MLS-SVR反演带来了新的挑战.

（1）本文将MLS-SVR应用于频率域一维航空电磁反演，对多个地电模型参数进行同时输出. 将参数的二次规划问题转化为线性方程，简化了运算过程. 理论数据、含噪数据及实测数据的反演都表明了该方法的有效性. 在反演精度上，MLS-SVR优于M-SVR，明显高于S-SVR，表明多输出学习是FAEM数据反演更优的策略.

（2）对于实测数据的反演，利用MLS-SVR能够预测地下电性的简单变化趋势，可以为实际工程应用提供一定的参考，或作为其他反演方法初始模型参数选择的依据. 为了节约计算成本、提高反演精度，应建立与勘探区域地质情况相符的样本数据集，设置MLS-SVR模型的相关参数.

（3）对于更复杂地电模型的反演，样本数据集的规模过于庞大，MLS-SVR对于小规模数据集的训练表现更好. 如何简化样本数据集以适应MLS-SVR的训练以及进一步优化MLS-SVR参数，有待进一步的研究.