机器学习中的地磁场

2025年3月17日 | 阅读 8 分钟

数千年来，科学家们一直对地球的地磁场着迷。它对各种地质和环境过程有着重大影响，从导航系统到气候模式都受到影响。随着机器学习技术的突破，研究人员现在可以以新的方式探索和理解地磁场的复杂性。

当地球磁场有效转移来自太阳风的能量时，就会发生地磁暴。随之而来的磁场变化会导致磁导航的精度增加。扰动-地暴时间指数（Dst）衡量地磁暴的强度。

Dst指数是磁层动力学的一个关键参数，它驱动着地磁扰动模型，如NOAA/NCEI的高清晰度地磁模型-实时（HDGM-RT）。此外，磁力测量员、政府组织、学术机构、卫星运营商和电网运营商也利用Dst指数来确定地磁暴的严重程度和持续时间。

自1975年以来，人们就提出了经验模型，仅根据拉格朗日（L1）点上的太阳风数据来预测Dst，这些数据来自NOAA的深空气候观测站（DSCOVR）和NASA的先进成分探测器（ACE）等卫星。在过去的三十年里，已经提出了各种预测太阳风Dst的模型，包括经验模型、基于物理的模型和机器学习技术。虽然机器学习模型优于传统技术，但总有改进的空间，尤其是在预测严重事件方面。更重要的是，我们想要能够处理原始的实时数据流并且不受传感器故障和噪声影响的系统。

机器学习在地磁场研究中的应用

地磁数据分析：机器学习技术用于分析来自地面观测站、卫星和其他来源收集的海量地磁数据。这些算法可以识别数据中的模式、异常和趋势，从而更好地理解地磁事件，包括磁暴、电离层扰动和地磁逆转。
预测建模：机器学习技术，如神经网络和支持向量机，被用于构建预测模型来预测地磁活动。这些模型可以提供空间天气事件的早期预警，这些事件可能干扰卫星通信、电网和导航系统，从而减少潜在的危险。
数据融合与集成：机器学习能够集成来自多个来源的异构数据集，例如地面测量、卫星观测和数值模拟。通过合并这些数据集，研究人员可以更全面地了解地磁场的动力学及其与电离层和磁层等其他地球系统的相互作用。
探索磁异常：机器学习方法有助于发现和表征磁异常，这是由于矿产资源、断层线和构造板块边界等地质特征引起的，与预测的地磁场行为的偏差。理解这些异常对于地质研究、资源测绘和地震监测至关重要。

机器学习在地磁场研究中面临的挑战

机器学习极大地增进了我们对地球地磁场的认识，使得预测更加准确，数据解释更加完善，可视化能力更强。然而，仍然存在重大的挑战，包括需要高质量的数据，开发能够处理复杂时空相互作用的强大算法，以及将机器学习模型整合到现有的地球物理框架中。

代码

我们现在的目标是创建能够突破预测性能界限的Dst预测模型，同时保持可操作的实用性，利用来自NOAA的DSCOVR和NASA的ACE卫星的实时太阳风（RTSW）数据流。改进的模型可以提供地磁暴的早期预警，并最大限度地减少磁导航系统中的不准确性。

TPU

我们将使用TensorFlow初始化一个TPU。

import tensorflow as tf
# detect and init the TPU
try:
    tpu = tf.distribute.cluster_resolver.TPUClusterResolver() 
    print('Running on TPU ', tpu.master())
except ValueError:
    tpu = None

if tpu:
    tf.config.experimental_connect_to_cluster(tpu)
    tf.tpu.experimental.initialize_tpu_system(tpu)
    tpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)
else:
    tpu_strategy = tf.distribute.get_strategy() 

print("REPLICAS: ", tpu_strategy.num_replicas_in_sync)

输出

Geomagnetic Field Using Machine Learning

导入库

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from tensorflow import keras 
from tensorflow.keras.layers import Dense, Conv2D, BatchNormalization, Activation 
from tensorflow.keras.layers import AveragePooling2D, Input, Flatten 
from tensorflow.keras.optimizers import Adam 
from tensorflow.keras.callbacks import ModelCheckpoint, LearningRateScheduler 
from tensorflow.keras.callbacks import ReduceLROnPlateau 
from tensorflow.keras.preprocessing.image import ImageDataGenerator 
from tensorflow.keras.regularizers import l2 
from tensorflow.keras import backend as K 
from tensorflow.keras.models import Model 

%matplotlib inline

读取数据集

DATA_PATH = Path("../input/magnet-nasa")

dst = pd.read_csv(DATA_PATH / "dst_labels.csv")
dst.timedelta = pd.to_timedelta(dst.timedelta)
dst.set_index(["period", "timedelta"], inplace=True)

sunspots = pd.read_csv(DATA_PATH / "sunspots.csv")
sunspots.timedelta = pd.to_timedelta(sunspots.timedelta)
sunspots.set_index(["period", "timedelta"], inplace=True)

solar_wind = pd.read_csv(DATA_PATH / "solar_wind.csv")
solar_wind.timedelta = pd.to_timedelta(solar_wind.timedelta)
solar_wind.set_index(["period", "timedelta"], inplace=True)

EDA

现在，我们来看看我们的数据集。

# first 5 rows of Dst
print("Dst shape: ", dst.shape)
dst.head()

输出

# grouping on the basis of period
dst.groupby("period").describe()

输出

# Shape of solar wind dataset
print("Solar wind shape: ", solar_wind.shape)
solar_wind.head()

输出

#  Shape of Sunspot dataset
print("Sunspot shape: ", sunspots.shape)
sunspots.head()

输出

# Grouping of solar_wind dataset on the basis of period
solar_wind.groupby("period").describe().T

输出

# Grouping of sunspots dataset on the basis of period
sunspots.groupby("period").describe().T

输出

数据集可视化

