起伏地形大地电磁二维反演
熊彬1, 罗天涯1, 蔡红柱2, 刘云龙1, 吴延强1, 郭胜男1
1.桂林理工大学 地球科学学院,广西 桂林 541006
2.美国 犹他大学 矿业与地球科学学院,犹他州 盐湖城 84112

作者简介: 熊彬(1974-),男,湖北仙桃人,博士,教授,第三届地球电磁专业委员会常务委员会委员,主要从事电磁场理论及反演成像方面的教学与科研工作。E-mail: xiongbin@msn.com

摘要

为了适应地形起伏的实际地质情况,开展带地形的最小二乘二维反演研究。鉴于大地电磁(MT)反演的不适定问题,引入Tikhonov的正则化方法,从而获得关于总目标函数的方程,利用光滑约束最小二乘法求解总目标函数方程。由于正则化因子值与反演精度以及稳定性相关,采用主动约束平衡方法获取最优化的正则化因子,以确保反演精度和稳定性都达到最佳。与此同时,利用电磁场互易定理以节省反演迭代过程求解雅可比矩阵的计算时间。构建了若干地质构造模型进行试算,分别讨论TE、TM模式以及二者联合模式的反演结果,并与前人研究工作对比以说明本文方法的反演效果。

关键词: 起伏地形; 大地电磁; 正则化; 最小二乘; 反演
中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)03-0587-07 doi: 10.11720/wtyht.2016.3.22
Two-dimensional magnetotelluric inversion of topography
XIONG Bin1, LUO Tian-Ya1, CAI Hong-Zhu2, LIU Yun-Long1, WU Yan-Qiang1, GUO Sheng-Nan1
1.College of Earth Sciences, Guilin University of Technology, Guilin 541006, China
2.College of Mines & Earth Sciences, University of Utah, Salt lake city, UT 84112, USA
Abstract

In order to simulate the actual geological conditions, the authors present the least squares inversion by incorporating topography into a forward model. In consideration of an ill-posed inverse problem of MT, the authors introduce Tikhonov regularization to obtain the equation of the total objective function and utilize smoothness-constrained least-squares inversion to solve the total objective function. As the regularized factor controls resolution and stability of the inverse problem, the authors put forward active constraint balancing (ACB) to obtain an optimized regularized factor that balances the resolution as well as the stability of the inversion process. Meanwhile, for the purpose of speeding up the calculation of the field Jacobian for 2-D magnetotelluric inversion, the principle of electromagnetic reciprocity is applied. Finally, the authors discuss the inversion results of TE mode, TM mode and joint inversion of TE and TM mode using some synthetic models, in comparison with some of previous work.

Keyword: topography; magnetotelluric; regularization; least-squares; inversion

大地电磁测深法以天然交变电磁场为场源, 当交变电磁场以波的形式在地下介质中传播时, 由于电磁感应作用, 地面电磁场的观测值将包含地下介质电阻率分布信息, 根据电磁场的集肤效应, 研究大地对天然电磁场的频率响应, 可获得地下不同深度介质电阻率信息[1]。在大地电磁数据解释中, 反演是不可或缺的一部分。随着正演理论方法的发展, 反演也不断取得新的进展, 经历吉洪诺夫和卡尼亚时期的一维反演, 到20世纪70年代开始的二维反演 25, 至今, 三维地质结构的反演 612成为研究的主流。然而三维反演非常耗时, 对计算机内存要求比较高, 此外, 大地电磁野外数据采集周期往往比较长, 勘探成本很高, 因而其在实际应用中并未普及。而对一些地质条件, 二维结构的假设也是合理的, 且二维反演具有计算速度快、数据存储量小等优势, 因此, 深入研究MT二维正反演仍不失为减小勘探成本的有效手段。

常用的大地电磁二维反演方法有OCCAM反演、REBOCC(Reduced basis OCCAM)反演、NLCG(nonlinear conjugate gradients)反演、快速松弛反演(rapid relaxation inversion, RRI)等。OCCAM二维反演算法是deGroot-Hedlin等[13]由OCCAM一维反演方法扩展而来的, 其中正演模拟方法为有限单元法, 在反演计算中, 每一次迭代都将求解模型参数的偏导数矩阵, 因此, 反演过程会随着测点数的增加以及网格剖分的加密而导致计算时间变长; REBOCC是Siripunvaraporn等[14]在OCCAM二维反演方法的基础上改进得来, 采用模型参数的协方差矩阵将反演问题从M维模型参数空间转变为N维数据空间, 由于NM, 使得反演计算时间大大缩减; SBI(sharp boundary inversion)反演方法是deGroot-Hedlin等[15]在原有的OCCAM反演算法的基础上, 构建新的模型粗糙度矩阵, 提出一种针对陡边界条件的二维反演算法; 为了避免OCCAM法直接线性搜索, Smith和Booker[16]提出RRI方法, 通过分析每个测深点之下的电阻率变化, 将二维反演问题降为一维反演问题, RRI法不直接求雅可比矩阵, 而是对正演求得的电场值做积分运算, 获得视电阻率对模型参数的偏导数, 每次迭代只要做一次正演, 提高了计算速度, 但在反演时若某些参数控制不好, 便得不到理想结果, 甚至不能收敛[17]。大地电磁正则化反演的传统做法是用光滑函数表示模型参数, 通过求解最小范数或者光滑稳定的函数, 即可得到光滑的反演模型, 但对于地质体界面分辨率较低。基于此, Portniaguine和Zhdanov[18]提出了聚焦反演方法, 并成功应用于重力反演中; 随后, Mehanee和Zhdanov[19]将聚焦反演方法引入大地电磁二维反演问题的求解, 且取得很高的反演精度; NLCG反演算法是Rodi等[20]提出的一种快速、稳定、收敛的二维反演方法, 其正演模拟方法是基于有限差分法, 反演结果具有较高的分辨率, 计算所需内存小, 另外, NLCG反演方法可以进行联合反演, 如TE和TM视电阻率和阻抗相位的联合。然而, NLCG对初始模型的选择存在依赖性, 正则化参数选择需要凭经验输入。

在前人研究的基础上, 针对反演分辨率、反演过程稳定性以及偏导数计算等问题, 将求解病态问题的正则化方法应用到最小二乘优化方法中; 同时, 用电磁场互易定理来计算雅可比矩阵, 最后, 对若干典型模型进行计算分析。

1 有限元正演

在直角坐标系中, 令地质构造走向沿y轴无限延伸, 假设电、磁场V(x, z)在求解区域Ω 内满足微分方程DV=f, 在Ω 边界∂ Ω上满足给定的边界条件V=f0; D为微分算符, f为关于场源的函数, f0为已知函数, 根据变分法可知, 求解偏微分方程相当于寻找使得拉格朗日函数l最小的场值V, 拉格朗日函数l表示为

l=ΩL(V, xV, zV)dxdz, (1)

L为拉格朗日密度。当场值V满足Euler-Lagrange方程

zL(zV)+xL(xV)=LV2

时, 拉格朗日函数取得最小值, 拉格朗日密度L可表示为

L=12η(zV)2+12η(xV)2+γ2V2, (3)

式中, η γ 的含义见文献[21]。将式(3)代入式(2)即可获得电、磁场满足的具体微分表达式。加入给定的边界条件获得边值问题以及用有限单元法求解边值问题的详细过程参见文献[21-24]。

2 计算辅助场

在获得V之后, 为了求解TM、TE模式的阻抗

Zxy=ExHy=IV,  Zyx=EyHx=-VI, (4),

还需求出辅助场I

zV=-τI, (5)

式中, I= ExTMI=- HxTE。沿用Lee等[24]的算法, 地表场值沿xz轴的导数可表示为其法向和切向分量导数与地表倾斜角度θ 的正、余弦函数的组合

VxVz=xnznxtzt-1VnVt=sinθ-cosθcosθsinθ-1VnVt, (6)

利用含有地表节点的单元构造方程

KsVs=Vns, (7)

式中, Vs=( V1s, V2s, …, VNxs)T, Vns= V1sn, V2sn, , VNxsnT, Ks是包含地表节点单元的系数矩阵, Vs是地表节点的场值, VnsVs的法向导数, Nxx方向上地表节点数, 从而可以求出场值的法向导数∂ V/∂ n。切向导数∂ V/∂ t可利用地表沿x方向相邻节点的差商求取。将∂ V/∂ n∂ V/∂ t代入式(6)即可解出∂ V/∂ z并代入式(5)可求出辅助场I, 进而可通过式(4)获得地下阻抗信息。

