CN110068857A - 基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 - Google Patents
基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 Download PDFInfo
- Publication number
- CN110068857A CN110068857A CN201910261107.4A CN201910261107A CN110068857A CN 110068857 A CN110068857 A CN 110068857A CN 201910261107 A CN201910261107 A CN 201910261107A CN 110068857 A CN110068857 A CN 110068857A
- Authority
- CN
- China
- Prior art keywords
- star
- principal component
- magnetic field
- data
- swarm
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000513 principal component analysis Methods 0.000 title claims abstract description 27
- 239000002243 precursor Substances 0.000 title claims abstract description 26
- 238000000034 method Methods 0.000 title claims abstract description 23
- 230000000694 effects Effects 0.000 claims abstract description 47
- 238000011160 research Methods 0.000 claims abstract description 30
- 239000004615 ingredient Substances 0.000 claims abstract description 17
- 239000000284 extract Substances 0.000 claims abstract description 11
- 239000011159 matrix material Substances 0.000 claims description 40
- 238000004458 analytical method Methods 0.000 claims description 18
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000002159 abnormal effect Effects 0.000 claims description 10
- 238000005259 measurement Methods 0.000 claims description 7
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 230000005856 abnormality Effects 0.000 claims description 4
- GNFTZDOKVXKIBK-UHFFFAOYSA-N 3-(2-methoxyethoxy)benzohydrazide Chemical compound COCCOC1=CC=CC(C(=O)NN)=C1 GNFTZDOKVXKIBK-UHFFFAOYSA-N 0.000 claims description 3
- FGUUSXIOTUKUDN-IBGZPJMESA-N C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 Chemical compound C1(=CC=CC=C1)N1C2=C(NC([C@H](C1)NC=1OC(=NN=1)C1=CC=CC=C1)=O)C=CC=C2 FGUUSXIOTUKUDN-IBGZPJMESA-N 0.000 claims description 3
- 230000005389 magnetism Effects 0.000 abstract description 5
- 239000005433 ionosphere Substances 0.000 description 8
- 238000000605 extraction Methods 0.000 description 6
- 230000008901 benefit Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 4
- 238000005516 engineering process Methods 0.000 description 4
- 238000000926 separation method Methods 0.000 description 4
- 241001493796 Paraboea crassifolia Species 0.000 description 3
- 238000012545 processing Methods 0.000 description 3
- 230000002547 anomalous effect Effects 0.000 description 2
- 230000007812 deficiency Effects 0.000 description 2
- 239000004744 fabric Substances 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 1
- 241000406668 Loxodonta cyclotis Species 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000005684 electric field Effects 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 239000002075 main ingredient Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/01—Measuring or predicting earthquakes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
- G01V3/08—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
- G01V3/087—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices the earth magnetic field being modified by the objects or geological structures
Landscapes
- Life Sciences & Earth Sciences (AREA)
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Physics & Mathematics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Electromagnetism (AREA)
- Geochemistry & Mineralogy (AREA)
- Acoustics & Sound (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及一种基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法，读取Swarm卫星A星和C星的磁场Y分量数据，并根据数据的标志位Flags_B剔除有误数据；根据地震选择研究的时间范围、研究区域；选取Swarm A星和C星的地方时为夜晚的轨道数据；利用IGRF模型去除地磁主磁场的影响；通过对Swarm A星和C星的磁场Y分量数据进行主成分分析，将原始的数据投影到一组新的空间正交基上，得到方差从大到小排列的各个主成分；通过各个主成分与地磁指数的对比，找到与地磁活动相关性高的成分并将其去除，只对剩余的主成分进行分析；通过偏度和峰度定义的偏峰度系数对剩余的主成分进行异常提取，提取地震前兆异常。本发明去除地磁活动干扰，精确提取地震前兆异常。
Description
技术领域
本发明涉及电离层磁场地震前兆异常提取领域，具体涉及一种基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法。
背景技术
越来越多的证据表明，电磁现象可用于对大地震进行短期预测。除了地面电磁现象的研究，基于LAIC(Lithospheric-Atmospheric-Ionospheric Coupling岩石圈-大气层-电离层)机制，地震电离层异常的研究依赖于卫星的发展也成为地震前兆异常的一个方向。DEMETER卫星是第一颗专门用于监测地震相关信号的卫星。大量学者通过震例分析，统计分析和多手段综合等分析方法对DEMETER卫星载荷测量的各种参数，如电子密度，离子密度/温度，电场，磁场，高能粒子和电离层总电子含量进行分析，结果表明，地震前后均存在异常现象。地震引起的电磁波频率范围相当宽(DC-VHF频段)，而ULF频段电磁干扰被认为是最有希望的地震前兆之一。Swarm是欧洲航天局的任务，旨在精确测量地球的磁场，测量采样率低(磁场1Hz/50Hz，电子密度2Hz)，但精度高，近年来，也有一些学者利用Swarm卫星测量的电离层参数进行震前异常提取。因此，基于Swarm卫星磁场数据的地震前兆异常提取是可行的。
主成分分析法能够将原始时间序列由强到弱投影到一组新的空间正交基上，达到降维和主特征分离目的，在本发明中，通过对Swarm卫星A星和C星磁场进行主成分分析，能够将地磁活动与其他影响因素分离到不同的成分中，通过去除与地磁指数相关性高的成分以去除地磁活动干扰，避免后续提取的异常由地磁活动引起；偏度和峰度是一种统计方法，分别表示数据分布的对称性和厚尾性，当数据中出现极端值时，偏度值和峰度值会相应的增大，因此利用偏度和峰度定义的偏峰度系数对去除地磁影响后的主成分进行分析，能够有效提取地震前兆异常。
CN101793972A公开了一种强地震短临预测卫星热红外和短波亮温异常技术，主要利用极轨卫星、静止卫星、气象卫星和小卫星座进行观测，并结合其他载有红外波段扫描仪的卫星，还包括卫星接收设备和图象处理设备，此时，卫星热红外云图的彩色密度分割档次，可对不同季节、不同纬度地区都适用；扫描仪所得灰度值可经大气模型校正获得实际温度值；可能在孕震区云层上方出现降温环现象；同时可能出现怪状云；发生在内陆高原区域的地震其应力热线交汇部位为未来震中，从而使其对地震短临预测的成功率提高到50％以上。
CN1047924公开了一种利用卫星热红外异常做中强以上地震三要素临震预报的方法。利用卫星、卫星接收的处理设备，利用卫星热红外异常和有关气象震兆作临震三要素预报。本发明资料准确可靠，覆盖面广，信息量大，速度快，能抓住地震前兆，并能掌握时空动态变化，因而能用于地震三要素的临震预报。
Akhoondzadeh等利用中值法分析了2016年4月16日的厄瓜多尔地震，对2015年11月1日至2016年4月30日的Swarm卫星群A星、B星和C星的多个载荷的数据进行异常提取；Santis等分析了2015年4月25日的尼泊尔地震，对地震前后一个月内的Swarm卫星A星测量的东向磁场分量提取到的异常按天进行累计，发现在震前和震后异常的增长呈现不同的速度，在异常增长速度最快的时候可能为地震异常。但到目前为止，尚未有基于Swarm双星磁场数据联合分析的地震前兆异常提取方法，即使是多参数多颗卫星的分析，也只是单独分析，对提取的异常结果进行对比，且只对静磁时间段的数据进行分析，没有实质上的双星联合进行异常提取的研究；也缺少将地磁活动的影响去除后进行异常提取的研究。
发明内容
本发明所要解决的技术问题在于提供一种基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法,充分利用Swarm A星和C星并列飞行且轨道高度几乎相同的特性，将主成分分析应用于两颗星Y分量的磁场数据进行特征分离，去除地磁活动的干扰后，利用偏度和峰度进行异常提取，弥补了上述技术中针对单颗卫星分析，或只是对多颗卫星数据的异常结果结合分析，且只对静磁时间段的数据进行分析，不能充分利用Swarm卫星的A、C两颗星几乎并列飞行、同时测量磁场的优势，去除地磁活动干扰后，提取地震前兆异常的不足。
本发明是这样实现的，
一种基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法,其特在于，该方法包括：
a读取Swarm卫星A星和C星的磁场Y分量数据，根据数据的标志位Flags_B剔除有误数据；
b根据地震选择研究的时间范围、研究区域；
c选取Swarm A星和C星的地方时为夜晚的轨道数据；
d将磁场Y分量数据减去IGRF(International Geomagnetic Reference Field)模型去除主磁场的影响；
e、对减去模型后的磁场Y分量残余值进行主成分分析，将不同的信号分离到不同的成分中；
f、计算各个主成分与地磁指数的相关系数，找到与地磁活动相关的主成分并去除；
g、对剩余的主成分求差分，求偏度和峰度，并定义偏峰度系数，利用偏峰度系数进行异常提取，输出异常结果。
进一步地，步骤b所述的根据地震选择研究的时间范围、研究区域，是选取地震前90天至地震后30天为研究时间范围，根据Dobrovolsky's公式R＝100.43M，其中M为地震的震级，计算地震的影响区域，并选择以地震为圆心，以R为半径的圆形区域的外接正方形区域为研究区域。
进一步地，步骤c所述的选取Swarm卫星A星和C星的地方时夜晚数据，是计算SwarmA星和C星的地方时，LT(地方时)＝UTC(世界时)+α/15(α为地理经度)，LT指的地方时，UTC指的是世界时，α为地理经度，夜晚指的是夜晚地方时为18：00-06：00。
进一步地，步骤d所述的将磁场Y分量数据减去IGRF模型去除主磁场的影响包括：通过国际给出的IGRF模型，分别计算Swarm A星和C星磁场Y分量对应的IGRF模型值，并分别从测量的磁场Y分量中减去模型值以去除主磁场的影响，得到Swarm A星和C星的Y分量磁场残余值。
进一步地，步骤e所述的对减去模型后的Swarm A星和C星的Y分量磁场残余值进行主成分分析，将不同的信号分离到不同的成分中，具体包括：
将减模型后的Swarm A星和C星的磁场残余值按照时间序列表示为：
Bi＝[bi(1),bi(2)...,bi(m)]，i＝1,2,...,n
其中m为轨道长度，n为磁场数据维度，得到矩阵Y的表达式为：
计算矩阵Y的协方差矩阵CY(n×n)的元素γuv，计算公式为：
其中，biu和biv分别为矩阵Y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值；
计算协方差矩阵的特征值和特征向量：
CY＝RΛRT
其中，Λ(n×n)为从大到小排列的特征值对角矩阵，R(1×n)为对应特征值的特征向量，将特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)；
利用矩阵R将矩阵Y的原始数据进行线性投影得到主成分Φ，
Φ＝R·Y＝[Φ1,Φ2,...,Φn]T
其中，Φ1,Φ2,...,Φn为第1至第n个主成分(1×m)，主成分是原始数据投影到新的空间的数据，按照方差从大到小的顺序排列，并且各个主成分互不相关。
进一步地，步骤f所述的计算各个主成分与地磁指数的相关系数，找到地磁活动相关的主成分并去除，包括将求得的所有轨道的各个主成分的特征值表示为：
Ei＝[λi1,λi2,...,λip]，i＝1,2,...,n
其中p为研究的轨道个数，n为数据维度，得到矩阵Γ＝[E1 T,E2 T,...,En T]T，Γ的表达式为：
矩阵Γ：
计算时间序列x(t)和y(t)相关系数的公式为：
其中σx和σy为x(t)和y(t)的方差，和分别是x(t)和y(t)的均值,Cxy(k)是时滞为k时的两时间序列的协方差；rxy(k)为两个时间序列在时滞为k时的相关系数；
根据公式计算所有轨道各个主成分的特征值Ei与对应的地磁指数ap在时滞为0时的相关系数：
相关系数高的说明相应的主成分主要反应地磁活动情况，相关系数低的成分表明几乎与地磁活动无关，通过去除与地磁指数相关性高的主成分，以去除地磁活动的影响。
进一步地，步骤g对剩余的主成分求差分，求偏度和峰度，包括：设上一步去除地磁活动影响后，与地磁指数相关性低的主成分为第q个主成分，被认为几乎不受地磁活动的干扰，则对第q个主成分Φq沿轨道进行差分得到：
dΦq＝[dΦq1,dΦq2,...,dΦqp]
其中p为研究的轨道个数。
进一步地，求差分结果的偏度和峰度，利用偏度和峰度定义偏峰度系数包括：
偏度的计算公式为：
峰度的计算公式为：
其中x为分析的数据，X为x的均值；
计算差分后的第q个主成分的偏度和峰度：
进一步地，结合偏度和峰度定义偏峰度系数γ：
其中p为研究的轨道的个数。
进一步地，求系数γ的均值μγ和标准差σγ，并设阈值为Thγ＝μγ+2×σγ，当轨道偏峰度系数大于设定阈值时，该轨道被认为是异常轨道。
本发明与现有技术相比，有益效果在于：
本发明充分利用Swarm A星和C星并列飞行且轨道高度几乎相同的特性，将主成分分析应用于两颗星Y分量的磁场数据进行特征分离，去除地磁活动的干扰后，利用偏度和峰度进行异常提取，弥补了上述技术中针对单颗卫星分析，或只是对多颗卫星数据的异常结果结合分析，且只对静磁时间段的数据进行分析，不能充分利用Swarm卫星的A、C两颗星几乎并列飞行、同时测量磁场的优势，去除地磁活动干扰后，提取地震前兆异常的不足。
附图说明
图1为基于主成分分析的Swarm双星磁场数据异常提取方法流程图；
图2为厄瓜多尔地震震中、影响区域和研究区域示意图；
图3为Swarm A星和C星磁场Y分量原始曲线；
图4为减模型后Swarm A星和C星磁场Y分量残余值曲线；
图5为Swarm A星和C星磁场Y分量残余值主成分分析结果曲线，图5(a)和5(b)分别为第一主成分曲线和第二主成分曲线；
图6为研究轨道的特征值与对应的地磁指数曲线，其中(a)第一特征值、(b)第二特征值与(c)地磁指数ap；
图7为研究轨道第二主成分差分后的偏度(a)和峰度曲线(b)；
图8为研究轨道偏峰度系数曲线。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本发明基于主成分分析的Swarm双星磁场数据地震前兆异常提取，首先读取Swarm卫星A星和C星的磁场Y分量数据，并根据数据的标志位Flags_B剔除有误数据；根据地震选择研究的时间范围、研究区域；选取Swarm A星和C星的地方时为夜晚的轨道数据；利用IGRF(International Geomagnetic Reference Field)模型去除地磁主磁场的影响；通过对Swarm A星和C星的磁场Y分量数据进行主成分分析，将原始的数据投影到一组新的空间正交基上，得到方差从大到小排列的各个主成分，以达到特征分离的目的；通过各个主成分与地磁指数的对比，找到与地磁活动相关性高的成分并将其去除，只对剩余的主成分进行分析；通过偏度和峰度定义的偏峰度系数对剩余的主成分进行异常提取，提取地震前兆异常。
具体包括：
a、读取Swarm A星和C星的磁场Y分量数据，并根据标志位Flags_B剔除有误数据选择有效数据；
b、根据地震选择研究的时间范围、研究区域；
c、选取Swarm卫星A星和C星的地方时夜晚数据；
d、将磁场Y分量数据减去IGRF模型去除主磁场的影响；
e、对减去模型后的磁场Y分量残余值进行主成分分析，将不同的信号分离到不同的成分中；
f、计算各个主成分与地磁指数的相关系数，找到与地磁活动相关的主成分并去除；
g、对剩余的主成分求差分，求偏度和峰度，并定义偏峰度系数，利用偏峰度系数进行异常提取，输出异常结果；
步骤a所述的读取Swarm卫星A星和C星的磁场Y分量数据，并根据标志位Flags_B剔除有误数据选择有效数据，是读取Swarm A星和C星的矢量磁力仪测得的磁场的Y分量，即NEC(North-East-Center)中指向东向的分量磁场数据，已经有相关文献表述，Y分量磁场能够更清晰的记录到电离层异常变化；卫星矢量磁场数据的质量主要通过标志位Flags_B进行筛选，剔除无效数据后，仅对有效数据进行分析；
步骤b所述的根据地震选择研究的时间范围、研究区域，是选取地震前90天至地震后30天为研究时间范围，根据Dobrovolsky's公式R＝100.43M(M为地震的震级)计算地震的影响区域，并选择以地震为圆心，以R为半径的圆形区域的外接正方形区域为研究区域；
步骤c所述的选取Swarm卫星A星和C星的地方时夜晚数据，是计算Swarm A星和C星的地方时，LT(地方时)＝UTC(世界时)+α/15(α为地理经度)，由于电离层在白天受到干扰会产生异常，为了避免受到白天电离层活动的干扰难以解释，因此选取夜晚地方时为18：00-06：00的轨道，由于卫星沿轨飞行，测得的数据是随时间和空间变化的量，为了后续便于分析，对数据按轨道进行分析和异常提取；
步骤d所述的将磁场Y分量数据减去IGRF模型去除主磁场的影响，是由于卫星测量磁场为地磁场，主要包含主磁场、岩石圈磁场、变化磁场和感应磁场四部分。地磁场幅值为104nT量级，而地震前兆引起的变化主要表现在岩石圈磁场中，只有几nT甚至更小，相对于总体幅值来说非常微弱，直接对原始测量磁场进行操作则无法提取异常，因此为了后续的异常提取，通过国际给出的IGRF模型，分别计算Swarm A星和C星磁场Y分量对应的IGRF模型值，并分别从测量的磁场Y分量中减去模型值以去除主磁场的影响，得到Swarm A星和C星的Y分量磁场残余值；
步骤e所述的对减去模型后的Swarm A星和C星的Y分量磁场残余值进行主成分分析，将不同的信号分离到不同的成分中，主要包括：
将减模型后的Swarm A星和C星的磁场残余值按照时间序列表示为：
Bi＝[bi(1),bi(2)...,bi(m)]，i＝1,2,...,n
其中m为轨道长度，n为磁场数据维度，得到矩阵Y的表达式为：
计算矩阵Y的协方差矩阵CY(n×n)的元素γuv，计算公式为：
其中，biu和biv分别为矩阵Y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值。
计算协方差矩阵的特征值和特征向量：
CY＝RΛRT
其中，Λ(n×n)为从大到小排列的特征值对角矩阵，R(1×n)为对应特征值的特征向量，将特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)；
利用矩阵R将矩阵Y的原始数据进行线性投影得到主成分Φ，
Φ＝R·Y＝[Φ1,Φ2,...,Φn]T
其中，Φ1,Φ2,...,Φn为第1至第n个主成分(1×m)，主成分是原始数据投影到新的空间的数据，按照方差从大到小的顺序排列，并且各个主成分互不相关。
步骤f所述的计算各个主成分与地磁指数的相关系数，找到地磁活动相关的主成分并去除，包括将求得的所有轨道的各个主成分的特征值表示为：
Ei＝[λi1,λi2,...,λip]，i＝1,2,...,n
其中p为研究的轨道个数，n为数据维度，得到矩阵Γ＝[E1 T,E2 T,...,En T]T，Γ的表达式为：
矩阵Γ：
计算时间序列x(t)和y(t)相关系数的公式为：
其中σx和σy为x(t)和y(t)的方差，和分别是x(t)和y(t)的均值,Cxy(k)是时滞为k时的两时间序列的协方差；rxy(k)为两个时间序列在时滞为k时的相关系数。
根据公式计算所有轨道各个主成分的特征值Ei与对应的地磁指数ap在时滞为0时的相关系数：
相关系数高的说明相应的主成分主要反应地磁活动情况，相关系数低的成分表明几乎与地磁活动无关，通过去除与地磁指数相关性高的主成分，以去除地磁活动的影响。
步骤g对剩余的主成分求差分，求偏度和峰度，并定义偏峰度系数，利用偏峰度系数进行异常提取，包括：设上一步去除地磁活动影响后，与地磁指数相关性低的主成分为第q个主成分，被认为几乎不受地磁活动的干扰，则对第q个主成分Φq沿轨道进行差分得到：
dΦq＝[dΦq1,dΦq2,...,dΦqp]
其中p为研究的轨道个数。
求差分结果的偏度和峰度，利用偏度和峰度定义偏峰度系数，利用偏峰度系数进行异常提取。
偏度和峰度是一种统计手段，偏度反应了数据分布的对称性，峰度反映了数据分布的厚尾性，当数据中出现极端值时，偏度和峰度会相应的增大。偏度的计算公式为：
峰度的计算公式为：
其中x为分析的数据，X为x的均值；
计算差分后的第q个主成分的偏度和峰度：
结合偏度和峰度定义偏峰度系数γ：
其中p为研究的轨道的个数。
求系数γ的均值μγ和标准差σγ，并设阈值为Thγ＝μγ+2×σγ，当轨道偏峰度系数大于设定阈值时，该轨道被认为是异常轨道。
以下结合实例对本发明进行进一步详细说明。
针对2016年4月16日23:58:36(世界时)发生的震级为7.8级的厄瓜多尔地震，以Swarm卫星群中的A星和C星的矢量磁场Y分量(1Hz)为例。
a、读取Swarm A星和C星的磁场Y分量数据，并根据标志位Flags_B剔除有误数据选择有效数据，是读取Swarm A星和C星的矢量磁力仪测得的磁场的Y分量，即NEC(North-East-Center)中指向东向分量磁场数据，已经有相关文献表述，Y分量磁场能够更清晰的记录到电离层异常变化；卫星矢量磁场数据的质量主要通过标志位Flags_B进行筛选，剔除无效数据后，仅对有效数据。
b、根据Dobrovolsky's公式R＝100.43M(M为地震的震级)选取厄瓜多尔地震的研究区域，计算得到R＝2259.4km，故选取以震中为中心的边长为4518.8km为边长的正方形区域为研究区域，选取震前90天至震后30天为研究的时间范围，即2016年1月17日至5月16日为研究的时间范围。厄瓜多尔地震震中，地震影响区域及研究区域如图2所示，其中星号为地震震中，圆形区域为地震影响区域，正方形为研究区域。
c、选取Swarm卫星A星和C星的地方时夜晚数据，是计算Swarm A星和C星的地方时，LT(地方时)＝UTC(世界时)+α/15(α为地理经度)，由于电离层在白天受到干扰会产生异常，为了避免受到白天电离层活动的干扰难以解释，因此选取夜晚地方时为18：00-06：00的轨道，由于卫星沿轨飞行，测得的数据是随时间和空间变化的量，为了后续便于分析，对数据按轨道进行分析和异常提取，通过筛选，得到224条有效轨道的磁场Y分量数据，以2016年3月27日轨道3为例，Swarm A星和C星经过研究区域，且世界时对应相同、地方时为夜晚的轨道的Y分量磁场曲线如图3所示：
由图3可以看到，Swarm A星和C星的磁场Y分量随纬度的变化趋势几乎一致，但是由于幅值太大，由于地震产生的微小变化被淹没，因此需要后续的处理，才能提取地震前兆异常。
d、将磁场Y分量数据减去IGRF模型去除主磁场的影响，是由于卫星测量磁场为地磁场，主要包含主磁场、岩石圈磁场、变化磁场和感应磁场四部分。地磁场幅值为104nT量级，而地震前兆引起的变化主要表现在岩石圈磁场中，只有几nT甚至更小，相对于总体幅值来说非常微弱，直接对原始测量磁场进行操作则无法提取异常，因此为了后续的异常提取，通过国际给出的IGRF模型，分别计算Swarm A星和C星磁场Y分量对应的IGRF模型值，并分别从测量的磁场Y分量中减去模型值以去除主磁场的影响，得到Swarm A星和C星的Y分量磁场残余值，以2016年3月27日轨道3为例，减模型后的A星和C星的磁场Y分量的残余值曲线如图4所示：
由图4可以看到，经过减IGRF模型，一些微小的变化被显现出来，而且由于IGRF模型是利用球谐公式计算的理论主磁场值，是平滑的曲线，与测量值相减不会在数据中产生假异常。两颗星的磁场残余值曲线随纬度的低频变化趋势十分相似，同时曲线中还有一些叠加在低频趋势上的高频的微弱的抖动。
e、对减去模型后的Swarm A星和C星的Y分量磁场残余值进行主成分分析，将不同的信号分离到不同的成分中，主要包括：
将减模型后的Swarm A星和C星的Y分量磁场残余值按照时间序列表示为：
Bi＝[bi(1),bi(2)...,bi(m)]，i＝1,2,...,n
其中m为轨道长度，n取2，得到矩阵Y的表达式为：
第一、二列分别表示A星和C星的磁场Y分量数据，计算矩阵Y的协方差矩阵CY(2×2)的元素γuv，计算公式为：
其中，biu和biv分别为矩阵Y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值；
计算协方差矩阵的特征值和特征向量：
CY＝RΛRT其中，Λ(2×2)为从大到小排列的特征值对角矩阵，R(1×2)为对应特征值的特征向量，将特征值表示为λ1,λ2(λ1＞λ2)；
利用矩阵R将矩阵Y的原始数据进行线性投影得到主成分Φ，
Φ＝R·Y＝[Φ1,Φ2]T
其中，Φ1,Φ2为第1和第2个主成分(1×m)，主成分是原始数据投影到新的空间的数据，按照方差从大到小的顺序排列，并且各个主成分互不相关。
以2016年3月27日轨道3为例，经过主成分分析，得到的两个主成分如图5(a)和5(b)所示：从图可以看到第一主成分的幅值大，且频率低，与减模型后的残余值的低频变化趋势相似，而第二主成分的幅值小，且频率较高，与叠加在原始数据低频趋势上的高频抖动相似；通过主成分分析，将A星和C星的磁场Y分量按照方差排列分离到两个互不相关的主成分中。
研究的224条轨道的各个主成分对应的特征值及对应的地磁指数ap指数曲线如图6所示：由图可以看到，第一特征值与地磁指数多次同时达到峰值，意味着，当地磁指数高、地磁活跃时，相应的第一特征值的幅值也大，猜测第一主成分可能与地磁活动相关；而第二特征值的变化则与地磁活动的变化无明显相似。为了定量分析主成分与地磁活动的关系，进行以下步骤。
f、所述的计算各个主成分与地磁指数的相关系数，找到地磁活动相关的主成分并去除，包括将求得的所有轨道的各个主成分特征值表示为：
Ei＝[λi1,λi2,...,λip]，i＝1,2
其中p为研究的轨道个数，得到矩阵Γ＝[E1 T,E2 T]T，Γ的表达式为：
矩阵Γ：
计算两个时间序列x(t)和y(t)相关系数的公式为：
其中σx和σy为x(t)和y(t)的方差，和分别是x(t)和y(t)的均值,Cxy(k)是时滞为k时的两时间序列的协方差；rxy(k)为两个时间序列在时滞为k时的相关系数。
计算所有轨道两个主成分的特征值Ei与对应的地磁指数ap在时滞为0时的相关系数：
相关系数高的说明响应的主成分反应地磁活动情况，相关系数低的成分表明几乎与地磁活动无关，其中可能包含地震相关的信号。
经过计算，第一特征值与相应的地磁指数的相关系数达到60.45％，相关程度非常高，说明第一主成分主要反应地磁活动情况，由于电离层参数的变化主要受到地磁活动，太阳活动的影响。而第二主成分与地磁指数的相关系数为19％，相关程度非常小，因此第二主成分几乎不受到地磁活动的影响。经过主成分分析，将受地磁活动大的成分分离到第一主成分中，因此接下来只对第二主成分进行分析。
g、通过与地磁指数的相关系数的对比，发现第二主成分几乎不受地磁活动影响，对第2个主成分Φ2沿轨道进行差分得到：
dΦ2＝[dΦ21,dΦ22,...,dΦ2p]
其中p为研究的轨道个数。
求差分结果的偏度和峰度，利用偏度和峰度进行异常提取。
偏度和峰度是一种统计手段，偏度反应了数据分布的对称性，峰度反映了数据分布的厚尾性，当数据中出现极端值时，偏度和峰度会相应的增大。偏度的计算公式为：
峰度的计算公式为：
其中x为分析的数据，X为x的均值；
计算差分后的第2个主成分的偏度和峰度：
研究的224条轨道的偏度和峰度曲线如图7所示：
进一步的还包括，结合偏度和峰度定义偏峰度系数γ：
其中p为研究的轨道的个数。
偏峰度系数由偏度与峰度占总值的比重求和得到，可以同时反应偏度和峰度的变化情况，对于偏度和峰度幅值大的在该系数中都能体现出来。
求系数γ的均值μγ和标准差σγ，并设阈值为Thγ＝μγ+2×σγ，当轨道偏峰度系数大于设定阈值时，该轨道被认为是异常轨道。偏峰度系数曲线如图8所示：图中，竖直的虚线为地震当天，水平直线为设定的阈值，从图中可以看到地震附近的偏峰度系数明显高于平时，且随着地震的发生增大，震后减小恢复到正常水平。故将图中提取出来的偏峰度超出阈值的轨道作为厄瓜多尔地震的前兆异常进行输出。
本发明基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法，利用IGRF模型去除地磁主磁场的影响；通过对Swarm A星和C星的磁场Y分量数据进行主成分分析，通过各个主成分与地磁指数的对比，找到与地磁活动相关性高的成分并将其去除，只对剩余的主成分进行分析；通过偏度和峰度定义偏峰度系数对剩余的主成分进行异常提取，提取地震前兆异常。弥补了上述技术中针对单颗卫星分析，或只是对多颗卫星数据的异常结果结合分析，且只对静磁时间段的数据进行分析，不能充分利用Swarm卫星的A、C两颗星几乎并列飞行、同时测量磁场的优势，去除地磁活动干扰，全面提取地震前兆异常的不足。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法,其特在于,该方法包括：
a读取Swarm卫星A星和C星的磁场Y分量数据，根据数据的标志位Flags_B剔除有误数据；
b根据地震选择研究的时间范围、研究区域；
c选取Swarm A星和C星的地方时为夜晚的轨道数据；
d将磁场Y分量数据减去IGRF模型去除主磁场的影响；
e、对减去模型后的磁场Y分量残余值进行主成分分析，将不同的信号分离到不同的成分中；
f、计算各个主成分与地磁指数的相关系数，找到与地磁活动相关的主成分并去除；
g、对剩余的主成分求差分，求偏度和峰度，并定义偏峰度系数，利用偏峰度系数进行异常提取，输出异常结果。
2.按照权利要求1所述的方法,其特在于,
步骤b所述的根据地震选择研究的时间范围、研究区域，是选取地震前90天至地震后30天为研究时间范围，根据Dobrovolsky's公式R＝100.43M，其中M为地震的震级，计算地震的影响区域，并选择以地震为圆心，以R为半径的圆形区域的外接正方形区域为研究区域。
3.按照权利要求1所述的方法,其特征在于,步骤c所述的选取Swarm卫星A星和C星的地方时夜晚数据，是计算Swarm A星和C星的地方时，LT＝UTC+α/15，LT指的地方时，UTC指的是世界时，α为地理经度，夜晚指的是夜晚地方时为18：00-06：00。
4.按照权利要求1所述的方法,其特征在于,步骤d所述的将磁场Y分量数据减去IGRF模型去除主磁场的影响包括：通过国际给出的IGRF模型，分别计算Swarm A星和C星磁场Y分量对应的IGRF模型值，并分别从测量的磁场Y分量中减去模型值以去除主磁场的影响，得到Swarm A星和C星的Y分量磁场残余值。
5.按照权利要求1所述的方法,其特征在于,
步骤e所述的对减去模型后的Swarm A星和C星的Y分量磁场残余值进行主成分分析，将不同的信号分离到不同的成分中，具体包括：
将减模型后的Swarm A星和C星的磁场残余值按照时间序列表示为：
Bi＝[bi(1),bi(2)...,bi(m)]，i＝1,2,...,n
其中m为轨道长度，n为磁场数据维度，得到矩阵Y的表达式为：
计算矩阵Y的协方差矩阵CY(n×n)的元素γuv，计算公式为：
其中，biu和biv分别为矩阵Y中的第i行u列元素和第i行v列元素；和分别是第u列和第v列元素的均值；
计算协方差矩阵的特征值和特征向量：
CY＝RΛRT
其中，Λ(n×n)为从大到小排列的特征值对角矩阵，R(1×n)为对应特征值的特征向量，将特征值表示为λ1,λ2,...,λn(λ1＞λ2＞...＞λn)；
利用矩阵R将矩阵Y的原始数据进行线性投影得到主成分Φ，
Φ＝R·Y＝[Φ1,Φ2,...,Φn]T
其中，Φ1,Φ2,...,Φn为第1至第n个主成分(1×m)，主成分是原始数据投影到新的空间的数据，按照方差从大到小的顺序排列，并且各个主成分互不相关。
6.按照权利要求1所述的方法,其特征在于,
步骤f所述的计算各个主成分与地磁指数的相关系数，找到地磁活动相关的主成分并去除，包括将求得的所有轨道的各个主成分的特征值表示为：
Ei＝[λi1,λi2,...,λip]，i＝1,2,...,n
其中p为研究的轨道个数，n为数据维度，得到矩阵Γ＝[E1 T,E2 T,...,En T]T，Γ的表达式为：
矩阵Γ：
计算时间序列x(t)和y(t)相关系数的公式为：
其中σx和σy为x(t)和y(t)的方差，和分别是x(t)和y(t)的均值,Cxy(k)是时滞为k时的两时间序列的协方差；rxy(k)为两个时间序列在时滞为k时的相关系数；
根据公式计算所有轨道各个主成分的特征值Ei与对应的地磁指数ap在时滞为0时的相关系数：
相关系数高的说明相应的主成分主要反应地磁活动情况，相关系数低的成分表明几乎与地磁活动无关，通过去除与地磁指数相关性高的主成分，以去除地磁活动的影响。
7.按照权利要求1所述的方法,其特征在于,
步骤g对剩余的主成分求差分，求偏度和峰度，包括：设上一步去除地磁活动影响后，与地磁指数相关性低的主成分为第q个主成分，被认为几乎不受地磁活动的干扰，则对第q个主成分Φq沿轨道进行差分得到：
dΦq＝[dΦq1,dΦq2,...,dΦqp]
其中p为研究的轨道个数。
8.按照权利要求1或7所述的方法,其特征在于，求差分结果的偏度和峰度，利用偏度和峰度定义偏峰度系数包括：
偏度的计算公式为：
峰度的计算公式为：
其中x为分析的数据，X为x的均值；
计算差分后的第q个主成分的偏度和峰度：
9.按照权利要求8所述的方法,其特征在于，结合偏度和峰度定义偏峰度系数γ：
其中p为研究的轨道的个数。
10.按照权利要求8所述的方法,其特征在于，求系数γ的均值μγ和标准差σγ，并设阈值为Thγ＝μγ+2×σγ，当轨道偏峰度系数大于设定阈值时，该轨道被认为是异常轨道。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910261107.4A CN110068857B (zh) | 2019-04-02 | 2019-04-02 | 基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910261107.4A CN110068857B (zh) | 2019-04-02 | 2019-04-02 | 基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110068857A true CN110068857A (zh) | 2019-07-30 |
| CN110068857B CN110068857B (zh) | 2020-02-04 |
Family
ID=67367015
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910261107.4A Active CN110068857B (zh) | 2019-04-02 | 2019-04-02 | 基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110068857B (zh) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111814699A (zh) * | 2020-07-13 | 2020-10-23 | 中国地震局地震预测研究所 | 一种面向swarm电磁卫星数据的深度学习地震预测方法 |
| CN112305606A (zh) * | 2020-10-16 | 2021-02-02 | 宁夏回族自治区地震局 | 一种基于自然正交函数展开的地震活动场分析方法 |
| CN112327371A (zh) * | 2020-11-06 | 2021-02-05 | 吉林大学 | 基于变分模态分解的卫星磁场数据时变背景场建立方法 |
| CN113361476A (zh) * | 2021-07-02 | 2021-09-07 | 中国地震局地震预测研究所 | 一种基于人工智能技术的张衡一号震前异常信号识别方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7277797B1 (en) * | 2005-03-29 | 2007-10-02 | Kunitsyn Viatcheslav E | Prediction system and method |
| CN106021710A (zh) * | 2016-05-13 | 2016-10-12 | 南京航空航天大学 | 基于大气电离层参数的震前卫星轨道异常识别方法 |
| CN106918836A (zh) * | 2017-03-31 | 2017-07-04 | 吉林大学 | 基于主成分分析的钻孔应变数据异常提取方法 |
| CN107356969A (zh) * | 2017-09-06 | 2017-11-17 | 四川易利数字城市科技有限公司 | 一种基于卫星热红外数据及gis的地震前兆分析方法 |
| CN109031403A (zh) * | 2018-08-20 | 2018-12-18 | 吉林大学 | 一种基于s-k特征的钻孔应变数据异常提取方法 |
-
2019
- 2019-04-02 CN CN201910261107.4A patent/CN110068857B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7277797B1 (en) * | 2005-03-29 | 2007-10-02 | Kunitsyn Viatcheslav E | Prediction system and method |
| CN106021710A (zh) * | 2016-05-13 | 2016-10-12 | 南京航空航天大学 | 基于大气电离层参数的震前卫星轨道异常识别方法 |
| CN106918836A (zh) * | 2017-03-31 | 2017-07-04 | 吉林大学 | 基于主成分分析的钻孔应变数据异常提取方法 |
| CN107356969A (zh) * | 2017-09-06 | 2017-11-17 | 四川易利数字城市科技有限公司 | 一种基于卫星热红外数据及gis的地震前兆分析方法 |
| CN109031403A (zh) * | 2018-08-20 | 2018-12-18 | 吉林大学 | 一种基于s-k特征的钻孔应变数据异常提取方法 |
Non-Patent Citations (2)
| Title |
|---|
| 宋奕瑶等: "基于SOM网的DEMETER卫星电场数据聚类异常分析", 《软件导刊》 * |
| 朱凯光等: "基于主成分分析的航空电磁数据噪声去除方法", 《中国有色金属学报》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111814699A (zh) * | 2020-07-13 | 2020-10-23 | 中国地震局地震预测研究所 | 一种面向swarm电磁卫星数据的深度学习地震预测方法 |
| CN111814699B (zh) * | 2020-07-13 | 2023-07-28 | 中国地震局地震预测研究所 | 一种面向swarm电磁卫星数据的深度学习地震预测方法 |
| CN112305606A (zh) * | 2020-10-16 | 2021-02-02 | 宁夏回族自治区地震局 | 一种基于自然正交函数展开的地震活动场分析方法 |
| CN112327371A (zh) * | 2020-11-06 | 2021-02-05 | 吉林大学 | 基于变分模态分解的卫星磁场数据时变背景场建立方法 |
| CN112327371B (zh) * | 2020-11-06 | 2021-07-30 | 吉林大学 | 基于变分模态分解的卫星磁场数据时变背景场建立方法 |
| CN113361476A (zh) * | 2021-07-02 | 2021-09-07 | 中国地震局地震预测研究所 | 一种基于人工智能技术的张衡一号震前异常信号识别方法 |
| CN113361476B (zh) * | 2021-07-02 | 2023-07-25 | 中国地震局地震预测研究所 | 一种基于人工智能技术的张衡一号震前异常信号识别方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110068857B (zh) | 2020-02-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110068857A (zh) | 基于主成分分析的Swarm双星磁场数据地震前兆异常提取方法 | |
| CN106021710B (zh) | 基于大气电离层参数的震前卫星轨道异常识别方法 | |
| CN109003422A (zh) | 用于山体滑坡的监测数据处理方法和山体滑坡预报方法 | |
| Richmond et al. | Mapping electrodynamic features of the high‐latitude ionosphere from localized observations: Combined incoherent‐scatter radar and magnetometer measurements for January 18–19, 1984 | |
| Knudsen et al. | Samples from the Lomonosov Ridge place new constraints on the geological evolution of the Arctic Ocean | |
| Hartog et al. | Subduction‐induced strain in the upper mantle east of the Mendocino triple junction, California | |
| CN109740453A (zh) | 一种基于小波变换的卫星磁场数据地震前兆异常提取方法 | |
| Buland et al. | Matched filtering for the seismic moment tensor | |
| He et al. | A study to investigate the relationship between ionospheric disturbance and seismic activity based on Swarm satellite data | |
| Darbyshire et al. | A first detailed look at the Greenland lithosphere and upper mantle, using Rayleigh wave tomography | |
| He et al. | Is there a one-to-one correspondence between ionospheric anomalies and large earthquakes along Longmenshan faults? | |
| Varga et al. | Identification of Saharan dust particles in Pleistocene dune sand-paleosol sequences of Fuerteventura (Canary Islands) | |
| Husebye et al. | Tomographical mapping of the lithosphere and asthenosphere beneath southern Scandinavia and adjacent areas | |
| Claoué-Long et al. | The duration of the Strangways Event in central Australia: Evidence for prolonged deep crust processes | |
| Haines-Young | The use of remotely-sensed satellite imagery for landscape classification in Wales (UK) | |
| Lan et al. | An analysis of mechanical constraints when using superconducting gravimeters for far-field pre-seismic anomaly detection | |
| Percival et al. | Detiding DART® buoy data for real-time extraction of source coefficients for operational tsunami forecasting | |
| Chen et al. | Lancang—Gengma earthquake: a preliminary report on the November 6, 1988, event and its aftershocks | |
| Langel et al. | A satellite magnetic anomaly map of Greenland | |
| Thorarinsson et al. | Directional spectral analysis and filtering of geophysical maps | |
| Eckstaller et al. | The geophysics observatory at Neumayer stations (GvN and NM-II) Antarctica | |
| Arias et al. | Rapid source characterization of the Maule earthquake using Prompt Elasto‐Gravity Signals | |
| Halkin et al. | Combined automated and off-line computer processing system for seismic monitoring with small aperture arrays | |
| Hobara et al. | Ionospheric perturbation in association with seismic activity. A statistical study | |
| Moldovan et al. | The geomagnetic method on precursory phenomena associated with 2004 significant intermediate-depth Vrancea seismic activity |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |