地球表面约70%被高导海水覆盖，海陆间明显的电导率差异会产生海洋效应(Geomagnetic Coast Effect，GCE)，从而影响地磁测深资料解释(Parkinson and Jones，1979；Fainberg and Singer, 1980).Kuvshinov等(1990)、Takeda(1993)、Tarits(1994)、Weiss和Everett(1998)、Kuvshinov等(1999)先后根据海陆导电性分布对沿海地磁测深响应的影响进行了模拟，表明海洋效应对周期1—15天之间的地磁响应的影响是不可忽略的.Kuvshinov等(2002)进一步的计算表明海洋效应甚至对沿海周期长达20天的地磁响应产生影响；另一方面，Everett等(2003)用有限元方法模拟了全球地表电导分布对卫星地磁数据的影响，认为海洋效应对2小时至2天的小周期卫星数据的影响也是不可以忽略的.

Kelbert等(2009)在利用全球电磁感应数据研究地幔转换带含水量时也指出，海洋效应会影响地磁数据的反演解释，沿海地区的C-响应需要进行校正.Kuvshinov和Utada(2010)在研究日本沿海地磁台站数据时，发现海洋效应而非地下三维结构是导致日本沿海磁场垂直分量反常的主要原因.基于同样的认识，Shimizu等(2010)利用进行过海洋效应校正的太平洋周围地磁台站数据，并联合海底电缆和长周期大地电磁测深反演获得了更为合理的北太平洋地幔一维电性结构.Semenov和Kushinov(2012)对全球119个地面地磁台站数据进行反演时指出部分台站不合理的C-响应主要由海洋效应造成的.Püthe等(2015)利用卫星数据反演时也考虑了海洋效应，并对C-响应做了海洋效应校正，得到最新的地球一维电性结构.事实上，除了从数值模拟方面模拟海洋效应，Huang和Ikeya(1998)在研究地震电磁信号的传播过程中，从物理模拟实验方面指出了海洋对电磁波信号的影响.

中国东部地区有漫长的海岸线，因此海洋效应同样可能对沿海的地磁观测产生影响.张帆等(2012)研究海岸效应对近海大地电磁数据的畸变作用时指出，当测区与海岸线的距离小于目标频率的趋肤深度时，海洋会严重影响测区电磁场分布；王桥和黄清华(2016)在分析华北13个地磁台站超过5年的地磁感应矢量时空特征时，推测靠近渤海的大连台站可能受到海洋效应的影响.而赵国泽等(2001)在利用大地电磁网法探测长春附近上地幔电性结构研究中，认为日本海和太平洋对远离海岸线的长春地区的长周期大地电磁观测影响可能不大.需要强调的是，Kushinov等(2002)分别使用大陆地幔模型和海洋地幔模型对C-响应进行了海洋效应，结果表明不同模型的校正量存在差异；Semenov等(2012)也提出海洋效应校正量依赖于海洋下方的1-D背景电导率模型.

本文将就海洋效应对中国地区地磁观测影响开展深入研究工作.海洋效应模拟采用“3-D地表电导+1-D”模型(Kelbert，2014).在数值算法上，通过球坐标系下交错网格有限差分技术进行海洋效应的数值模拟.然后利用不同的1-D背景模型，计算了海洋效应对中国地区地磁测深的影响.最后，以广州台站为例，获得了海洋效应校正系数，对受海洋效应影响较大的广州台的观测数据进行比值法校正，对校正后C-响应进行反演，并对反演模型进行了讨论.本文工作的意义在于，当针对全国地磁测深数据利用1-D反演技术研究中国地区地幔电性结构时，本文的工作不仅大致给出了海洋效应的影响范围，还提供了计算海洋效应校正技术手段.

取正谐变场的时间因子e^iωt，地磁测深三维导电模型产生的磁场H满足的亥姆霍兹方程为

(1)

式中，H为磁场强度；ω、ρ和μ₀分别是角频率、电阻率和真空磁导率.假设电导率为均匀各向同性，真空磁导率为4π×10^-7H·m^-1.

地磁测深三维模拟基于球坐标下的交错网格有限差分策略(Uyeshima and Schultz, 2000)，李建平等(2017)对算法的思路和精度进行了验证.在图 1球坐标下，采用正交的曲面坐标对地球进行离散，其中南方向为θ分量，东方向为ϕ分量，地球半径方向r分量.因此，单元的节点编号(i，j，k)对应球坐标三个分量(ϕ，θ，r).磁场切向分量H_ϕ和H_θ分别定义为经度和纬度方向且分别指向东和南，垂直分量H_r定义在半径方向.