3 最小二乘光滑约束反演
3.1 基本公式

针对大地电磁反演的不适定性, 引入Tikhonov的正则化方法

Pα(m)=φ(m)+αs(m), (8)

其中, Pa(m)为总目标函数; φ (m)为观测数据与预测数据之差的平方和, α 为正则化因子, s(m)为模型约束目标函数。根据不同的先验约束条件, 模型约束目标函数的构建方法有多种, 常用的3种是最小模型约束、最平缓模型约束、最光滑模型约束[25], 文中采用基于先验模型的最光滑模型约束。因此, 反演问题的总目标函数可表示为[26]

Pα(m)=Wd[dobs-F(m)]2+αWm(m-mref)2, (9)

式中, Wd为观测数据权系数矩阵, dobs为观测数据, F为正演响应函数, Wm为光滑度矩阵, m为模型参数, mref为先验模型。经过适当的运算, 可得到模型修改量Δ m的表达式

Δm=(JTkWTdWdJk+αWTmWm)-1·[JTkWTdWd(dobs-dk)+αWTmWm(mref-mk)], (10)

式中, J为雅克比矩阵, mk为模型的第k次迭代值, dk=F(mk)。

3.2 雅可比矩阵的计算

考虑到雅可比矩阵求解需要耗费大量时间的问题, 文中利用电磁场的互换定理[27]来求解雅可比矩阵。此外, 反演过程中, 针对电阻率和视电阻率很大的动态范围和便于在目标函数中采用均方误差, 将模型参数表示为对数形式[24]:

m=(logρ1, logρ2, , logρM)T, (11)

同时也将视电阻率转换成其对数形式, 在二维地质结构中, TM模式阻抗相位角0° < ϕ < 90° , TE模式阻抗相位角-180° < ϕ < -90° , 为了消除两种极化模式联合反演过程中数值的差异性, 将TE模式的阻抗相位角加上180° 。

将观测数据表示为

d=(logρa, 1, , logρa, N, ϕ1, , ϕN)T, (12)

其中, ρ a, iϕ i分别表示第i个视电阻率和相位, 则第i个视电阻率和第i个相位关于第k个模型参数的偏导数分别表示为

logρa, ilogρk=2ωμReZiρkZi* ρkρa, i, (13)ϕilogρk=ρkcos2ϕiImZiρk-tanϕiReZiρkRe(Zi), (14)

式中Zi是第i个阻抗, * 表示共轭复数。同理可得

Ziρk=Eiρk-ZiHiρk/Hi, (15)

式中, EiHi分别是第i个电场和磁场。

3.3 正则化因子选取

正则化因子α 的取值将会决定反演过程的稳定性和解的精度, α 很大主要拟合先验约束信息, 反演过程稳定性得到提高但是反演结果分辨率降低; α 很小则主要拟合实际观测数据, 但可能会引起反演过程发生不稳定现象。为了权衡反演稳定性及精度, 采用主动约束平衡(active constraint balance, ACB)方法确定正则化因子, ACB方法主要是根据参数分辨率矩阵R和扩展函数SP自动选取最优化正则化因子[24, 28]。参数分辨率矩阵R由伪逆矩阵J-g和雅可比矩阵J的乘积获得

R=J-gJ, (16)J-g=(JTJ+λCTC)-1JT, (17)

根据式(16), 扩展函数可写成

其中:M是反演参数总个数; wij是权重因子, 由第i个和第j个模型参数的距离来定义; Sij是决定是否把约束条件或者正则化因子强加于第i个模型参数和其临近的模型参数的因子, 当Cij=0, Sij=0, Cij≠ 0, Sij=1, 扩展函数取得一个较大的值, 说明此参数的分辨率低, 则应该分配较大的正则化因子[28]

根据扩展函数和模型参数分辨率的反比关系, 且由于扩展函数随模型参数的坐标按指数形式变化, 因而对正则化因子取对数, 第i个模型参数的正则化因子λ i(x, z)的对数值表示为[24]

log(λi)=log(λmin)+log(λmax)-log(λmin)log(SPmax)-log(SPmin)×{log(SPi)-log(SPmin)}, (19)

