J. Plasma Fusion Res. Vol.97, No.2 (2021)６
４‐６５

小特集

磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展
Progress of Physical Modeling by Data-Driven Approachin Magnetically Confined Fusion Plasmas

１．はじめに
1. Introduction
今寺賢志
IMADERA Kenji
京都大学大学院エネルギー科学研究科（原稿受付：２
０
２
０年９月３
０日）

近年，計算機を用いた研究アプローチとして，従来のシションである．磁場閉じ込め核融合プラズマ分野では，第ミュレーションに加え，データ駆動科学ともいうべき新し一原理であるジャイロ運動論に準拠した大規模シミュレー

い方法が台頭している［１］．本学会誌においても近年，『核ションが精力的に行われているが，時間的にも経済的にも融合プラズマにおけるデータマイニングの活用』（２０１６年

高コストである．例えば，筆者が行っている大域的ジャイ５月号）［２］や，『プラズマ・インフォマティクス−データ

ロ運動論シミュレーションの場合，用いる計算機やプラズ

駆動科学のプラズマへの応用−』（２０１９年１１月号）［３］など

マサイズに拠るが，大型計算機を用いて数日から数週間程の小特集が取り上げられており，ここ数年のプラズマ・核度の計算時間が必要であることが標準的である．それ故

融合学会でも多くの関連講演が行われている．

に，物理パラメータの探索についても多くの時間を費やす

シミュレーションは，「何らかの数値計算手法を用いて，必要がある．そこで考えられるのが，第一原理シミュレー

与えられた支配方程式の近似解を求める」という演繹的方

ションや実験で得られた大規模データをデータ駆動科学に法論と定義できる一方，データ駆動科学は，「何らかの統基づいて処理することで，簡約化された物理モデルの開発

計数理的手法を用いて，得られたデータの背後にある法則を行い，高精度を維持しつつシミュレーションのコスト削性を抽出・推定する」という帰納的方法論である．また物減を行う方法である．この点については，文献［３］のまとめ理学では，要素還元主義に代表されるように，多くの事象でも「本小特集で紹介できなかった例で，今後重要になるを少数の法則で数学的に矛盾なく説明することが可能であと考えられのが，大規模数値シミュレーション（HPC）のるという考えに立脚する場合が多く，シミュレーションも低コスト化である」と言及されている．これは，シミュ

歴史的にこの価値観を共有してきた．一方，データ駆動科

レーションとデータ駆動科学が融合した良い事例である．

学は，そのような汎化性はもちろんめざしながらも，あく

また，統計解析の観点からもデータ駆動科学を考察してまで与えられたデータの内挿を高精度に行うことに主眼をみよう．古典的な統計解析は，「何らかの統計手法を用いおいており，少なくともモデルの単純さは多くの場合に優てモデルの特性を理解する」方法論である．例えば線形回先度が低いと考えられてきた．このような点が，データ駆帰を行えば，２つの変数の線形応答性を理解することがで

動科学を導入する際に従来の科学的アプローチとの差異をきるであろう．一方でデータ駆動科学では，モデルは（あ感じる要因ではないかと筆者は考えている．

くまで人間が直感的に理解困難という意味で）ブラック

一方で，そのような異なった特性を有したシミュレー

ボックス的な関数であって，モデルそれ自体を理解するこションとデータ駆動科学を併用することで，相乗的な結果とが必ずしも主目的にはならない場合が多い．いわゆる機を得られる場合がある．その一例が低コストシミュレー

械学習はこのケースに該当する．つまり端的には，統計解

Graduated School of Energy Science, Kyoto University, Gokasho, Uji, KYOTO, 611-0011, Japanauthor’s e-mail: imadera@energy.kyoto-u.ac.jp

６４!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

1. Introduction

K. Imadera

析は近似解をめざすという意味でシミュレーション的なア関係を推定する方法を紹介する．大域的流体シミュレー

プローチであり，機械学習はその再現をめざすという意味

ションによって得られたケルビンヘルムホルツ乱流にでエミュレーション的なアプローチに分類される．

DMD 法を適用した結果，乱流場は背景分布の変動，帯状

しかしながら，そのような説明性の観点からもデータ駆流，コヒーレント揺動，スパイラル構造に分解され，得ら

動科学の分野で様々な試みが行われており［４］，その１つ

れた乱流構造の強度を示す指標を導入しその時間発展を見が次元縮約手法の一種である主成分分析［５］の応用である．

ることで，構造間の因果関係が同定可能であることを本章

機械学習の分野では，計算コスト削減の観点からしばしば

でまとめる．また，実験データから乱流の時空間発展の予主成分分析が用いられてきたが，本質的に重要な特徴量を測に応用した試みについても紹介する．

抽出して次元を下げることで分析や解釈が容易となり，特第５章 『データ駆動アプローチを用いた雪崩的乱流輸に３次元以下に次元削減を行えば，データセットの様子を送現象の解析』では，データ駆動科学的手法である主成分

直接可視化することができるため，説明可能な新たな物理分析を５次元位相空間構造の時系列データに適用すること

モデルの開発に貢献することも期待される．これは統計解

で，ジャイロ運動論シミュレーションにおける突発的輸送

析とデータ駆動科学によって相補的に手法が進展した良い現象と関連したモードの位相空間構造を同定する方法を紹事例である．

介する．多次元データ内の乱雑さを，特異値分解を用いて

本小特集は，磁場閉じ込め核融合プラズマの分野で，そ定量化することで位相空間ダイナミクスを整理した結果，のようなシミュレーション，統計解析，データ駆動科学の突発的輸送の前に線型モードの形成による乱雑さの減少が融合によって行われた物理モデリングに関する興味深い４

起き，突発的輸送の直後に位相混合による乱雑さの増加が編の論文から構成されている．第２章 『機械学習による起きることを本章でまとめる．また，主成分分析による位

輸送モデリング ‐JT-60U における実践例‐』では，データ

相空間構造データの圧縮技術を開発した点についても言及

駆動科学的手法により，物理現象の再現性を保ちつつ大幅する．

に計算コストを削減した輸送モデリングとして，第一原理

しばしばデータ駆動科学は，データの「補間」，つまり高シミュレーションの結果を学習して得られたニューラル

次元空間におけるデータ点を滑らかに繋ぐ曲線を求める作

ネットワーク（NN）モデルと，大域的最適化手法で生成さ業に帰結されると言われている．本小特集では更に踏み込れた膨大なデータを学習して得られた NN 代理モデルを紹んで，シミュレーション，統計解析，データ駆動科学が介する．また，これらのモデルを実装し，JT-60U 実験を対各々の特性を相互に「補完」することで，物理モデリング

象とした輸送シミュレーションの結果についても言及すの進展に資することに焦点を当てている．データ源としてる．

は，シミュレーションを実験に置き換えても場合によって第３章 『データ同化による輸送モデリング ‐LHD におは差し支えないし，今回は磁場閉じ込め核融合プラズマのける実践例‐』では，データ駆動科学の一手法であるデータ

事例に限定しているが，その方法論自体は超領域性を有し同化を用いた輸送モデリングに関する研究を紹介する．こており，広くプラズマ科学の分野にも応用可能であると考の手法を用いて，ある輸送モデルに基づいたシミュレー

えられる．こうした観点を含めて，本小特集に掲載した４ションで得られた計算値と実際の観測値を結びつけ，その編の論文を読んでいただければ幸いである．

関係が統計的に妥当になるように輸送モデルを最適化する
方法について言及する．併せて，LHD 実験を対象とした

参考文献

データ同化システムによる輸送シミュレーションの結果［１］T. Hey et al., The Foruth Paradigm Data-Intensive Scientific
Discovery, (Microsoft Research, Redmond, Washington,2009).
［２］山本 聡 他：プラズマ・核融合学会誌 92, 332 (2016).
［３］浜口智志 他：プラズマ・核融合学会誌 95, 535 (2019).
［４］R. Guidotti, ACM Computing Surveys 51, 93 (2018).
［５］I.T. Jolliffe, Principal Component Analysis (Springer, 2002).

や，熱拡散係数の回帰モデルの妥当性を検証した結果についても言及する．
第４章 『データ駆動アプローチを用いた動的乱流現象の解析』では，データ駆動科学的手法である動的モード分
解（Dynamic Model Decomposition; DMD 法）を用いて乱流の時空間構造の抽出を行い，抽出された乱流構造間の因果

６５

J. Plasma Fusion Res. Vol.97, No.2 (2021)６
６‐７１

小特集 磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展

２．機械学習による輸送モデリング‐JT-60U における実践例‐
2. Machine-Learning Assisted Transport Modeling-Practical Cases in JT-60U成 田 絵 美，本 多充

NARITA Emi and HONDA Mitsuru
量子科学技術研究開発機構（原稿受付：２
０
２
０年９月３
０日）

核融合プラズマは幅広い時空間スケールの物理現象に支配されるため，その性能予測や物理解明には，各時空間スケールの物理現象を扱うモジュール群から構成される統合モデルが必要不可欠である．構成要素の一つである輸送モジュールはプラズマの巨視的な時間発展を扱う．ここでは，乱流輸送モデルを用いて乱流流束を予測するが，物理モデルの精緻化により計算コストが増大したため，統合モデル全体の計算速度を低下させていた．
そこで，データ駆動科学的手法により，物理現象の再現性を保ちつつ大幅に計算コストを削減した輸送モデルが国内外で開発されている．本章では，第一原理計算と JT-60U 実験で得られたデータに基づく人工ニューラルネッ
トワーク（NN）モデル及び，大域的最適化手法の計算途上で生成される膨大な輸送計算データを学習することで構築される，輸送モデルを模擬する NN 代理モデルを紹介するとともに，JT-60U 実験に対する輸送シミュレー
ションの結果を示す．
Keywords:
machine learning, neural network, gyrokinetic simulation, stiff transport model, integrated code, global optimization,
genetic algorithm, Nelder-Mead method

１０３回以上繰り返すことは現実的ではなく，ましてや統合

２．
１ はじめに核融合プラズマの性能予測や物理研究に用いられる統合

コードにおける輸送モデルとして用いることはできない．

モデルは，輸送ソルバを核として，MHD 平衡，ダイバー

そこで輸送モジュールには理論モデルを数値計算でフィッ

タ，加熱などの異なる時空間スケールの物理を扱う複数のティングした関数形を用いる CDBM［２］や，ジャイロ運動モジュール（コード）群から成る［１］．核融合出力を決定論モデルを簡約化した TGLF［３］や QuaLiKiz［４］などのづけるプラズマのコア領域における密度と温度を予測する

高速に計算できる輸送モデルが導入されている．TGLF やという観点からは，粒子と熱が磁気面を横切りプラズマ外QuaLiKiz はジャイロ運動論コードと比べて１０６倍程度高速部へと運び出される正味の量である粒子束と熱流束を正確化されているものの，統合コードの計算速度を低下させる

かつ高速に求めることが重要となる．プラズマの定常状態

主要因となっており，定常状態の予測に並列計算で数日かにおける温度と密度の予測には，設定する時間刻み幅にもかる場合がある．限られた放電の実験解析などには有用で３６よるが，輸送方程式を１０ から１０ 回繰り返し解き定常状態

あるが，多くの試行回数を要するトカマク運転シナリオのへとプラズマを緩和させる必要があり，必然的に同じ回数

開発には不向きであるため，CDBM などのより簡略化されもしくはそれ以上に粒子束と熱流束を算出しなくてはならた輸送モデルが用いられることが多く，第一原理計算に基ない．JT-60U などのトカマクプラズマでは，乱流が輸送をづく乱流輸送物理は統合モデルに十分に取り込まれている

支配しているため，密度と温度の分布予測精度は乱流に起とは言い難い状況であった．

因する粒子束と熱流束の予測精度に強く依存する．粒子束

統合コードにおける輸送モジュールの高速化のたと熱流束を乱流輸送物理に最も忠実に計算する手段とし

め，２０１４年頃から機械学習モデルの一種である人工ニューて，ジャイロ運動論コードを用いた第一原理計算が挙げら

ラルネットワーク（NN）モデルの導入が欧米で始まった．

れるが，一般にジャイロ運動論コードの実行はスーパーコ

輸送モジュールの高速化に用いられる NN は回帰問題をンピュータを用いた大規模並列計算が必要である．更に，扱っており，乱流輸送物理を特徴付ける複数の物理量を入一度の計算に数日から一週間程度を要するため，これを力変数として，粒子束や熱流束を出力するモデルとなって

National Institutes for Quantum and Radiological Science and Technology, Naka, IBARAKI 319-1112, Japancorresponding author’s e-mail: narita.emi@qst.go.jp

６６!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

2. Machine-Learning Assisted Transport Modeling -Practical Cases in JT-60U‐

E. Narita and M. Honda

いる．米国の研究グループは，OMFIT［５］と呼ばれる統合

で，統合シミュレーションにおいては核融合出力を左右す

モデルフレームワークの中でDIII-Dの実験データから得らる密度や温度の径方向分布を上手く制御して理想的な状態

れた熱流束を再現する NN の開発を皮切りに［６］，現在ではを得ることが求められるため，フラックスを生み出す背景

TGLF が予測する熱流束・粒子束・運動量流束を模擬する

物理と分布形成の関係を理解することは有益である．そこ

NN 代理モデル（surrogate model）や，ペデスタル安定性

で，粒子束と熱流束が準線形理論で表現できることを仮定を評価する EPED1 モデルの挙動を模擬する NN 代理モデすることでフラックスと熱力学的力の関係を陽に書き下しルを用いるまでに拡張を進めている［７］．欧州の研究グ

た，電 子 の 粒 子 束 と 熱 流 束 を 予 測 す る 輸 送 モ デ ル

ループでは，QuaLiKiz の挙動を模擬する NN 代理モデルがDeKANIS が構築された［１３］．

開発されている．最初に発表されたモデルは QuaLiKiz が#& !!$$%"#' !!%%"
Γ &""#（１）$&"#&##' !!"$$%"#& !!#%"（２）&

予測する熱流束のみを再現する NN であったが［８］，米国&&&のグループと同様に粒子束と運動量流束まで予測範囲が拡

ここで，Γ &と $&は規格化された電子の粒子束と熱流束，張され，JET の実験データを基に NN の学習に用いるデー
タを生成する統合モデリングの枠組みの構築と共に，開発#'" と #&は規格化された電子の粒子・熱拡散係数，%"が精力的に進められている［９，１０］．学習に用いられるデーと %"#& は電子の規格化密 度・温 度 勾 配 で あ る．（１），タに違いはあるが，いずれのモデルでもTGLFやQuaLiKiz（２）式の右辺の括弧内における第一項は対角項，第二・三と比べて１０５倍程度の高速化に成功している．特に欧州で項は非対角項であり，それぞれ拡散・対流という輸送過程&&

は，従来型の大規模な統合コードよりも簡素化することでに対応している．非対角項における係数 !%!$!#!"$ は対流の高速化された統合コード RAPTOR［１１］に NN 輸送モデル

大きさと向きを表しており，ジャイロ運動論シミュレーを導入し，実時間シミュレーションをめざしている［９］．

ション に よ っ て 求 め る こ と が で き る［１４，１５］．こ こ で，

実時間シミュレーションについては，［１２］も参考のこと．

#& は Γ &と $&を駆動する熱力学的な力である%"#' と %"&&

本章では国内で開発されている２種の NN 輸送モデルをが，Γ &と $&の一方の対角項を駆動する熱力学的 な 力紹介する．これらの特徴を欧米で開発されているモデルとは，他方の非対角項の駆動力となっている．このとき，互共に表１にまとめる．２．
２節で紹介するモデルは他のモデいの非対角項は等しくなるという Onsager 対称性が成り立ルと異なり，第一原理計算と実験データの組み合わせを学

つ．この対称性を崩さないよう，!%!$!#!"$，"，#&の計６つ

習に用いている．２．
３節で紹介するモデルは米国の研究グの係数を求めることで，拡散・ピンチに起因する Γ &とループと同様に TGLF を模擬する NN 代理モデルである$&を無矛盾に予測している．電子の粒子束と熱流束に加が，大域最適化手法に基づく定常輸送コードが実行時に大えて，DeKANIS はイオンの熱流束を次のように表現して量に生成するデータを用いて訓練を行った，予測精度に重予測する．

