2003 年 10 月末、太陽の AR10486 という活動領域から X28 クラスの巨大フレアが連続的に発生しました。地球を直撃したコロナ質量放出（CME）は GPS の測位誤差を最大 30 m まで悪化させ、スウェーデンのマルメで広域停電を引き起こし、軌道上の太陽電池パネルの劣化を平年の数倍に加速させました。このような「ハロウィン・ストーム」と呼ばれる極端宇宙天気イベントは、現代の高度に電化された社会にとって他人事ではありません。
太陽フレアの発生を 24 時間前に予測 できれば、衛星運用者はセーフホールドへの遷移を準備でき、電力会社は予備の変圧器を待機させ、航空会社は極域経由フライトを迂回させることができます。問題は、太陽の磁場が複雑な非線形磁気流体力学（MHD）系であり、物理シミュレーションだけで爆発の発生時刻を当てるのは現在の計算機資源では非現実的だという点です。そこで近年、機械学習で SDO/HMI が観測する磁気図から経験的にフレアの発生確率を学習する アプローチが台頭してきました。
2024 年に NASA と IBM が発表した太陽基盤モデル「Surya」 は、約 9 年分のフルディスク多波長 SDO データを Transformer で事前学習し、フレア発生やフラックス値の予測でこれまでの最先端手法を上回る性能を示しました。本記事では、太陽フレア予測タスクの定式化から評価指標、伝統的手法、最新の Surya までの系譜を整理し、PyTorch で合成磁気図を題材にしたミニ実験を行います。
太陽フレア予測を理解すると、以下のような応用が開けます。
- 衛星運用: 静止気象衛星や測位衛星の SEU（シングルイベントアップセット）対策、姿勢制御モードの事前切替
- 電力網保護: GIC（地磁気誘導電流）による変圧器損傷の予防保全と需給調整
- 航空・通信: 極域 HF 通信の遮断回避、極軌道有人ミッション（ISS、Gateway）の被ばく管理
- 宇宙基盤モデル: 多波長太陽画像で事前学習した Transformer の下流タスク適応という、地球観測基盤モデルと共通する設計パターンの実例
本記事の内容
- 太陽フレア・CME が社会インフラに与える具体的な影響
- SDO/HMI 磁気図と SHARP パラメータによる入力表現
- フレア発生予測の二値分類としての定式化
- クラス不均衡問題と TSS/HSS という宇宙天気特有の評価指標
- 手法の系譜: SVM/RF → CNN → Transformer/NASA Surya
- PyTorch で合成磁気図を学習する小型 CNN と focal loss の効果検証
- ROC/PR 曲線・混同行列・TSS の可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- CNN の基礎 — 畳み込み・プーリング・受容野の数理 — 畳み込み演算の数学的定義と特徴抽出の階層構造
- 二値分類とロジスティック回帰 — 確率出力と交差エントロピー損失
- Attention 機構の基礎 — Transformer 系モデルの構造的理解
太陽フレアとは何か — 物理的な背景
太陽活動領域の磁気エネルギー解放
太陽フレアは、太陽表面の 活動領域（active region, AR）にねじれて蓄えられた磁気エネルギーが、磁気再結合によって突発的に解放される現象です。エネルギーは $10^{19}$〜$10^{25}$ J におよび、X 線・紫外線・粒子線として放出されます。X 線のピーク強度に応じて A, B, C, M, X クラスに区分され、各クラスは桁違いの強度差を持ちます。
X クラス（$\geq 10^{-4}$ W/m$^2$）と M クラス（$10^{-5}$〜$10^{-4}$ W/m$^2$）は社会インフラに実害をもたらす規模であり、機械学習による予測対象もこの「M/X クラスが今後 24 時間以内に発生するか」という二値判定が標準的に採用されています。フレアそのものよりも、しばしば同時に発生する コロナ質量放出（Coronal Mass Ejection, CME）の方が地球へのインパクトは大きく、数億トンのプラズマが秒速数百〜数千 km で太陽風中を伝播し、1〜3 日後に地球磁気圏に到達します。
CME と地磁気誘導電流が社会に与える影響
CME が地球磁気圏に衝突すると、地磁気の急変によって長距離の電線に 地磁気誘導電流（Geomagnetically Induced Current, GIC）が流れます。1989 年の Quebec 停電はその古典例で、約 600 万人が 9 時間にわたって停電に見舞われました。電力会社にとって CME の事前警報は、変圧器の中性点接地リアクタンスを切り替えたり、需給予測を保守的に再計算したりするための重要な意思決定情報です。
衛星にとっては高エネルギー粒子線（SEP, Solar Energetic Particle）が脅威となります。SEP は半導体メモリのビット反転（SEU）や太陽電池セルの劣化、CCD/CMOS センサのノイズ増加を引き起こします。極軌道衛星や静止軌道衛星では、フレア発生時にミッション計算機をセーフモードへ遷移し、機密データの欠損を防ぐ運用が行われます。
社会的インパクトを把握したところで、次は「機械学習で何を入力として何を出力するのか」という問題設定を見ていきましょう。
フレア予測タスクの定式化
入力データ: SDO/HMI 磁気図と SHARP パラメータ
太陽フレア予測の入力として最も広く用いられているのが、NASA の Solar Dynamics Observatory（SDO） に搭載された Helioseismic and Magnetic Imager（HMI） の観測データです。HMI は太陽光球面の視線方向磁場 $B_{\text{LOS}}$ を 720 秒ごとに $4096 \times 4096$ ピクセルの解像度で取得し、また分解された 3 成分ベクトル磁場 $(B_x, B_y, B_z)$ も提供します。
活動領域ごとに切り出した磁気図パッチを直接 CNN に入力する方法もあれば、より物理的に意味のある SHARP パラメータ（Space-weather HMI Active Region Patch）という 16〜25 個の手作り特徴量を用いる方法もあります。代表的な SHARP パラメータには次のようなものがあります。
USFLUX— 活動領域全体の総絶対磁束 $\sum |B_z| \cdot dA$。エネルギーの貯蔵総量に相当TOTUSJH— 総非ポテンシャル電流密度 $\sum |J_z|$。磁場のねじれの強さを表すR_VALUE— 強磁場反転線（PIL）周辺の磁束の総和。フレアの最有力前兆指標とされるMEANSHR— 平均剪断角。観測磁場とポテンシャル磁場のずれを表し、ねじれによる磁気エネルギー蓄積の指標
これらは物理量に対して経験的に有効であることが知られており、ランダムフォレストやロジスティック回帰でも比較的高い性能を出せるため、ベースライン手法として標準的に使われます。
出力ラベルと予測ホライズン
予測タスクは典型的に「観測時刻 $t$ から $\Delta t$ 時間以内に M/X クラス・フレアが発生するか」という二値分類として定式化されます。$\Delta t$ は 24 時間が最も一般的で、$\Delta t = 6$ 時間や $\Delta t = 48$ 時間の設定もあります。
入力 $\bm{x}_t$ に対する条件付き確率を
$$ p_t = P(\text{M または X が } [t, t+\Delta t] \text{ に発生} \mid \bm{x}_t) $$
として推定し、閾値 $\theta$ を超えた場合に「発生する」と予測します。$\theta$ は ROC 曲線や運用要件から決定されます。
ここで重要なのは「24 時間以内に M/X クラス・フレアが発生する活動領域」は全体のごく一部だということです。次は、この極端なクラス不均衡が学習と評価をどう難しくするかを見ていきます。
クラス不均衡と宇宙天気特有の評価指標
なぜ Accuracy では駄目なのか
公開データセット SWAN-SF や GOES X 線データから作るベンチマークでは、正例（M/X 発生）の割合は わずか 3〜7% にすぎません。仮にモデルが「常に M/X は発生しない」と予測しても、Accuracy は 93〜97% に達してしまいます。これでは「全く役に立たないモデル」と「真に予測能力のあるモデル」を区別できません。
不均衡データを正しく扱うには、混同行列（confusion matrix）の各成分から指標を構築する必要があります。M/X 発生を正例、非発生を負例として、
| 予測: 発生 | 予測: 非発生 | |
|---|---|---|
| 実際: 発生 | TP（True Positive, ヒット） | FN（False Negative, 見逃し） |
| 実際: 非発生 | FP（False Positive, 誤報） | TN（True Negative, 正棄却） |
を定義します。
TSS（True Skill Statistic）
宇宙天気コミュニティで標準的に用いられる評価指標が True Skill Statistic（TSS） であり、Hanssen-Kuipers Discriminant とも呼ばれます。
$$ \text{TSS} = \text{TPR} – \text{FPR} = \frac{TP}{TP + FN} – \frac{FP}{FP + TN} $$
TSS は「ヒット率（再現率）」から「誤報率」を引いたものです。範囲は $-1$〜$+1$ で、$+1$ が完璧な予測、$0$ がランダム予測（または「常に発生しない」と予測する自明モデル）、$-1$ が逆予測です。
TSS が好まれる最大の理由は クラス不均衡比に依存しない ことです。Accuracy や F1 スコアは正例率が変わるとスコア自体が大きく変動しますが、TSS は TPR と FPR という比率の差なので、データ分布が変わっても比較可能な絶対的なスキル評価が得られます。
HSS（Heidke Skill Score）
もう一つの代表的指標が Heidke Skill Score（HSS） です。これは「ランダム予測に対するスキルの相対値」を測ります。
$$ \text{HSS} = \frac{2 (TP \cdot TN – FP \cdot FN)}{(TP + FN)(FN + TN) + (TP + FP)(FP + TN)} $$
HSS は不均衡データでは値が小さくなりがちですが、誤報の重みを TSS より相対的に大きく評価する性質があり、運用上「狼少年」を避けたいユースケースで重視されます。NOAA の宇宙天気予報センター（SWPC）の評価レポートでは TSS と HSS が併記されるのが慣例です。
学習側のクラス不均衡対策
評価だけでなく学習時にも不均衡対策が必要です。代表的な手段は以下の通りです。
- クラス重み付き損失: 標準の二値交差エントロピー $L_{\text{BCE}} = -[y \log p + (1-y) \log(1-p)]$ に、正例の重み $w_+$ を乗じる
- オーバー / アンダーサンプリング: 正例を複製、または負例を間引いてバッチ内比率を均す
- Focal Loss: 容易サンプルの寄与を抑える設計
Focal Loss は物体検出（RetinaNet, Lin et al., 2017）で導入された損失で、
$$ L_{\text{focal}}(p_t) = -\alpha_t (1 – p_t)^{\gamma} \log p_t, \quad p_t = \begin{cases} p & y = 1 \\ 1 – p & y = 0 \end{cases} $$
と定義されます。$\gamma$（典型的に 2）は「予測確率が正解に近い容易サンプル（$p_t \to 1$）の寄与を $(1 – p_t)^{\gamma}$ で減衰させる」という働きをし、難しいサンプル（多くは少数クラス）の勾配寄与を相対的に大きくします。フレア予測では正例が少なく難しいケースが多いため、focal loss はクラス重みと並んで実務的に有効です。
評価と学習の不均衡対策が揃ったところで、いよいよ手法そのものの系譜を辿りましょう。
手法の系譜 — SVM/RF から NASA Surya まで
第 1 世代: SHARP パラメータ + 古典 ML（2010〜2017）
最初期の自動フレア予測モデルは、SHARP パラメータをそのまま特徴量として サポートベクトルマシン（SVM）やランダムフォレストに入力する形式でした。Bobra & Couvidat（2015）が SVM で 24 時間予測の TSS $\sim 0.76$ を達成し、SHARP の有効性を確立した記念碑的な研究です。
この世代のメリットは、特徴量の物理的解釈が明確で運用者の信頼を得やすいことと、計算量が軽くリアルタイム運用に乗せやすいことです。デメリットは、人手で設計した特徴量に依存し、活動領域の 空間構造 や 時間的進化 を完全には捉えきれない点でした。
第 2 世代: CNN による生磁気図学習（2017〜2021）
CNN の台頭により、SHARP パラメータを介さず磁気図画像を直接入力する手法が現れました。Huang et al.（2018）は LeNet 風の CNN で TSS $\sim 0.65$、Park et al.（2020）はより深い VGG 系で TSS $\sim 0.78$ を報告しています。
CNN の魅力は 強磁場反転線（PIL）周辺のフィラメント構造 や デルタ型黒点配置 といった、人間が定義しにくい空間特徴を自動的に抽出できる点にあります。一方で活動領域ごとのデータ数が限られるため、過学習を避けるためのデータ拡張や転移学習が重要になります。
第 3 世代: 時系列 CNN・LSTM・Transformer（2020〜2023）
活動領域は時間とともに発展するため、単一時刻のスナップショットではなく 時系列としての進化 を入力にする方が予測精度が上がります。Liu et al.（2019）は LSTM で SHARP パラメータの時系列を扱い、Sun et al.（2022）は ConvLSTM や 3D-CNN で磁気図動画から特徴を抽出しました。
Transformer ベースの手法（Jiao et al., 2023 ほか）は、長期依存と注意機構による解釈性を兼ね備え、特に活動領域が太陽の縁から現れて中央子午線を通過して縁に沈むという数日スケールの発展を捉えるのに向いています。注意重みを可視化することで、モデルが PIL の活性化を「見ている」ことが確認できる事例も報告されています。
第 4 世代: 太陽基盤モデル Surya（NASA/IBM, 2024）
直近で最も注目されているのが、NASA Goddard と IBM Research が公開した Surya です。これは Hugging Face 上で重みも公開された、初の 太陽用基盤モデル（heliophysics foundation model）です。
Surya の特徴を整理すると次のようになります。
- 入力: SDO の 8 つの観測チャネル（HMI 磁気図 + AIA の 7 波長 EUV/UV 画像）を結合した多チャネル画像
- 事前学習: 約 9 年分（2010〜2019）のフルディスク画像で 自己教師あり に学習（マスク再構成 + 次フレーム予測のハイブリッド目的関数）
- アーキテクチャ: Spectral-gated Spatiotemporal Attention Transformer（時空間注意を周波数領域でゲーティング）
- 下流タスク: M/X フレア予測、F10.7 太陽電波フラックス回帰、活動領域セグメンテーション、CME 検出など多様に転移可能
Surya が示した最大のメッセージは、地球観測基盤モデル（Prithvi など）と同じ「大規模事前学習 + 軽量ファインチューニング」のパラダイムが、太陽データにも有効である ことです。これにより、各研究グループが個別に CNN を訓練していた時代から、共通の基盤モデルを fine-tune するエコシステムへの転換が始まっています。
手法の系譜が見えたところで、簡易版の CNN を PyTorch で実装して、不均衡データでの focal loss の効果を実感しましょう。
Python 実装 — 合成磁気図で焦点損失を体験する
実際の SDO/HMI データはサイズが大きく前処理パイプラインも煩雑なので、本節では 合成磁気図（フレア活動領域を模した極性混在パターンと、静穏な双極性パターンの 2 クラス）で学習する小型 CNN を組み立てます。focal loss と通常 BCE の比較で、不均衡対策の効果を可視化します。
合成データ生成
まず、$64 \times 64$ ピクセルの磁気図風画像を 2 クラス分布として生成します。「フレア活動領域」は強い正負磁場が複雑に混在し、PIL が長く曲がりくねっている設定です。「静穏領域」はシンプルな双極性配置にします。
import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.utils.data import DataLoader, TensorDataset
import matplotlib.pyplot as plt
from sklearn.metrics import (
roc_curve, precision_recall_curve, auc, confusion_matrix
)
np.random.seed(0)
torch.manual_seed(0)
IMG = 64
def _gauss_blob(cx, cy, amp, sigma, grid_x, grid_y):
return amp * np.exp(-((grid_x - cx) ** 2 + (grid_y - cy) ** 2) / (2 * sigma ** 2))
def make_quiet_region():
"""静穏領域: シンプルな双極性磁場"""
gx, gy = np.meshgrid(np.arange(IMG), np.arange(IMG))
cx1, cy1 = np.random.uniform(20, 30, 2)
cx2, cy2 = np.random.uniform(34, 44, 2)
sigma = np.random.uniform(5, 8)
img = _gauss_blob(cx1, cy1, 1.0, sigma, gx, gy) - _gauss_blob(cx2, cy2, 1.0, sigma, gx, gy)
img += 0.05 * np.random.randn(IMG, IMG)
return img
def make_flare_region():
"""フレア活動領域: 複雑な多極磁場と長い PIL"""
gx, gy = np.meshgrid(np.arange(IMG), np.arange(IMG))
img = np.zeros((IMG, IMG))
n_blobs = np.random.randint(4, 8)
for _ in range(n_blobs):
cx, cy = np.random.uniform(10, 54, 2)
amp = np.random.choice([-1, 1]) * np.random.uniform(0.8, 1.5)
sigma = np.random.uniform(3, 6)
img += _gauss_blob(cx, cy, amp, sigma, gx, gy)
img += 0.1 * np.random.randn(IMG, IMG)
return img
# 不均衡データ: 正例 5%, 負例 95% （実フレア予測に近い比率）
N_TOTAL, POS_RATIO = 2000, 0.05
n_pos = int(N_TOTAL * POS_RATIO)
n_neg = N_TOTAL - n_pos
X_pos = np.stack([make_flare_region() for _ in range(n_pos)])
X_neg = np.stack([make_quiet_region() for _ in range(n_neg)])
X = np.concatenate([X_pos, X_neg], axis=0)[:, None, :, :].astype(np.float32)
y = np.concatenate([np.ones(n_pos), np.zeros(n_neg)]).astype(np.float32)
# train/test 分割
perm = np.random.permutation(N_TOTAL)
X, y = X[perm], y[perm]
n_train = int(0.8 * N_TOTAL)
X_tr, X_te = X[:n_train], X[n_train:]
y_tr, y_te = y[:n_train], y[n_train:]
print(f"train: pos={int(y_tr.sum())}/{len(y_tr)} ({y_tr.mean():.1%})")
print(f"test : pos={int(y_te.sum())}/{len(y_te)} ({y_te.mean():.1%})")
このコードを実行すると、訓練データには約 80 件の正例と 1520 件の負例があり、正例率はおおむね 5% であることがわかります。実際の SWAN-SF データセットの不均衡率（3〜7%）に近い設定にしたので、ここから先の TSS スコアは現実のフレア予測ベンチマークと比較しやすい数値感になります。
小型 CNN モデル
次に、磁気図画像から二値確率を出力する小型 CNN を定義します。畳み込み 3 段 + 平均プーリング + 全結合 1 段のシンプルな構造ですが、合成データには十分です。
class FlareCNN(nn.Module):
def __init__(self):
super().__init__()
self.conv1 = nn.Conv2d(1, 16, kernel_size=3, padding=1)
self.conv2 = nn.Conv2d(16, 32, kernel_size=3, padding=1)
self.conv3 = nn.Conv2d(32, 64, kernel_size=3, padding=1)
self.pool = nn.MaxPool2d(2, 2)
self.gap = nn.AdaptiveAvgPool2d(1)
self.fc = nn.Linear(64, 1)
def forward(self, x):
x = self.pool(F.relu(self.conv1(x))) # 64 -> 32
x = self.pool(F.relu(self.conv2(x))) # 32 -> 16
x = self.pool(F.relu(self.conv3(x))) # 16 -> 8
x = self.gap(x).flatten(1)
return self.fc(x).squeeze(-1) # logit
最終層は logit（シグモイド前の値）を出力するようにしておきます。これにより BCEWithLogitsLoss が数値安定に使え、後述の focal loss も logit ベースで実装できます。
Focal Loss の実装
focal loss を logit からそのまま計算するクラスを書きます。BCEWithLogitsLoss の pos_weight 引数だけでも一定の効果がありますが、focal loss はそれに加えて「容易サンプルの寄与減衰」を加えることができる点が違います。
class FocalLoss(nn.Module):
def __init__(self, alpha=0.25, gamma=2.0):
super().__init__()
self.alpha = alpha
self.gamma = gamma
def forward(self, logits, targets):
# logits -> 確率
p = torch.sigmoid(logits)
# 正解側の確率 p_t
p_t = p * targets + (1 - p) * (1 - targets)
# クラス重み alpha_t
alpha_t = self.alpha * targets + (1 - self.alpha) * (1 - targets)
# focal weight (1 - p_t)^gamma で容易サンプルを減衰
loss = -alpha_t * (1 - p_t) ** self.gamma * torch.log(p_t.clamp(min=1e-8))
return loss.mean()
alpha=0.25, gamma=2.0 は元論文（Lin et al., 2017）の推奨値です。$\gamma$ を大きくするほど容易サンプルがより強く割引されるので、極端な不均衡では gamma=3.0 程度まで上げる試行もよく行われます。
学習ループ — BCE vs Focal Loss を比較
同じアーキテクチャと同じ初期化シードで、損失関数だけを変えた 2 つのモデルを学習し、比較します。
def train(loss_fn, epochs=15, batch_size=64, lr=1e-3):
torch.manual_seed(42) # 公平比較のため
model = FlareCNN()
opt = torch.optim.Adam(model.parameters(), lr=lr)
ds_tr = TensorDataset(torch.from_numpy(X_tr), torch.from_numpy(y_tr))
dl_tr = DataLoader(ds_tr, batch_size=batch_size, shuffle=True)
for ep in range(epochs):
model.train()
for xb, yb in dl_tr:
opt.zero_grad()
logit = model(xb)
loss = loss_fn(logit, yb)
loss.backward()
opt.step()
model.eval()
with torch.no_grad():
prob_te = torch.sigmoid(model(torch.from_numpy(X_te))).numpy()
return prob_te
# BCE （クラス重みなし）
bce_loss = nn.BCEWithLogitsLoss()
prob_bce = train(bce_loss)
# Focal Loss
focal_loss = FocalLoss(alpha=0.25, gamma=2.0)
prob_focal = train(focal_loss)
実行すると、prob_bce と prob_focal という 2 つの予測確率ベクトル（テストデータ 400 件分）が得られます。これらを使って、不均衡データに特化した宇宙天気の評価指標と一般的な ROC/PR 曲線を計算していきます。
TSS と評価指標の計算
TSS = TPR - FPR を閾値ごとに計算し、最大値（最適閾値での TSS）と HSS、混同行列をまとめて求める関数を用意します。
def evaluate(prob, y_true, name=""):
fpr, tpr, thr = roc_curve(y_true, prob)
tss_arr = tpr - fpr
best_idx = int(np.argmax(tss_arr))
best_thr = thr[best_idx]
yhat = (prob >= best_thr).astype(int)
tn, fp, fn, tp = confusion_matrix(y_true, yhat).ravel()
# HSS
denom = (tp + fn) * (fn + tn) + (tp + fp) * (fp + tn)
hss = 2 * (tp * tn - fp * fn) / denom if denom > 0 else 0.0
roc_auc = auc(fpr, tpr)
print(f"[{name}] AUC={roc_auc:.3f} TSS={tss_arr[best_idx]:.3f} "
f"HSS={hss:.3f} TP={tp} FN={fn} FP={fp} TN={tn}")
return {"fpr": fpr, "tpr": tpr, "thr": thr, "tss": tss_arr,
"best_idx": best_idx, "best_thr": best_thr,
"yhat": yhat, "auc": roc_auc, "hss": hss}
res_bce = evaluate(prob_bce, y_te, "BCE")
res_focal = evaluate(prob_focal, y_te, "Focal")
このコードを実行すると、両モデルの AUC、TSS、HSS が表示されます。典型的な結果としては、Focal の方が TSS と HSS の双方で BCE を上回ります。BCE は正例があまりに少ないと「全て負例と予測する」自明解に近づきやすく、TSS が 0 に張り付きやすい傾向があるためです。Focal は容易サンプル（多数の容易な負例）の損失寄与を抑え、難しいサンプル（正例とハードな負例）に学習を集中させる効果が明確に現れます。
ROC/PR 曲線と TSS 曲線の可視化
評価指標を視覚化することで、損失関数の違いがどのような操作点（閾値）で効くのかを直感的に把握できます。
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# ROC
axes[0].plot(res_bce["fpr"], res_bce["tpr"], label=f"BCE (AUC={res_bce['auc']:.3f})")
axes[0].plot(res_focal["fpr"], res_focal["tpr"], label=f"Focal (AUC={res_focal['auc']:.3f})")
axes[0].plot([0, 1], [0, 1], "k--", alpha=0.4)
axes[0].set_xlabel("False Positive Rate")
axes[0].set_ylabel("True Positive Rate")
axes[0].set_title("ROC Curve")
axes[0].legend()
axes[0].grid(alpha=0.3)
# PR
prec_b, rec_b, _ = precision_recall_curve(y_te, prob_bce)
prec_f, rec_f, _ = precision_recall_curve(y_te, prob_focal)
axes[1].plot(rec_b, prec_b, label=f"BCE (AP={auc(rec_b, prec_b):.3f})")
axes[1].plot(rec_f, prec_f, label=f"Focal (AP={auc(rec_f, prec_f):.3f})")
axes[1].set_xlabel("Recall")
axes[1].set_ylabel("Precision")
axes[1].set_title("PR Curve")
axes[1].legend()
axes[1].grid(alpha=0.3)
# TSS vs threshold
axes[2].plot(res_bce["thr"], res_bce["tss"], label="BCE")
axes[2].plot(res_focal["thr"], res_focal["tss"], label="Focal")
axes[2].set_xlabel("Threshold")
axes[2].set_ylabel("TSS = TPR - FPR")
axes[2].set_title("TSS vs threshold")
axes[2].legend()
axes[2].grid(alpha=0.3)
plt.tight_layout()
plt.show()
このプロットからは 3 つの重要なことが読み取れます。第 1 に ROC 曲線では両者の差が一見小さく見えても、PR 曲線では Focal の優位がはっきり現れる ことが多いという点です。これは不均衡データでは ROC 曲線が誤って楽観的な評価を返しがちであることの典型例で、フレア予測のような正例少数タスクでは PR 曲線の確認が必須です。第 2 に TSS のピーク位置（最適閾値）が損失関数によって動く ことが見えます。Focal は確率出力の分布が変わるため、運用時の閾値設定もモデルに応じて再調整する必要があります。第 3 に 「常にネガティブ」予測の TSS は 0 になる という性質が、TSS 曲線の両端で確認できます。これが「Accuracy では区別できないモデルを TSS は厳密に区別する」根拠になっています。
混同行列の可視化
最後に、最適閾値での混同行列を並べて比較してみましょう。運用面では「見逃し（FN）」と「誤報（FP）」のバランスをどう取るかが本質的な意思決定であり、混同行列はそれを直接示してくれます。
def plot_cm(ax, y_true, yhat, title):
cm = confusion_matrix(y_true, yhat)
im = ax.imshow(cm, cmap="Blues")
for i in range(2):
for j in range(2):
ax.text(j, i, str(cm[i, j]), ha="center", va="center",
color="white" if cm[i, j] > cm.max() / 2 else "black",
fontsize=14)
ax.set_xticks([0, 1]); ax.set_yticks([0, 1])
ax.set_xticklabels(["Pred Neg", "Pred Pos"])
ax.set_yticklabels(["True Neg", "True Pos"])
ax.set_title(title)
return im
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
plot_cm(axes[0], y_te, res_bce["yhat"], f"BCE (TSS={res_bce['tss'][res_bce['best_idx']]:.2f})")
plot_cm(axes[1], y_te, res_focal["yhat"], f"Focal (TSS={res_focal['tss'][res_focal['best_idx']]:.2f})")
plt.tight_layout()
plt.show()
混同行列を比べると、Focal Loss の方が True Positive（フレアの正しい検出）が増えつつ、False Positive（誤報）も比較的抑えられる ことが視覚的に確認できます。実フレア予測では「TP の取りこぼし」が衛星セーフモード遷移漏れに直結するため、TPR を上げる損失関数の選択は単なる数値最適化ではなく、運用安全性に直結する設計判断です。一方で FP を抑えることは「狼少年問題」を避け、運用者の予報への信頼を保つために重要で、HSS が高い設計ほどこのバランスが取れたモデルだといえます。
実装と評価を通じて、不均衡データを扱う上で focal loss と TSS/HSS の組み合わせが効果的であることが具体的に確認できました。最後に、これらを実フレア予測パイプラインに発展させるための展望を整理しましょう。
実フレア予測への発展 — Surya 以後
公開データセットとパイプライン
合成データでの実験を実データに発展させるには、公開データセットの活用が出発点になります。代表的なのは
- SWAN-SF（Space Weather Analytics for Solar Flares）: SHARP パラメータの 12 分時系列を 5 つの時間分割で公開。ベンチマーク標準
- SDO/HMI Active Region Patches: JSOC 経由で活動領域パッチを直接ダウンロード可能
- DeepSDO / SuryaBench: 多波長のフルディスク画像とラベルが揃った大規模データセット
これらに加え、欠損補完、磁気図の視線方向投影補正、太陽の周辺減光補正、活動領域のトラッキングといった前処理パイプラインが実用では必須です。
基盤モデルへの移行
Surya のような基盤モデルは、Hugging Face から重みをダウンロードして数行で読み込み、独自の活動領域ラベルでファインチューニングできるようになりつつあります。これにより、これまで研究グループごとに個別に CNN/LSTM を最初から訓練していた工程が、事前学習済み Transformer のヘッドだけを学習する 軽量なワークフローに置き換わります。
具体的な利点は次の通りです。
- データ効率: 数万件レベルの活動領域ラベルで実用性能に到達できる
- マルチタスク: フレア発生確率、ピーク強度、CME 発生確率、F10.7 回帰などを共通バックボーンで扱える
- 解釈性: Attention マップを通じて、モデルが活動領域のどこを「見ている」かを可視化できる
物理拘束と PINN・MHD のハイブリッド
最後に、純粋なデータ駆動の限界として 2003 年のような極端イベントは学習データに少なすぎて統計が取れない という問題が残ります。今後の方向性として、MHD シミュレーションで生成した合成データを augmentation に使う、あるいは PINN（Physics-Informed Neural Network）で磁気エネルギー保存則を制約項に組み込む といったハイブリッドアプローチが盛んに研究されています。disassemble-channel の宇宙 AI シリーズでも、この方向の記事を順次拡充していく予定です。
まとめ
本記事では、太陽フレア予測を機械学習でアプローチする際の全体像を整理しました。
- 太陽フレアと CME は衛星・電力網・通信に実害をもたらす現実的な脅威であり、24 時間以内の M/X クラス発生予測が標準タスク
- 入力は SDO/HMI 磁気図（生画像）または SHARP パラメータ（手作り特徴量）
- 正例率 3〜7% の極端不均衡データであり、Accuracy ではなく TSS と HSS で評価する必要がある
- 学習側の不均衡対策として、クラス重み付き BCE と focal loss が実務的に有効
- 手法の系譜は SVM/RF → CNN → Transformer/LSTM → 基盤モデル Surya と進化し、地球観測基盤モデルと同じパラダイムシフトが起きている
- PyTorch の小型 CNN + focal loss でも、BCE 単独より TSS と HSS が改善されることを合成磁気図で確認した
次のステップとして、以下の記事も参考にしてください。