式中, SPminSPmax分别是扩展函数的最小值和最大值, λ minλ max分别是正则化因子的最小值和最大值, 由此方法可自动给分辨率高的模型参数分配较小的正则化因子, 这相当于在反演过程中, 扩展函数SPi取得较小的值, 反之, 则给分辨率低的模型参数分配较大的正则化因子。

4 算例分析
4.1 正演算法验证

为了验证本文正演算法的有效性, 首先计算层状介质模型响应并与解析解对比。构建一个H型地电模型:h1=1.2 km, ρ 1=1 000 Ω · m, h2=1.41 km, ρ 2=100 Ω · m, ρ 3=1 000 Ω · m。采用四边形网格剖分, 正演结果如图1所示。从图中可知, 数值解和解析解吻合情况令人满意, 平均相对误差均小于0.5%。

图1 数值解与解析解对比

4.2 模型2

参照Mehanee和Zhdanov[19]所用模型(图2)。在均匀半空间中赋存两个异常体, 沿测线以240 m的间距布设51个测站, 正演频率取0.01, 0.03, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 100, 300, 1 000 Hz共14个频点。初始模型取电阻率为50 Ω · m的均匀半空间, 先验信息为零。从反演结果(图3~图5)可以得出, 两异常体顶部、底部埋深, 以及沿测线方向上的展布与理论模型基本一致。同时三种反演效果还存在一些细微的差异, TE模式对高阻体的反映在纵向上分辨率较TM模式和TE+TM联合模式稍低, 但对低阻体而言, 三者的反演精度均较高。由此可知, 相对于高阻体, 低阻体的反演分辨率较高; TE+TM联合反演精度比TE或者TM模式单独反演精度高。在笔记本计算机(IBM T60, Genuine Intel(R) CPU T2400@1.83GHz)上, TM模式反演时间不足5 min, 与TE及TE+TM联合反演所需时间相当, 数据拟合差如图6所示。

图3 模型2TE模式反演结果

图4 模型2TM模式反演结果

图5 模型2TE+TM联合反演结果

图6 模型2数据拟合差

Mehanee和Zhdanov[19]利用的是聚焦正则化反演方法, 他们引入最小支撑泛函的稳定器并与惩罚泛函相结合, 对模型参数变化大和不连续的区域用一个新的稳定泛函来描述, 使得目标函数集中到最小的面积, 更好的判断地质体界面的存在。图7为他们的TE+TM联合反演结果, 异常体与围岩之间几乎没有过渡带; 而在本文中, 从异常体到围岩, 电阻率值是一个渐变过程, 这是由于这里采用的是光滑函数表示模型参数的变化。因此, 我们应该结合先验信息, 选择适当的函数来描述模型参数的变化, 来获得更加符合实际情况的反演结果。

图7 模型2TE+TM联合聚焦反演

4.3 模型3

如图8所示的地质构造, 岩石电阻率取测试样本中的均值[29], 模型的构建及网格化由COMSOL获得, 沿测线布设41个测点, 点距为500 m; 计算频率范围取10-2~103 Hz, 按对数等间隔共截取11个频点, 将基于非结构化有限元算法的正演响应作为本文的拟合数据; 初始模型为100 Ω · m的均匀半空间, 先验信息取零。从图9~图11可得, 三种模式反演结果较为一致, 均反映出浅部低阻砂岩覆盖层, 在其之下, 可以看到明显的断裂带, 断裂带西侧为低阻铁质页岩, 东侧为高阻石灰岩; 起伏盐类岩、硅质页岩和基底花岗岩也得到较明显的反映。同时, 基底花岗岩中出现一些虚假异常; 反演数据拟合差如图12所示。

图9 模型3TE模式反演结果

图10 模型3TM模式反演结果

图11 模型3TE+TM联合反演结果

图12 模型3数据拟合差

4.4 模型4

如图13所示, 模型4为明显的沉积环境, 地层分界线清晰可见; 按照地层运动规律推断, 地表之下三层砂岩为晚期形成的覆盖层, 它们之下的地层可能遭受横向挤压, 引发褶皱构造, 向斜被角砾岩填充, 背斜遭遇风化剥蚀削平。除频点数为5外, 获取拟合数据的其他信息与上例相同; 初始模型为100 Ω · m的均匀半空间, 先验信息取零。从其反演结果(图14~图16)来看, 异常的分布具有显著的分层现象, 由地表高阻砂岩层向低阻砂岩层转换, 再过渡到高阻基底, 由于两低阻角砾岩填充体的存在, 映现出两低阻异常, 同时可以观察到TM模式及TE+TM联合模式反演结果与实际电性结构更为吻合, TE模式横向分辨率不足致使两低阻异常衔接为一体。数据拟合差如图17所示。

