张衡一号卫星感应式磁力仪数据处理方法及系统与流程

1.本发明涉及卫星遥感技术领域，特别是涉及一种张衡一号卫星感应式磁力仪数据处理方法及系统。


背景技术：

[0002]“张衡一号”卫星是中国第一颗自主研制的地球物理场卫星，感应式磁力仪是卫星上搭载的八种载荷之一。


技术实现要素：

[0003]
本发明要解决的技术问题是提供一种张衡一号卫星感应式磁力仪数据处理方法及系统，能够将卫星上感应式磁力仪的原始探测数据进行处理，生成实际可用的数据产品。
[0004]
为解决上述技术问题，本发明提供了一种张衡一号卫星感应式磁力仪数据处理方法，所述方法包括：感应式磁力仪0级数据进行数据格式转换及参数校正，生成1级数据产品；根据星下点的位置、卫星姿态等信息，利用坐标变换生成带有地理和地磁坐标信息的三分量波形和功率谱物理量数据，并标注地震、磁情指数等辅助信息，生成2级数据产品；选取2级数据产品中的特征频段波形或功率谱数据作为输入，对变化磁场数据进行重采样生成全球范围及中国区域的重访轨道观测数据的时间序列产品，标注地震及磁情指数信息，获得 3级数据产品；选取2级数据产品中的特征频段波形或功率谱数据作为输入，滑动选取全球范围和中国区域上空变化磁场波形和特定频段功率谱数据，得到空间分布背景场及动态变化幅度，生成4级数据产品；其中，1级数据产品包括：每条轨道记录的磁力仪三分量在ulf、elf、vlf频段内波形和功率谱物理量数据随时间的变化；2级数据产品包括：每条轨道记录的感应式磁力仪在 ulf、elf、vlf频段内三分量波形和功率谱物理量数据随时间和空间的变化趋势；3级数据产品包括：在2级数据的基础上，对变化磁场数据进行重采样生成全球范围及中国区域的重访轨道观测数据的时间序列产品，标注地震及磁情指数信息；4级数据产品包括：全球磁场空间分布的图像的文件名及其包含的子图。
[0005]
在一些实施方式中，1级、2级、3级及4级数据产品均包括：科学数据、图像产品和数据处理报告。
[0006]
在一些实施方式中，感应式磁力仪0级数据进行数据格式转换及参数校正，生成1级数据产品，包括：对二进制的0级数据进行格式转换，根据科学数据的类型判断是否需要进行傅里叶变换，并在频域范围逐条谱线进行正交性校正，利用磁传感器各频段的传递函数进行校正和物理量转换，如果是波形数据还要经过傅里叶逆变换转到时域范围内，频谱数据转换为功率谱密度数据，最终生成校正后的三分量波形和功率谱物量数据，及相应的处理报告和快视图像产品。
[0007]
在一些实施方式中，在校正过程中，提取在轨定标数据与地面试验获取的定标信号校准参数进行误差分析判断感应式磁力仪的工作状态。
[0008]
在一些实施方式中，根据星下点的位置、卫星姿态等信息，利用坐标变换生成带有
地理和地磁坐标信息的三分量波形和功率谱物理量数据，并标注地震、磁情指数等辅助信息，生成2级数据产品，包括：利用传感器坐标系和卫星坐标系之间的关系，将三分量磁场数据转换到卫星坐标系下；根据转换矩阵，由卫星坐标系sbc转换到地理坐标系geo。
[0009]
在一些实施方式中，选取2级数据产品中的特征频段波形或功率谱数据作为输入，对变化磁场数据进行重采样生成全球范围及中国区域的重访轨道观测数据的时间序列产品，标注地震及磁情指数信息，获得3级数据产品，包括：根据cses卫星轨道重访周期，利用当前轨道号，可获取前30天相应的重访轨道数据；针对半轨道的2a级波形和功率谱数据，按纬度间隔0.05
°
和0.1
°
进行划分网格；通过时序分析，求取每个网格里的中值；根据求取的中值，绘图与产出数据处理报告。
[0010]
在一些实施方式中，根据求取的中值，绘图与产出数据处理报告，包括：将当前5天的重访轨道数据，以及由前30天的重访轨道数据计算得到的、序列，即重采样的时间序列，区分三分量和三频段绘制成图，同时，在图标注当前轨道时间的磁情指数；绘出每条轨道正负5个单位的经度区域范围的前30 内天发生的震级大于5级、6级及7级以上的地震事件序列；绘出当前轨道所在全球位置，以及标注前7内天发生的震级大于6级的地震事件；生成数据处理报告。
[0011]
在一些实施方式中，选取2级数据产品中的特征频段波形或功率谱数据作为输入，滑动选取全球范围和中国区域上空变化磁场波形和特定频段功率谱数据，得到空间分布背景场及动态变化幅度，生成4级数据产品，包括：按照 2.5
°
(纬度)
×5°
(经度)的空间间隔将全球划分为若干网格；执行空间分布分析；加入地震事件和磁情指数，生成4级数据处理报告。
[0012]
在一些实施方式中，执行空间分布分析，包括：按照30天时间间隔及5
ꢀ°
(经度)
×
2.5
°
(纬度)空间间隔，滑动选取全球和中国区域上空的磁场2级数据，计算每个空间间隔内所有轨道的中值b
m
、四分位点(及)及分位差iqr，以及异常判别值b
c
；并通过对各个间隔内中值b
m
进行插值，得到全球及全国范围空间分布背景场；对于当前5天变化磁场数据，计算每一个网格中变化磁场的中值其及四分位差；计算每天每个间隔内中值与之前30天的背景中值的差值，并使用每个格子内的最大值进行归一化，并通过对各个间隔归一化后的值进行插值，得到全球及全国的每天动态变化空间分布图。
[0013]
此外，本发明还提供了一种张衡一号卫星感应式磁力仪数据处理系统，所述系统包括：一个或多个处理器；存储装置，用于存储一个或多个程序，当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现根据前文所述的张衡一号卫星感应式磁力仪数据处理方法。
[0014]
采用这样的技术方案后，本发明至少具有以下优点：
[0015]
本发明提供的张衡一号卫星感应式磁力仪数据处理方法及系统，通过对感应式磁力仪原始数据的格式转换、数据定标、物理量校正、坐标系变换、反演等一系列运算，生成数据精度及格式均符合要求的数据产品。
附图说明
[0016]
上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。
[0017]
图1是1～4级数据命名规则的原理示意图；
[0018]
图2是科学数据处理报告的命名示意图；
[0019]
图3是科学数据产品图像的命名示意图；
[0020]
图4是感应式磁力仪1级数据hdf5文件格式示意图；
[0021]
图5a是感应式磁力仪vlf波形示意图；
[0022]
图5b是感应式磁力仪vlf功率谱数据示意图；
[0023]
图6是感应式磁力仪2级科学数据hdf5格式示意图；
[0024]
图7a是感应式磁力仪2级数据波形示意图；
[0025]
图7b是感应式磁力仪2级数据功率谱数据示意图；
[0026]
图8是4级数据产品示意图；
[0027]
图9是感应式磁力仪0
‑
4级数据处理流程示意图；
[0028]
图10是感应式磁力仪1级数据处理流程图；
[0029]
图11是感应式磁力仪2级数据处理流程图；
[0030]
图12是3级数据处理流程图；
[0031]
图13是4级数据处理流程图。
具体实施方式
[0032]
以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。
[0033]
首先介绍下文中需要使用的缩略语。这些缩略语的含义如表1所示：
[0034]
表1缩略语
[0035]
[0036][0037]
1感应式磁力仪简介与数据产出
[0038]
感应式磁力仪(search coil magnetometer，scm)是一种基于法拉第电磁感应定
律的星载仪器，主要用于获取cses卫星轨道位置处10hz～20khz频率范围内变化磁场矢量的波形和频谱信息。
[0039]
感应式磁力仪可以依据星务管理分系统注入指令可进行工作模式切换。不同工作模式下感应式磁力仪产出的数据产品如表2所示。
[0040]
表2感应式磁力仪数据产出表
[0041][0042]
2标准数据产品
[0043]
2.1分级定义
[0044]
按照《卫星对地观测数据产品分类分级规则》(gb/t 32453
‑
2015)，以及《电磁监测卫星数据产品分级及定义》，感应式磁力仪观测数据包括了0级，1 级、2级，3级和4级数据产品，感应式磁力仪标准数据产品定义见表3。
[0045]
表3 scm数据产品分级定义及数据产品一览表
[0046]
[0047][0048]
2.2命名规则
[0049]
为了使各级数据产品检索查询方便，各级数据产品的命名应包含卫星名称、载荷名称、轨道号、数据起止时间等其它必要的标识。感应式磁力仪0～4 级数据产品、图像和处理报告的命名规则示例如图1所示。
[0050]
其中：
[0051]
(1)卫星名称(4位字符)：用cses表示；
[0052]
(2)卫星编号(1位数字)：从1开始，依次递增；
[0053]
(3)载荷编码(3个字符)：分别为hpm、scm、efd、lap、rpa、gro、tbb、 hep分别代表8个载荷；
[0054]
(4)载荷序号(1位数字)：用于区分一个载荷装载多个同类测项探头的情况，从1开始，依次递增，1代表编号为“1”的探头，2代表编号为“2”的探头，依次递推，如果不区分探头用0代表；除高能粒子载荷外，其他载荷均为0；对高能粒子载荷，1～4依次代表低能探头、高能探头、意大利载荷和 x射线；
[0055]
(5)数据分级编码(1位符号2位数字)：从左到右，第一位是“l”，右边两位表示数据级别，分别用00、01、02/2a、03和04表示0
‑
4级；
[0056]
(6)观测对象编码(2位数字)：依据《地震电磁卫星测项分类与代码》(送审稿)，设置观测对象编码。对0级数据，观测对象编码为00；
[0057]
(7)轨道号(5位数字)：从星上时间为准，00001开始，用于按轨道组织数据文件，无法标注轨道号的数据产品以“00000”表示；
[0058]
(8)升降轨标志(1位数字)：升轨为1，降轨为0；
[0059]
(9)数据起始时间，采用14位数字表示，其中年(4位)、月(2位)、日 (2位)、时(2位)、分(2位)、秒(2位)；
[0060]
(10)数据结束时间，格式同(9)；
[0061]
(11)接收站编码(3位数字)：此项专门为三频信标接收机预留，对于感应式磁力仪不涉及地面接收站信息，因此标记为000；
[0062]
(12)文件扩展名：当文件扩展名为dat时，代表以二进制格式存储的数据文件；当文件扩展名为h5时，表示数据以hdf格式存储。
[0063]
对于科学数据处理报告命名，处理报告以txt为文件后缀，以ascii码格式存储，其命名是在科学数据命名的基础上，通过修改部分字段并添加相应的“_rp”后缀组成，如图2所示。
[0064]
对于科学数据产品图像命名与处理报告文件命名类似，它是在科学数据命名的基础上，通过修改部分字段并添加相应的“_xx.png”后缀组成，示例如图3所示。
[0065]
(13)表示观测对象拓展代码，由两位字符组成，第一位用于区分不同区域：1表示中国地区，2表示全球，0表示不区分区域；第二位用于区分同一载荷的多张图像，从1开始依次递增。
[0066]
2.3各级数据产品介绍
[0067]
2.3.11级数据产品
[0068]
将0级数据进行格式转换、正交性校正及传递函数校正后生成按时间排列的感应式磁力仪三分量波形和频谱的物理量数据，处理过程中所需辅助数据为标定校正数据。1级数据产品包括：科学数据、图像产品和数据处理报告。
[0069]
(1)1级科学数据
[0070]
●
ulf：详查、巡查和定标工作模式下10～200hz频带内三分量波形数据；
[0071]
●
elf：详查、巡查和定标工作模式下200～2.2khz频带内三分量波形数据；
[0072]
●
vlf：详查、巡查和定标工作模式下1.8k～20khz频带内三分量波形数据；巡查和定标工作模式下1.8k～20khz频带内三分量功率谱数据。
[0073]
表4感应式磁力仪1级数据文件属性说明表
[0074][0075][0076]
表5感应式磁力仪1级数据表格式说明表
[0077]
[0078][0079]
(2)1级图像产品
[0080]
每条轨道记录的磁力仪三分量在ulf、elf、vlf频段内波形和功率谱物理量数据随
时间的变化；各个参量展示如图5所示。
[0081]
(3)1级数据处理报告
[0082]
主要记录1级科学数据的在处理工程中的相关信息，包括数据处理时间、程序版本信息、数据的连续性描述、载荷工作状态描述等信息。组成要素：
[0083]
①
处理软件版本号：ver0.1
[0084]
②
起始时间：
[0085]
yyyymmdd hh:mm:ss.zzz
[0086]
③
输入数据：0级数据产品文件名
[0087]
辅助数据
[0088]
④
校正参数：正交性校正数据；三分量三频段频率响应电压
‑
磁场转换参数；在轨标定信号校准参数；
[0089]
⑤
载荷工作状态描述
[0090]
●
电子学箱温度：正常
[0091]
●
监控电压数据
[0092]
●
3.3v电压：数据包号
‑
错误情况，0
‑
＜阈值范围1
‑
＞阈值范围2
‑
无效数值；
[0093]
●
2.5v电压：数据包号
‑
错误情况，0
‑
＜阈值范围1
‑
＞阈值范围2
‑
无效数值；
[0094]
●
1.9v电压：数据包号
‑
错误情况，0
‑
＜阈值范围1
‑
＞阈值范围2
‑
无效数值；
[0095]
●
625hz信号检测电压：数据包号
‑
错误情况，0
‑
＜阈值范围1
‑
＞阈值范围2
‑
无效数值；
[0096]
●
10khz信号检测电压：数据包号
‑
错误情况，0
‑
＜阈值范围1
‑
＞阈值范围2
‑
无效数值；
[0097]
⑥
数据处理过程
[0098]
●
正常
[0099]
●
缺数：数据包号
‑
出错项目；
[0100]
●
数据文件损坏：包序号
‑
时间；
[0101]
●
处理过程用到的方法：dft计算，idft变换；fft计算，ifft变换；正交性校正；传递函数校正；
[0102]
⑦
在轨定标数据判定：正常/异常；
[0103]
●
625hz 时间 vx δx vy δy vz δz；
[0104]
●
10khz 时间 vx δx vy δy vz δz；
[0105]
(*vx、vy、vz为三分量定标信号数据；δx、δy、δz为三分量定标信号数据误差值；)
[0106]
⑧
结束时间
[0107]
yyyymmdd hh:mm:ss.zzz
[0108]
⑨
输出数据：1级科学数据文件名
[0109]
2.3.2 2级数据产品
[0110]
2级数据是利用坐标转换矩阵将1级数据从传感器坐标系变换到地理坐标系，生成带有地理、地磁坐标和姿态信息的三分量波形和功率谱物理量数据。处理过程中所需辅助数据有传感器
‑
卫星坐标系转换矩阵、星地坐标转换矩阵等。2级数据产品包括：科学数据、
图像产品和数据处理报告。
[0111]
(1)2级科学数据
[0112]
●
ulf：详查、巡查和定标工作模式下10～200hz频带内三分量波形数据；
[0113]
●
elf：详查、巡查和定标工作模式下200～2.2khz频带内三分量波形数据
[0114]
●
vlf：详查、巡查和定标工作模式下1.8k～20khz频带内三分量波形数据；巡查和定标工作模式下1.8k～20khz频带内三分量功率谱数据。
[0115]
表6感应式磁力仪2级数据文件属性说明表
[0116][0117][0118]
表7感应式磁力仪2级数据表格式说明表
[0119]
[0120][0121]
(2)2级图像产品
[0122]
每条轨道记录的感应式磁力仪在ulf、elf、vlf频段内三分量波形和功率谱物理量数据随时间和空间的变化趋势。展示图如图7a及图7b所示，包括 ulf、elf和vlf x、y、z波形以及vlf功率频谱。
[0123]
(3)2级数据处理报告
[0124]
主要记录2级科学数据的在处理工程中的相关信息，包括数据处理时间、程序版本信息、坐标转换矩阵、数据质量等信息。组成要素：
[0125]
①
处理软件版本号:ver0.1
[0126]
②
起始时间：yyyymmdd hh:mm:ss.zzz
[0127]
输入数据：列出1级数据产品文件名。
[0128]
③
辅助数据信息
[0129]
●
校准参数：传感器卫星坐标转换矩阵；星地转换矩阵；
[0130]
④
处理过程
[0131]
●
处理过程正常
[0132]
●
异常现象及可能的原因
[0133]
●
数据文件损坏：包序号
‑
时间；
[0134]
●
处理过程用到的方法：传感器
‑
卫星坐标转换；卫星
‑
地理坐标转换；地理
‑
地磁坐标转换；
[0135]
⑤
输函数是否修正：0
‑
否；
[0136]1‑
是：修正开始时间、传递函数修正系数；
[0137]
⑥
结束时间:yyyymmdd hh:mm:ss.zzz
[0138]
⑦
输出数据:列出2级数据产品文件名。
[0139]
2.3.3 3级数据产品
[0140]
感应式磁力仪3级数据是在2级数据的基础上，对变化磁场数据进行重采样生成全球范围及中国区域的重访轨道观测数据的时间序列产品，标注地震及磁情指数信息。3级数据产品包括：科学数据、图像产品和处理报告。
[0141]
(1)2级科学数据
[0142]
3级科学数据文件内容包括：
[0143]
当前轨道在网格区域内的中值、均值、标准差等统计量
[0144]
一定周期内重访数据在网格区域内的中值，上下四分位值及其二者差值 (即四分位差)、标准差等统计量
[0145]
超出上下界限的标识(0为正常；超出边界则标注出当前观测值与上下边界的差值)
[0146]
这些数据的格式见表8和表9。
[0147]
表8感应式磁力仪3级数据文件属性说明
[0148]
序号属性名称属性内容备注1payloadid仪器代码 2orbitnum轨道号 3orbitflag升降轨标识升轨1、降轨04psd_tresvlf功率谱的时间分辨率2秒5psd_fresvlf功率谱的频率分辨率12.5hz6softversion程序版本号0.1
[0149]
表9感应式磁力仪3级数据表格式说明(以四分位算法为例)
[0150]
[0151][0152]
(2)3级图像产品
[0153]
(3)3级数据处理报告
[0154]
感应式磁力仪3级数据处理报告组成要素如下：
[0155]
①
处理软件版本号：ver0.1
[0156]
②
起始时间：yyyymmdd hh:mm:ss.zzz
[0157]
③
输入数据：(a)当前轨道数据，列出其文件名；(b)前5条重访轨道的 3级数据(以卫星30天探测时间跨度为例)，列出其文件名；缺失数据情况。
[0158]
④
辅助数据:
[0159]
(a)震目录：yyyymmdd.earthquakes.txt(b)磁情指数：
[0160]
(b)当前5天的kp(3小时)、dst(1小时)、f107(1天)、ae(1分钟沿纬度按照0.1
°
间隔采样，进而采样逐日滑动的四分位差方法；
[0161]
⑤
处理方法
[0162]
⑥
结束时间：yyyymmdd hh:mm:ss.zzz
[0163]
⑦
输出产品：列出3级数据产品文件名
[0164]
⑧
处理结果
[0165]
(a)处理过程正常
[0166]
(b)处理过程中每个格子缺数情况；反演中存在问题；
[0167]
⑨
异常情况说明
[0168]
(a)连续3个点，在2020年3月23日01时34分至2020年3月23日01 时36分，星下点经纬度范围为(23
°‑
33
°
n，110
°‑
112
°
e)出现超界现象。 (b)建议进一步分析。
[0169]
2.3.4 4级数据产品
[0170]
感应式磁力仪4级数据按若干个重访周期时间跨度，以一定的经纬度空间间隔，滑动选取全球范围和中国区域上空变化磁场波形和特定频段功率谱数据，得到空间分布背景场；以当前5天的观测得到变化磁场波形和特定频段功率谱的空间分布，减去之前一定时间段内(下文以30天为例)空间分布背景场，获得动态变化幅度。
[0171]
(1)4级数据产品
[0172]
表10感应式磁力仪4级数据文件属性说明
[0173]
序号属性名称属性内容备注1payloadid仪器代码 2softversion程序版本号0.13orbittype轨道类别 [0174]
表11感应式磁力仪4级数据表格式说明
[0175]
[0176]
[0177][0178]
以x分量，625hz频点，降轨数据为例，说明4级数据产品——全球磁场空间分布的图像的文件名及其包含的子图如图8所示，具体内容介绍如下。
[0179]
子图1：x分量当前5天，625hz功率谱密度幅度空间分布图；
[0180]
子图2：x分量前30天，625hz功率谱密度幅度空间分布图；
[0181]
子图3：子图1减去子图2之后的差值,并标注30天内大于6级以上的地震事件；
[0182]
子图4：y分量当前5天，625hz功率谱密度幅度空间分布图；
[0183]
子图5：y分量前30天，625hz功率谱密度幅度空间分布图；
[0184]
子图6：子图4减去子图5之后的差值,并标注30天内大于6级以上的地震事件；
[0185]
子图7：z分量当前5天，625hz功率谱密度幅度空间分布图；
[0186]
子图8：z分量前30天，625hz功率谱密度幅度空间分布图；
[0187]
子图9：子图7减去子图8之后的差值,并标注30天内大于6级以上的地震事件。
[0188]
对于中国区域的图像与全球图像内容类似。对于ulf和elf重采样后的波形幅值全球分布图，与特定频点的功率谱密度全球分布图类似。
[0189]
(3)4级数据处理报告
[0190]
感应式磁力仪4级数据处理报告组成要素如下：
[0191]
①
处理软件版本号：ver0.1
[0192]
②
起始时间:yyyymmdd hh:mm:ss.zzz
[0193]
③
输入数据
[0194]
(a)前35天内的所有2级数据文件
[0195]
(b)缺失数据情况
[0196]
④
辅助数据
[0197]
震目录：yyyymmdd_earthquakes.txt
[0198]
磁情指数：当前5天的kp(3小时)、dst(1小时)、f107(1天)、ae(1 分钟).
[0199]
⑤
处理方法
[0200]
⑥
结束时间:yyyymmdd hh:mm:ss.zzz
[0201]
⑦
输出数据:4级数据文件名
[0202]
⑧
处理结果
[0203]
(a)处理过程正常；
[0204]
(b)处理过程发生数据冗余溢出，程序终止因发生冗余错误导致程序中断，未生成相关图片和分析结果
[0205]
⑨
异常情况说明:举例如
①
在2020年3月23日01时34分至2020年3 月23日01时36分，星下点经纬度范围为(23
°‑
33
°
n，110
°‑
112
°
e)出现超界现象，超出值为：。此时dst＝
‑
60nt，kp＝4，f107最大值＝150。
[0206]
3数据处理流程
[0207]
感应式磁力仪0级数据首先进行数据格式转换，如果是波形数据需经过傅里叶变换到频域内进行定标参数校正再进行傅里叶逆变换生成校正后的波形数据，而频谱数据直接在频域内进行校正生成校正后进行功率谱计算生成功率谱密度数据，同时提取在轨定标数据进行误差分析判断感应式磁力仪的工作状态，最终生成1级数据产品；根据星下点的位置、卫星姿态等信息，利用坐标变换生成带有地理和地磁坐标信息的三分量波形和功率谱物理量数据，并标注地震、磁情指数等辅助信息，生成2级数据产品。选取的特征频段波形或功率谱数据作为输入，分别获得对应的3级和4级数据产品。
[0208]
4数据处理方法
[0209]
4.1 0级数据生成1级数据
[0210]
图10是感应式磁力仪0级数据生成1级数据的处理流程，具体为：对二进制的0级数据进行格式转换，根据科学数据的类型判断是否需要进行傅里叶变换，并在频域范围逐条谱线进行正交性校正，利用磁传感器各频段的传递函数进行校正和物理量转换，如果是波形数据还要经过傅里叶逆变换转到时域范围内，频谱数据转换为功率谱密度数据，最终生成校正后的三分量波形和功率谱物量数据，及相应的处理报告和快视图像产品。在校正过程中，提取在轨定标数据与地面试验获取的定标信号校准参数进行误差分析判断感应式磁力仪的工作状态。
[0211]
4.1.1二进制转换成十进制数据
[0212]
将0级处理生成的二进制科学数据(ulf、elf的波形数据，vlf的波形和频谱相位数据)和自采集的工程参数换为十进制数据。
[0213]
(1)感应式磁力仪科学数据格式转换
[0214]
感应式磁力仪科学数据由ad芯片采集，数据存储为16位带符号二进制数，最高位为符号位，0表示正数，1表示负数，它能够表示
‑
(2
n
‑1‑
1)
…
(2
n
‑1‑
1)范围 (即
‑
32767～32767)的十进制数值，其中n＝16。根据ad采集芯片所设定的测量电压量程
±
10v，那么二进制数值所代表的电量测量精度为10v/32767。因为ad芯片采集的物理量量程已设定为
±
10v，所以为了将其采集的二进制数值 (为了表达简洁使用下文使用16进制写出二进制数值)与10进制电压值对应起来，那么16进制的原码0000对应的电压值为0v，而16进制的7fff对应的电压值为 10v。负电压采用补码表示，即16进制的8001对应
‑
10v。根据此对应关系，将二进制数值hh转为十进制的电压v的转换公式如式(4
‑
1)和式 (4
‑
2)。
[0215]
正数转换公式为：
[0216][0217]
单位为伏特(v)。
[0218]
负数转换公式为：
[0219][0220]
单位为伏特(v)。
[0221]
感应式磁力仪的科学数据分为三种类型：波形数据、频谱幅度数据和相位数据。其数据格式转换方法分别如下：
[0222]
●
波形数据
[0223]
波形数据为16位有符号整数，最高位为符号位：“0”表示正；“1”表示负。则正数转换公式为式(4
‑
1)，负数转换公式为式(4
‑
2)。
[0224]
●
频谱幅度数据
[0225]
幅度数据为16位有符号整数，最高位为符号位，同时用于表示幅度数据是否经过仪器内部放大：“0”表示正，无放大；“1”表示负，放大
‑
100倍。则正数无放大，转换公式为(4
‑
1)；负数为经过仪器放大
‑
100后补码表示的数据，转换公式为式(4
‑
3)。
[0226][0227]
单位为伏特(v)。如果数据较小，数据保存时程序会进行100倍放大(多取数据位数)以提高数据精确度，并在最高位以”1”标示，因为是16为带符号数据，转换为十进制后数据值除以
‑
100。
[0228]
●
频谱相位数据
[0229]
相位数据为16位有符号数据，最高位为符号位：“0”表示正，“1”表示负，该数值经过仪器内部1000倍的放大。
[0230]
正数转换公式为：
[0231][0232]
单位为度(
°
)。
[0233]
负数转换公式为：
[0234][0235]
单位为度(
°
)。其中，为圆周率，取值为3.14159。因为计算出的相位数据比较小，在数据保存时程序进行1000倍放大(多取数据位数)以提高数据精确度。
[0236]
4.1.2傅里叶变换
[0237]
波形数据的校正主要是在频域内进行的，因此需要将波形数据进行傅里叶变换到频域内，再逐个谱线利进行幅度校正和物理量转换，最后利用傅里叶逆变换将信号从频域转换到时域中，从而生成校正后的波形数据。
[0238]
(1)ulf、elf波形数据的离散傅里叶变换(dft)
[0239][0240]
其中，x(n)为ulf、elf时域的离散信号数据；n为连续数据的个数，根据数据包中科学数据的个数，ulf波形数据取n＝82，elf波形数据取n＝820；
[0241]
x(k)为dft计算获得的频域内数据。
[0242]
(2)vlf波形数据的快速傅里叶变换(fft)
[0243]
首先使用fft将信号从时域变换到频域。fft是一种线性积分变换，其原理是将n点x数据分解为两个n/2点的x1(n)和x2(n)数据，对x1(n)和x2(n) 进行n/2点的dft傅里叶变换。具体公式见式(4
‑
7)和式(4
‑
8)。
[0244][0245][0246]
其中，x(n)为时域离散的信号；n为连续数据的个数，vlf波形数据取 n＝4096；x1(k)、x2(k)分别为x1(n)和x2(n)的n/2点离散傅里叶转换，见式 (4
‑
6)；x1(n)、x2(n)为x(n)按奇偶拆分获得的数组。
[0247]
4.1.3正交性校正
[0248]
因加工、装配过程引起的误差，感应式磁力仪传感器三分量的磁轴、几何轴不可能完全正交。需要通过磁正交性标定试验和几何形位测量试验获得传感器磁正交坐标系和结构正交坐标系之间的关系，对测量的空间磁场变化数据进行正交性校正。因此，正交性校正过程需要的辅助数据为正交校正矩阵m
orth
，其正交性校正公式为：
[0249][0250]
其中，为三分量频谱幅度数据；为三分量正交校正后幅度数据； m
orth
是正交性校正矩阵。
[0251]
4.1.4传递函数校正
[0252]
在频域内将电压量的科学数据，通过磁传感器的电磁转换因子逐条谱线进行幅度校正和物理量转换，生成校正后的三分量磁场物理量数据。此过程中需要的辅助数据为传感器包膜内温度数据和带有温度信息的三分量电压
‑
磁场转换参数(即磁传感器的传递函数)：
[0253]
●
传感器监控温度数据；
[0254]
●
ulf电压
‑
磁场转换频率响应参数；
[0255]
●
elf电压
‑
磁场转换频率响应参数；
[0256]
●
vlf电压
‑
磁场转换频率响应参数；
[0257]
其中，磁传感器温度数据从感应式磁力仪自采集温度工程参数中获取；三分量传电压
‑
磁场转换频率响应参数在实验室不同温度下条件下磁传感器的频带响应实验获取。根据磁传感器实时温度监控数据，选取对应温度区间内的电压
‑
磁场转换因子数据，将谱线相应频点的电压量转换为磁场b1，转换公式如下：
[0258]
[0259]
4.1.5傅里叶逆变换
[0260]
科学数据在频域内完成幅度校正后利用傅里叶逆变换将波形数据转换到时域中生成波形物理量数据。
[0261]
(1)ulf、elf的离散傅里叶逆变换(idft)
[0262]
通过idft将信号从频域变换到时域：
[0263][0264]
其中x(k)为ulf和elf波形数据的频域幅度值；x(n)为idft转换后的波形数据；n为参与计算的数据个数，对于ulf，n＝82，对于elf，n＝820。
[0265]
(2)vlf快速傅里叶逆变换ifft
[0266]
通过ifft将信号从频域变换到时域，将n点傅里叶逆变换idft分解为 n/2个两点idft进行计算，其变换公式见式(4
‑
12)和式(4
‑
13)。
[0267][0268][0269]
其中，x1(n)、x2(n)为奇偶拆分获得的频域数据利用式(4
‑
8)进行idft 转换的时域数据；x(n)为时域的离散信号；n为连续数据的个数vlf波形数据取n＝4096。
[0270]
4.1.6功率谱密度计算
[0271]
将vlf频段的频谱数据计算为功率谱密度(psd)数据，功率谱密度计算采用以下公式：
[0272][0273]
其中，psd
ch1
(kδf)、psd
ch2
(kδf)、psd
ch3
(kδf)为vlf各分量的功率谱密度，单位为nt
·
hz
‑
1/2；k为离散频域数据序号，k＝0,1,2，
···
n/2；fvlf为vlf 频段的采样频率，常量，fvlf＝51.2khz。n为vlf频段的fft计算样本量，取值为4096；δf为vlf频谱数据的频率间隔，δf＝f
band3
/n＝12.5，单位为hz； b1
x
(k
·
δf)、b1
y
(k
·
δf)、b1
z
(k
·
δf)为vlf各分量频谱数据的幅度值，单位为nt。
[0274]
4.1.7在轨定标数据提取
[0275]
在定标工作模式下，感应式磁力仪自身具有产生625hz和10khz定标信号的功能，定标信号将磁通负反馈线圈作为标定线圈，在磁传感器内部产生标定磁场信号，通过对磁传感器对此标定磁场感应电动势的采集。根据感应式磁力仪的工作状态，提取在轨定标模式下的数据，并在在频域范围内提取vlf频带内625hz和10khz两频点的电压量数据与地面
测得的在轨定标信号校准参数进行误差分析，判断感应式磁力仪的工作状态。
[0276]
4.2 1级数据生成2级数据
[0277]
将星务数据、磁情指数等信息与每轨数据时空匹配，对1级数据进行传感器到卫星坐标系和星地坐标系的转换，并标注地理和地磁坐标和姿轨控数据。处理流程见图11。
[0278]
4.2.1传感器坐标系转换到卫星坐标系
[0279]
利用传感器几何坐标系探测坐标系(传感器坐标系)和卫星结构基准立方镜坐标系(卫星坐标系)之间的关系，将三分量磁场数据转换到卫星坐标系下。需要的辅助参数有传感器到卫星坐标系转换矩阵m1，其为3*3的坐标转换矩阵，在传感器装配和卫星总装过程中测量获得，由载荷和卫星厂家提供；
[0280]
假设卫星平台坐标系的三个轴分别为x，y和z，感应式磁力仪传感器的三个轴分别为x1，y1和z1(传感器的x1轴与平台x轴平行)，传感器输出信号分别为b
1x
，b
1y
和b
1z
，则观测数据由传感器坐标系(探测坐标系)转换到卫星坐标系的转换关系为：
[0281][0282]
其中，b
2x
，b
2y
和b
2z
为卫星坐标系下三分量磁场数据，单位，nt；m1在卫星总体装配后由载荷厂家提供。
[0283]
4.2.2星地坐标转换
[0284]
由于卫星姿态的变化，设卫星坐标系相对于地理坐标系(地心地固系)之间的三个欧拉角度是φ、ψ、θ。由卫星坐标系sbc到地理坐标系geo之间的转换矩阵设为m2，则m2可表示为：
[0285][0286]
设卫星坐标系下的三分量磁场数据用(b
2x
，b
2y
和b
2z
)表示，地理坐标系中的交流磁场用(b
3x
，b
3y
和b
3z
)表示，则有：
[0287][0288]
设地理坐标系geo到地磁坐标系mag之间的转换关系为：
[0289][0290]
其中，θ、分别为偶极轴的地理余纬度、经度。
[0291]
4.3 2级数据生成3级数据
[0292]
根据感应式磁力仪的2数据产品，按照以下步骤得到3级数据产品及处理报告，流
程如图12所示。其中背景值的求取方法，以30天重访轨道数据的中值为例进行描述。
[0293]
第一步：搜索重访轨道。
[0294]
搜索重访轨道的具体算法：
[0295]
根据cses卫星轨道重访周期，利用当前轨道号，可获取前30天相应的重访轨道数据，具体公式为：
[0296]
s＝l
‑
n
×
m,n＝1,
…
,6
ꢀꢀ
(4
‑
19)
[0297]
其中，为重访轨道号；为当前轨道号；为cses重访轨道固定差。
[0298]
第二步：划分网格。
[0299]
针对半轨道的2a级波形和功率谱数据，按纬度间隔0.05
°
(中国区域：) 和0.1
°
(其它区域)进行划分网格。
[0300]
第三步：时序分析，求取每个网格里的中值。
[0301]
中值的计算方法采用下列算法描述。假设每个网格内有从小到大排列的n 个磁场数据b1,
…
,b
n
。当n 1为4的整数倍数时，中值b
m
及第一、第三四分位数计算公式为式：
[0302][0303]
当不为4的整数倍数时，通过上式计算出来的四分位数位置带有小数，这时，对该小数四舍五入取整。
[0304]
求取背景值的方法：对落在每个网格内的数据，即当前轨道的前30天重访轨道的波形和(或)频谱数据，使用式(4
‑
20)求取中值及四分位值。将此中值定义为当前轨道所在网格的背景值。
[0305]
重采样的方法：对落在每个网格内的当前轨道的数据(波形和(或)频谱数据)，使用公式(4
‑
20)求取中值，并将其作为此网格区域内的重采样值。
[0306]
用于判定异常的阈值的计算方法，见式(4
‑
21)：使用和，计算四分位差 iqr(inter
‑
quartile range，式(4
‑
22))。
[0307]
b
c
＝b
m
±
iqr
ꢀꢀ
(4
‑
21)
[0308][0309]
第四步：绘图与产出数据处理报告。
[0310]
(1)将当前5天的重访轨道数据(分升轨和降轨)，以及由前30天的重访轨道数据计算得到的b
m
、b
c
序列，即重采样的时间序列，区分三分量和三频段绘制成图。同时，在图标注当前轨道时间的磁情指数。
[0311]
(2)绘出每条轨道正负5个单位的经度区域范围的前30内天发生的震级大于5级、6级及7级以上的地震事件序列(具体图标的选取地震事件原则，参考demeter卫星快视图中对地震事件序列的绘图方式)。
[0312]
(3)绘出当前轨道所在全球位置，以及标注前7内天发生的震级大于6级的地震事
件。
[0313]
(4)生成数据处理报告。
[0314]
4.4 2级数据生成4级数据
[0315]
根据感应式磁力仪的2数据产品，按照以下步骤得到4级数据产品及处理报告，流程如图13所示。其中背景值的求取方法，以30天重访轨道数据的中值为例进行描述。
[0316]
图13中的绿色网络区域示意的是全球区域的每个网格大小：2.55。蓝色线表示卫星降轨轨迹，红色表示卫星升轨轨迹。在5天的回归周期内，迟到上两条相邻升(降)轨轨迹间距约为530km。
[0317]
第一步：划分网格。
[0318]
按照2.5
°
(纬度)
×5°
(经度)的空间间隔将全球划分为若干网格。
[0319]
第二步：空间分布分析。
[0320]
(1)求背景值
[0321]
按照30天时间间隔及5
°
(经度)
×
2.5
°
(纬度)空间间隔，滑动选取全球和中国区域(0
‑
60
°
纬度；60
‑
140
°
经度)上空的磁场2级数据，计算每个空间间隔内所有轨道的中值、四分位点()及分位差，以及异常判别值；并通过对各个间隔内中值进行插值，得到全球及全国范围空间分布背景场；
[0322]
(2)当前5天变化磁场数据
[0323]
对于当前5天变化磁场数据，参考式(4
‑
20)计算每一个网格中变化磁场的中值其及四分位差。
[0324]
以当前5天时间间隔及5
°
(经度)
×
2.5
°
(纬度)空间间隔，计算观测得到磁场在每个空间间隔内所有轨道的中值；并通过对各个间隔内中值进行插值，得到全球及全国的空间分布图，滑动更新每天的全球及全国范围空间分布图；
[0325]
(3)当前5天与背景场的差值，并归一化
[0326]
计算每天每个间隔内中值与之前30天的背景中值的差值，并使用每个格子内的最大值进行归一化，并通过对各个间隔归一化后的值进行插值，得到全球及全国的每天动态变化空间分布图。
[0327]
第三步：加入地震事件和磁情指数，生成4级数据处理报告。
[0328]
在第二步绘图结果的基础上，标注30天内大于7级地震事件记录。将以上处理过程中的数据错误等异常，写入数据处理报告中，生成4级数据产品处理报告。
[0329]
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。