きを置いた代理モデル開発という点に特徴がある．これら$("のモデルを JT‐６０U 実験に対する輸送シミュレーションの結果と共に紹介する．

#(!&''
'&
# ( (
# %"
#&!&'' &!&'' & '& &&(（３）

ここで，#&''は $ を対角項のみで表現したときの実効熱拡

２．
２ 第一原理計算と実験データに基づくニュー
ラルネットワーク輸送モデリング'&と散係数であり，%"#& はイオンの規格化温度勾配，'("(&("&&はイオンと電子の密度・温度比である．実効熱拡散

既存の輸送モデルは拡散係数もしくは輸送フラックスと#&!&''を求めることでイオンの熱流束を予測す
係数比 #(!&''"

して出力するため，そのフラックスを算出した背後にあるる こ と が で き る た め，上 記 と 合 わ せ て７つ の 係 数

様々な熱力学的力に対する依存性を個別具体的に知ること"#&!&''，"，#&）を決定する必要がある．
（!%!$!
#!
"$，#
(
!
&'
'

はパラメータサーベイでもしない限り困難である．一方７つの係数は NN を用いて求められるが，その学習に用いられるデータはジャイロ運動論コード GKW［１６］の数値
計算結果と JT-60U の実験データによって構築されている．

表１

国内外の NN モデルを用いた輸送モデル開発の特徴．

#&!&''は，計算コストが低い GKW の線形計!%!$!#!"$ と #(!&''"
算によって得ることができるが，"や #&の決定には乱流揺動の飽和レベルを評価するために非線形計算を必要とす
る．前述のようにその計算には数日以上の時間を要するた
め，学習データとして十分な回数を実施することは現実的に不可能である．そこで，JT‐６０U の実験値から " を半経験的に求めることで，輸送過程の理解に重要な非対角項の係数はジャイロ運動論コードによる第一原理計算に依拠しつつ，十分な学習データを得ることを可能にしている［１５］．残る #&は Onsager 対称性から一意に決まる．
JT-60U の H モード放電の実験データに基づき構成された７千点程度の学習データを訓練・検証用データに分け，
６７

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

NN を構築した（図１）．ここで，NN の出力変数は６つのることとなる．十分広い変数領域に訓練データが及んでい!
"'!'((，"）とし，"'は含めていない．
係数（!&!%!$!#%，")!
'((

ない場合は，このような外挿が起こりやすく，特に本例の"'はOnsager対称性を崩さないようにこれらの係数から求ように NN を輸送モデルとして用いる際は，訓練データのめている．一方，入力変数には温度や密度の勾配など，範囲内から大幅に外れた位置における過度な外挿が輸送方

GKW の入力変数群に近い１２の物理量としている．NN は全程式を安定に解くことを妨げる可能性がある．そこで Soft-

結合型のモデルであり，Python ライブラリである Tensor-

sign 関数の採用によって出力値の幅を制限することによ

Flow［１７］をバックエンドとして，Keras［１８］を用いて構築

り，安定した輸送計算を実現している．構築した NN を用された．NN の表現力（複雑さ）を決める隠れ層の数や，各いて DeKANIS は電子の粒子束及び電子・イオンの熱流束

隠れ層が持つユニットの数などのハイパーパラメータは，を予測する．学習データの構築に含まれないテスト用のベ イ ズ 最 適 化 を 利 用 す る Python ラ イ ブ ラ リ Hyperopt

JT-60U の H モード放電に対する DeKANIS と TGLF の予［１９，２０］を Keras のラッパーである Hyperas［２１］を介して測誤差を図２に示す．JT-60U の H モード放電に 対して用いることで調整している．隠れ層の数を１から４の範囲

TGLF よりも精度が良いことが確認できる．

で，各隠れ層のユニット数は１０から３００の範囲で最適化を

DeKANIS は TOPICS［２２］と GOTRESS［２３，２４］に導入さ実施した結果，２層の隠れ層を有し，ユニット数はそれぞれ，密度や温度の分布予測に用いられている．これらのれ２５０と１７０の NN が得られた．ここでは，ドロップアウト

コードは Fortran で書かれているため，Python を用いて構率の最適化も行っており，各隠れ層での平均は約 0.02 であ築された NN を高速化のために Fortran で書き換えて輸送

る．ドロップアウト率とは，モデルの訓練時にランダムにモデルとして利用している．図３はテスト用の JT-60UのH学習に参加させないユニット数の割合を示しており，推論

モード放電３９０６１に対して密度及び温度の定常分布予測を時には影響を及ぼさないパラメータである．ドロップアウ行なった例である．TOPICS では密度と温度を同時に解いトによって過学習のリスクを軽減している．活性化関数にており，ここで得られた密度を用いて，GOTRESS では温は，出力層に Softsign 関数を，出力層以外に Rectified Lin-

度のみを解いている．TOPICS の輸送ソルバは輸送方程式

ear Unit（ReLU）を用いている．回帰問題を解く NN の出を時空間に離散化して解く従来型の輸送コードである一力層の活性化関数は線形関数とする例が多いが，Softsign

方，GOTRESS は遺伝的アルゴリズムとネルダー・ミード関数を用いることで上限・下限値を設けることができる．

法による大域最適化手法を用いて輸送方程式の解を求めてここで，上限・下限値は訓練データの各出力変数の平均値いる．このように異なる解法でも，両輸送コードで一致す#の値で設定している．ここで #は訓練データを標から !"る分布が得られることが確認できる．分布予測に要する時本としたときの標本標準偏差である．線形関数を用いる

間は，TOPICS を用いた場合，１コアで計算して３時間程と，上限・下限値に近くデータ点数の少ない変数領域の訓度であり，従来のTGLFやQuaLiKizなどを用いた予測より練データに強く依拠して訓練データの範囲外の値を予測すも十分早い．一方，GOTRESS を用いた場合，図３のケー
スでは２６コアを利用し１分以下で収束解が得られている．
DeKANIS は乱流フラックスが密度や温度の空間勾配に強く依存するジャイロ運動論シミュレーションの結果を基に構築されているため，硬い輸送モデルに分類される．硬い輸送モデルを用いて TOPICS のような従来型の輸送コード
で分布予測を安定して行うには，時間刻み幅を非常に細かく取る必要があり，計算時間が増大する傾向にあるが，GOTRESS では分布の勾配値を大域的最適化手法で直接求図１めるために数値振動が起きづらく，より高速な分布予測が

DeKANIS が用いる NN モデルの構造．出力された変数から
電子粒子束及び電子・イオン熱流束が得られる．

可能となっている．次節では，この GOTRESS の枠組みの中で構築する NN 代理モデルを紹介する．

図２

テスト用の１
４の JT-60U H モードプラズマに対して DeKANIS と TGLF を用いて（a）
電子粒子束及び（b）
電子・（c）
イオン熱流束を予測したときの予測誤差（二乗平均平方根誤差 #）の径方向分布．

６８

Special Topic Article

2. Machine-Learning Assisted Transport Modeling -Practical Cases in JT-60U‐

E. Narita and M. Hondaの多く の 試 行 が な さ れ た と い う こ と は，NN モ デ ル で
TGLF の振る舞いを模擬するに足る多数のデータが収集できたということに他ならない．さらに，ランダムや一様なサンプリングと異なり，GOTRESS はあくまで解を求めよ
うとしており解の周辺に多くのデータを生成させるため，
そのデータから NN モデルを作れば，現在扱っているケースにおける正しい解の近くの TGLF の振る舞いをよく知っていると期待できる．そのため，TGLF の NN 代理モデルの構築によるシミュレーションの高速化を試みる．図４に
NN 代理モデルの構築から定常解を得るまでの過程を示す．データ収集のための GOTRESS シミュレーションにおいて，遺伝的アルゴリズムにおける個体数を１００，世代数上限を通常シミュレーション時の上 限 の１０分 の１で あ る
１，０００世代とし，異種粒子種間の等温化過程計算のための繰り返し計算を３回で終了させることとした．径方向位置
を空間８点とすると，100×1,000×8×3＝2.4×106 のデータが取得できる．前処理として，有効数字５桁で四捨五入し
図３

テスト用の JT‐
６
０U H モードプラズマに対する（a）
電子密度
及び（b）
電子・（c）
イオン温度の DeKANIS による予測と計測値の比較．温度の予測には統合コード TOPICS と GOTRESS を用いた２種の結果を示す．

た数値データが全て一致したデータセットは同じであるとして捨てることで，データ量を百万程度にまで減らせる．
このシミュレーションでは温度や温度勾配の可動域をかなり広く取っているため，現実にはあり得ない値も評価してし ま う（例 え ば 電 子 温 度 "$!0.1 keV，イ オ ン 温 度"%!20 keV という組み合わせも数値的に作り得るため，温２．
３ 遺伝的アルゴリズムを利用したニューラル
ネットワーク代理モデルの構築

度比は核融合プラズマでは非現実的な "$"!
"%!!
!!#となっ輸送モデルは統合コードの中で定常解を得るまでに何万てしまう）ことが起こりうる．このような外れ値は前提条

回，時には何百万回以上も呼び出されることを考えれば，件が非物理的であるだけでなく，そもそもの輸送モデルのジャイロ運動論シミュレーションからは大幅に簡約化され適用範囲外である可能性が高いため，得られる輸送フラッたとは言え，TGLF モデルをそのまま統合コードの中で運クスも意味をなさないおそれがある．これらは NN モデル

用するには時間的な制約が大きい．その数値コストからの訓練においても有害であるため，NN モデルを訓練する

TGLF は内部で波数ごとに MPI 並列化されているが，標準#に入ったために用いるデータは，各入力変数の分布の ""的な２３並列（１波数１コア）として JFRS-1 上で計算してものだけに限定する．刈り込んだ後のデータ数はおよそ８０

も，評価１回当たり約 0.15−0.2 秒掛かる．統合シミュレー

万となる．なお，TGLF の NN 代理モデルの入力変数は２
３，

ション終了まで TGLF を百万回呼ぶ必要があると仮定する

出力変数は１２である．出力変数は TGLF に倣い，電子，主と，TGLF の評価だけで約５０時間以上も掛かってしまうこ

イオン，不純物それぞれに対する粒子束，熱流束，低波数とになり，他のモジュールに必要な計算時間と照らし合わからの寄与のみの熱流束，運動量流束となっているが，以せても改善の余地は大きい．

下では熱流束の情報しか用いない．入出力変数の詳細はGOTRESS は他の輸送コードと大きく異なり，遺伝的ア［７，２３］参照のこと．NN モデルは３層の隠れ層を有する全ルゴリズムなどの大域的最適化手法を用いて定常輸送方程結合 NN であり，各層に２００ユニットを配しており，ドロッ

式を解いている．GOTRESS の数値アルゴリズムの詳細はプアウト率は 0.25 とした．

［２３，２４］に譲るが，端的に言えば，各径方向位置における温このようにして訓練された TGLF の NN 代理モデルが予度と温度勾配を１つの組として輸送方程式を満たす解の組

測した電子とイオンの熱流束で決定係数 !"を算出するとをパラメータ空間内で直接探索している．GOTRESS で採用している大域的最適化手法はニュートン法といった局所最適化手法と異なり内部で微分計算を必要としないため，勾配値に敏感な硬い輸送モデルを扱う上で有利な点を有している一方で，解に到達するためには多くの試行（この場
合は TGLF の計算）を積み重ねなくてはいけない欠点もあ

図４る＊．しかし，この欠点はデータ駆動科学の観点からは有利と取れる面もあり，広いパラメータ空間において TGLF

GOTRESS によって NN 代理モデルを構築し，温度の定常分布を予測する流れ図．NN 代理モデルは図１に示すモデルと同様の構造である．

＊：但し，GOTRESSが安定に求解できる点を考慮すれば，数値安定性のために時間刻み幅を大幅に細かくしなければいけない従来型の輸送コードと比べて結果的に計算時間に遜色ないかむしろ高速である可能性がある点には留意が必要である．前述のDeKANIS を使用した例では１桁程度高速である．

６９

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

それぞれ0.962，0.897となり，かなり良い再現性を示していミュレーションによって明らかにする．図５は GOTRESS

ることがわかる．NN 代理モデルを１回評価するのに必要に よ る 定 常 時 の（a）温 度 と（b）拡 散 係 数 の 予 測 で あ り，な計算時間は１コアで 0.1 ミリ秒であり，１コアで計算しTGLF を用いた場合と最適化 NN 代理モデルを用いた場合

た TGLF と比較すると 3×104 倍程度高速であることが分の結果を比較している．同じ初期条件から始め，収束に至かった．NN 代理モデルを作成するまでに必要な時間は，るまでを計算している．平衡や加熱のデータには，JT-60U

データ生成に約３時間，前処理に２１５秒，NVIDIA GeForceの H モード放電３９１１７の解析結果を用いた．最適化 NN 代理

GTX 1080 Ti を用いた NN モデルの訓練に１８１秒となってい

モデルは拡散係数の分布をほぼ完全に再現し，その帰結とる．学習済みの NN 代理モデルを用いた GOTRESS シミュ

して温度分布も極めて良好に再現していることがわかる．

レーションで定常分布が得られるまでに掛かった時間は６４

このようにして GOTRESS データを元に作られた NN

秒であり，TGLF をそのまま用いた状況と比べて５００倍ほ代理モデルは，他の輸送コードでも利用可能である．実際，

ど高速である．

TOPICS に NN 代理モデルの重みデータを渡すことで同様

次に，前述の Hyperas を用いてモデルの決定係数を向上のシミュレーションが可能であることが明らかになっていさせるべく NN モデルの構造とハイパーパラメータの最適

る．すなわち，データ生成が得意な GOTRESS によってモ化を行った．先程のモデルでは隠れ層は３層で固定であっデルの訓練を行い，TOPICS でそのモデルを用いて時間発たが，必要に応じて４層目を付け加えて良いとし，さらに展シミュレーションを行う，という運用も可能となってく

各層のユニット数（離散的に変化）とドロップアウト率

る．そのためにはモデルの汎化性能を向上させるなど工夫（連続的）を可変とした．また，バッチサイズも可変としすべき点が残っているが，原理的には NN 代理モデルを用た．最適化の結果，４層目は不要と判定され，各層のユいた高速シミュレーションは実証されたと言える．

ニット数の平均はおよそ４００，ドロップアウト率は平均で２．
４ おわりに

約 0.09，バッチサイズは４倍の８，
１９２となった．この結果はベイズ最適化の結果として得られたものであるが，この選NN 輸送モデルは欧米で数年先行して導入されたが，国択の意味を解釈することができる．ユニット数の増加は，内でも欧米のモデルとは異なる特徴を生かした開発が進ん用いられた訓練データに比して元のモデルの表現力（複雑さ）が足りていなかったことを意味している．ドロップア
ウト率の減少は，訓練データが十分な量でかつ多様性に富んでいるため，ドロップアウトによって未知のデータに対する汎化性能を上昇させる必要があまり無いと判断したと思われる．また，バッチサイズの増加は，訓練データ内に含まれる相対的外れ値に対する耐性を獲得するためと推定される．NN モデルはバッチサイズでまとめられたサブ
セットデータを入力として流し込んでからネットワークの重みを更新する．もしバッチサイズが小さいと，サブセッ
トデータに含まれた外れ値によってネットワークが影響を受ける割合が高まり，その結果最適でない重み値でモデルが更新されてしまう可能性が高まる．そうなるとなかなか最適なモデルに収束しなくなるため，小さなバッチサイズ
は好ましくない＊＊．構築されたモデルの性能を測る指標として損失関数があるが，それによるとモデルは倍以上性能が向上していることがわかる．実際に熱流束に対する決定係数を算出すると，それぞれ 0.987，0.943 となり，大幅に性能が向上していることがわかる．なお，最適化には約８時間を要したが，２時間を経過した段階でほぼ最適に近い性能が得られていたため，そこで打ち切っても実用上問題はない．また，最適化によって得られたハイパーパラ
メータ情報は，同種のデータセットに対しても同様の性能向上が見込めることが経験的にわかっているため，これらの情報を初めから他のケースのハイパーパラメータとして
図５

設定しておけば，最適化の手間無く良好な性能を持つ NN
モデルを構築することができる．
最適化された NN 代理モデルの再現性を GOTRESS シ

GOTRESS によって行われた，JT-60U H モードプラズマに対する NN 代理モデルと TGLF による（a）
温度と（b）
拡散係数の定常分布の予測．

＊＊：訓練に掛かる計算時間の観点からも，一般に訓練データが大きくなるとバッチサイズは大きくした方が良い．但し，今回の最適化では計算時間を指標としていないため，計算時間を短くするために大きなバッチサイズを選択したわけではない．

７０

Special Topic Article

2. Machine-Learning Assisted Transport Modeling -Practical Cases in JT-60U‐

E. Narita and M. Honda［４］C. Bourdelle et al., Plasma Phys. Control. Fusion 58,
014036 (2016).
［５］O. Meneghini et al., Nucl. Fusion 55, 083008 (2015).
［６］O. Meneghini et al., Phys. Plasmas 21, 060702 (2014).
［７］O. Meneghini et al., Nucl. Fusion 57, 086034 (2017).
［８］J. Citrin et al., Nucl. Fusion 55, 092001 (2015).
［９］K. L. van de Plassche et al., Phys. Plasmas 27, 022310(2020).
［１０］A. Ho et al., Nucl. Fusion 59, 056007 (2019).
［１１］F. Felici et al., Nucl. Fusion 58, 096006 (2018).
［１２］本多 充 他：プラズマ・核融合学会誌 95, 453 (2019).
［１３］E. Narita et al., Nucl. Fusion 59, 106018 (2019).
［１４］E. Fable et al., Plasma Phys. Control. Fusion 52, 015007(2010).
［１５］E. Narita et al., Plasma Phys. Control. Fusion 60, 025027(2018).
［１６］A.G. Peeters et al., Comput. Phys. Commun. 180, 2650
(2009).
［１７］M. Abadi et al., "TensorFlow: Large-Scale Machine
Learning on Heterogeneous Systems", (2015) http://tensorflow.org/［１８］https://keras.io［１９］J. Bergstra et al., in Proc. 30th International Conference
on Machine Learning, Atlanta, Georgia, USA (2013), Vol.
28.
［２０］http://jaberg.github.io/hyperopt/［２１］http://maxpumperla.com/hyperas/［２２］N. Hayashi and JT-60 Team, Phys. Plasmas 17, 056112(2010).
［２３］M. Honda and E. Narita, Phys. Plasmas, 26, 102307 (2019).
［２４］佐竹真介，菅野龍太郎，本多 充：プラズマ・核融合学
会誌 96, 405 (2020).
［２５］F. Imbeaux et al., Nucl. Fusion 55, 123006 (2015).

でいる．機械学習の分野は日進月歩であるため，国際会議を始めとした分野横断の活動で最新の技術・知見を取り入れながら，NN 輸送モデルの改良に留まらずプラズマ分野において広範な応用が進んでいる．本章では NN 輸送モデ
ルの開発に重点を置いた解説を行ったが，実験解析や運転
シナリオ開発などでの活用は実用に近いレベルに達している．そ の 成 果 は 国 際 ト カ マ ク 物 理 活 動（ITPA: International Tokamak Physics Activity）における輸送と閉じ込め
トピカルグループ及び統合運転シナリオトピカルグループ
において継続的に報告されており，ITER 等の将来装置に適用可能なモデルとして位置付けられている．
今後の進展のために不可欠な点として，機械学習に用いるデータの点数と変数領域の拡充が挙げられる．これにより，対象とする実験装置に依らない安定した輸送シミュ
レーションが実現することが期待される．既述のように，
欧米では DIII-D や JET の実験データを取り込んだ統合フレームワークの中で NN 輸送モデルを構築する環境が整っており［７，
１０］，より強力に開発が進むと考えられる．国内
でも，独自の特徴を生かしつつ，例えば ITER プラズマ実
験のモデリングと実験解析に向けた統合フレームワーク
IMAS（Integrated Modelling and Analysis Suite）［１，２５］な
どを利用しながら開発環境を整え，NN 輸送モデルの開発を加速させる必要があるだろう．

参考文献［１］林 伸彦 他：プラズマ・核融合学会誌 95, 423 (2019).
［２］A. Fukuyama et al., Plasma Phys. Control. Fusion 37, 611(1995).
［３］G.M. Staebler et al., Phys. Plasmas 12, 102508 (2005).

７１

J. Plasma Fusion Res. Vol.97, No.2 (2021)７２‐７８

小特集 磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展

３．データ同化による輸送モデリング‐LHD における実践例‐
3. Transport Modeling Applying Data Assimilation Techniques−Practical Cases in LHD−森 下 侑 哉，村 上 定 義，横 山 雅 之１，３），上 野 玄 太２，３，４）
MORISHITA Yuya, MURAKAMI Sadayoshi, YOKOYAMA Masayuki１，３）and UENO Genta２，３，４）
京都大学大学院工学研究科，１）自然科学研究機構 核融合科学研究所，２）情報・システム研究機構 統計数理研究所，
３）

総合研究大学院大学，４）情報・システム研究機構 データサイエンス共同利用基盤施設（原稿受付：２
０
２
０年９月３
０日）

核融合プラズマの挙動を高速かつ高精度に予測・制御できるシステムを実現するため，統合輸送シミュレー
ションコードをベースとしたデータ同化システムを構築している．データ同化とは，シミュレーションと観測を結びつけ，観測によるシミュレーションモデルの最適化を実現する技術である．データ同化は，核融合の分野においても重要な役割を果たすことが期待される．本章では，データ同化の概要と大型ヘリカル装置（LHD）における実践例を紹介するとともに，現状のデータ同化システムの課題と今後の展望について述べる．
Keywords:
data assimilation, ASTI, EnKF, integrated transport simulation, turbulent transport model, LHD

３．
１ はじめにデータ同化とは，観測を用いてシミュレーションモデル

複雑な現象が混在する核融合プラズマの挙動を予測，解を最適化する統計的手法であり，大規模な統合シミュレー

析する有効な手法の一つとして，統合シミュレーション

ションを必要とする気象予報や海洋解析の分野で発達してコードの開発が世界的に進められている［１］．統合シミュ

きた．データ同化は，主に以下のような目的で用いられる．

レーションでは，各物理現象をモジュールに分けて記述

・予測のための初期値，モデルパラメータ推定

し，それらを統合することで核融合プラズマ全体の挙動を・非観測量を含む状態の推定計算する．しかしながら，核融合プラズマの挙動を完全に・観測値の補正，時空間補間再現できるようなモデルは，事実上存在しないと言わざる・観測システムの効率化，感度解析を得ず，時間発展を追うごとに計算結果の持つ不確実性は

地球磁気圏のプラズマ密度分布および電場分布の推定［３］増大する．また，大きな不確実性をもつモジュールがあるや，台風の予測・解析［４］，ゲリラ豪雨予測［５］などさまと，モジュール間で相互作用し，計算結果全体が大きな不ざまな分野でデータ同化が活用され，成果を挙げている．

確実性を持つことになる．京都大学と核融合科学研究所

核融合の分野でも，欧州で開発が進められている統合輸送（NIFS）で開発を進めている統合輸送シミュレーション

コード RAPTOR に比較的簡易なデータ同化手法である拡

コード TASK3D［２］においても，乱流輸送モデルが大きな

張カルマンフィルタが導入されており，安全係数分布の推不確実性を持ち，計算結果に大きく影響する．このような定などに用いられている［６］．

不確実性の問題に対して，従来では，シミュレーションモ

筆者らは，核融合プラズマのリアルタイム高精度予測・

デルの精緻化が行われてきたが，近年では，観測とシミュ制御を目標として，統合輸送コード TASK3D に基づく

レーションを結びつけ，最適な状態を推定するデータ同化

データ同化システム ASTI（Assimilation System for Toroi-と呼ばれる統計的手法が注目を集めている．データ同化を

dal plasma Integrated simulation）の構築を行っている［７］．

用いることで，逐次的に得られる観測情報をもとにシミュ

ASTI は，乱流輸送モデルの推定をはじめとする様々な物

レーションモデルを最適な状態に更新することができ，よ理現象の強力な解析ツールとしての利用も期待される．本り高精度な予測を実現することができる．また，時系列観章では，データ同化手法の概要を説明するとともに，構築

測データからシミュレーションモデルを用いて情報を抽出中の ASTI とその実ショットへの適用例を紹介する．最後

するような場合には，モデルが持つ複数の不確実性を同時に今後の展望について述べる．なお，本章での「輸送モデに考慮に入れた解析が可能となる．

リング」は，輸送係数のモデルに限らず，統合輸送シミュ

Kyoto University, KYOTO 615-8530, Japancorresponding author’s e-mail: morishita.yuya.47r@st.kyoto-u.ac.jp

７２!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

Y. Morishita et al.

3. Transport Modeling Applying Data Assimilation Techniques -Practical Cases in LHD-

レーションモデル全体の構築や最適化を指すことに注意されたい．

３．
２ データ同化３．
２．
１ データ同化の手法
データ同化の手法は，大きく二つに分けられる．一つは，カルマンフィルタ（KF）［８］に代表される，逐次ベイズフィ
ルタに基づく手法［９］である．本研究で用いているアンサ
ンブルカルマンフィルタ（EnKF）［１０］もこちらの手法に属する．逐次ベイズフィルタは，シミュレーションモデルに

逐次ベイズフィルタの概略図（イオン温度 Ti での例）．

図１

よる予測とベイズの定理に基づく最適化（フィルタ）を繰り返して，逐次的に観測データをシミュレーションに同化（図１中の #*'*）が所与であるとすると，式（１）のシステムさせていく手法である．近年のデータ同化の広がりは，逐モデルによって次の観測時点 *"の予測分布（一期先予測分!

!

次ベイズフィルタに基づく手法の導入の簡便さと推定の強

布 #*'*）を得ることができる．観測時点 *"において，一期

力さ，計算機性能の向上によるところが大きいと考えられ

先予測分布を事前分布とし，式（２）の観測モデルを用いてる．本章では，３．
２．
２，３．
２．
３節において，その計算を簡単計算される尤度分布を用いると，ベイズの定理が適用できに説明する．

る．こうして状態ベクトルの確率分布は，観測の情報（$
*）

"

!

"

もう一方のデータ同化手法は，４次元変分法（アジョイを取り込んだ事後分布（フィルタ分布 #*'*）へと更新されント法）と呼ばれる手法［１１］である．４次元変分法は，ある．この最適化操作がフィルタであり，この事後分布を用る時間区間（同化ウィンドウ）に含まれる観測データに基いてさらに予測とフィルタを繰り返す（逐次ベイズフィル

づき，状態の同時確率を最大化する初期値セット（事後確タ）．この操作によって，シミュレーションモデルが真の率最大解）を反復計算によって求める手法である．大規模

時間発展へと近づくことが期待される．次節では，逐次ベなシミュレーションに強く，多くの気象予報システムで採

イズフィルタに基づくデータ同化手法の一つであるアンサ

用されている手法である．本章では詳細な説明は行わな

ンブルカルマンフィルタについて説明する．

い．詳細は，［１１，１２］等を参照されたい．

３．
２．
３ アンサンブルカルマンフィルタ""

システムモデルの '*と観測モデルの (*がともに線形で，

３．
２．
２ 状態空間モデルと逐次ベイズフィルタ

確率分布としてガウス分布を仮定すると，状態空間モデル

データ同化では，シミュレーションと観測を結びつけ，は以下のようになる．

シミュレーションモデルの最適化ができるようにするため，状態空間モデルを考える．状態空間モデルでは，状態!!
#*#"*#*!!"!*!!*$$%
%*&
$
!
#
"
!
&
!
##
""
$$%
*
**
*
*
*&

を表す変数を成分とするベクトル（状態ベクトル）により，系のある時刻における状態を表現する．状態空間モデル

（３）
（４）

は，以下のシステムモデル（１）と観測モデル（２）によって構

システムモデルと観測モデルが線形であることから，式成される．

（１），（２）における '*，(*を行列 "*，#*で書いている．ま

#*#'
#*!!!
!*&
*%

（１）

$
#*!
"*&
*#(
*%

（２）

た，!$$%#!
"&は !が平均 #，分散共分散行列 "のガウス分布に従うことを意味する．この状態空間モデルの仮定
における逐次ベイズフィルタはカルマンフィルタ（KF）とここで，#*は時刻 *における状態ベクトルであり，$*は時呼ばれる．この仮定の下では，予測分布とフィルタ分布は

刻 *における観測ベクトルである．システムモデル（１）共にガウス分布に保たれる．よって，確率分布の時間発展

は， '*によって 状態 ベ ク ト ル を 時 間 発 展 さ せ るとともやフィルタは，平均と分散共分散行列のみで議論することに，時間発展に伴うノイズ（システムノイズ）!*の付与をができる．

行う．すなわち，'*はシミュレーションモデルに相当する．

カルマンフィルタは線形なシステムに対するデータ同化

観測モデル（２）は，(*によって状態空間と観測空間の間の変手法であるが，一般にデータ同化を必要とするような複雑換を行い，観測ノイズ "*を付与する．この観測ノイズになシステムは，非線形なシステムである．非線形なシステ

は，観測機器に起因する測定誤差の他にシミュレーション

ムでも実行できるデータ同化手法として，アンサンブルカ

モデルで表現しきれない（興味のない）誤差分も含まれる．

ルマンフィルタ（EnKF）［１０］
や粒子フィルタ（PF）［９，１３］

この二つのモデルを連立することで，観測データからの情がある．これらの手法では，確率分布をアンサンブル近似

報をシミュレーションモデルに取り込むことができ，状態

し，各アンサンブルメンバーを以下の様に個々に時間発展ベクトルを最適化する土台ができる．

させることで，確率分布の時間発展を追う．

図１は逐次ベイズフィルタの概略図である．一般的なシ

%
&
%
&
%
&
)
)
)
!*
#*
#'#
!*!!'*!!!
"'!!
*

ミュレーションが一つの状態ベクトルの時間発展を追うの（５）に対して，データ同化では，状態空間における確率分布の!!までの観測情報を同化した，時
ここで，#*'*!!は，時刻 *時間発展を追う．ある時刻 *! の状態ベクトルの確率分布

刻 *の状態ベクトルを表す．また，%)&はインデックス )７３

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