图13 模型4示意

图14 模型4TE模式反演结果

图15 模型4TM模式反演结果

图16 模型4TE+TM联合反演结果

图17 模型4数据拟合差

4.5 模型5

如图18所示, 在0 m处有一座按照半周期余弦函数变化的山脊, 山脊顶点海拔为1 km, 山脊跨度 5 km, 距水平地表1.4 km处发育一低阻块, 宽2.5 km, 高1 km, 电阻率为10 Ω · m, 围岩电阻率为100 Ω · m。测点、频率信息与模型3相同, 反演结果如图19~图21所示。由于地形起伏, 大地电磁场受到极大扰动[30], 而此处反演结果基本消除了地形的影响, 围岩中虚假异常较少, 尤其是TE+TM联合反演, 仅仅在异常体正下方有一不明显的异常。三种模式反演所得异常体的电性以及空间展布与真实模型几乎一致, 数据拟合差如图22所示。

图18 模型5示意

图19 模型5TE模式反演结果

图20 模型5TM模式反演结果

图21 模型5TE+TM联合反演结果

图22 模型5数据拟合差

5 结论

1) 本文方法对地形的适应性较好, 反演结果基本消除了地形的影响, 异常体的信息明显, 在围岩中虚假的异常较少。

2) 在反演分辨率方面, TE+TM联合反演的分辨集合了TE模式和TM模式的优点, 弥补二者的不足之处, 这也反映了具备更丰富的地球物理信息, 可以使反演结果与真实模型更为接近。

3) 对模型进行光滑约束时, 在电性突变区域, 反演结果表现出渐变电阻率过渡带, 这导致在电性结构比较复杂的地质条件中, 一些比较小的异常体, 譬如模型3中尖灭的盐类岩, 其异常信息没有清晰呈现。

The authors have declared that no competing interests exist.

