木星磁気圏の光学観測データに対するLassoを利用した変動周期のスパース推定
鈴木文晴1、吉岡和夫1、吉川一朗1
1. 東京大学大学院 新領域創成科学研究科 複雑理工学専攻

2018/9/13 科学とデータ研究集会 @NICT

木星磁気圏とイオプラズマトーラス（ＩＰＴ）

図：イオ衛星の噴火の様子(c) NASA/Johns Hopkins University Applied Physics
Laboratory/ Southwest Research Institute

図：木星とIPTの模式図[after Audouze et al., 1988]• 木星の衛星イオは火山活動が活発で、その規模は太陽系の全ての天体の中でも最大級である。
• 噴火によってイオの公転軌道付近に供給されたガスが電離し、円環状に分布している（イオプラズマトーラス）。
• 硫黄、酸素の原子やイオンが電子衝突励起に起因して可視域、紫外域で発光している（発光強度は～1012 W）。

IPTの観測意義

図：木星磁気圏における粒子輸送の概念図[Yoshikawa et al. 2017]• 木星半径の約100倍に広がる磁気圏において、IPTは木星から木星半径の6倍と比較的内側に位置している。
• 粒子、エネルギーの輸送の始点は尾部のリコネクションと考えられているため、IPTの光学観測結果を外部磁気圏に関する観測結果と比較することで、内向き輸送の到達点が分かる。

観測データ※Rj ： 木星半径

◎ひさき衛星• 2013年にイプシロンロケットによって打ち上げられた地球周回衛星。
• 極端紫外光（波長52～148 nm）を観測波長とする分光器を搭載しており、木星をはじめとして、太陽系の惑星（時折彗星や月）を観測している。

朝夕

図2.1 EXCEED観測機による木星観測の概要
EXCEED観測機が木星観測時に使用するスリットはダンベル型の形状をしており、スリットの太い部分にIPT、細いオーロラ全体を入れて観測を行っている。

ひさき衛星によって観測されたIPTの発光強度
200

IPT emission [GW]

◎IPTの発光強度には顕著な周期変動が見られる
・周期：10時間、42時間◎関連すると考えられる周期・イオ衛星の公転周期：42.46 時間
・木星の自転周期：9.925 時間

IPT Dawn
朝側100300~40 hours~40 hours~10 hours

Dusk
IPT 夕側2001

定常的に存在するこれらの周期成分を推定しないと、突発的な増光現象を特定しにくい~40 hours

234567Day of Year in 2014 (2014/1/1 0:00 = 1)

図；IPTの朝側、夕側の発光強度の時間変動。

8

光学観測データの周期解析◎ 惑星科学における周期解析：（目的）系の定常状態に影響する物理現象を明らかにするため、あるいは定常状態と過渡的な変動を区別するために、特徴的な周期変動を抽出する。
（定式化）時系列データを有限個の周期関数の線形和で表す。

42時間
10時間

◎ イオトーラスに存在する周期変動• 42時間周期変動：イオ衛星の公転周期に同期→ イオ衛星と磁気圏プラズマ間の相互作用を示唆• 10時間周期変動：木星自転周期に同期→ IPTの磁気経度異方性を示唆※周期変動は、系に内在する多様な現象を表す証拠となる

図：Lomb Normalized periodogram法によって導出されたイオトーラスの発光強度の周期成分。 [Tsuchiyaet al., 2015]

しかし、ノイズによって生じた偽のピークを見分けるのは困難である。

モデル評価の基準◎ 統計解析の最終目的 → 「真のモデル」の推定※ 考えるべき数理モデルの要素• 関数系： 多項式、正弦関数、ガウス関数、指数関数etc.
• フィッティングの自由度（次元）• フィッティングパラメータの値、精度◎ 真のモデルに「近い」モデルを選択する

評
価
基
準