３．
２．
４ アンサンブルカルマンスムーザの ア ン サ ン ブ ル メ ン バ ー で あ る こ と を 表 す．式
（５）は
$
%
'
で決まる初期値を物理モデルで時間発展させれば
")
!!&
!!
)
よいことを示しており，手を加えることなく EnKF や PF

ため，フィルタ時点よりも未来に観測された情報は，当然の枠組みに物理モデルを取り込むことができる．

ながら，その推定には影響しない．予測を行う場面では，

EnKF のフィルタ操作は，KF に基づいて，過去の情報しか用いることはできないが，解析的な場面で

$
%
$%
$
%
$
%
'
'
'
#)$
#
")
#")&'
""
%
)"!) !!)"
&
!!
&
!!
)
)
)
)

（６）

!!
# !"$
#)%
# !""##)#$"&!! ) !)$)
&
!! )
)
)
)

（７）

フィルタによる推定は，過去の観測のみに基づいている

フィルタ時点よりも未来の観測情報も含めた時空間的に一貫性のある推定を行いたい場合には，ある時点の観測情報をそれよりも過去のフィルタ推定に反映させる必要があ#はアンサンブルから計算さで表される［９，１０］．ここで，$る．この時間を遡る最適化操作は平滑化（スムーザ）と呼れる分散共分散行列である．図２に，EnKF（ASTI）の計

ばれる．EnKF の場合には，フィルタの際にアンサンブル

算の流れの概略図を示す．EnKF は，KF のアルゴリズムをを再構成する行列 &を保存しておくことで，平滑化を容易

非線形なシステムに対応させた手法であるため，システムに実現できる（アンサンブルカルマンスムーザ，EnKS）．

の非線形性が強く変数間の関係を線形の範囲で捉えきれな

EnKF におけるフィルタ分布のアンサブルメンバーは，一い場合（分布がガウス分布で近似できない場合）や観測モ

期先予測分布のアンサンブルメンバーの線形結合で表現でデルが非線形な場合には，必ずしも推定がうまくいくとはきる．すなわち，各アンサンブルメンバーを列に持つ行列限らない．こうした状況には，より一般的（非線形・非ガ%を用いると，フィルタ操作は，行列 &を用いて，ウス）なデ ー タ 同 化 手 法 で あ る PF が 用 い ら れ る．PF%)&)#%)&)!!&)

は，尤度関数によるアンサンブルメンバーの重み付けに（８）