参考文献
[1] 石应骏, 刘国栋, 吴广耀, . 大地电磁测深法教程[M]. 北京: 地震出版社, 1985. [本文引用:1]
[2] Jupp D L B, Vozoff K. Two-dimensional magnetotelluric inversion[J]. Geophys. J. R. astr. SOC. , 1977, 50: 333-352. [本文引用:1]
[3] Ogawa Y, Uchida T. A two-dimensional magnetotelluric inversion assuming Gaussian static shift[J]. Geophysical Journal International, 1996, 126: 69-76. [本文引用:1]
[4] Didana Y L, Thiel S, Heinson G. Magnetotelluric imaging of upper crustal partial melt at Tendaho graben in Afar, Ethiopia[J]. Geophysical Research Letters, 2014, 41(9): 3089-3095. [本文引用:1]
[5] Martí A. The role of electrical anisotropy in magnetotelluric responses: from modelling and dimensionality analysis to inversion and interpretation[J]. Surv. Geophys. , 2014, 35(1): 179-218. [本文引用:1]
[6] Newman G A, Alumbaugh D L. Three-dimensional magnetotelluric inversion using non-linearconjugate gradients[J]. Geophys. J. Int. , 2000, 140: 410-424. [本文引用:1]
[7] Siripunvaraporn W, Egbert G, Lenbury Y, et al. Three dimensional magnetotelluric inversion: data-space method[J]. Physics of the Earth and Planetary Interios, 2005, 150: 3-14. [本文引用:1]
[8] Siripunvaraporn W. Three-Dimensional magnetotelluric inversion: An introductory guide for developers and users[J]: Surv. Geophys. , 2012, 33(1): 5-27. [本文引用:1]
[9] Zhang L L, Koyama T, Utada H, et al. A regularized three-dimensional magnetotelluric inversion with a minimum gradient support constraint[J]. Geophys. J. Int. , 2012, 189(1): 296-316. [本文引用:1]
[10] Grayver A V. Parallel three-dimensional magnetotelluric inversion using adaptive finite-element method. Part I: theory and synthetic study[J]. Geophys. J. Int. , 2015, 202(1): 584-603. [本文引用:1]
[11] Lindsey N J, Newman G A. Improved workflow for 3D inverse modeling of magnetotelluric data: Examples from five geothermal systems[J]. Geothermics, 2015, 53: 527-532. [本文引用:1]
[12] Kordy M, Wannamaker P, Maris V, et al. 3-D magnetotelluric inversion including topography using deformed hexahedral edge finite elements and direct solvers parallelized on SMP computers — Part I: forward problem and parameter Jacobians[J]. Geophys. J. Int. , 2016, 204(1): 74-93. [本文引用:1]
[13] deGroot-Hedlin C, Constable S. Occam's inversion to generate smooth, two-dimensional models from magnetotelluric data[J]. Geophysics, 1990, 55(12): 1613-1624. [本文引用:1]
[14] Siripunvaraporn W, Egbert G. An efficient data-subspace inversion method for 2-D magnetotelluric data[J]. Geophysics, 2000, 65(3): 791-803. [本文引用:1]
[15] deGroot-Hedlin C, Constable S. Inversion of magnetotelluric data for 2D structure with sharp resistivity contrasts[J]. Geophysics, 2004, 69: 78-86. [本文引用:1]
[16] Smith J T, Booker J R. Rapid inversion of two and three dimensional magnetotelluric data[J]. Geophys Res, 1991, 96(B3): 3905-3922. [本文引用:1]
[17] 陈向斌, 吕庆田, 张昆. 大地电磁测深反演方法现状与评述[J]. 地球物理学进展, 2011, 26(5): 1607-1619. [本文引用:1]
[18] Portniaguine O, Zhdanov M S. Focusing geophysical inversion images[J]. Geophysics, 1999a, 64: 874-887. [本文引用:1]
[19] Mehanee S, Zhdanov M. Two-dimensional magnetotelluric inversion of blocky geoelectrical structures[J]. Journal of Geophysical Research, 2002, 107(B4): EPM 2-1-EPM 2-11. [本文引用:3]
[20] Rodi W L, Mackie R L. Nonlinear conjugate gradient algorithm for 2-D magnetotelluric inversion[J]. Geophysics, 2001, 66(1): 174-187. [本文引用:1]
[21] Rodi W L. A technique for improving the accuracy of finite element solutions for magnetotelluric data[J]. Geophys. J. Int. , 1976, 44(2): 483-506. [本文引用:1]
[22] 陈乐寿. 有限元法在大地电磁场正演计算中的应用及改进[J]. 石油勘探, 1981, (3): 84-103. [本文引用:1]
[23] Zienkiewicz O C. The finite element method in engineering science[M]. New York: McGraw-Hill, 1971. [本文引用:1]
[24] Lee S K, Kim H J, Song Y, et al. MT2DInvMatlab-A program in MATLAB and FORTRAN for two-dimensional magnetotelluric inversion[J]. Computers & Geosciences, 2009, 35(8): 1722-1734. [本文引用:4]
[25] 陈小斌, 赵国泽, 汤吉, . 大地电磁自适应正则化反演算法[J]. 地球物理学报, 2005, 48(4): 937-946. [本文引用:1]
[26] 柳建新, 童孝忠, 郭荣文, . 大地电磁测深勘探——资料处理、反演与解释[M]. 北京: 科学出版社, 2012. [本文引用:1]
[27] de Lugão P P, Wannamaker P E. Calculating the two-dimensional magnetotelluric Jacobian in finite elements using reciprocity[J]. Geophys. J. Int. , 1996, 127(3): 806-810. [本文引用:1]
[28] Yi M J, Kim J H, Chung S H. Enhancing the resolving power of least-squares inversion with active constraint balancing[J]. Geophysics, 2003, 68(3): 931-941. [本文引用:2]
[29] 李磊. 湘南骑田岭锡铅锌多金属矿区岩矿石电性研究[J]. 物探与化探, 2007, 31(S1): 77-80. [本文引用:1]
[30] Wannamaker P E, Stodt J A, Rijo L. Two-Dimensional Topographic Responses in Magnetotellurics Modeled Using Finite Elements[J]. Geophysics, 1986, 51(11): 2131-2144. [本文引用:1]