def show_raw_visualization(data):
    fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 15), dpi=80)
    for i, key in enumerate(data.columns):
        t_data = data[key]
        ax = t_data.plot(
            ax=axes[i // 2, i % 2],
            title=f"{key.capitalize()}",
            rot=25,
        )

    fig.subplots_adjust(hspace=0.8)
    plt.tight_layout()


cols_to_plot = ["bx_gse", "bx_gsm", "bt", "density", "speed", "temperature"]
show_raw_visualization(solar_wind[cols_to_plot].iloc[:1000])

输出

我们将为solar_wind、sunspots和dst中的变量生成一个相关性热力图。

corr = solar_wind.join(sunspots).join(dst).fillna(method="ffill").corr()

plt.figure(figsize=(10, 5))
plt.matshow(corr, fignum=1)
plt.xticks(range(corr.shape[1]), corr.columns, fontsize=14, rotation=90)
plt.gca().xaxis.tick_bottom()
plt.yticks(range(corr.shape[1]), corr.columns, fontsize=14)


cb = plt.colorbar()
cb.ax.tick_params(labelsize=14)
plt.title("Feature Correlation Heatmap", fontsize=16)
plt.show()

输出

from numpy.random import seed
from tensorflow.random import set_seed

seed(2020)
set_seed(2021)

标准缩放器

我们将创建一个预处理管道，以确保在用于建模之前，数据得到妥善的清理、聚合和缩放。

from sklearn.preprocessing import StandardScaler, MinMaxScaler

# subset of solar wind features to use for modeling
SOLAR_WIND_FEATURES = [
    "bt",
    "temperature",
    "bx_gse",
    "by_gse",
    "bz_gse",
    "speed",
    "density",
]

# All of the features we'll use, including sunspot numbers
XCOLS = (
    [col + "_mean" for col in SOLAR_WIND_FEATURES]
    + [col + "_std" for col in SOLAR_WIND_FEATURES]
    + ["smoothed_ssn"]
)


def impute_features(feature_df):
    """Imputes data using the following methods:
    - `smoothed_ssn`: forward fill
    - `solar_wind`: interpolation
    """
    # forward fill sunspot data for the rest of the month
    feature_df.smoothed_ssn = feature_df.smoothed_ssn.fillna(method="ffill")
    # interpolate between missing solar wind values
    feature_df = feature_df.interpolate()
    return feature_df


def aggregate_hourly(feature_df, aggs=["mean", "std"]):
    """Aggregates features to the floor of each hour using mean and standard deviation.
    e.g. All values from "11:00:00" to "11:59:00" will be aggregated to "11:00:00".
    """
    # Group by the floor of each hour using timedelta index
    agged = feature_df.groupby(
        ["period", feature_df.index.get_level_values(1).floor("H")]
    ).agg(aggs)
    # flatten hierarchical column index
    agged.columns = ["_".join(x) for x in agged.columns]
    return agged


def preprocess_features(solar_wind, sunspots, scaler=None, subset=None):
    """
    Preprocessing steps:
        - Subset the data
        - Aggregate hourly
        - Join solar wind and sunspot data
        - Scale using standard scaler
        - Impute missing values
    """
    # select features we want to use
    if subset:
        solar_wind = solar_wind[subset]

    # Aggregate solar wind data and join with sunspots
    hourly_features = aggregate_hourly(solar_wind).join(sunspots)

    # subtract mean and divide by standard deviation
    if scaler is None:
        scaler = StandardScaler()
        scaler.fit(hourly_features)

    normalized = pd.DataFrame(
        scaler.transform(hourly_features),
        index=hourly_features.index,
        columns=hourly_features.columns,
    )

    # impute missing values
    imputed = impute_features(normalized)

    # We want to return the scalar object as well to use later during the prediction
    return imputed, scaler

features, scaler = preprocess_features(solar_wind, sunspots, subset=SOLAR_WIND_FEATURES)
print(features.shape)
features.head()

输出

YCOLS = ["t0", "t1"]


def process_labels(dst):
    y = dst.copy()
    y["t1"] = y.groupby("period").dst.shift(-1)
    y.columns = YCOLS
    return y


labels = process_labels(dst)
labels.head()

输出

data = labels.join(features)
data.head()

输出

分割数据集

现在我们将数据集分割成训练集、测试集和验证集。

def get_train_test_val(data, test_per_period, val_per_period):
    """Splits data across periods into train, test, and validation"""
    # assign the last `test_per_period` rows from each period to test
    test = data.groupby("period").tail(test_per_period)
    interim = data[~data.index.isin(test.index)]
    # assign the last `val_per_period` from the remaining rows to validation
    val = interim.groupby("period").tail(val_per_period)
    # The remaining rows are assigned to train
    train = interim[~interim.index.isin(val.index)]
    return train, test, val

train, test, val = get_train_test_val(data, test_per_period=6_000, val_per_period=3_000)

ind = [0, 1, 2]
names = ["train_a", "train_b", "train_c"]
width = 0.75
train_cnts = [len(df) for _, df in train.groupby("period")]
val_cnts = [len(df) for _, df in val.groupby("period")]
test_cnts = [len(df) for _, df in test.groupby("period")]

p1 = plt.barh(ind, train_cnts, width)
p2 = plt.barh(ind, val_cnts, width, left=train_cnts)
p3 = plt.barh(ind, test_cnts, width, left=np.add(val_cnts, train_cnts).tolist())

plt.yticks(ind, names)
plt.ylabel("Period")
plt.xlabel("Hourly Timesteps")
plt.title("Train/Validation/Test Splits over Periods", fontsize=16)
plt.legend(["Train", "Validation", "Test"])

输出

print(train.shape)
train.head()

输出

print(test.shape)
test.head()

输出

print(val.shape)
val.head()

输出

批次大小和时间步长

from keras import preprocessing


data_config = {
    "timesteps": 34,
    "batch_size": 34,
}


def timeseries_dataset_from_df(df, batch_size):
    dataset = None
    timesteps = data_config["timesteps"]

    # iterate through periods
    for _, period_df in df.groupby("period"):
        # realign features and labels so that the first sequence of 32 is aligned with the 33rd target
        inputs = period_df[XCOLS][:-timesteps]
        outputs = period_df[YCOLS][timesteps:]

        period_ds = preprocessing.timeseries_dataset_from_array(
            inputs,
            outputs,
            timesteps,
            batch_size=batch_size,
        )

        if dataset is None:
            dataset = period_ds
        else:
            dataset = dataset.concatenate(period_ds)

    return dataset


train_ds = timeseries_dataset_from_df(train, data_config["batch_size"])
val_ds = timeseries_dataset_from_df(val, data_config["batch_size"])

print(f"Number of train batches: {len(train_ds)}")
print(f"Number of val batches: {len(val_ds)}")

输出

模型

我们现在将创建一个包含多个旨在处理序列数据的层的模型架构。

from tensorflow.keras.layers import Dense, LSTM,Bidirectional,Dropout,Embedding
from tensorflow.keras.models import Sequential
model_config = {"n_epochs": 15, "n_neurons": 256, "dropout": 0.1, "stateful": False}
# define our model
def creat_model():
    model = Sequential()
    model.add(Input(shape=(data_config["timesteps"],
                                 len(XCOLS))))
    model.add(tf.keras.layers.RNN(tf.keras.layers.LSTMCell(256,dropout=0.3), return_sequences=True,
                                         name="lstm1"))
#     model.add(tf.keras.layers.RNN(tf.keras.layers.LSTMCell(256), return_sequences=True,dropout=0.1,
#                                          name="lstm1.5"))
    model.add(tf.keras.layers.RNN(tf.keras.layers.LSTMCell(256), name="lstm2",))
    model.add(tf.keras.layers.Dense(units=2))

    return model

with tpu_strategy.scope():
        model=creat_model()
        model.compile(
        loss=tf.keras.losses.MeanSquaredError(),
        optimizer=tf.keras.optimizers.Adam(learning_rate=0.0005),
         )

model.summary()

输出

训练模型

我们现在将使用训练数据集（train_ds）训练模型，并在验证数据集（val_ds）上评估其性能。

history = model.fit(
        train_ds,
        batch_size=data_config["batch_size"],
        epochs=model_config["n_epochs"],
        verbose=1,
        shuffle=False,
        validation_data=val_ds,
    )

输出

求值

现在，我们将评估模型。

for name, values in history.history.items():
    plt.plot(values)

输出

测试模型

# Let us check the accuracy of the model 
test_ds = timeseries_dataset_from_df(test, data_config["batch_size"])
mse = model.evaluate(test_ds)
print(f"Test RMSE: {mse**.5:.2f}")

输出

在评估过程中计算出的192.2285的损失值是模型预测值与测试数据集上实际目标值之间的平均差值。虽然这个数字对于评估模型在训练期间的性能很有用，但必须根据数据集的大小和复杂性来解释。

13.86的RMSE表明，模型的预测值与实际值平均相差约13.86个单位。虽然这个数字提供了模型性能的信息，但关键是要考虑问题域的背景和具体标准，以评估模型的性能是否足够。

随着机器学习的进步，未来的地磁场研究将受益于深度学习、强化学习和可解释人工智能方法的进步。此外，地球物理学家、数据科学家和领域专家之间的多学科合作将是处理复杂的地磁现象和在实际应用中使用机器学习驱动的见解所必需的。

下一话题使用机器学习进行图像生成

机器学习中的地磁场

机器学习在地磁场研究中的应用

机器学习在地磁场研究中面临的挑战

TPU

导入库

读取数据集

EDA

数据集可视化

标准缩放器

分割数据集

批次大小和时间步长

模型

训练模型

求值

测试模型

联系信息

关注我们

教程

面试题

在线编译器

Python

Java

.Net Framework

AI, ML and Data Science

Cloud Technology

B.Tech and MCA

Web Technology

PHP

Software Testing

Technical Interview

Java Interview

Python

Web Interview

Database Interview

B.Tech / MCA

Important Interview

Software Testing Interview

Company Interviews

Online Compilers

Multiple Choice Questions

机器学习

监督式学习

分类

杂项

相关教程

面试题

机器学习中的地磁场

机器学习在地磁场研究中的应用

机器学习在地磁场研究中面临的挑战

TPU

导入库

读取数据集

EDA

数据集可视化

标准缩放器

分割数据集

批次大小和时间步长

模型

训练模型

求值

测试模型

相关帖子

机器学习中的 BLEU 分数

ML | 主动学习

机器学习在交易中的应用

什么是人机协作机器学习

5G 技术对物联网有哪些好处？

使用神经网络进行讽刺检测

机器学习中的 IPL 预测

使用梯度下降进行线性回归

机器学习中的感知器

VGGNet-16 架构

订阅 Tpoint Tech

联系信息

关注我们

教程

面试题

在线编译器