よって，フィルタ分布のアンサンブルを推定する．どんなと書ける．ある時点 (のフィルタアンサンブルに，次時刻モデルにも対応できる柔軟性を持つ反面，分布の情報とし("!の観測を反映させた分布，すなわち，時刻 (におけるてアンサンブ ル の２次 の モ ー メ ン ト ま で し か 用 いない"!までの観測に基づく平滑化アンサンブルは，時刻
時刻 (EnKF に対して，より大きなアンサンブルメンバー数を必("!の行列 &を用いて以下のように表せる［９，１０］．

要とする．本章で紹介する適用例では，リアルタイムでの%(&("!#%(&(&("!

予測を念頭において，より計算コストの低い EnKF を用い
ている．

（９）

この式から平滑化アンサンブルに関する漸化式，

データ同化の大きな利点として，観測のない状態変数を%(&)"!#%(&)&)"!

も最適化できる点が挙げられる．EnKF において，この能（１０）

力を担うのは，変数間の相関関係である．個々のアンサンを導くことができる．ただし，)!(である．つまり，フィブルメンバーをシステムモデルを用いて時間発展させるこ

ルタ毎に得られる行列 &さえ保持していれば，任意の時刻とにより，変数間の相関関係が自動的に予測アンサンブル

までの観測を用いた平滑化分布の推定が可能である．

の共分散に現れる．観測量がある変数とない変数は，この３．
３ 統合輸送コード TASK３D を用いたデータ同化システム

共分散による相関関係でつながれており，観測のある状態変数の最適化による修正が観測のない変数まで伝搬する．

本節では，統合輸送シミュレーションコード TASK3D

この仕組みによって，シミュレーションモデルの内部までに基づくデータ同化システム ASTI［７］について説明す

最適化が行き届く．

る．ASTI は EnKF による核融合プラズマの予測・解析を目的としたコードであり，Fortran で実装され て い る．
Fortran によるシミュレーションコード（TASK3D）の時間発展計算を ASTI から呼び出せるようにモジュール化す
ることで，既存のシミュレーションコードをシステムモデ
ルとしてデータ同化（状態空間モデル）の枠組みに取り込
むことができる．本節では，状態ベクトルに温度，密度，乱流モデル係数，中性粒子ビーム入射（NBI）加熱分布をとった場合の EnKF による推定の例を示す．図２は ASTIの概略図である．
３．
３．
１ 状態空間モデルと観測データ
EnKF のシステムモデルとして用いている TASK3D は，熱と粒子の一次元輸送方程式をベースとする統合輸送シ
ミュレーションコードである［２］．現在，乱流熱拡散係数を与えるモデルとして，電子に gyro-Bohm モデル，イオンに gyro-Bohm gradT モデルを用いている．
図２

ASTI の概略図．

７４

Special Topic Article

' +
-*
$.
$".! 0"! ."-Y. Morishita et al.

3. Transport Modeling Applying Data Assimilation Techniques -Practical Cases in LHD-

/
!

' +
&+
-*
$0
$"0! 0"! 0"! 0-"
- /
!
+0

EnKF がうまく機能することを確認した．
まず，システムノイズは，標準偏差（)4の対角成分の平（１１）方根）をその時点の状態変数の平均値に比例する形で推定した．この比例定数の組み合わせは対数尤度に基づいて決

（１２）

定した
［７］．観測列 &4)4 に対する EnKF のパラメータ &"%ここで，!は磁場強度，'0はイオンの Larmor 半径，-はプの対数尤度 2&&)'4)4 (は，以下の式で近似することができ

ラズマ小半径である．".，"0は実験結果に基づいて調整さる［１９］．

"

れる係数である［１４］．乱流モデルの係数を実験に合わせて

4
%

&
2
&)
'
&
&
&
"
2/3'
$' 1
4)4 (
4)
4)4 (

調整することで，温度分布の再現が比較的よくできている"［１５］．NBI 加熱による加熱分布は GNET-TD コード
［１６］に

%

%'
'(
#410 $

（１３）

1!"

（１７）
&

'(
2/&
1
2/%' #410&!% 1

4
"
1 $4

状態ベクトルを以下の TASK3D 内の変数で構成する．

"

1

4
"
1 $4
4
%

よって予め計算されたものを用いている．
%$'
".-"
"0-"
$.
$0
!.-"
!0-"
#.-"
#0-(
"
"

%

0$"

"
.
*4)
)
#
%(
('
1

'
( - !"'(
0
0& !#%4
&
".43#
!"'
(*4 '
(
$
4 !#%
)
4)
4
4
# 4

添字 /は電子，0はイオンを表す．それぞれの変数は，小半

1

径方向１次元の空間グリッド点上で定義される．本研究で

1

1!"

1

1

1!"

（１８）

は，空間を６０グリッドで分割しているため，各々の変数は，

ここで，& はアンサンブルメンバー数，.は観測ベクトル

以下のような構造を持つ．
".$'
+.""
+.#"
+.&!(
*"

1

&の次元（本設定では２４０）である．)4を対角行列（分散の（１４）

み，共分散０）で仮定した場合，ノイズを加えた後のアン

上付きの数字は，対応するグリッド点を表す．ここで，"サンブルメンバーは，空間的に不連続になる．大きな不連は状態ベクトルの温度部分，$は密度部分である．! は乱続があると，TASK3D で計算ができな く な る．そ の た流モデル係数の部分であり，状態ベクトルに含めるにあた

め，システムノイズは空間的に滑らかになるように加え方

り，空間分布を持たせている．また，#は，GNET-TD によを工夫する必要がある．本研究では，グリッドにおけるノ

るNBI加熱分布の計算結果に不確実性を導入するための因イズが近隣のグリッドに指数関数的に減衰するように，分子であり，これにも空間分布を許す．GNET-TD により計

散共分散行列の共分散成分を決定した［９］．

算された加熱分布 (,*+は0(,#*$10(,0*++

一方の観測ノイズは，フィルタ毎に，標準偏差（*4の対角成分の平方根）を予測誤差（観測値−予測値）に比例す（１５）る 形 で 推 定 し た．比 例 定 数 は，尤 度 を 最 大 に す る 値

として，TASK3D に読み込まれる．これらの状態変数の初（0.8）を採用した．この仮定は観測値に応じた適切な観測

期値は，温度と密度に関しては初期観測値を用い，"，1

ノイズ R を設定していることに相当する．予測誤差が大きに関しては，従来のシミュレーションで仮定していた値を用いときには観測ノイズを大きくし，過適合にならないよういる．今回の場合では，".$#$!
!，"0$(!
!'，1.$10$"!
!

に働く．予測誤差が大きい場面でフィルターによる最適化である．

が緩和され，結果としてフィルタ後の背景誤差 ,4)4が大き

ここに，温度と密度の観測データを同化する．本研究でく評価される．これによってより柔軟な空間分布形状が調は，空間的に補間したデータを用いるため，径方向６０グ

べられる．また，予測がうまく行っている状態では，背景リッド上に対応する観測値をそのまま同化することができ誤差を広げないように働き，最適な状態に保つ．

る．状態ベクトルは４８０次元（６０グリッド×８変数），観測

３．
３．
３ 実ショットへの適用例
ここでは，LHD における NBI ショット（１
１４０５
３）を実際

ベクトルは２４０次元（６０グリッド×４変数）である．
観測モデルに関しては，本研究における状態ベクトルのに TASK3D の シ ミ ュ レ ー シ ョ ン に 同 化 し た 例［７］を 示観測有り部分と観測ベクトルとが一対一に対応するため，

す．密度変化は緩やかであるため，TASK3D では熱輸送の以下の観測行列を導入する．

みを解くこととし，密度はフィルタによる最適化のみで時$#4$'#%!"#%! '#%!"#%!(

間変化を追った．図３は，データ同化を行わない TASK3D（１６）による従来型のシミュレーションの結果とデータ同化を$は単位行列，' は零行列である．

80 ms おきに行った場合の電子温度とイオン温度の結果で

３．
３．
２ ノイズ

ある．なお，アンサンブルメンバー数は &$#!!!とした．

EnKF の性能を左右するハイパーパラメータとして，シデータ同化を行わない単純なシミュレーションの結果に対ステムノイズと観測ノイズの分散共分散行列がある．これ

して，データ同化による予測値（塗りつぶし点）の方が，観らの推定には，観測との誤差から統計的に推定する方法

測値の時空間分布を高精度で再現しているのがわかる．今［１７］やベイズ推定に基づく方法［１８］など，状況に応じた

回は同化周期 80 ms の 結 果 の み を 示 す が，同 化 周 期 を色々な方法が提案されている．本研究では，以下に述べる

240 ms 程度まで大きくしても，予測精度は大きく下がらなようにノイズの分散共分散行列を調整していくことで，いことを確認している．温度の時間発展を主に乱流モデル
７５

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

図３

EnKF による温度の予測結果（上段：TASK3D のみ，下段：80 ms おきにデータ同化）．規格化小半径 "= 0.0 と 0.6 について示す．下の棒グラフは，"= 0.0 における背景誤差と観測ノイズの標準偏差を表す．また，イオン温度の灰色のハッチは，観測データが無いた
め，代わりに電子温度を同化させた部分である．

係数で調整していることを考えると，本ショットにおいて，! は 240 ms 程度，あるいは，それよりも遅い時間ス
ケールで変化していると考えられる．
図４は，EnKF によって推定された乱流モデル係数 !
の時空間分布（上段）と EnKS によって推定された時空間分布（下段）である．グラフの値は，推定値を従来仮定していた係数の値で割った比であり，１が従来仮定していた値に相当する．図４の EnKF による結果を見ると，!$，!%
ともに，空間的には大きく変化しているが，時間的には緩やかに変化しているのがわかる．これが，EnKF で同化周期をある程度大きくとっても予測がうまくいった理由である．ま た，図４の EnKS に よ る 結 果 を 見 る と，!$の"!!
!
#，1.2 s 付近におけるピークや，!%の "!!
!
"，1.4 s付近におけるピークが少し過去方向に伸びており，未来の情報が過去のフィルタ分布を修正しているのがわかる．!

図４と "の平滑化時空間分布を乱流熱拡散係数のモデルに取り

乱流熱拡散モデル係数の推定分布．

込んだ TASK3D によるシミュレーション（データ同化を度，密度の観測データを同化した結果を紹介した．これに行わない）結果を図５に示す．シミュレーション結果と観より，温度の時間発展の高精度な再現とそれを実現する乱測値が良い一致を示していることから，この EnKS による流モデル定数が得られた．

推定は妥当であると言える．

リアルタイム予測に向けて，現在の ASTI が持つ大きな課題は，計算時間である．EnKF の計算自体は，B-spline

３．
４ まとめと今後の課題及び展望

補間等を利用して，同化する径方向のグリッド点数を５〜

本章では，LHD におけるデータ同化の取り組みを紹介し１０に絞ることで，状態ベクトルの次元を節約することがでた．核融合プラズマの高精度な予測・解析を実現するたきるため（多くとも５０次元程度），現在の計算では問題にめ，データ同化システム ASTI を構築中である［７］．そのならない．計算時間のほとんどは，TASK3D 自体によるも適用例として，温度，密度，乱流モデル係数，加熱分布をのであり，タイムステップ１ms，径方向６０グリッドで計算

状態ベクトルにとり，TASK3D のシミュレーションに温する場合，1 s 分の時間発展を計算するのに約１０分かか
７６

Special Topic Article

3. Transport Modeling Applying Data Assimilation Techniques -Practical Cases in LHD-

図５

Y. Morishita et al.

EnKS による推定結果を用いた TASK3D によるシミュレーション結果．

る．そのため，TASK3D 自体に効率的な並列化を施すととイオン温度の結果である．式（１５）の加熱分布と同様に，因もに，データ同化によるサポートを前提とした，モジュー

子 "（空間分布を持つことを仮定）によって，回帰モデル

ル（モデル）の簡約化を行う必要がある．こうしたモデルに よ る 熱 拡 散 係 数 に 不 確 実 性 を 導 入 し，こ の "を 式の簡約化の過程にもデータ同化を用いることができる．

（１３）の乱流モデル係数 ! の代わりに状態ベクトルに含め

データ同化によって推定されたモデルのパラメータ等は，た．"の空間分布が一様に１の場合が回帰モデルそのもの

現状のモデルに欠けている部分の情報を持っている．そのに対応する．ショットは３．
３．
３節と同様に１
１４０５３を用いてため，データ同化によって推定された値を上手くモデルにいる．このグラフを見ると，最初から 0.4 s にかけて予測誤取り込むことができれば，より高精度なモデルを作り出す差が大きくなっており，それと同時に背景誤差の標準偏差ことができる．モデルの簡約化とモデルの構築を繰り返すも大きくなっているのがわかる．一方で 0.4 s を超えたあた

ことで，目的の精度，速度に向けてモデルを鍛錬することりから誤差は減少し始め，0.8 s 以降はかなり精度が良くができる．また，時間刻みや空間グリッドを粗くすることなっている．0.8 s 以降から "の空間分布が１でほぼ一様でに伴う不確実性の増大（精度の悪化）もデータ同化によっあったことから，このモデルは，温度上昇部（0.0 s−0.4 s）て補うことができる．現在，こうした簡約化（粗視化）のでは，モデルの予測能力は低いが，安定部（0.6 s）以降で方針に従って，ASTI の高速化を進めている．

は，予測能力が高いことがわかる．こうした性質は，回帰

本章で紹介した例の様に，１ショットに関してデータ同モデルが対象としたデータセットに依存すると考えられ化を行った場合は，その推定結果は，そのショット専用のる．一般的なシミュレーションによる検証では，序盤の予ものとなる．そのため，普遍的な物理や構造を１ショット

測に引っ張られて，こうした性質を把握することは難しだけから抽出することはデータ同化では難しい．しかし，い．データ同化を用いることで，モデルの能力や妥当性の状態ベクトルを複数ショットについて拡張する（１つの状

検証を時空間的に区切って行うことができる．この点もモ態ベクトルが，複数のショットの情報を持つ様にする）と

デルの精度向上に有用である．

複数のショットに共通する性質や構造を抽出できるように

ASTI を実用的なものにするためには，計算の高速化となる．また，シミュレーションモデルを複数のショットに多ショット解析への拡張が必要であり，まだまだ課題が多同時に当てはめたり，複数のショットで同時に評価するこい状況である．しかしながら，統合コードは，データ同化とができる．こうした多ショットに対する同時解析もデーによってより一層，強力な予測・制御システムかつ解析

タ同化に期待されることの一つであり，早期実現のためにツールになり得る．データ同化は，核融合炉の実現におい

ASTI の開発を進めている．

て要となる技術であると考えている．

上記の他に，データ同化はモデルの能力，妥当性の詳細

謝な検証に用いることができる．図６は，LHD における熱拡

辞

本章に記載した研究進展について，データサイエンス共散係数の回帰モデル［２０］を用いて，データ同化した場合の同利用基盤施設データ同化研究支援センターにおける技術
指導，統計数理研究所共同研究 2019-ISMCRP-2027，2020ISMCRP-2026，核融合科学研究所共同研究 NIFS20KLPT
007，および，原型炉研究開発共同研究「核融合の大規模データを活用するデータ駆動型モデリング手法の研究」による支援に感謝いたします．

参考文献

図６

［１］林 伸彦 他：プラズマ・核融合学会誌 95, 423 (2019).
［２］S. Murakami et al., Plasma Phys. Control. Fusion 57,
119601 (2015).
［３］S. Nakano et al., J. Geophysical Res. : Space Phys. 113, 11(2008).
［４］T. Honda et al., Mon. Wea. Rev. 146, 213 (2018).

回帰モデルを用いたイオン温度のデータ同化結果．

７７

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

［５］T. Miyoshi et al., Bull. Amer. Meteor. Soc. 97, 1347 (2016).
［６］M.C.C. Messmer et al., Plasma Phys. Control. Fusion 61,
035011 (2019).
［７］Y. Morishita et al., Nucl. Fusion 60, 056001 (2020).
［８］R.E. Kalman, J. Basic Engineering 82, 1 (1960).
［９］樋口知之 他：データ同化入門−次世代のシミュレー
ション技術−（朝倉書店，２
０
１
１）
．
［１
０］G. Evensen, Ocean Dyn. 53, 343 (2003).
［１１］O. Talagrand and P. Courtier, Q.J.R. Meteorol. Soc. 113,478 (1987).
［１２］淡路敏之 他：データ同化 観測・実験とモデルを融合するイノベーション（京都大学学術出版会，２
０
０
９）
．

［１
３］上野玄太：統計数理 67, 241 (2019).
［１４］A. Wakasa et al., Proc. 39th EPS Conf. and 16th Int. Conf.
Plasma Physics P2.028 (2012)［１５］H. Yamaguchi et al., JPS Conf. Proc. 1, 015045 (2014).
［１６］H. Yamaguchi et al., Plasma Fusion Res. 8, 2403099 (2013).
［１７］G. Desroziers et al., Q.J.R. Meteorol. Soc. 131, 3385 (2005).
［１８］G. Ueno and N. Nakamura, Q.J.R. Meteorol. Soc. 142, 2055
(2016).
００５）．
［１９］北川源四郎：時系列解析入門（岩波書店，２［２０］M. Yokoyama and H. Yamaguchi, Nucl. Fusion 60, 106024(2020).

７８

J. Plasma Fusion Res. Vol.97, No.2 (2021)７
９‐８５

小特集 磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展

４．データ駆動アプローチを用いた動的乱流現象の解析
4. Data-Driven Analyses of Dynamical Turbulence Phenomena
佐 々 木 真１，２），河 原 吉 伸３，４），草 場
１，
２）

SASAKI Makoto

彰１）

３，
４）

, KAWAHARA Yoshinobu

and KUSABA Akira１）

１）

九州大学応用力学研究所，２）九州大学極限プラズマ連携研究センター，

３）

九州大学マス・フォア・インダストリ研究所，４）理化学研究所 革新知能統合研究センター（原稿受付：２
０
２
０年１
０月１
４日）

本章ではデータ駆動科学的手法として近年注目されている動的モード分解について，その理論的背景からプ
ラズマ乱流への応用例まで紹介する．動的モード分解は流体現象の時空間構造を観測／計測データから抽出するための数値解析手法の一つとして提案されたものであり，現在では様々なシステムへの応用がなされている．動的モード分解を用いると，時空間発展データから統計的に重要な空間構造とその構造に対応する複素周波数を得ることができ，より少数の自由度でシステムを理解することが可能となる．ここでは，（１）物理機構の理解をめざす例として，動的モード分解をリミットサイクル振動を伴う乱流に適用し，構造間の因果関係を抽出する試みについて紹介する．また，（２）乱流の時空間発展の予測に応用する試みについても紹介する．
Keywords:
plasma turbulence, abrupt phenomena, mode decomposition, DMD, causality, prediction

４．
１ 背景

観測データに基づき，システムの時空間発展を特徴付ける

磁場閉じ込めプラズマや宇宙のプラズマ現象はエネル

ことができる新たな方法である［５］．システムの時間発展

ギーの流入出がある非平衡開放系であり，自律的なリズムを考慮したモード分解であるため，系の動的挙動の理解やや自己組織的な空間パターンの形成が起こるため，システ

予測に適している．近年，流体力学分野をはじめとして，

ムが本質的に動的である．特に，極限的に空間不均一な系

様々な分野における時空間ダイナミクスに適用され，そのである磁場閉じ込めプラズマでは，プラズマ装置に対し有用性が認識されてきている．プラズマ研究への適用も始て，ミクロスケールの乱流が発達し，その非線形ダイナミまっており，例えば［６‐１０］などがある．ここでは，動的クスによってプラズマ閉じ込めが規定されている．ミクロ

モード分解の基礎を概観し，プラズマ乱流研究への応用例

スケール乱流が駆動する装置スケールのマクロ構造や中間を紹介し，物理的理解に向けた試みと定量的予測の試みにスケールのメゾスケール構造の重要性が認識されているついて紹介する．

［１］．そのため，プラズマの動的挙動を理解するためには，

４．
２ 動的モード分解

乱流とその非線形性によって駆動される構造間の相互作用を理解することが必要である．

動的モード分解は，当初，流体現象の時空間構造を観

磁場閉じ込めプラズマにおける重要な動的挙動には周辺

測／計測データから抽出するための数値解析手法の一つと局在モード（edge localized mode，ELM）［２］や閉じ込め遷して提案されたものである［１１，１２］．しかしその後，応用数移［３］等の自発的なものから，燃料供給時の外部摂動に対学分野などで議論される力学系のクープマン解析との関連する応答［４］などがある．これらの挙動の物理機構の理解も知られるようになり，流体分野に限らず広く科学・工学

やその予測は重要な問題である．これらの物理機構の解明

分野へ応用されるようになっている．

をめざした研究が理論・シミュレーション・実験の側面か流体分野における類似したモード分解としては，固有直ら精力的に進められており，理解が大きく進んできている交分解（POD）も古くからよく用いられている［１３］．固有

ものの，これらの現象の定量的な予測や制御には至ってい直交分解では，空間方向を一つの列としてまとめた行列ないのが現状である．近年，データ駆動科学的手法が急速（つまり行方向は時間に対応）に対して特異値分解を計算に発達してきており，この新たな手法を第４の科学的手法することで，時間発展に寄与する特徴的な空間構造の抽出として利用する方向性の研究も必要となっている．

を行う．統計的にはこの計算は，時間方向の相関を無視す
ればいわゆる主成分分析と同様であり，変動の大きな空間

本章で焦点を当てる動的モード分解は，データ駆動的にResearch Institute for Applied Mechanics Kyushu University, Kasuga, FUKUOKA 816-8580,Japancorresponding author’s e-mail: sasaki@riam.kyushu-u.ac.jp

７９!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

構造，つまり物理的な解釈としてはエネルギーが大きい空間構造を取り出していることになる．そのように，固有直交分解は時間方向の周期性を陽に抽出している訳ではな
く，多くの場合に特徴的な周期性を持つ空間構造の抽出が困難であることが知られていた．しかし後に見るように，
動的モード分解は，時間方向の関係を陽に用いたモード分解を計算する方法であり対応する周波数（とその減衰・発展）も同時に得られる．

図１

クープマン作用素を用いた力学系の解析．

ここでは，まず動的モード分解の理論的背景について概説し，そしてその代表的な具体的手順を示しつつ，動的#'%%&'"*"'%&"*%+&%&-%&"*"" #&#"!$!"% &"!

モード分解の解釈について説明する．
４．
２．
１ クープマン解析（１）'個

ここでは動的モード分解の背景として，力学系のクープマン解析について概説する．動的モード分解は，クープマ

ただし，#&（$$）は固有関数 %&に対応する固有値を表す．

ン解析へのデータを用いた数値解析手法の一つとしても位

式（１）中の各コンポーネント中で時間発展に寄与する部分置付けられる．

は#&のみであるため，観測量の時間発展が#&で決まる固有

まず，状態空間"（一般には実ベクトル空間の部分集合）の周期性（#&の位相）と減衰率（#&の絶対値）を持つダイ$における（決定的な）力学系 ,（,$"）を考え
*
,&
*
&
'
"$
&
'
'

ナミクスに分解されていることがわかる．

る１．$: "→"は，系の軌跡を定める写像である．ここではクープマン作用素のスペクトル分解を用いた非線形力学

データを用いた解析を考えるため，次の離散時間での表現系の解析に関する応用として重要なものの一つに，縮約理も適宜用いる．

論との関係があげられる．特に，リミットサイクル振動子などの周期的な運動に近い挙動を示す系で用いられる位相

&
'!!'
"*

,
,
"
",
'!!"&
''
'!#
"*&

'
"*

$
,&
$
#$
&
'
'

縮約は，クープマン作用素の固有関数と直接的な関係を持つ［１７］．またもう一つの重要な応用としては，力学系から

ただし，,
'
"*
&
'であり，&"*: "→"は対応するフロー' ",得られる観測が流体現象のように空間的に広がりを持つ場である．これらの表現は状態の軌跡を記述することでダイ

合には，分解（１）は，時空間的なコヒーレンス構造の抽出とナミクスを表している．

して有用であることが知られている．つまり，式（１）の各一方で，ダイナミクスのもう一つの表現として観測量のコンポーネント中の固有値以外の部分，つまり +&%&にあた

時間発展を表す方法として，クープマン作用素を用いた記る部分は，各固有値で表される時間的変動に関する %への述が知られている［１５］．状態空間"上に定義されたなんら寄与として解釈できる．各観測量が各空間上の点におけるかの観測空間!内の関数 ': "→%（$!）（観測量）を考えなんらかの量（例えば，流体現象であれば流体速度や渦度たとき，クープマン作用素は次式により定義される［１６］など）に対応している場合には，この寄与の相対的な違い（図１も参照）．

を見ることで，各時間変動に対応する空間的な同調性を知ることができる．流体現象の解析においては，このような#%"%%&"*+#%&
,
,
"%&
''
'!!'

空間的なコヒーレンス構造（組織構造）を知ることは重要

ただし %は関数の合成である．定義から#の線形性は容易な課題として従来から盛んに研究されており，これが動的に示される．今，作用素#が点スペクトルのみを持ち固有

モード分解への注目のきっかけにもなっている．

値分解可能であるとし，さらに観測量 'は作用素#の固有

４．
２．
２ 基本的な手順とその意味

関 数 %&: "→$を 用 い て 一 次 展 開 で き る，つ ま り

ここでは，一般によく用いられる動的モード分解（厳密

'"! #
&%
&の よ う に 表 す こ と が で き る と 仮 定 す る
&"!+

DMD
［１８］）を取り上げ，その手 順 を 説 明 し つ つ，動 的（+&$$は展開係数）．このとき，観測量 'の時間発展は次モード分解の性質について述べる．
動的モード分解は，未知の系に関するクープマン分解式のように展開することができる（クープマン・モード分（１）を，その系から得られたデータから近似的に推定する解とも呼ばれる）．

手順である．この際に，データは，状態変数を入力とする( 個の観測関数 '&: "→$（&"!
!
(）を通して得られて,!
１ 決定的でない場合は，例えば，%$"#（"）
を観測量としてとれば，それは確率的な力学系を扱うことに相当して，任意の$,に対して次のようにクープマン作用素を定義することができる
［１
４］
．
"
'
"#)"*&
#,
(
,
,
%&
,
%&
#'!!'
#' ",
,
#%&
''
'!!(
''
''
'!!"!)'*
ただし，)"*&"*における状態 ,,,'!!'"*における状' が与えられた際の時刻 &'!!(''は遷移状態確率である．つまり，ある時刻 '態の期待値を与える．

８０

Special Topic Article

4. Data-Driven Analyses of Dynamical Turbulence Phenomena

M. Sasaki et al.

いると仮定する．そのように，まずデータとしては，解析

コヒーレンス構造の情報を与えている．流体現象における，##!'したい未知の系 $"/に関して 0
0
0(として
##$
"(
"/'

適用では，こういった情報は，データから空間的な組織構生成された時系列データが与えられたとする造を抽出する目的などに利用できる［１３，１５，２４］．この動的T）．例えば，この系から生成したと想定さ（!: #*%"!
%, +,!

モードに関するコヒーレンスという視点は応用的にも重要

#
#
#
,!
れる有限長%の時系列データ #!!
"!
#!
% が与えられるであり，流体のみでなく，脳科学や感染症伝播，経済学な#場合には，（#（/#"!
%）のようにペアを構成し
,!
/）
!"!
/

ど，多くの分野でこの特徴を用いた解析が研究されているて用いる．このとき，厳密 DMD の手順は次のようになる．

［２５‐２８］．また，このコヒーレンス構造をデータの特徴として予測（教師あり／教師なし学習）に用いるための枠組と(#:#*
１．２つの行列(":#*
#
#
#
#
#
#
,!
,!
!!
"!
"!
#!
%!"+
%+

みも提案されており［２９，３０］，いくつかの応用場面においを構成する．

ても有用であることが確認されている［３１，３２］．

２．("の（打ち切り）特異値分解 ("%&$'H を計算する（!H はエルミート転置，採用ランクを .とする）．
/: #&H(#'$!"の固有値 #.
３．行列 !
+と固有ベクト ル

４．
３ プラズマ乱流の時空間ダイナミクス解析への応用/+を計算する（+#"
"
!
!
.）．
#
,!

４．
３．
１ 乱流構造間の因果関係同定への応用.
.
/+（動的モード）と対応する固有値 #
４．"
$!""
#'
+#(
+

乱流と背景場との動的な相互作用の理解に向けた試みとを得る．

して，動的モードを利用した構造間因果関係抽出方法の紹

この手順からわかるように，動的モード分解は比較的容易に実介をする．動的な相互作用の例として，大域的乱流シミュ装が可能である．なおこの手順は，)(*!&" )(#!!(")".
の最小二乗解 ! の固有値と固有ベクトルを計算している

レーションで得られたリミットサイクル振動に動的モード,",

分解を適用した成果について述べる［６］．
直線磁化プラズマにおける長谷川・若谷方程式の大域的ことに相当するが，主要な固有値のみが得られれば十分で.の固有値分解を計算する代わりに，これあるので，直接 !

乱流シミュレーションで得られた時系列データを解析に使/#&H !& の固有値分解を計
を & を用いて射影した行列 !

用した．本ミュレーションでは，密度・静電ポテンシャ

算している．なお，& は直交行列であるので，これらの固ル・平行速度場の空間３次元的分布の時間発展が出力さ有値分解で得られる最初の .個の固有値は同じものであれ，乱流データは１ケースおおよそ数十ギガバイト程度とる．一般に，DMD は高次元なデータに対して適用する必.も高次元になり固有値計算な
要があり２，そのため行列 !

なる．基礎プラズマ実験でしばしば観測されてる強い周方

どが困難になるため，上記のような計算がよく用いられる

流れの空間不均一性によって駆動されるケルビンヘルムホ

訳である．

ルツ不安定性の非線形発展を得た．渦度ソースの強度を増向流れ［３３］を渦度ソースを導入することで模擬し，周方向そして計算された動的モードとその固有値を用いて，例

加させると，定常的な非線形飽和状態からリミットサイク

えば時系列データの場合には，各サンプルは次のように分ル振動を伴うダイナミックな状態へ変化した［３４］
．ここで

解される．

は，リミットサイクル振動を伴うケースに動的モード分解.

を適用した結果に焦点を当てる．図２に乱流シミュレー

!"
/
.
.
#
+ 1
/%! #
+"
+

（２）

ションで得られた非線形飽和状態でのフーリエモード+#",!
-(の時間発展（,，-はそれぞれ周方向・軸方向モー
'

.
.: #*
T は # #'.
.
.
）
ただし，$:#*
（'
"
"
"
1
1
1
,!
,!
$
.+
.+
"!
#!
"!
#!
!

ド数）を示す．ここで背景場は'!
!
!(に対応している．軸方を満たす定数である．, 個の観測関数 )
（*
!
,）が，
,!
* #"

向に均一で周方向に有限の波数を持つ揺動が発達してい４．
２．
１節で述べた関数空間!内の要素であると仮定すると，る．これはケルビンヘルムホルツ揺動の特徴である．揺動

この分解（２）は，分解（１）の有限データによる近似を与えてエネルギーが大きく変動し，それと同期して背景場も変動いることが確認できる．実際，ある数学的仮定の下で，上している様子がわかる．また，図２下段で示しているス

記の動的モード分解により推定される固有値が，内在するナップショットを見ると，揺動の空間パターンも時間的に真のクープマン作用素のスペクトルの情報を捉えることも変動していることがわかる．リミットサイクルの周期は示されている［２３］．

%%$& $"!!程度であり，これは乱流の周期 %,-+'$"!に比

なお，動的モード分解で得られる分解（２）の各コンポー.
ネントは，推定した固有値 #+で決まる固有の周期性と減衰べ格段に長い（時間はイオンジャイロ周波数で規格化）．

率を持つ要素に分解されていることは，クープマン・モー

る．時刻 /における観測ベクトル #/を，そ の 時 刻 に お

ド分解（１）と同様である．そして，分解（２）における動的

け る 各 空 間 点 に お け る 密 度 #/'!
+(か ら 作 成 し，

このような動的乱流状態に動的モード分解を適用す

モード "+が，これらのダイナミクスがどのように各観測量

#
!
!
!
#/'
!
#/'
!
#/'
-!
+を用意する．そして，時
,(
/#*
"(
#(

%*へ相対的な重みで寄与するか，つまりある種の空間的な空間情報マトリックス ("#*
#
#
#
-!
!!
"!
%!"+と，観測時刻

２ 厳密 DMD において用いるデータの次元は，本質的には観測量 %*の数に相当する．分解（２）（または分解（１））では，%*の張る空間が十分大きいことを仮定しており（例えば［１８］），低次元なデータに対して適用する場合には，数理的に空間を広げる必要がある．具体的な手法もいくつか提案されており，例えば，基底関数を事前に準備する方法［１９］や，再生カーネルを用いた方法［２０］，ニューラルネットをデータから学習して用いる方法
［２
１，
２
２］
などが提案されている．

８１

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

図２

上段：乱流シミュレーションで得られた非線形飽和状態でのフーリエモード（m, n）の時間発展（m, n はそれぞれ周方向・軸方向
モード数）．背景流れ場は（0, 0）に対応している．下段：t ＝ 3050，3150，3200 における密度揺動のスナップショット［６］．

ルツ揺動，Mode4，5 はスパイラル構造に対応している．上述のように，解析している乱流状態は，プラズマの背景分
布，乱流揺動が乱流スケールに比べ長い時間でリミットサ
イクル振動しており，乱流の空間パターンも大きく変化する．この場合，各 DMD モードの時間発展を決める複素周
波数は指数関数的増加・減少しか記述することができない
ため，各モードの増減を表現するには工夫が必要である．
そこで，畳み込み積分を用いた固有関数の強度を評価する方法を提案する．この方法では，各時間ステップにおいて，得られた固有関数とオリジナルの乱流データとの間で畳み込み積分を行い，その空間最大値を固有関数（特徴的構造）の 強度 とする．すなわち，Mode-&の 強度 の瞬時
値 !&#($は，次のように計算される．
図３

支配的固有関数の空間パターン（Mode1~5）
［６］
．
!!

!
!'$
!!!'$
!'$
##(#
"&#
!
(
"&#
'!'" ##(#
'!'，
$
"!
（３）

を１ステップずらした %""%
"
"
"
(!
!!
"!
$ &とを準備し，前節で述べた特異値分解による低次元化した動的モード分解

!
!&#
(
"&#
!
(
)%
$
"(&
$
&
.

を行なった．ここで低次元化に用いたランクは '"%とし!

（４）

た．こうして，系の特徴量である固有値 #&とそれに対応す得られた !&#($を図４に示す．!&はうまく DMD モードのる固有関数 "&#!$を得ることができる．以下では，リミッ

増減を表現しており，プラズマ背景分布の変形とスパイラ

トサイクル振動とリミットサイクル周期に比べ高周波の乱ル構造の出現・消失まで含めた乱流パターンの変化を明確

流の挙動について，動的モード分解より得られた固有関数に捉えることができている．
DMD モードの 強度 !&を用いて，プラズマ分布や乱流（特徴的空間構造）間の因果関係を抽出する試みを紹介する．

構造間の因果関係を評価することができる．一例として，

図３に得られた固有関数（寄与率の大きな上位５つ）を

コ ヒ ー レ ン ト 揺 動（Mode3）と ス パ イ ラ ル 構 造

示す．ここでは，寄与率の最も大きなモードから順番に（Mode4）について，図４下段に示すように !#，!$の時間

Mode1 から Mode5 と呼ぶことにする．Mode１は背景密度発展からリサージュを書くことができる．コヒーレント揺分布，Mode2 は径方向に有限の波数を持ち周方向に一様な動の振幅が小さな時はコヒーレント揺動とスパイラル構造

帯状流成分，Mode３はコヒーレントなケルビンヘルムホの強度は共に増加していくが，ある程度コヒーレント揺動
８２

Special Topic Article

図４4. Data-Driven Analyses of Dynamical Turbulence Phenomena

M. Sasaki et al.

上段：各DMDモードの相関係数Cj の時間発展．下段：各リミットサイクル周期で描いた（C3, C4）間のリサージュ．ここでC3 はコヒーレントな KH 揺動，C4 はスパイラル構造に対応［６］．

が大きくなると，スパイラル揺動の強度が急増し，コヒーの次元%は時間ステップ数"よりも極端に小さくなること

レント揺動の強度が大きく減少することがわかる．この増が多い．これは設置できるセンサー数に限りがあるためで減を繰り返し，リミットサイクル振動が実現していることあり，数値シミュレーションとの大きな違いである．固有がわかる．また，スパイラル構造は，コヒーレント揺動の値分解に基づく標準的動的モード分解アルゴリズムでは，非線形飽和振幅を規定していることを示唆している．

モード数は最大でも % しか得られない．仮に，モード数が以上のように，DMD モードを用いて乱流構造間の因果""!%#以上あれば，観測ノイズまで含めてデータを正確関係を同定する方法を紹介した．リミットサイクル振動をに再構成できる．換言すると，回帰モデル #"!!#! は，縦伴うケルビンヘルムホルツ乱流に適用した結果を示した．

長のデータ行列には厳密にフィットし得るのに対し，極端リミットサイクルのように乱流の増減を繰り返す場合は，に横長のデータ行列には表現力が十分でない．別の見方を

畳み込み積分による構造 強度 の評価が有効である．得らすると，動的システムを構成する本質的なモードの数よりれた構造強度について，モード間のリサージュを描くことも少ないモード数でフィッティングすると，挙動の近しいで因果関係を抽出できることを示した．コヒーレント揺動

複数の本質的なモードが近似的に１つのモードでまとめての強度がスパイラル構造の出現に規定されていることを示表現される．さらに観測ノイズがあれば，ノイズに対応し唆する結果を得た．この方法は，他の様々な揺動観測に適て付加されるモードもまとめて表現される．こうして得ら用することが容易であり，幅広い波及効果を持つと期待しれたモードの重ね合わせでは，十分正確にデータを再構成ている．また，さらに拡張された発展的な動的モード分解

することができない．

手法もあり（例えば recursive DMD［３５］），今回の解析でこのような問題を解決するために，動的システムが有すは使用していない複素周波数の時間変化や位相ダイナミクるモード（本質的なモードとノイズ由来のモード）の数以スを抽出することで，乱流ダイナミクスの理解から輸送へ

上のモードを用いてフィッティングすることで必要なモーの理解への展開も期待される．また，ニューラルネット

ドを正確に推定したのちに，不要なモードを取り除くといワークを応用した構造間相互作用への発展も開発が進めらうアプローチを提案する．それぞれ，ハンケル行列作成とれている．

スパース動的モード分解によって実装することができる．

４．
３．
２ 乱流現象の定量的予測への応用

ハンケル行列は，元の時系列データから窓幅 'とスライド
幅 &で切り出したデータを $回積み上げて，次のように定動的モード分解をプラズマ乱流の観測データに適用する義される．

ことで，乱流挙動を予測する試みについて紹介する．ここでは，ハンケル行列とスパース動的モード分解を用いた方法について述べる．
一般に，実験から時系列データを取得する場合，データ
８３

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

"
&&
!
%!('
!
&
%
%&&
!
%!("'&&% & &+"(+'#&%!-(""'%&#: $%
%!
&
&&
!
%!("#
#
'
，
&
%
&
%
.
$
#
&&
(!
%'
!

:$+
.
!
!
!
-!
&&
)
!
)
"*
'
""!
"#!
"*,
)!
)
)
)

（５）

（６）

このハンケル行列に対して動的モード分解を適用すること
で，モード数は (+ だけ増加する．標準的動的モード分解
では，時刻 .における観測ベクトル !
.は行列 !の固有値分解から，式（２）のように展開された．スパース動的モード
分解では以下の手続きに よ り 寄 与 率 0*を 決 め 直す．まず，観測した時空間行列と再構成した時空間行列の差のFrobenius ノルムを最小化するための目的関数を $
"'と
&

図５

DMD モードの複素周波数．青がスパース DMD で除去された DMD モード，赤の DMD モードのみを解析に使用した．

図６

乱流時空間発展の学習と予測：上段から順にオリジナル計
測データ，ハンケル・スパース DMD を用いた再構成と予測，標準的 DMD を用いた再構成と予測［４０］．t < 3000［μs］
は学習領域，t > 3000［μs］
は予測領域．

し，L１正則化項を加えた次の最適化問題を解く．
"'
0
+)
,)
+)
0
'
&
"$' *
.
" $
**

（７）

*

最適化でゼロとなった寄与率 0*を拘束した条件下で，再び正則化項を外した $"'のみ最小化して，最終的な寄与率&0*を決定する．以上の手続きで，重要度の低いモードが取り除かれる．アルゴリズムの詳細は，文献［３６］を参照されたい．
我々は，動的モード分解による乱流データの解析を進めてきた［７］．解析に使用した乱流データは，九州大学にある直線磁化プラズマ PANTA 装置で観測されたイオン飽和
電流を用いた．周方向に６４チャンネル設置されているラン
グミュアプローブアレイによって１マイクロ秒ステップで計測された時空間発展を使用した．PANTA 装置では，放電条件によって帯状流やストリーマ，孤立波状態等の様々な乱流状態が存在するが［３７］，ストリーマ構造が支配的な条件で解析を行った［３８］．この放電では，5 kHz 程度のドリフト波乱流と 1 kHz 程度のストリーマ構造が共存している．このデータからハンケル行列を作成することは，上流側にプローブアレイを擬似的に増設していることに相当する．ここでは，-$"!!と ($%の設定で，時間区間%$3 ms測である［４０］．L1 正則化のパラメータ $は，DMD モードのデータで学習し，その先の予測を試みた．

の数が計測点数６４以下となるように調節した．標準的動的

この計測データにおいて，スパース動的モード分解の適モード分解では，時空間パターンは時間とともに減衰し用は，時空間データからのノイズ除去に相当する．ノイズて，細部も正確に再構成できていないのに対し，提案アプ除去の様子を，DMD モードの複素周波数の観点から図５

ローチではこの問題が大幅に改善されている．このようにに示す．ここで，複素平面における単位円付近の成分（赤精度良く再構成する動的モード分解を用いると，乱流の空丸）は動的システムの本質的なモードであり，単位円から間構造を含めた時間発展をより正確に予測できる．予測性離れた成分（青丸）はノイズ由来のモードと考えられる．

能を定量的に議論するため，以下の誤差指標"&.
'を導入す

実際，青丸で示した DMD モードは減衰が大きく，周波数る．

が高く，寄与率が小さい．ただし，ここで述べるノイズ除&$

去の方法は高周波なランダムノイズに対するもののみに適

"&
.
'
$
(

用可能なものであるため，計測器の基本的なノイズ対策は必要である［３９］．以降の結果は，赤丸で示した DMD モー#'' )$"(
/
.
.
'
!/
'
)
)&
)&
&$

#
' )$"/
.
'
)&

（８）

'
ここで，/
，計測
.
.
"%#.
'と /
'は，それぞれ，時刻 .
&
'
)&
)&

ドのみから得られるものである．

点 )における予測データと検証データである．誤差指標に

図６は計測データ（訓練データと検証データ），提案ア
プローチおよび標準的動的モード分解における再構成と予

対して，
８４

Special Topic Article

4. Data-Driven Analyses of Dynamical Turbulence Phenomena

M. Sasaki et al.

17H06089，16H02442，15H02155，15H02335，18K03578，
19J00871，九州大学応用力学研究所共同利用研究および
JST CREST JPMJCR1913 の助成をいただきました．R.O.
Dendy 教授，稲垣滋教授との有益な議論に感謝いたしま
す．

参考文献

図７

［１］P.H. Diamond et al., Plasma Phys. Control. Fusion 47, R35(2005).
［２］J.W. Connor et al., Phys. Plasmas 5, 2687 (1998).
［３］F. Wagner, Plasma Phys. Control. Fusion 49, B1 (2007).
［４］J.D. Callen and M.W. Kissick, Plasma Phys. Control. Fusion 39, B173 (1997).
［５］P.J. Schmid, J. Fluid Mech. 656, 5 (2010).
［６］M. Sasaki et al., Plasma Phys. Control. Fusion 61, 112001(2019).
［７］A. Kusaba et al., Plasma Fusion Res. 15, 1301001 (2020).
［８］A.A Kaptanoglu et al., Phys. Plasmas 27, 032108 (2020).
［９］R. Taylor et al., Rev. Sci. Instrum. 89, 053501 (2018).
［１０］H. Natsume et al., Phys. Plasmas 27, 042301 (2020).
［１１］C.W. Rowley et al., J. Fluid Mech. 641, 115 (2009).
［１２］P.J. Schmid, J. Fluid Mech. 656, 5 (2010).
［１３］K. Taira et al., Aiaa Journal 55, 4013 (2017).
［１４］S. Klus et al., J. Nonlinear Sci. 30, 283 (2020).
［１５］I. Mezi!, Annu. Rev. Fluid Mech. 45, 357 (2013).
［１６］B. Koopman, Proc. National Academy of Sciences of theUnited States of America 17, 315 (1931).
［１７］A. Mauroy et al., Physica D 261, 19 (2013).
［１８］J.H. Tu et al., J. Comput. Dyn. 1, 391 (2014).
［１９］M.O. Williams et al., J. Nonlinear Sci. 25, 1307 (2015).
［２０］Y. Kawahara, Neural Inf. Process. Syst. 29, 911 (2016).
［２１］N. Takeishi et al., Adv. Neural Inf. Process. Syst. 30, 1130
(2017).
［２２］B. Lusch et al., Nat. Commun 9, 4950 (2018).
［２３］M. Korda and I. Mezic, J. Nonlinear Sci. 28, 687 (2018).
［２４］B.J. Cantwell, Annu. Rev. Fluid Mech. 13, 457 (1981).
［２５］B.W. Brunton et al., J. Neurosci. Methods 258, 1 (2016).
［２６］J.L. Proctor and P.A. Eckhoff, Int. Health 7, 139 (2015).
［２７］J. Mann and J.N. Kutz, Finance 16, 643 (2016).
［２８］K. Fujii et al., Sci. Rep. 9, 16755 (2019).
［２９］K. Fujii et al., Proc. 2017 European Conf. on Machine
Learning and Principles and Practice of Knowledge Discovery in Databases (ECML-PKDD'17), pages 127 (2017).
［３０］T. Bito et al., Proc. 2019 Int'l Joint Conf. on Neural Networks (IJCNN'19), number 19278 (2019).
［３１］K. Fujii et al., Sci. Rep. 10, 3005 (2020).
［３２］H. Shiraishi et al., J. Neural Eng. 17, 036009 (2020).
［３３］B.M. Annaratone et al., Phys. Plasmas 18, 032108 (2011).
［３４］M. Sasaki et al., Phys. Plasmas 26, 042305 (2019).
［３５］B.R. Noack et al., J. Fluid Mech. 809, 843 (2016).
［３６］M. R. Jovanovi!et al., Phys. Fluids 26, 024103 (2014).
［３７］T. Kobayashi et al., Plasma Fusion Res. 12:1401019 (2017).
［３８］T. Yamada et al., Nat. Phys. 4, 721 (2008).
［３９］永 島 芳 彦，稲 垣 滋：プ ラ ズ マ・核 融 合 学 会 誌 96, 2(2020).
［４０］A. Kusaba et al., submitted to Plasma Fusion Res. (2020).
［４１］F. Kin et al., Phys. Plasmas 26, 042306 (2019).

予測誤差指標 E（t）の時間発展．t ＞ 3000［μs］
が予測領域．

!#
"
"
!
"
!$
#%$#
$
!!#
$

（９）

というフィッティングを行うことで，予測可能時間 "を評価する．図７は !#"$の対数プロットであり，"""!!!マイクロ秒程度であると評価できる．
現状では支配的乱流揺動（周期２００マイクロ秒程度）の数
周期程度の予測が空間構造を含めてできている．解析に用いたショットでは，ストリーマに起因するミリ秒スケールの空間局在した輸送が起こっている［４１］．今後，より長時間の予測が可能となれば，このような突発的・間欠的な輸送現象の空間構造を含めた予測にも応用できる可能性があ
る．

４．
４ まとめ
本章ではデータ駆動科学的手法として近年注目されている動的モード分解について，理論的基礎から応用例まで紹介した．動的モード分解は流体現象の時空間構造を観測／
計測データから抽出するための数値解析手法の一つとして提案されたものであり，ここではプラズマ乱流への応用について紹介した．動的モード分解を用いると，時空間発展
データから統計的に重要な空間構造とその構造に対応する
複素周波数を得ることができ，より少数の自由度でシステ
ムを理解することが可能となる．物理機構の理解をめざした応用例として，大域的乱流シミュレーションで得られた乱流場に適用した例を紹介した．ここでは動的モード分解
によって抽出された特徴的空間構造を元データとの畳み込
み積分を行うことで構造間の動的因果関係を抽出する方法を提案した．また，時空間発展予測へ向けた応用例として，
実験で観測された乱流の時空間発展について，ハンケル行
列とスパース動的モード分解を適用し，空間構造を含む乱流の発展を予測した試みを紹介した．今後，突発的輸送等の空間発展予測を可能にすると期待される．

謝辞

本研究は九州大学若手研究者グローバルリーダー育成型
Progress 100（NB８
０６４５０２８）および九州大学数理・データ
サイエンスに 関 する 教 育・研 究 支 援 プ ロ グ ラ ム の支援を受けました．また JSPS 科研費 JP18H03287，16K18335，

８５

J. Plasma Fusion Res. Vol.97, No.2 (2021)８
６‐９２

小特集 磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展

５．データ駆動アプローチを用いた雪崩的乱流輸送現象の解析
5. Data-Driven Analyses of Avalanche Like Turbulent Transport Phenomena
朝 比 祐 一，藤 井 恵 介１）
ASAHI Yuuichi and FUJII Keisuke１）
日本原子力研究所，１）京都大学（原稿受付：２
０
２
０年１
０月１
２日）

本研究では，第一原理的ジャイロ運動論的シミュレーションによる大規模データを，データ駆動科学的手法
により解析した．まず，少数の波が支配的なコヒーレントな状態と様々な波が入り乱れる乱雑な状態の判別を，特異値分解を用いて行った．これにより突発的に起こる熱輸送現象のあとプラズマは乱雑な状態になること，乱雑さはその後自発的に減少すること，次の突発現象はそのような自己組織化の後に起きることが明らかになった．この過程は Landau 減衰を始めとする速度空間構造の変化と密接に変化していると考えられる．しかし，従来手法では５次元位相空間構造の時系列解析は不可能であった．そこでさらに主成分分析による位相空間構造デー
タの圧縮技術を開発した．圧縮されたデータを利用しても突発的輸送が表現できることや，どのような位相空間構造が突発的輸送と関連しているかを論じる．
Keywords:
gyrokinetic simulation, plasma turbulence, singular value decomposition, dimensionality reduction,
avalanche like transport

５．
１ はじめにターン形成と輸送の関連性が調べられてきた．例えば，文

磁場閉じ込め核融合装置の炉心プラズマの閉じ込め性能

献［３］では突発的熱輸送が起きているときに，静電ポテン

は乱流輸送によって支配される．ジャイロ運動論的シミュ

シャル内のコヒーレントなモード構造が見られることが指レーション［１］は理論，実験解析問わず，プラズマ乱流輸送

摘されている．一方で高速粒子輸送研究などからは位相空解析に幅広く利用される第一原理的解析手法である．計算間内で形成されるパターンと突発的な輸送現象との関連性

は空間３次元，速度２次元，時間１次元の合計６次元空間が示唆されている［４］．本研究ではこの突発的熱輸送現象

内で行われ，膨大なシミュレーションデータが生成されに着目する．特に，この現象がプラズマ内のどのような波る．このデータを解析する上では，大きく分けて三つの困のダイナミクスによって起こっているかを調べる．

難が存在する：

まずはじめに大規模シミュレーションデータが，少数の波が支配的に励起されているコヒーレントな状態と，様々

１．非線形相互作用による乱流自体の解釈困難さな波が入り乱れている状態のどちらに位置するかを定量化

２．データの高次元性による可視化の困難さすることを考えた．具体的には，特異値分解を用いてプラ３．データの大規模性によるデータ処理の困難さ

ズマの空間分布から空間的な乱雑さ，つまり擬似的なエン

本研究ではこのような大規模，高次元データから物理的知

トロピーを定義した．これにより突発的輸送現象によりプ見を抽出することをめざし，データ駆動科学的手法によりラズマの乱雑さが増大すること，プラズマの乱雑さはその上記のジャイロ運動論的シミュレーションで得られたデー

後自発的に減少すること，次の突発現象はそのような自己タの解析を行った．

組織化の後に起きることが明らかになった．この過程は本研究で用いる full-F ジャイロ運動論的シミュレーショ

Landau 減衰を始めとする速度空間構造の変化と密接に関ンでは，分布と揺動の時間発展を自己無頓着に解くため，連していると考えられるが，従来手法では５次元位相空間

自発的に励起される突発的な輸送現象を扱うことができ構造の時系列解析は不可能である．そこで主成分分析による．この突発的熱輸送現象は空間的に伝播すること，またる位相空間構造データの圧縮技術を開発した．圧縮された熱流束の大きさの周波数スペクトルが !
!
!分布に近くなる

データを利用しても突発的輸送が表現できることや，どのことから，自己組織化臨界現象（self-organized criticality，

ような位相空間構造が突発的輸送と関連しているかを論じ

SOC）［２］との関連が議論されてきた．従来研究では，このる．これらの解析と関連する先行研究については５．
３節と突発的熱輸送現象に関して流体モーメントの３次元空間パ５．
４節で別途論じる．

Japan Atomic Energy Agency, CCSE, 178-4-4 Wakashiba, Kashiwa, CHIBA, Japancorresponding author’s e-mail: asahi.yuichi@jaea.go.jp

８６!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Special Topic Article

5. Data-Driven Analyses of Avalanche Like Turbulent Transport Phenomena

５．
２ 特異値分解とデータセット

Y. Asahi and K. Fujii

場合は固定されたフーリエ基底に対応する係数を求める

本節では５．
３節および５．
４節における共通の手法である特が，特異値分解は式（２）を最小化するように基底も同時に異値分解（Singular Value Decomposition，SVD）とデータ求めていることになる．このように特異値分解はフーリエセットについて説明する．

変換と共通する性質がある．そのため例えばプラズマ中で５．
２．
１ 特異値分解（SVD）

ある波が卓越しているときの主な主成分は，その波の sin，

本研究では，ある物理量に関する時系列データを解析対

cos 成分として現れることが多い．

象とし，これを ) 行，*列の行列 $で表現する（ここでは一方で，フーリエ基底は一様空間でのラプラシアン固有)$*と仮定する）．行列 $の特異値分解は以下で表され関数である．そのため一様乱流データを表現するには，る．

フーリエ基底がもっとも性能がよいと期待される．しかし$#%Σ&T，

トーラス形状をしているプラズマのデータを最もよく表す（１）

基底がどんなものかは自明ではない．

ここで %を左特異値行列（) 行，) 列），&を右特異値行

また特異値分解では，Frobenius norm の意味で最も元の列（*行，*列）と呼ぶ．どちらも正規直交行列である．ま

情報をよく表した低ランク行列が得られる．特異値の分布

た，Σ は ) 行 *列の非負値対角行列であり，特異値行列とがデルタ関数的であるとき，つまり少数の特異値が大きく呼ばれている．Σ の対角成分 $"%$#%(%$) &%!'を特異てそれ以外の特異値が小さい場合，大きな特異値に対応す値と呼ぶ．

る右特異値行列・左特異値行列を使うことで元データを効

特異値の大きさは $に対するFrobenius normの意味での率よく圧縮できることになる．

寄与の大きさと対応する．そのため，寄与の小さな特異値

５．
２．
２ シミュレーションデータ

に対応する成分は無視しても $を近似的に表現できる．よ本研究で用いるジャイロ運動論コードGT5D
［５，
６］は，５
次元分布関数 %の時間発展を計算する．分布と揺動の時間

り正確に書くと，

発展を自己無頓着に解くfull-Fモデルである．図１は，磁気#!&$&'!$&''

面平均された物理量の時空間発展を示す［７］．図１（a）は，

（２）

!
&
'

""! ドリフトに起因するイオン乱流エネルギーフラック
ス #&" を示しており，径方向へ伝搬する突発的な輸送が起を最小化するランク +") の行列 $は，寄与の小さな特異値（(#+）を無視したこっていることがわかる．
図２（a）と（b）に示すように雪崩フェイズと静的フェイ+$&'#! ,&($(-'(

ズの静電ポテンシャルや密度の揺動成分の構造を見比べる（３）と，雪崩フェイズにおいてはコヒーレントなバルーニング(#!

に一致するということである．なお，ここで .
&'は行列 '

構造が見られるのに対し，静的フェイズでは明瞭なモードの &，'成分を表す．

構造が見られない．また，図２（e）と（f）を比較するとわか例として，行列 $として時系列シミュレーションデータ

るように，密度と温度では構造が異なる．これは速度空間を使うことを考えよう．時間方向を $の行方向，空間方向においても構造形成が起こっていることを示唆している．

を $の列方向に並べるとする．つまり，ある時間スライス&，ある空間点 'でのデータを $&'とする．特異値分解は，この行列を時間に依存しない正規直交の基底 &T と空間に依

５．
３ 乱流シミュレーションデータにおける乱雑さの定量化と突発的現象の関連存しない係数 %Σ に分解する．これは正規直交基底展開の

図２（a）と（d）を見比べると雪崩フェイズにおいてはコ一つであるフーリエ変換とよく似ている．フーリエ変換のヒーレントなモード構造が存在するのに対し，静的フェイ

図１ （a）
乱流エネルギーフラックスQEi ，（b）
規格化イオン温度勾配R0/Lti，（c）
径電場シアdEr/dr の時系列データ．
５．
４節では赤い点線に囲まれた部分について解析を行った［７］．

８７

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

図２

雪崩フェイズ（tR0/vti ＝ 250）における（a）
静電ポテンシャル $，（b）イオン密度 #ni，（c）
イオン温度 #Ti と静的フェイズ（tR0/vti ＝ 320）
イオン温度 #Ti のポロイダル断面図［７］．
における（d）
静電ポテンシャル $，（e）イオン密度 #ni，（f）

ズでは構造が乱雑化していることが定性的に理解できる．

が存在し，静的フェイズでは構造が乱雑化しているという少数の波が卓越しているか，多くの波が乱立しているかの観察結果と一致する．

違いを擬エントロピーで表すことを考える．つまり，単一

このようにして求めた擬エントロピーの時間変化を図４の波しか存在しないようなコヒーレントな状態は擬エント（c）に，熱流束の時間変化図４（a）とともに示す．図中縦のロピーが小さく，様々な波が入り乱れている状況は擬エン

点線で，半径位置 0.4 にて熱流束が極大値をとった時刻をトロピーが大きいと表すわけである．これまで，トーラス

示す．擬エントロピーは，熱流束が極大となる付近で最も形状など任意の形状のプラズマに対して，そのエントロ

上昇し，熱流束が小さくなった時刻付近で減少に転じていピーの議論はあまりされてこなかった．これは，トーラス

ることがわかる．

プラズマでは１つのモードを単一のフーリエモードで記述

不安定バルーニングモードの一つ（トロイダルモード&""#）の振幅の時間変化を図４（b）に示す．この不安定

できないことに関連すると考えられる．
トーラスプラズマに対してこの擬エントロピーを定量的モードの振幅は熱流束の大きさとほとんど同じ変化をしてに評価することを考える．まず，図３（a）に示すように，各いる．このことから，熱流束はこの不安定モードが成長す

時刻において小さな窓空間における 揺 動 の イ オ ン密度ることで誘起されていることがわかる．また，この不安定#&"&!&%"!!&"!の空間データを様々な窓位置に対して取モードが成長するに従って擬エントロピーも大きくなって得し，それを行列 $"#!"!
!#!
!! $で表現する．ただし，%!

いく．擬エントロピーがある程度大きくなると，この不安!は窓の総数であり，!#はそれぞれの窓内部のデータと対定モードは急激に減衰する．擬エントロピーは，どの程度応する．式（１）で定義される特異値分解の結果から以下の様々な波が系に存在するかを示すものであることから，大&'%式で示される擬エントロピー "を評価した．

きくなった不安定モードが非線形相互作用により他の波を励起していること，また励起した他の波へのエネルギー輸!

"&'% "!" %#)*(#%#，

送により不安定モードの振幅が小さくなっているのだと推

（４）

察される．

#

ここで，%#"%#"!!
$は規格化された特異値である．擬エ$%&'%は，入力データ $を記述するのに実効的にさを定量化することができた．二次元プラズマ乱流に対し必要な基底の数を求めることに対応する．例えば，単一のて 固 有 直 交 分 解（Proper Orthogonal Decomposition，

不安定モードだけが存在するときは小さい値が得られ，多POD）を適用した例でも同様のエントロピーが定義さ数のモードがほぼランダムに分布するときは大きな値が得れ，時系列データ全体の乱雑さとの関連が述べられている

ントロピー "

このように，特異値の分布を用いることでデータの乱雑られる．実際，図３（b）を見ると，雪崩フェイズでは少数の［８］．本研究はこれと異なり，各時刻の乱雑さを定量化す

基底でデータが表現できているのに対し，静的フェイズでることで，系の乱雑さの変化と雪崩的輸送を関連づけるこはデータを表現するのに多数の基底が必要となっている．

とができた．この手法の特徴は，データの次元数に制限が

これは雪崩フェイズにおいてはコヒーレントなモード構造ないことである．今回はプラズマ密度の３次元分布を用い
８８

Special Topic Article

5. Data-Driven Analyses of Avalanche Like Turbulent Transport Phenomena

Y. Asahi and K. Fujii

図３ （a）
突発的輸送現象がおこっている時およびその直後におけるプラズマ密度ゆらぎの分布．この３次元データから複数の小さな窓領域を取得したのち，特異値分解を適用する．（b）
突発的輸送現象の最中およびその直後のデータから得た特異値の分布．輸送現象の最中では特異値の分布がより尖っている（よりコヒーレント）であるのに対し，直後では特異値の分布がより平坦（より乱雑）である．挿入図で，SVD で得られた基底の典型的なものを示す．

１）
SVD を元にした擬エントロピーの時空間分布．（a‐２），
図４ （a‐
１）
熱流束の時空間分布，（b‐１）
n ＝ 12 モード強度の時空間分布，（c‐
（b‐
２）
，（c‐
２）
にはそれぞれ，#＝ 0.4 の位置での時間発展を示す．

たが，速度空間を含めることもできる．これまで特に速度縮性に関する議論が行われているのに対し，本研究では，

空間でのエントロピーが精力的に調べられてきた
［９，
１０］．

より大規模（〜10 TB）の大域的乱流のデータ圧縮およびそれらの研究との比較が今後の課題である．

データ可視化手法に関して論じる．分布関数の揺動部分の時間発展 "!
"!
!!
#"!!
$"!に対して主成分分析を適用する

５．
４ 主成分分析による５次元時系列データの分解とその解釈と，##%Σ&$

本節ではデータ解析において幅広く使われる主成分分析%（PCA）を用いて分布関数時系列データの解析を行う．関

& $
&"
"! &"

連 研 究 と し て は，ジ ャ イ ロ 運 動 論 の 時 系 列 デ ー タ を"""

Tucker 分解によって圧縮した例［１１］が知られる．この研と 近 似 で き る．こ こ で %は 用 い た 主 成 分 の 数 で あ り，

究では，比較的小規模（〜10 GB）の局所的乱流のデータ圧

&
&"
%Σ%
"$
" を主成分 "の係数（%Σ の "列目），&" を主成

８９

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

分 (の基底（-の (列目）と呼ぶこととする．+%"
，
&
&
!"

析した結果の右特異値ベクトル（以下基底と呼ぶ）を示し，&の平均であり，主成分分析ではデータから平均を引"&は "

図７は，左特異値ベクトル（厳密には左特異値ベクトルといた上で特異 値 分 解 を 行 う．主 成 分 分 析 を 行 う た めに

特異値行列の積で，以下係数と呼ぶ）の雪崩フェイズと静は，５次元分布関数の時系列データ "&を行列の形状へ変換

的フェイズの空間構造を示す．ここで，図６と図７を見比する必要がある．この変換には自由度が存在し "は例えば，&べることで可視化困難であった５次元データを３次元以下， '，$.$/ (
$,!
$&!
$+!
$$!
$&(
$+!
$$(
$.$/ ($,!
#'#''の画像の組み合わせとして表現できていることに注意され

$,(
$+!
$$!
$&!
$.$/ (などと表せるとわかる．ここで'#')(*(は)行，*列の行列を意味し括弧内の次元につい'#'たい．プラズマの実験解析においては，同様の手法によりては一次元データとして平坦化することとする．図５は，ばしば行われる［１２‐１４］．本解析では同様の考えを進め，同

$,!
!
!
!
!
!
$+!
$$!
$&!
$.!
$/ (
)!!
"!(
#&'
$#
*'
#!(の時系'%'様の解析を高次元性に起因する可視化の困難さを克服する

列データを異なる形状の行列として行列化し，主成分分析ために利用している．

)

)

)

データの時間情報と空間情報を分解して解析することがし)した際の累積寄与率とデータ圧縮率の関係を示す．なお，

図７において，第０主成分は磁場強度分布，第３主成分

電子とイオンの分布関数に関しては独立の行列として処理

は対流セルと呼ばれる構造であることがわかる．第１，第した．ただし，主成分 +（ここでは +%'%）までの累積寄与

２主成分は，*%"#のバルーニングモードと対応してい率は以下で定義した．

る．図７の（a）に示す雪崩フェイズと（b）に示す静的フェイズを比較すると，第０主成分と第３主成分は定常的であり+## (%"%(.
*## (%"%(

ほぼ変化しないのに対し，第１，第２主成分は雪崩フェイズで卓越していることがわかる．これは流体モーメントにおいてコヒーレント構造が卓越していることと対応していデータ圧縮率は，係数と基底のバイト数を全データのバイる．各主成分が雪崩的熱輸送とどのように関連しているか

ト数で規格化したものである．

を調べるため，以下の分布関数再構成の式を用いて各主成
分のエネルギー輸送に対する寄与を調べた． &自体の再構

図５（a）より，電子に関しては分解方法によらずデータ

成結果については文献［７］を参照されたい．

を表現できていることがわかる．一方図５（b）に示すようにイオンに関しては，速度空間を基底とする行列の方が

#
)'.)
3"
3
%
!"#! $
.)/##
&+
(
"#!"&
!
.'
'%$/
#

データの表現力が高いとわかる．Tucker 分解により５次元分布関数時系列データの圧縮を行った先行研究［１１］にお3(，
310.2"# %3!!"%%%

いても，速度空間に関しては少数の成分でデータを表現で(きると指摘されている．
ここで，以下ではこれらの特性を加味し，３桁のデータ圧縮率と高い累積寄与率を両立している'$&!
$,!$+!$$(
$.)!
$/ (
#'

* ,
3
&
&
-(
"# -(
%&
!!""

の分解でデータ解析を行う．

(

図６は，'
$&!
$,!
$+!
$$(
$.!
$/ (の行列を主成分分
#'
)

図５

異なる行列形状に対する（a）
電子分布関数と（b）
イオン分布関数の累積寄与率の圧縮率依存性．

９０

（５）

Special Topic Article

5. Data-Driven Analyses of Avalanche Like Turbulent Transport Phenomena

図６

図７

Y. Asahi and K. Fujii$＝ 0 断面［７］第１６主成分までの位相空間基底 VTk ($, v‖, w)の（a）
w ＝ 0.12 断面と（b）．

第１６主成分までの（a）
雪崩フェイズおよび（b）
静的フェイズの空間係数 uḱ（r, #（）ポロイダル断面）
［７］
．

#
'%+)
-!!$#)+)""
!%"$#
#
&)
(
!"!"$
!
+'
!!
)
#

５．
５ まとめ

#
'%+)
-+*(,$#)+)""
!%"$#
#
$
&)
(
!"!"!
!
+'
)
#

ミュレーションの時系列データの解析手法について示し

#
'%+)
* &
-& $#)+)""
!
#
)
'&，
'
#
&
(
!"!"*&!
+ %"$)#た．このシミュレーションでは，突発的な輸送現象を扱うことができる．シミュレーションデータは，乱流の非線形である．図８は，式（５）の右辺最終項に現れる各主成分の解析が困難となっており，従来は平均操作などを用いた上エネルギーフラックスへの寄与である．図から熱輸送をで解析が行われてきた．本章では，データ駆動科学的手法担っているのは主成分１と２であることがわかる．図６とにより，シミュレーションデータの乱雑さの定量化と高次

本章では，特異値分解を利用したジャイロ運動論的シ性，データの高次元性，大規模性などの要素を有するため

７からこれらは ($"#のバルーニングモードと対応し，そ

元データの次元圧縮による可視化および解析を行った．この速度空間構造は雪崩フェイズにおいて #%!の悪い曲率

れらはともに特異値分解を利用しており，複雑な現象から側で形成される速度空間構造と対応していることが明らか

より重要な要素を抽出することを目的としている．

となった．つまり，主成分分析によって雪崩的輸送を引き

まず，プラズマ乱流ではコヒーレントな状態と様々な波起こす空間，速度空間構造が物理的な知見を利用せず抽出が入り乱れている状態が繰り返されていることに着目し，できたと言える．

その乱雑さを定量的に評価することを試みた．密度の時空間データを様々な窓に対して取得し，それを行列で表し
９１

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

つ，８３％の累積寄与率を保持できた．これは，データを３次元の位相空間基底と２次元の空間係数の時系列データとして表現したことに対応する．位相空間基底と係数を用いて，各主成分のエネルギーフラックスへの寄与を調べ，!!!"のバルーニングモードのみが雪崩的輸送を駆動していることを特定した．対応する位相空間基底を調べることで，これらの速度空間構造は悪い曲率側で形成される速度
空間構造と対応していることが明らかとなった．主成分分析によって雪崩的輸送を引き起こす空間，速度空間構造が物理的な知見を利用せず抽出できたと言える．

参考文献

図８

［１］X. Garbet, Y. Idomura, Nucl. Fusion 50, 043002 (2010).
［２］P. Bak, C. Tang and K.Wiesenfeld, Phys. Rev. Lett. 59, 381(1987).
［３］W. Wang et al., Nucl. Fusion 60, 066010 (2020).
［４］H.L. Berk et al., Phys. Plasmas 6, 3102 (1999),
https://doi.org/10.1063/1.873550.
［５］Y. Idomura et al., Comput Phys. Commun. 179, 391 (2008).
［６］Y. Idomura, J. Comput. Phys. 313, 511 (2016).
［７］Y. Asahi et al., Phys. Plasmas 28, 012304 (2021).
［８］S. Futatani et al., Phys. Plasmas 16, 042506 (2009),
https://doi.org/10.1063/1.3095865.
［９］T.-H. Watanabe et al., J. Comput. Phys.: Conference Series
510, 012045 (2014).
［１０］T. Tatsuno et al., Phys. Rev. Lett. 103, 015003 (2009).
［１１］D.R. Hatch et al., J. Comput. Phys. 231 4234 (2012).
［１２］T.D. deWitt, Plasma Phys. Control. Fusion 37, 117 (1995).
［１３］J. Levesque et al., Nucl. Fusion 53, 073037 (2013).
［１４］B.S. Victor et al., Plasma Phys. Control. Fusion 57, 045010(2015).

各主成分のイオンエネルギーフラックスへの寄与［７］．

た．この行列を特異値分解し，その特異値から擬エントロ
ピーを定義した．この擬エントロピーは先ほどの時空間
データを表現するのに必要な基底の数を求めることに対応
し，単一の不安定モードだけが存在するときは小さい値が得られ，多数のモードがほぼランダムに分布するときは大きな値が得られる．この手法によりデータの乱雑さを定量化し，雪崩的輸送の前後での不安定モードの挙動を抽出することができたと言える．
次に，主成分分析による５次元時系列データの解析を行った．５次元時系列データを位相空間基底と空間係数へと分 解 す る こ とで，３桁 に 及 ぶ デ ー タ 圧 縮 を 実 現しつ

９２

J. Plasma Fusion Res. Vol.97, No.2 (2021)９３‐９４

小特集 磁場閉じ込め核融合プラズマにおけるデータ駆動アプローチによる物理モデリングの進展

６．まとめ
6. Summary
今寺賢志
IMADERA Kenji
京都大学大学院エネルギー科学研究科（原稿受付：２
０
２
０年９月３
０日）

本小特集では，磁場閉じ込め核融合プラズマの分野で，る融合研究の進展が挙げられるであろう．理化学研究所のシミュレーション，統計解析，データ駆動科学の融合に

スーパーコンピュータ「富岳」は，２０２０年６月期と１１月期

よって行われた４編の論文を紹介したが，大別すると統合に TOP５００で世界第１位を獲得したが，人工知能（AI）でシミュレーションで用いる輸送モデリングの高速化と高精主に用いられる単精度や半精度演算処理に関する性能ベン度化を行った前半（２章，３章）と，主成分分析によってチマーク HPL-AI においても世界第１位を獲得している

説明可能な物理モデリングの検討を行った後半（４章，５［１］．また研究テーマとしても，様々なシミュレーション章）に分かれる．更に前半を分類すると，２章ではジャイ研究に加え，シミュレーションと AI の融合研究が重要

ロ運動論／流体シミュレーションから得られたデータをテーマとして位置付けられている［２］．海外においても，

ニューラルネットワークで学習し，計算コストが高い元の

例えばアメリカのオークリッジ国立研究所の「SUMMIT」

モデルに取って代わる「新たな輸送モデルを開発」したの（TOP５００で世界第２位）は，アーキテクチャ的にもAIに特に対して，３章では実験が再現されるように，統合シミュ化したシステムであり，その点を重要視していることは間レーションを行っていく中で「既存の輸送モデルを最適違いないであろう［３］．これらのハードウェアを有効利用化」していったという点でアプローチが異なっているが，することによって，今後も深層学習等の計算の高速化と高共に実験を高精度に再現している点は興味深い．また今回精度化が相補的に進展し，新たな物理モデリングの開発に得られた輸送モデルが高精度であることは，経験則から求応用されることが期待される．

めた単純なモデルで輸送現象を表現することが難しいこと

多階層複雑系であるプラズマを研究する際，扱う物理現を示唆しており，多次元パラメータ空間で非滑らかな関数

象やその時空間スケールに応じてある特定の要素に還元しを表現可能なデータ駆動科学の特徴が発揮された結果であて問題の解明に取り組む方法が王道であるが，近年ではそるとも言える．そのようなモデリングの柔軟性は，逆説的の要素を全体として統合した試みが多く行われているには人間が現象を理解する上で障壁となる可能性がある［４，５］．一方で，それが進むにつれて，現象の非線形性，非が，高い再現性は背景にある物理を推論する上での保障と局所性，多階層性，複合性は増していき，より柔軟なモデなるという意味で，輸送モデリングが一歩進展したことはリングが求められるであろう．そのような問題に取り組む

間違いないであろう．

ためにも，シミュレーション（実験），統計解析，データ駆

一方，後半では，４章で動的モード分解によって抽出さ動科学の相互補完が重要であり，演繹（deduction）と帰納れたケルビンヘルムホルツ揺動とスパイラル構造の動的因（induction）を併用する必要性は高まっていくと考えられ果関係を同定しており，５章で雪崩的突発輸送を引き起こる．現在，データ駆動科学と核融合・プラズマ科学の融合

す主モードを実空間・速度空間双方から同定している．い

研究を促進させる機会創出の試みも進んでおり［６］，これずれも従来の統計解析では困難であった物理現象を明らかからの計算機性能の向上と並行する形で，このような融合にしており，説明性の観点から物理モデリングの進展に大研究は進展するであろう．ただし，論理的推論のもう１つきく資するものである．また４章では実験データから乱流の方法である仮説形成（abductive reasoning）を含め，プの時空間発展を予測する試みが，５章では位相空間構造

ラズマ科学の知識（ドメイン知識）が研究の根幹であるこ

データの圧縮技術開発がそれぞれ紹介されており，適切なとは明白であり，プラズマ科学者が主となってこのような次元削減を行うことで，少数モードで主要な物理現象が表研究が進展することを期待したい．

現できる主成分分析が持つ様々な可能性を提示している．
今後の展望としては，大型計算機の発展と連動した更なGraduated School of Energy Science, Kyoto University, Gokasho, Uji, KYOTO, 611-0011, Japanauthor’s e-mail: imadera@energy.kyoto-u.ac.jp

９３!2021 The Japan Society of Plasma
Science and Nuclear Fusion Research

Journal of Plasma and Fusion Research Vol.97, No.2 February 2021

［５］プラズマ・核融合学会 編：プラズマシミュレーショ
ン―多階層複雑現象の解明へ―（京都大 学 学 術 出 版
会，２
０
１
８）
．
［６］統計数理研究所 共同研究集会「諸科学における大規
模データと統計数理モデリング」
（２
０
２
０）
．

参考文献
［１］https://www.top500.org/［２］https://www.r-ccs.riken.jp/jp/［３］https://www.olcf.ornl.gov/summit/［４］岸本泰明 他：プラズマ・核融合学会誌 79, 460 (2003).

!!
いま

でら

けん

小特集執筆者紹介

じ

始や何やらと多忙な日々が続いた一方で，運動不足による体重増．これは何とかせねばと最新計測機器を導入し，倍返しに成功，これもデータ駆動？

今 寺 賢 志京都大学大学院 エネルギー科学研究科 プラ
ズマ・核融合基礎学分野 准教授（エネルギー博士）．専門は大域的ジャイロ運動論シミュ
レーション研究ですが，最近はその「大域性」
をデータ駆動科学で解析できないか模索しています．Zoomに映った自分の顔を見て，６歳の可愛い息子に似ていると思い AI で診断したら，ほぼ他人でした．AI のことが少し嫌いになりました．
なり

た

え

よこ

み

だ

みつる

充

した

ゆう

うえ

さ

さだ

げん

た

き

まこと

佐々木

さ

真

九州大学応用力学研究所助教．２
０
０
９年東京大
学大学院理学系研究科・博士（理学）．日本電気株式会社中央研究所を経て２０
１
０年１月より現職．プラズマ乱流の非線形過程に関する理論・シミュレーション研究に従事．最近はデータ駆動科学的手法を用いた研究を通じて，プラズマ研究が幅広い非平衡システムの時空間ダイナミクスの研究にも有益だと認識してきた．

や

京都大学大学院工学研究科原子核工学専攻，博士課程１年．統合輸送シミュレーションに基づくデータ同化システム ASTI の研究・開発を行っています．ASTI の実用化を目指して，日夜 PC に向かっています．理論が理解できないときやプ
ログラムが言うことを聞かないとき，結果が思い通りにならないときは，鴨川で現実逃避をしています．
かみの情報・システム研究機構統計数理研究所モデリング研究系教授，同機構データサイエンス
共 同 利 用 基 盤 施 設 教 授．京 都 大 学 博 士（理学，地球惑星科学専攻）．学生のころはプラズ
マ不安定性の研究をしていましたが，修了後は統計科学，特にデータ同化の方法の研究に移りました．速度分布関数の代わりに状態ベクトルの確率分布が現れます．MHD の代わりにカルマンフィルタ，PIC に代わりに粒子フィルタ，……．手計算と近似計算とスパコンで戦っています．

森 下 侑 哉むらゆき

上 野 玄 太量子科学技術研究開発機構 那珂核融合研究所
先進プラズマ研究部 上席研究員．トカマク中の輸送現象や輸送シミュレーションを研究対象としていたような気もしますが，機械学習や最適化問題にも手を出すなど，自分が何をしているのか自分でもよく把握していません．ご多分に漏れず，子どもたち（男３人）が鬼滅の刃にどっぷりで，日々日輪刀で斬られています．人を斬るのは隊律違反ではないかと指摘しても一向に止まないところを見ると，私は鬼なのかもしれません．
もりまさ

自然科学研究機 核融合科学研究所 ヘリカル研究部 核融合理論シミュレーション研究系教授．工学博士（京都大学大学院原子核工学専攻）LHD 実験解析型統合輸送解析スイート
TASK3D-a の運用で蓄積した熱輸送データベースを用いて，統計手法やデータ駆動手法との連携を実践，模索しています．さらに，プラズマ・核融合をはじめとする諸科学における大規模データと統計数理分野との連携機会の創出に楽しく取り組んでいます．お声がけさせていただいた際には，楽しく巻
き込まれてください！

量子科学技術研究開発機構 那珂核融合研究所
先進プラズマ研究部 研究員．２
０１
５年大阪大学
大学院・博士（工学）．トカマクプラズマにおける乱流輸送を対象とした実験解析やモデル
開発，統合コードの開発に従事している．根っからのインドア派のため，コロナ禍になって以来，後ろめたさもなく，休日は堂々と家に引きこもれています．
ほんやま

横 山 雅 之成 田 絵 美

本 多

""

かわ

はら

よし

のぶ

河 原 吉 伸２００
８年３月東京大学大学院工学系研究科博士課程修了．博士（工学）．２０１
９年４月より九州
大学マス・フォア・インダストリ研究所 教授．また，２
０
１
６年９月より理化学研究所革新知能統合研究センターにおいてチームリーダーを兼任．令和２
年度科学技術分野の文部科学大臣表彰若手科学者賞などを受賞．統計的機械学習の基礎・応用に関する研究や，大学・大学
院における教育に従事．

よし

村 上 定 義京都大学大学院工学研究科 教授，理学博士（広島大学）．非軸対称トカマクを含む広い３次元系トーラスにおける加熱や新古典輸送など運動論的物理現象について研究を行っている．また最近は，核融合プラズマの実時間制御を目指して，
データ同化を用いた統合輸送シミュレーションコードの開発なども進めている．春頃のコロナ騒ぎでオンライン授業の開

９４

ば

あきら

ふじ

草 場

くさ

彰

藤 井 恵 介２
０
１
９年九州大学大学院工学府航空宇宙工学専
攻博士後期課程修了．博士（工学）．現在，九州大学応用力学研究所助教．専門は半導体結晶成長の理論解析．最近は材料プロセス設計
への機械学習応用を研究しております．また，異分野連携としてプラズマ乱流のデータ解析を行っておりますので，どうぞよろしくお願いいたします．
あさ

ひ

ゆう

い

けい

すけ

京都大学工学研究科 助教．２
０１
２年 同研究科
大学院・博士（工学）．確率的な描像で様々な現象を理解・説明することをめざして，プラズ マ 物 理・原 子 物 理・量 子 物 理 の 間 を さ ま
よっています．統計的な機械学習も好きですが，（残念ながら）まだうまく使いこなせていない気がします．最近は色々なところに出てくるべき則まわりの物理を好んで調べています．

いち

朝 比 祐 一日本原子力研究開発機構 システム計算科学センター 高度計算機技術開発室 研究
員．２
０１
５年東京工業大学博士後期課程修了．博士（工学）．専門は，プラズマ乱流シミュレー
ションと High performance computing（HPC）．最近 は 主 にHPC 関連の仕事（https://github.com/yasahi-hpc など）をしています．

９５