图 1

Fig. 1

图 1 球坐标系及磁场分量的定义 节点编号(i，j，k)对应球坐标(ϕ，θ，r)三个方向.磁场切向分量Hϕ和Hθ分别定义为东和南方向.垂直分量Hr定义在半径方向. Fig. 1 The magnetic field defined in a spherical coordinate system Node index (i, j, k) corresponds to a spherical co-ordinate (ϕ, θ, r). The tangential components Hϕ and Hθ of magnetic field are in east and south directions, respectively. The vertical component Hr points to radial direction.

边界条件的设置对于保证计算精度非常重要.空气层的上边界应离地球足够远，以保证地球内部产生的二次磁场在源的位置能被衰减.本文将计算域的外边界设置在距地心10倍的地球半径R_e(R_e=6371 km)的空气层处；计算域的底边界为核幔边界(CMB)，距离地表 2871 km，并且假设地核的电导率为无穷大.在此假设条件下，地核内磁场切向分量为零，但是核幔边界处磁场切向分量不为零.

地磁测深外源的加载非常关键(Sun and Egbert, 2012).在中、低纬度地区，对于中长周期的地磁响应，常常使用P₁⁰项来模拟磁层的等效环形电流(Schultz，1990).本文将环形电流在计算域上边界的值，赋值给磁场切向分量的环形积分，从而模拟合适的源，也就是求解方程(1)的Dirichlet边界问题.

在本文的数值计算中，均采用180×90×98的网格规模，即地球在经纬度方向被剖分为2°×2°网格、半径方向剖分为98层，其中空气层为17层.