真のモデルへの「近さ」をどう表現するか① 手持ちのデータを最も良く表現する → 最尤法→ 仮定するモデル関数の種類が正しければこれでよい② 未知のデータに対する適合度（予測精度）を最大にする
基準量AIC= -2(最大対数尤度）+2（自由度）を最小とするモデルを選択（赤池、1973）→ 「複雑さ」がモデルの評価基準として肝要 （単純なほど良いモデル）

データへの適合度

予測精度

疎性（スパース）モデリング◎高次元データに内在するスパース(疎)性を利用することで、データから最大限の情報を効率よく抽出できる技術の総称。近年、機械学習や情報数理学の分野で注目されている。

自然科学の統計解析でのモデリングに求められる事項①データへの適合度を高める②必要十分な数の変数によってデータを適切に説明する◎ 多くの変数から少数の変数を抽出することは、データの背後にある本質的な構造の理解につながる。

Lasso [Tibshirani, 1992](Least Absolute Selection and Shrinkage Operator)「L1ノルム正規化最小二乗回帰」→重要なパラメータの抽出が可能に{𝛷𝑗 𝑥 } = 1, 𝑥, 𝑥 2 , … … , sin 𝑥 , sin 2𝑥 , sin 3𝑥 , … … 等◎ 線形回帰最小二乗法：

𝑦𝑖 = 𝑓𝛽 𝑥𝑖 ; 𝑓𝛽 𝑥𝑖 =𝛽𝑗 × 𝛷𝑗 (𝑥𝑖 )𝐸 𝛽 =𝑦𝑖 − 𝑓𝛽 𝑥𝑖

2𝑖→ Lasso 𝐸𝐿𝑎𝑠𝑠𝑜における定式化

𝑗𝛽 =𝑦𝑖 − 𝑓𝛽 𝑥𝑖𝑖
データへの適合度を示す項の最小化
2+λ×|𝛽𝑗 | の最小化

𝑗
係数の大きさに制限を設けるための「罰則項」※ λ：罰則項の重みを調節し、モデルの複雑度を決定するパラメータ

Lassoの性質
𝐸𝐿𝑎𝑠𝑠𝑜𝛽 =𝑦𝑖 − 𝑓𝛽 𝑥𝑖2

数学的に等価+λ×𝑖|𝛽𝑗 |𝑗• 罰則項によって不要なパラメータの値が ’0’と推定され、パラメータの抽出が可能• データへの適合度の最適化と、パラメータの刈り込みによる推定精度の最適化を両立可能• 線形回帰や線形変換に基づいた様々な手法に対して、Lassoを用いた拡張は容易に可能𝑁

例：離散フーリエ変換

𝑦𝑖 =
𝑗 =1

2𝜋𝑓 𝑗 𝑥 𝑖
−𝛽𝑗 𝑒 𝑁

罰則項による制約
図：Lassoによるパラメータ抽出の様子の模式図。
最小二乗回帰の解βから離れるとデータへの適合度が減少する。

𝑁→𝐸𝛽 =

𝑦𝑖 −𝑖𝑗 =122𝜋𝑓 𝑗 𝑥 𝑖−𝛽𝑗 𝑒 𝑁+λ|𝛽𝑗 |𝑗Lassoを周期解析に用いた先行研究◎対象：KIC J1924変光星• 発光強度の変動周期の時間変化を調査• 従来手法（離散フーリエ変換）とLassoによる修正を施した手法を比較[Kato and Maehara, 2013]• Lassoを用いる子とで変光星の変動周期の時間変化を明確に示すことに成功している。

図：KIC J1924に対する周期解析の結果[Kato and Maehara, 2013 より改変]
縦軸が周波数、横軸が時間、グレースケールが強度
左図が短時間フーリエ変換、右図がそれにLassoによる修正を加えた手法の結果

• 対象：IPTの朝側から放出されるS+ 76.5 nm
輝線の発光強度の時系列データ• 期間：2013年12月21日～2014年1月15日

IPT emission [GW]

解析対象
3020100-9-6-3036912Day of Year in 2014 (2014/01/01 00:00 = 1)

図：ひさき衛星によって得られたIPT朝側から放出される76.5 nmの輝線の発光強度の時間変動

• ひさき衛星の取得したデータにはデータ欠損が存在するため、補間が必要である。
• 今回は、一次関数による内挿を行なった。

3020100-9-6-3036912

Day of Year in 2014 (2014/01/01 00:00 = 1)

欠損データの補間
IPT emission [GW]

• 離散フーリエ変換は、サンプリン
グ間隔が均一で欠落のないデータに対して使用可能な手法である。

IPT emission [GW]

欠損データの補間3020100-9-6-30369Day of Year in 2014 (2014/01/01 00:00 = 1)

12

結果

𝐸λ𝐿𝑎𝑠𝑠𝑜

𝛽 =𝑦𝑖 − 𝑓𝛽 𝑥𝑖𝑖

図：IPTの発光強度の観測データ（黒点）及びLassoによって得られたフィッティング関数（緑線）。

2+λ×|𝛽𝑗 |𝑗λ

図：Lasso（緑線）及び離散フーリエ変換（黒線）
によって得られたIPTのパワースペクトル。

• 罰則項の重みを決定するパラメータ λ を変えることで、フィッティング関数の複雑さとデータに対する適合度を調節できる。
• 最適な λ の決定が肝要な問題となっている。

10分割交差検証によるλの決定◎ 目的
推定精度が最適となるような λ の値を導出する。
• データセットを二つに分割し、その片方（訓練データ）でフィッティングし、残りのデータ（検証データ）に対する誤差を計算する。
• データをランダムに10分割し、分割された標本群それぞれを検証データとする10回の検証を行い、平均値と標準偏差を求めた。

訓練データ

検証
訓練データ

訓練検証データ◎ 結果• 検証データに対する誤差は、λ=0.6で最小となる。
• λ=0～1.1では誤差(1σ)の範囲で同程度の推定精度である。
01（複雑度大）（one standard errorルール[Tibshirani, 2001]）

………

訓練データ

検証データ

予測誤差
Square ErrorMean◎ 手法：交差検証

全データ

訓練データ

23

λ

456（複雑度小）結果◎ 交差検証の結果より、λ=0 （離散フーリエ変換）とλ=1.1とおいたモデルでは、誤差(1σ)の範囲で同程度の推定精度を持つといえる。
◎ λ=1.1とおいたモデルでは、4種類の周期成分が残っている。
 １０時間周期 → 木星の自転周期 [Sandeland Dessler 1978 など]
 ４２時間周期 → イオの公転周期 [Sandeland Dessler 1978 など]
 ２０時間周期、１５時間周期 → 過去の研究によって報告されていない周期変動。

42時間

15時間 10時間
20時間

図：Lasso（緑線）及び離散フーリエ変換（黒線）によって得られたIPTのパワースペクトル。

今回の手法は、一次変換である離散フーリエ変換をフィッティングの問題に落とし込んでいるため、欠損データを補間せずに解析を行う事が可能である。
結果：• 42時間、10時間、20時間、15時間周期に加え、5時間周期の存在が示唆される結果が得られた。
• 5時間周期成分の存在は、過去の研究によって既に示されたものである[Tsuchiya et al., 2015]

IPT emission [GW]

欠損データの補間による影響

30
20
100-9-6-3036912

Day of Year in 2014 (2014/01/01 00:00 = 1)

欠損データを補間せずに周期解析を行った
42時間

15時間
10時間
20時間

5時間

図：Lasso（緑線）によって得られたIPTのパワースペクトルまとめ• 離散フーリエ変換にLassoを適用することで、観測データにとって本質的な周期成分を同定する新規手法を示した。
• データに対する予測精度が離散フーリエ変換の結果と同等で、より変数の少ない単純な数理モデルによってIPTの周期変動をモデル化することに成功した。
• 提案手法は欠損の存在するデータに対しても妥当な結果を与えており、観測に時間的制約のある場合にも適用可能であることが示唆される。
• Lassoは周期解析に限らず数理モデルのフィッティングで普遍的に適用可能であり、多様な観測データへの効果が期待される。