2. 中国河北 054000 河北省地震局邢台中心台
2. Xintai Seismic Station, Hebei Earthquake Agency, Heibei Province 054000, China
傅里叶变换是一种重要的信号分析方法,在诸多领域涌现了大量研究成果。1965年,Cooley和Tukey提出快速傅里叶变换(FFT)算法(Cooley et al,1965),使得原有乘法数目由N 2次减少为(N/2)log2N次。FFT算法根据离散傅里叶变换的奇、偶、虚、实等特性进行改进,将原长序列组成若干短序列,利用旋转因子特性减少运算量(程佩青,2010),具有运算量相对减少、处理技术灵活、速度较快等优点。
在电磁学领域,伪谱时域方法(PSTD)(Liu,1997)、共轭梯度快速傅里叶变化算法(CGFFT)(Sarkar et al,1986)、自适应积分方程方法(AIM)(Bleszynski,1996)、预校正快速傅里叶变换方法(PFFT)(Nie et al,2005)、高阶预校正快速傅里叶变换方法(HOPFFT)(Gedney et al,2003)等应用广泛。张建国等(2012, 2013)利用傅里叶频谱分析方法,提取震前电磁异常;何康等(2014)利用此方法在不同台站大地电场数据频谱处理中,得到较显著的研究结果,为判断地电场数据的可靠性提供了依据。文中利用Matlab软件提供的FFT工具箱进行仿真,对中国地磁台网和地电台网部分台站记录的电磁信号做频谱周期分析,以便进一步认识电磁环境变化规律。
1 台网建设及台站选取地电场包括大地电场和自然电场,其中:大地电场是由固体地球外部,特别是电离层中的各种电流体系与地球介质相互作用,在地球内部产生的感应电场,具有区域性变化特征,如日变化、地电暴区域性同步变化等;自然电场是由地球介质局部的物理化学条件形成的局部性电场,例如地形的差异、物质成分的不同以及局部电化学作用等(钱家栋,2010)。一般,由地下场源产生的自然电场,具有相对稳定的变化特性,也会由于局部地下水系、裂隙等作用,在时间和空间上产生剧烈变化(莫承彬等,1995)。由于地电场对地下介质和结构的不均匀性反映敏感,国内外广泛用于地震与火山等地质灾害的监测预警。
地磁场是地球重要的基本物理场之一,由基本磁场、异常场、变化场和感应场组成,包含宇宙空间及地球内核的丰富信息,是地球科学、资源探测、航空航天、航海、通讯、国防、地震预测、空间天气等领域研究的基础。除地面观测,目前还对地磁场开展卫星、航空和海洋测量。但在卫星、飞机及测量船上观测的地磁场变化同时是时间和空间的函数,而地面观测点位固定,观测的地磁场变化仅为时间的函数(李琪,2007)。
1.1 台网建设(1)地电台网建设。1966年河北省邢台地震后,中国开始建设地电观测台站,目前已发展成为国际上独具特色的大规模、规范化的地电观测网,囊括102个地电场固定观测台,集中分布在中国大陆地震断裂带及附近。
(2)地磁台网建设。中国国家地磁台网中心管理了37个在线运行的国家级地磁台,包括97个数字化地磁观测台站及甘肃天祝、四川西昌2个磁通门磁力仪台阵,涉及记录仪器总计约180套,其中秒采样仪器约100套,分钟采样仪器约80套。
1.2 台站选取为全面而有效地开展地磁、地电台网观测数据分析,按纬度分布,由北向南选取所需地电、地磁台站。
地电场台站分布在48.5°—25.4°N范围内,包括德都、三岗、四平、(锦州)义县、温泉、昌黎、安丘、高邮、江油、嵩明10个地电台(表 1),均配备ZD9A-Ⅱ型地电仪进行观测,观测数据连续、真实,具有反映局部环境变化和大空间地球物理场的监测能力。
地磁台站分布在48.5°—25.0° N范围内,包括德都、三岗、乌鲁木齐、大连、昌黎、泰安、马陵山、蒙城、成都、南安10个地磁台(表 2),配备GM4磁通门磁力仪或FHDZ-M15型磁力仪进行观测,观测数据稳定可靠,月连续性均达到99%以上。采用相对观测仪(如FHD系列核旋观测仪、GM系列磁通门磁力仪、FGM系列磁通门磁力仪)与绝对观测仪(如DI无磁经纬仪、G856质子旋进磁力仪、Mag-01地磁经纬仪)结合方式,观测并计算得到分均值、日均值、月均值等产品数据。
地磁场日变化主要有静日和扰日2种,静日变化时不仅地磁观测分量形态清晰,地电场观测日变形态也规律明显,部分台站固体潮清晰可见。磁暴是磁扰中具有代表性的空间电离层活动,在地磁和地电场监测中具有同步性。基于地电场监测受局部环境变化灵敏的特点,各台站扰动程度不尽相同,且日常监测中多有复杂的干扰因素掩盖其中,对各台频谱特征的计算结果,可能不一致。
国际上,地磁活动情况指数(简称磁情指数)种类较多,中国国家地磁台网中心采用Kp指数分析磁扰活动,每月在官方网站公布磁情活动结果(http://10.2.210.50/)。Kp指数是用分级方法描述地磁活动性的一种数字指标,从0—9分为10级,数值越大扰动越剧烈,Kp≤2时为磁静日。磁暴分为3类:中常磁暴m(Kp = 5—6)、中烈磁暴ms(Kp = 7—8)、强烈磁暴s。
地电学科监测记录中的地电暴,与磁暴同步出现,二者具有共同外源场,受全球性空间电磁场扰动影响。目前,在地电学科观测日志中有关于地电暴事件的统一描述,但对地电暴的分类、扰动程度的量化计算等定义尚不明确,故文中采用地磁学科产出的磁情报告选取地电暴日。
选取磁静日记录的地磁、地电观测数据和地电暴日(磁暴时段)沿海及内陆地区部分台站记录的地电数据,进行傅里叶频谱分析。所选数据应具有长期稳定性、可靠性,无明显环境干扰和异常变化。
2.1 磁静日地电、地磁数据傅里叶频谱特征 2.1.1 地电场数据多数台站地电场观测包含NS、EW测向,且均采用长、短极距并行观测,若各观测电极附近无场地环境干扰,则同一测向长、短极距数据变化形态基本同步,故文中只采用2个测向的长极距数据进行分析。
在其他相关文献中,数据样本时间为单天或1个月,文中取10天为一个研究时段,以丰富相关研究内容。选取10个地电台站2018年6月3日—12日记录的分钟值数据,分析地电数据频谱特征,结果见图 1、表 3。
由图 1和表 3可知:①12 h周期成分显著,其次为8 h、24 h,部分台站存在周期成分为6 h、16 h、18 h、30 h、40 h、48 h等弱频信号;②对比每个台站NS、EW测向最高频幅值,其中7个台站EW向幅值大于NS向,个别高纬度台站(约39° N以上)存在相反现象;③温泉台NS测向、昌黎台EW和NS测向、安丘台EW和NS测向、义县台EW测向,最大频幅值相对偏高。
综合分析发现:①地电场观测中12 h周期成分发育显著,有8个地电台站的2个测向均出现12 h周期成分,24 h、8 h周期成分则相对显著;②EW测向最高谱值整体大于NS测向。
2.1.2 地磁场数据地磁场包含北向分量X、东向分量Y、垂向分量Z、水平强度H、磁偏角D、磁倾角I、总场强度F。根据麦克斯韦方程,地磁场各分量的一阶差分变化与地电场同向分量变化形态一致。为与地电场记录的NS、EW分量做同向分析,地磁数据采用X、Y分量做频谱分析。以10个地磁台站2018年5月18日—28日记录的分钟值数据为例,选取其中10天数据分析地磁频谱特征。因篇幅所限,只绘制主要周期成分百分比图件,并按照地磁台站由北向南的次序(纬度递减顺序),绘制X、Y分量最大谱值曲线,见图 2。
由图 2可见:①X、Y分量显著频谱周期成分为24 h、12 h、8 h(按谱值大小暂称为最大谱、中谱、最小谱),还包含20 h、36 h、54 h、72 h等非显著周期成分,Y分量未显示其他周期成分;②X分量最大谱值显示7个台站具有12 h周期成分,最小谱值显示6个台站具有8 h周期成分。Y分量最大谱值显示8个台站具有24 h周期成分,最小谱值显示10个台站均具有8 h周期成分。X、Y分量的中谱显示8个台站具有12 h较显著周期成分;③Y分量最大谱值为7.92×104 — 1.08×105 mV/km,X分量最大谱值为3.23×104 — 5.64×104 mV/km,且各地磁台Y分量均大于X分量,二者差值约3.12×104 — 6.33×104 mV/km;④地磁台站分布由北及南,Y分量最大谱值逐渐减小,X分量则两头略高,在40.0° — 35.0° N范围内变化平稳。
综合分析发现:①地磁场观测中X分量的12 h周期成分显著发育,Y分量24 h周期成分发育明显;②8 h周期成分在2个分量中均为最小谱值(相对来讲);③各地磁台Y分量最大谱值均大于X分量,与地电场计算结果基本一致。
2.2 地电暴日沿海及内陆地电场傅里叶频谱特征地磁台网观测不易受地表以下环境干扰,各台数据变化具有同步性、广域性,且对磁暴日地磁数据傅里叶频谱特征分析,徐文耀等(徐文耀,1987;徐文耀等, 1994, 2009)学者早有研究,本文不再赘述。地电台网观测受台址条件和测区环境变化影响较大,具有局域性和地域性特征,故只分析地电观测数据。
2017年9月26日23时47分至29日20时出现SC型磁暴,活动程度m,抽取该时段地电台网中沿海地区的大连、昌黎、嵩明地电台及内陆地区的红山、天水、乾陵台分钟值记录,分析地电暴日地电场NS、EW向傅里叶频谱变化特征,结果见图 3。
由图 3可见:①沿海地区地电台站:大连台NS测向主要周期成分依次为20 h、24 h、12 h、8 h,最大谱值3.372×104 mV/ km;EW测向周期成分周期依次为12 h、8 h,最大谱值2.552×104 mV/ km。昌黎台NS测向主要周期成分依次为36 h、24 h、12 h、8 h,最大谱值5.09×104 mV/ km;Y分量主要周期成分依次为12 h、8 h,最大谱值3.049×104 mV/ km。嵩明台NS测向主要周期成分依次为24 h、12 h、18 h、8 h、36 h,最大谱值3.245×104 mV/ km;EW测向主要周期成分依次为12 h、8 h、18 h、6 h,最大谱值2.344×104 mV/ km;②内陆地区地电台站:红山台NS测向主要周期成分依次为18 h、24 h、6 h、12 h,最大谱值4.253×104 mV/ km;EW测向主要成分周期依次为12 h、8 h、18 h、14.4 h,最大谱值2.597×104 mV/ km。天水台NS测向主要周期成分依次为18 h、12 h、6 h,最大谱值3.788×104 mV/ km;EW测向主要周期成分依次为12 h、8 h、18 h,最大谱值2.474×104 mV/ km。乾陵台NS测向主要周期成分依次为18 h、12 h,最大谱值3.501×104 mV/ km,EW测向主要周期成分依次为12 h、8 h、18 h,最大谱值2.411×104 mV/ km。
综合分析发现:①沿海地区NS测向FFT频谱变化主要周期成分比内陆地区复杂且偏大,大于18 h的周期成分相对活跃,最大振幅谱对应周期成分为20 h、24 h、36 h,而内陆地区最大周期成分为18 h,其他周期成分相对稳定;②各地电台EW测向谱值变化较一致,均体现了以12 h、8 h周期成分为主的变化特征,其他周期成分相对同步且稳定;③各地电台NS测向最大谱值均大于EW测向。
3 讨论电场和磁场是相互依存的关系,在连续地震监测记录中,二者均受地球内、外源场的影响,同时又互相影响。地电场与地磁场在磁静日期间具有12 h、24 h、8 h显著周期成分,且地电场EW向(或地磁场Y分量)频谱幅度大于NS向(或地磁场X分量),说明二者具有同源性。因地电场观测受局域性影响明显,各台捕捉的地球物理信号“纯净度”参差不齐,周期成分较复杂,而地磁场频谱成分则相对“纯净”,且Y分量比X分量更“纯净”,显著周期成分只有12 h、24 h、8 h。
3.1 地电场变化地电场变化与地下介质有一定关系,谱值大小与观测台址的浅层电阻率有关,电阻率越高,大地电场各种周期成分谱值越大(杜学彬等,2007;叶青等,2007)。地电场观测信号的复杂性与介质的非均匀性关系甚密(杜学彬等,2007;叶青等,2007;马钦忠等,2014),存在同台不同向信号幅度不一致现象,1—2天短时域FFT频谱各台也不甚一致,计算得到的长周期频谱较差(何康等,2014)。据实验研究,在对青岛换流站接地极向地下注入2 100—3 004 A大电流时,华北东部冀鲁豫苏皖地区21个地电场台站中,沿郯庐断裂带、跨郯庐断裂带西南、跨郯庐断裂带北西向部分台站(包括昌黎、安丘台),记录到的地电流信号特征变化非常大,明显优于其他台站,存在明显方向性特征和“穴位”(马钦忠等,2014)。由表 3可知,昌黎、安丘台的NS和EW测向、温泉台NS测向和义县台EW测向最大频谱值均偏高。根据表 1—表 2中各台基岩资料,昌黎和义县台近邻郯庐断裂带西北测,基岩以花岗岩为主,属高阻台站;安丘台位于郯庐断裂带,基岩以砂砾岩为主,电极埋设较深(20 m,表 1),位于基岩,介质各向异性度、电阻率均比黄土层高,造成傅里叶频谱值偏大;同理,基岩地处高阻层的温泉台傅里叶频谱值偏大。
3.2 地磁场变化据赵旭东等(2008)的分析,地磁场在磁静日的变化主要由1—5次谐波组成,周期分别为24 h、12 h、8 h、6 h、4.8 h,振幅依次减小。此外,存在27日谱线和其他次要谱线。全球Sq场日变化基本取决于纬度和地方时,纬度在± 60°内,X分量在周期成分为12 h时出现3个极区,分别位于赤道、60° N、60° S附近,且在纬度± 30°处,X值发生反向,说明X分量(或地电场EW分量)受地方时和纬度影响明显。在纬度± 60°内,Y分量等值线存在4个极区,以赤道为对称轴,南北半球各2个,其中北半球2个极区的极值出现时间分别约8 h、14 h。依据劳埃德季节划分,6月为夏至点月份(J月),此时北半球中纬度Sq等效电流体系的白天电流涡中心时间为11:30(徐文耀等,2009)。实际上Sq变化包含纬度、地方时、经度3部分,表达式可写为
$ S_{\mathrm{q}}(\theta, \lambda, t)=S_{t}(\theta, t)+S_{\lambda}(\theta, \lambda)+S_{t \lambda}(\theta, \lambda, t) $ |
依据全球台网资料,可进行傅里叶谐波分析,提取St (θ,t)、Sλ(θ,λ)、Stλ(θ,λ,t)(徐文耀等,2009)。地方时St (θ,t)成分比Sλ(θ,λ)、Stλ(θ,λ,t)大得多,其傅里叶变换谐波振幅中,中低纬度地区Y分量明显大于X分量(祁贵仲,1975;徐文耀,1987;徐文耀等,1994),由此验证了地电场和地磁场EW(或Y)分量频谱值大于NS(或X)分量。
本文只绘制傅里叶谐波分析中地磁X、Y分量地方时St (θ,t)成分的一、二次谐波,谐波振幅曲线见图 4,图中方框区域为文中涉及台站的纬度范围。由图 4可见:在研究区域内,X分量地磁谐波在30° N附近出现低值,向高低纬度增大;Y分量在44° N附近出现高值,向高低纬度减小。由此可知,在研究区内,地磁Y分量频谱值随纬度变化,由北向南依次减弱;纬度较高的德都和较低的南安地磁台X分量频谱值略高。
选取中国大陆25.0°—48.5°N范围内10个地电场台站和10个地磁台站记录,分别对地电场NS、EW分量和地磁场X、Y分量进行FFT频谱分析,可以得到以下结论。
(1)磁静日:①地电场和地磁场均以12 h周期成分为主,24 h、8 h周期成分为次,其他非显著周期成分在地电场观测EW向和地磁X分量中所占频次较大,地电场NS向和地磁Y分量较“纯净”,只存在12 h、24 h、8 h周期成分;②地电场台站EW向最大频谱值均大于NS向,地磁台站Y分量最大频谱值均大于X分量;③各地磁台X分量最大频谱值多以12 h周期成分表现显著,Y分量则24 h周期成分显著;④台站由北向南,Y分量最大谱值逐渐减小,X分量则两头略高,中间35.0°—40.0° N位置变化平稳。
(2)地电暴日:①沿海地区地电台站NS测向主要周期成分复杂,且频谱周期偏长,大于18 h的周期成分相对活跃,如20 h、24 h、36 h,而内陆地区最大周期成分为18 h,其次为12 h。EW测向谱值在各台变化较一致,均体现了以12 h、8 h周期成分为主的变化特征;②各台NS测向最大谱值均大于Y分量。
(3)在地电场观测中,由于地下介质不同或所处断裂带位置不同,傅里叶频谱变化具有局域特征,如昌黎和安丘台的NS和EW测向、义县台EW向、温泉台NS向最大频谱值偏高,可能与台站所处位置有关。
程佩青. 数字信号处理教程[M]. 北京: 清华大学出版社, 2010: 45-90. | |
杜学彬, 叶青, 赵杰, 等. 地电场日变化研究[J]. 地震, 2007, 27(Z1): 121-130. | |
何康, 洪德全, 李军辉, 等. 安徽省数字化地电场FFT频谱分析研究[J]. 防灾科技学院学报, 2014, 16(2): 63-67. DOI:10.3969/j.issn.1673-8047.2014.02.011 | |
李琪. 国内外地磁台网观测能力评估[J]. 国际地震动态, 2007(9): 20-28. DOI:10.3969/j.issn.0253-4975.2007.09.004 | |
马钦忠, 李伟, 张继红, 等. 与大电流信号有关的华北东部地区地电场空间变化特征的研究[J]. 地球物理学报, 2014, 57(2): 518-530. | |
莫承彬, 陈忠献, 陆怀成. 自然电场法剧变场的起因初探及其应用[J]. 物探与化探, 1995, 19(4): 315-318. | |
祁贵仲. 局部地区地磁日变分析方法及中国地区Sq场的经度效应[J]. 地球物理学报, 1975, 18(2): 104-117. | |
钱家栋. 地震电磁学理论基础与观测技术[M]. 北京: 地震出版社, 2010: 25-95. | |
徐文耀. 用中低纬度地磁资料确定Seq电流体系的焦点[J]. 空间科学学报, 1987, 7(2): 117-127. | |
徐文耀, 李卫东. Sq外源和内源电流体系的经度效应和UT变化[J]. 地球物理学报, 1994, 37(4): 440-447. DOI:10.3321/j.issn:0001-5733.1994.04.004 | |
徐文耀. 地球电磁现象物理学[M]. 合肥: 中国科学技术大学出版社, 2009: 260-310. | |
叶青, 杜学彬, 周克昌, 等. 大地电场变化的频谱特征[J]. 地震学报, 2007, 29(4): 382-390. DOI:10.3321/j.issn:0253-3782.2007.04.005 | |
张建国, 姚丽, 刘晓灿, 等. 地震电离层VLF电磁场频谱特征研究[J]. 大地测量与地球动力学, 2012, 32(3): 110-115. | |
张建国, 焦立果, 刘晓灿, 等. 汶川MS 8.0级地震前后ULF电磁辐射频谱特征研究[J]. 地球物理学报, 2013, 56(4): 1253-1261. | |
赵旭东, 杜爱民, 徐文耀, 等. Sq电流系午前午后不对称性现象的来源[J]. 地球物理学报, 2008, 51(3): 643-649. DOI:10.3321/j.issn:0001-5733.2008.03.005 | |
Cooley J W, Tukey J W. An algorithm for the machine calculation of complex Fourier series[J]. Math Comput, 1965, 19(90): 297-301. DOI:10.1090/S0025-5718-1965-0178586-1 | |
Bleszynski E, Bleszynski M, Jaroszewicz T. AIM:Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems[J]. Radio Science, 1996, 31(5): 1225-1251. DOI:10.1029/96RS02504 | |
Gedney S, Zhu A M, Tang W H, et al. A fast, high-order quadrature sampled pre-corrected fast-Fourier transform for electromagnetic scattering[J]. Microw Opt Technol Lett, 2003, 36(5): 343-349. DOI:10.1002/(ISSN)1098-2760 | |
Liu Q H. The PSTD algorithm:A time-domain method requiring only two cells per wavelength[J]. Microw Opt Technol Lett, 1997, 15(3): 158-165. DOI:10.1002/(ISSN)1098-2760 | |
Nie X C, Li L W, Yuan N, et al. Precorrected-FFT solution of the volume integral equation for 3-D inhomogeneous dielectric objects[J]. IEEE Trans Antennas Propag, 2005, 53(1): 313-320. DOI:10.1109/TAP.2004.838803 | |
Sarkar T, Arvas E, Rao S. Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies[J]. IEEE Trans Antennas Propag, 1986, 34(5): 635-640. DOI:10.1109/TAP.1986.1143871 |