为了模拟海洋效应，本文采用“3-D地表电导+1-D”模型，即在1-D模型的顶层采用全球地表三维电导薄层模型.3-D地表海洋陆地电导分布如图 2所示，海陆的平均厚度采用50 km，其中海洋部分由海水和下方沉积层组成.该电导率分布模型和Everett等(2003)及Kelbert等(2014)采用的模型一致.关于全球地表电导薄层的数据来源和详细介绍可参考网站(http://marineemlab.ucsd.edu/~steve/projects/MDAT/sedmap.html).

图 2

Fig. 2

图 2 全球地表 3D模型电导分布图 横坐标是经度，纵坐标为余纬度，色标单位为S. Fig. 2 Map of the global surface conductance map The horizontal axis is longitude and the vertical axis is co-latitude. The color scale is in units of S.

Kushinov等(2002)指出C-响应海洋效应时，高阻岩石圈背景模型的校正结果更好，并且需要2°×2°或者更细的网格剖分.Semenov等(2012)研究表明不同1-D背景电导率模型的海洋校正量差异很大，最后采用了位于全球不同区域的6个台站的电导率模型对不同台站进行海洋校正.

因此，为了研究中国沿海地磁台站的海洋效应，本文分别采用全球1-D电导率模型、中国地区平均1-D电导率模型、兰州地区1-D电导率模型进行海洋校正.全球1-D电导率模型由经过海洋效应校正后的10年卫星数据获得(Püthe et al., 2015).中国地区平均电导率模型由9个非沿海地区且地域分布均匀的地磁台站C-响应数据平均值反演得到，使用到的具体台站信息如表 1，地磁数据观测时间长度从15年到50年不等，但9个地磁台站数据质量较好.C-响应的计算由修改后的BIRRP代码(Chave and Thomson, 2004)计算得到，计算得到的具体的C-响应如图 3所示，其中全球平均C-响应采用Püthe等(2015)的结果.可以看出，全球平均C-响应在周期小于10天时，和中国平均以及兰州台差别较大.中国平均和LZH的C-响应的虚部非常接近，按照ρ^*-z^*原理(Schmucker，1987)，兰州台站和全国平均的电阻率随深度变化应该具有一致的变化规律；相反，两个台站计算的实部存在一定的差异，这种差异反映了中国地区上地幔的横向不均匀性(Semenov and Kushinov, 2012).

表 1 Table 1

表 1 中国平均电导率模型计算中使用的中国地磁台站信息 Table 1 Information of the geomagnetic stations in China 台站名 台站代号 数据采集时间长度/年 满洲里 MZL 1995—2015 长春 CNH 2001—2016 北京 BJI 1957—2013 武汉 WHN 1980—2015 乾县 QIX 1995—2015 兰州 LZH 1980—2012 格尔木 GLM 1995—2015 成都 CDP 1995—2015 通海 THJ 1995—2015

表 1 中国平均电导率模型计算中使用的中国地磁台站信息 Table 1 Information of the geomagnetic stations in China

图 3

Fig. 3

图 3 中国地区平均、全球平均和兰州地区C-响应 Fig. 3 C-responses of China, global and LZH

C-响应反演利用L-BFGS进行1-D反演(李世文等，2016).图 4给出了具体的1-D电导率模型，其中全球1-D模型直接采用Püthe等(2015)的结果.从图可以看出，无论中国地区平均模型还是兰州单站的导电模型，除了在800~1100 km差异较小外，其余深度的导电性明显好于全球平均模型，整体上高于全球模型电导率的2~10倍，这个结论和Karato (2011)关于中国东部电导率变化的规律基本一样.就中国平均模型和兰州模型相比而言，总体的电导率剖面基本规律一样，随深度增加，导电性均匀增加；但在转换带下，导电性呈现不连续的台阶状增加.尽管如此，两个模型最大的差别发生在地幔转换带附近，兰州模型导电性明显好于平均模型，反映出在该深度中国地区地幔转换带导电性的横向不均匀性.此外，由于地磁测深在600~1200 km深度范围内灵敏度较高(Khan and Shankland, 2012)，因此，推测兰州地幔平均1-D电导率高于中国平均电导率是可靠的.

图 4

Fig. 4

图 4 中国地区、兰州和全球1-D地幔电性结构模型 Fig. 4 1-D mantle conductivity profiles of China, LZH and global

对大地电磁而言，当高导海洋存在时，如果测点与海岸线的距离小于低频大地电磁的趋肤深度，海洋的存在会影响海岸附近电磁场分布，具体表现为感应磁场的垂直分量急剧增大，而且感应矢量向海洋方向倾斜(Yang et al., 2010；张帆等，2012).与此类似，对地磁测深C-响应的小周期而言，高导海水中东西向的电流遇到高阻的大陆时会发生严重畸变，其中南北向的畸变电流会产生感应磁场沿经度方向的B_ϕ分量；而沿着大陆海岸线偏南或者偏北弯曲的畸变电流则产生感应磁场沿纬度方向的B_θ分量和沿地球半径方向的B_r分量，从而使C-响应的实虚部发生变化(Ribaudo，2011).图 5显示了考虑海洋效应时，“3-D地表电导+全球1-D”模型周期为1.7天中国地区3-D模型响应C与1-D模型响应C₀的实部和虚部比值的空间分布特征.数据误差对反演的结果影响很大，李墩柱等(2009)指出加入5%的相对误差对大地电磁反演是可以接受的，随着误差增大，反演结果模型的细节差异越明显，尤其是模型越复杂，误差的影响越显著.如果以20%作为地磁数据可接受的数据误差，由图 5看出，实部和虚部20%的相对变化界限基本平行中国东部海岸线自哈尔滨向西南到贵阳一线.1.7天基本是地磁测深的最小周期，同时，在后面的讨论中发现随着周期的增加，海洋效应将逐渐减小.因此，该界限基本界定了海洋效应对中国地磁观测的影响范围.从图 5还可以看见，海洋效应的影响，对于C-响应的实部和虚部的影响程度基本一样.此外，在靠近海洋的广州台、福建的泉州台站和海南的琼中台站处于海洋效应影响的梯度带上，受到的影响尤其明显，说明沿海地磁数据的C-响应必须进行海洋效应校正.

图 5

Fig. 5

图 5 研究区周期1.7天的C-响应的空间分布 (a) C-响应的实部；(b) C-响应的虚部. Fig. 5 Spatial distribution of the C-response at T=1.7 days in China (a) Real parts of the C-response; (b) Imaginary parts of the C-response.

对于受海洋效应影响的C-响应数据，可以采用比值法进行海洋效应校正(Kuvshinov et al., 2002)，具体的校正为公式

(2)

(3)

其中，C^{obs, corr}是校正后的响应，C^obs是观测数据.k称为校正系数，其可以由一维模型响应C^1D和包含地壳电导3D分布的一维模型响应C^1D+sh得到，其中后者利用前面的三维模拟算法得到.

广州台站是中国南方受海洋影响较大的地磁台站之一.图 6a给出了采用图 4中3种不同1-D模型计算得到的广州台海洋效应校正系数.从图可见，总体上实部校正系数大于1，虚部校正系数小于1；就实部而言，中国地区1-D模型校正系数大于兰州1-D模型，并大于全球1-D模型校正系数；对于虚部，中国地区1-D模型校正系数大于兰州1-D模型，并大于全球1-D模型校正系数.就实虚部而言，虚部校正系数明显大于实部，表明虚部受海洋效应的影响更大；同时，就影响周期而言，在周期大于20天后，所有模型的校正系数都接近1，表明中国沿海地区的海洋效应影响周期可能达到20天.

图 6

Fig. 6

图 6 广州地磁台站海洋影响校正 (a)不同模型校正系数k；(b)实测数据不同模型海洋校正. Fig. 6 Ocean effect correction of C-responses for GZH (a) Correction coefficient k of different model; (b) Correction of the observed data.

广州台地磁数据采集时间范围为1960—2015年.用BIRRP程序处理后获得的C-响应绘制在图 6b中.从响应曲线看，周期大于40天数据误差大；而在小周期时，响应的实部明显减小，而虚部受感应作用明显变大.根据电阻率和C-响应的对应关系(Schmucker，1987)，虚部和电导率成正比，显然，如此强的响应，对应于地幔浅部极其高的导电性，就1-D模型而言这不符合实际情况.因此，广州台站的C-响应必须进行海洋效应的校正.

为此，采用图 6a中的三种校正系数对原始的C-响应进行校正，校正后响应曲线示于图 6b.从图可见，一个明显的特点是，使用全球1-D模型校正的C-响应虚部在小周期明显小于中国地区1-D模型和兰州1-D模型的校正结果；而使用中国地区1-D平均模型的C-响应校正结果和使用兰州1-D模型差别不大.在细节上，就实部来说，1-D模型的选择对响应的海洋效应校正影响不大，但虚部相对而言，在更小的周期上，受校正模型的选择影响较大.因此，就中国沿海地磁观测而言，可以选择中国平均模型或者兰州模型作为1-D模型进行海洋效应校正.表 2给出了广州台站原始C-响应和用中国平均1-D模型海洋效应校正后的具体C-响应数据.

表 2 Table 2

表 2 广州地磁台站(GZH)1960—2015年观测数据海洋校正前和校正后的C-响应 Table 2 The uncorrected and corrected C-responses of GZH from 1960—2015 Period/days uncorrected Re(C)/km uncorrected Im(C)/km δC/km coh2 corrected Re(C)/km corrected Im(C)/km 2.667 866.8164 -723.867 48.56 0.27 1092.488 -428.619 3.556 1044.483 -566.462 44.84 0.40 1240.262 -360.204 5.333 1220.044 -438.618 50.89 0.48 1366.486 -303.77 7.111 1268.078 -312.825 46.16 0.48 1377.174 -227.21 10.667 1279.916 -303.305 54.37 0.50 1349.185 -233.727 14.222 1301.194 -281.104 59.32 0.57 1352.288 -224.962 21.333 1339.532 -343.828 85.49 0.62 1372.718 -289.266 28.444 1351.619 -372.168 81.26 0.52 1376.853 -323.591 42.667 1375.564 -327.441 150.63 0.54 1392.429 -296.423 56.889 1399.341 -426.236 168.38 0.58 1413.284 -395.116 85.333 1450.794 -403.941 206.95 0.54 1451.772 -383.968 113.778 1502.952 -434.615 308.53 0.32 1511.196 -419.002

表 2 广州地磁台站(GZH)1960—2015年观测数据海洋校正前和校正后的C-响应 Table 2 The uncorrected and corrected C-responses of GZH from 1960—2015

海洋效应影响校正前后的C-响应采用有限内存拟牛顿方法(L-BFGS)进行一维反演(李世文等，2017).在进行多次反演试验后，图 7给出了广州地磁台站C-响应最终反演1-D模型及拟合情况.从图 7a的数据拟合情况看，未进行海洋校正反演模型的C-响应的实部和虚部在小于10天的周期上基本无法拟合原始数据，而几乎和校正后的C-响应重合；另一方面，海洋效应校正后的数据无论是实部还是虚部都拟合很好.

图 7

Fig. 7

图 7 广州地磁台站海洋校正前后的一维反演 (a)模型的C-响应；(b)一维电导率模型；(c)拟合差；(d)和北方大陆平均电导率模型的对比. Fig. 7 The 1-D inversion results for GZH with and without correction (a) C-response of the inversion model; (b) 1-D electrical conductivity model; (c) Rms; (d) Compared with the average conductivity model of north China.

图 7b给出了300~1200 km深度范围的反演得到的一维导电模型.比较可以看出，校正前后的反演结果在转换带410~670 km范围的电导率差异很大，也反映了海洋效应主要影响地幔转换带及上地幔电性结构的研究.由于原始C-响应虚部相对较大，对应于更好的导电性，即使反演得到如图 7b中红色曲线的低阻结构，其导电性仍然不够高，还是无法拟合观测响应.因此，无论就响应拟合情况还是模型的导电性结构，小周期的原始C-响应是不合理的.校正后响应的反演拟合明显变好，并且看出地幔转换带及以上的电阻率明显增加，平均约为100 Ωm.从图 7c可以看出，由于小周期数据的无法拟合，没有校正的C-响应拟合差明显较大；而海洋效应校正后反演的拟合收敛情况明显得到改善.此外，选择中国平均和兰州电导率模型作为背景进行海洋效应校正得到的反演结果差异不大，结合图 4，这进一步说明海洋效应校正的效果主要由背景模型中上地幔及过渡带的电导率决定.

图 7d为广州地区1-D电导率模型和中国北方大陆平均电导率模型的对比情况，其中北方大陆平均电导率模型采用北京(BMT)、兰州(LZH)、乾县(QIX)和武汉(WHN)4个地磁台站1-D反演得到的电导率模型平均后得到.由图 7d可见，中国北方大陆地幔转换带的电阻率变化范围约为10~70 Ωm，而且随着深度的增加电阻率迅速减小，该反演结果和徐光晶等(2015)研究华北地区地幔转换带的电导率结果是一致的.广州地区地幔转换带的电阻率随深度变化平缓，并集中在约80~100 Ωm范围，这比北方大陆过渡带电阻率高接近一个量级.广州地处华南大陆的南缘.前人的全球地磁测深三维反演结果也表明，广州地区的地幔转换带表现为相对高阻特征(Kelbert et al., 2009；Semenov and Kushinov, 2012)，本次反演得到的模型也验证了该特点.根据Zhao等(2017)的研究，菲律宾板块的俯冲并滞留在华南地区南部的深部地幔转换带中，导致广州地区地幔转换带的为高速异常.因此，推测校正后C-响应反演得到的相对高电阻率，可能是俯冲滞留的菲律宾海板块的电性反映.

本文以广州台站为例，讨论了海洋效应对中国地磁观测测深C-响应的影响.海洋效应三维模型采用“3-D+1-D”的地表电导率分布模型，数值模拟采用球坐标系下的交错网格有限差分技术.经过数值模拟，可以得到如下的结论和认识：

(1) 基于全球海陆电导分布模型，模拟了存在海洋效应时中国地区C-响应的空间分布规律.发现受海洋效应的影响，短周期的C-响应朝海洋方向实部减小，虚部增加，响应等值线和海岸走向基本平行；远离海洋的内陆距离越远，海洋效应影响减小；按照20%的误差水平估计，海洋效应对中国地区地磁观测的影响可达哈尔滨—贵阳一线.

(2) 广州台站的不同1-D模型下的海洋效应校正系数计算表明，海洋效应对C-响应的影响可达20天左右；同时，不同模型的校正系数不同，但采用不受海洋效应影响的中国大陆平均1-D模型或内地的兰州地区1-D模型均可以得到较为可靠的C-响应.针对中国沿海地磁台站的海洋效应校正，利用全球1-D模型进行校正并不合适，建议选择中国平均1-D电导率模型进行校正.

(3) 对广州台站海洋效应校正前后的C-响应进行了1-D反演，反演模型表明周期在20天以内的C-响应对过渡带的电导率敏感度很大.如果需要对地幔转换带电性结构解释和认识，在C-响应反演前必须对沿海地磁台站进行海洋效应校正.

(4) 广州地区上地幔及地幔转换带，电阻率相对中国平均和北方大陆深部的电阻率模型，高大约一个量级，反映了中国华南地区南部下方的地幔，特别是地幔转换带附近，可能处于相对冷的环境，并可能成为菲律宾海板块西向俯冲并滞留到华南板块下方地幔转换带内的电性证据.

值得指出的是，虽然本文只是以典型的广州台站为例，计算了海洋效应对中国地磁观测C-响应的影响.但当针对全国地磁测深数据利用1-D反演技术研究中国地区地幔电性结构时，本文的工作不仅大致给出了海洋效应的影响范围，还提供了海洋效应校正的技术手段.

特别感谢美国地质调查局的Anna Kelbert在三维数值模拟工作中就球坐标网格剖分等细节问题给出的宝贵建议和热情帮助.Chave提供了BIRRP数据处理程序.国家地磁台网中心提供了地磁数据，两位匿名评审专家提出了宝贵意见和建议，在此一并致谢.