J. Plasma Fusion Res. Vol.96, No.5 (2020)２６２‐２６６

研究・技術ノート

特異スペクトル解析法を用いた揺動成分の分離
Separation of Fluctuation Component with Singular Spectrum Analysis
坂 東 隆 宏＊，†，大 舘暁１）

BANDO Takahiro＊，† and OHDACHI Satoshi１）
総合研究大学院大学，１）核融合科学研究所（原稿受付：２０２０年２月３日／原稿受理：２０２０年３月１１日）

実験データの解析においては時系列データからトレンド成分と揺動成分を分離し解析する必要性に出くわすことがある．分離手法として広く用いられているものとしては Fourier 変換があげられる．Fourier 変換を用いる場合では，任意の時系列データが異なる振幅・周波数・位相を持つ正弦波の無限和によって表現できることを利用し，注目する周波数成分のみ分離する．一方で，Fourier 変換は定常信号の解析に有効であり，非定常信号には用いることができない．本稿では，非定常信号に適用でき，かつ，基底をあらかじめ指定しない特異スペクトル解析法（Singular Spectrum Analysis）を用いた揺動成分の分離手法を紹介する．特異スぺクトル解析法では，１次元の時系列データを１対１写像を用いて行列に変換し，その行列に対して特異値分解（Singular ValueDecomposition）を行う．あらかじめ基底を指定しないので，データ点が少ない場合などでも揺動成分の分離を柔軟に行うことができる．
Keywords:
Singular Spectrum Analysis (SSA), Singular Value Decomposition(SVD), signal processing, experimental data,
separation of fluctuationcomponent, noise reduction

１．はじめに大型ヘリカル装置 LHD において水素ペレット入射直後の実験データの解析においては時系列データからトレンド

軟 X 線信号に観測された Snake-like 構造［２］や，Helically-成分と揺動成分を分離し解析する必要性に出くわすことがtrapped-energetic-ion driven resistive interchange modeある．分離手法として広く用いられているものとしては（EIC）励起時に観測された電子温度の揺動成分にみられるFourier 変換があげられる．Fourier 変換を用いる場合でTearing parity を持つモード構造［３］の解析においては，特は，任意の時系列データが異なる振幅・周波数・位相を持異スペクトル解析法の有用性が示されている．他の分野につ正弦波の無限和によって表現できることを利用し，注目おける活用例については文献［４］などを御覧いただきたい．

する周波数成分のみ分離する．一方で，Fourier 変換は定常

本稿の第２節では，特異スぺクトル解析法の具体的な手信号の解析に有効であり，非定常信号には用いることがで順について述べる．第３節において，テスト信号およびプきない．非定常信号に適用可能な手法としては，Fourierラズマ実験で得られた軟 X 線信号に対する特異スぺクトル変換と同様に基底を元に分解する離散 Wavelet 変換を用い解析法の適用例を示す．第４節でまとめを述べると共に関た多重解像度解析があげられる．

連する文献を紹介する．

本稿では，非定常信号に適用でき，かつ，基底をあらか２．特異スペクトル解析法

じめ指定しない特異スペクトル解析法（Singular SpectrumAnalysis）を用いた揺動成分の分離手法を紹介する．特異特異スペクトル解析法では，時系列データの部分時系列

スぺクトル解析法では，１次元の時系列データを１対１写からなる Trajectory matrix を構成し，その Trajectory ma-像を用いて行列に変換し，その行列に対して特異値分解

trix に特異値分解を適用することで時系列データを分解す（Singular Value Decomposition）を行う．あらかじめ基底る．ここではまず特異値分解の性質を述べた後に，特異スを指定しないので，データ点が少ない場合などでも揺動成

ペクトル解析法を用いた分離手法の詳細について述べる．

分の分離を柔軟に行うことができる．一方で，特異値分解２．
１ 特異値分解を行うため，長大なデータ列を計算機で処理すると膨大な特異値分解とは行列に対する次の分解（式１）のことを言う．すなわち，任意の !!"実行列 !は適当な ! 次実直

計算コストがかかる点には注意が必要である．

交行列 " および "次実直交行列 # を取れば，プラズマ・核融合分野での実験データ解析への特異スペ
クトル解析法の適応例は少ない（たとえば［１］）ものの，＊†現所属は量子科学技術研究開発機構那珂核融合研究所
ORCID ID ：００００‐
０００３‐
１４９３‐
９１８５

corresponding author’s e-mail: bando.takahiro@qst.go.jp

２６２!2020 The Japan Society of Plasma
Science and Nuclear Fusion ResearchResearch and Technology Note

Separation of Fluctuation Component with Singular Spectrum Analysis!$#!$'

味する．Hankel 行列への理解を深めるため，Trajectory（１）matrix の具体的な様子を付録２に示した．
次に，Trajectory matrix %を特異値分解を用いての 形 に 分 解 で き る．!$(*%""%#"%!"')'/(に 含 ま れ る!(%!"%""%#"%/*を !の特異値と呼ぶ．ここで /$',+')/"である．#$，$$"!"
"""
"- .
#""
#..とすれば，#!"/"/"''!$%!"!#!
"%"""#""/"%/"/#'/ $!!"!""/"!/，!)$%)")#')%$%!"%""/"%/（１０）%)$%)")#')（１１）と 展 開 す る．/$%(で あ る．こ こ で，集 合
',+'

（２）

"""
/!
"
#)
/"
*か ら 任 意 の *個 の 整 数 を 取 り 出 し て 集 合)
)#$))
/"!"
""
**を作り，と展開できる．ここで，)#)*は正規直交基底と")*および )できる．

%#$%)"%)"/"%)!

特異値分解は行列 に 対 す る「自 然 な」１ノ ル ム である
Frobenius ノルム-

T. Bando and S. Ohdachi

"

（１２）

*を構成することを考える．%#が Hankel 行列であれば，ここから新たな時系列データ"#を構成することができる．

.

"' ' ()"*!,,"$( )$!*$!

しかしながら，一般に %#や )%)*は Hankel 行列ではな（３）*

い．%#$行列 %#から長さ &の時系列データを再構成する
ためには，%#を Hankel 行列に変換する必要があるため，と密接な関係がある．すなわち，%#を Hankel 行列に変換する演算子!をここで次のように.および &それぞれの成分を.$! &とし，&定義する．&.
2)"*として，)"*，2/""!,,"$' %)（４）""!),,"$%)（５）)$!

.
2!0+++",$)"*$ ' 2という関係がある．したがって，"'+,(&!0/""""!),!,$' %),,"$"$%)
)$!

（１３）

!0$)
+",+",$)&'("*$0*（６）となり，!)の !に対する寄与の大きさが %)の大きさに.はHankelは集合 !0の元の個数である．この&である．+!0+よって評価できることがわかる．式（４）の証明は文献［５］行列であり，,"!&," を最小化する Hankel 行列 " である．などを参照いただきたい．式（５）の証明は本稿の付録１に

この証明は本稿の付録３に記す．なお，演算子!は線形写記した．

像であるので，式（１２）で!を %#に作用させることは，!を２．
２ 特異スペクトル解析法による主要項と非主要項の分個々の %) に作用させて総和を取ることと同値である．
.#$! %#なる Hankel 行列化した %.#を用いれば，新たな%*離の手続き
まず，分解したい時系列データを"& $'1111!"""#/"&(

時系列データ"#を構成することができる．と置く．ここで，Lagged vector を次のように定義する．
'!)$'
11!
!
!
"1
（!%)%$）"!""%!!()"))

さて，先に述べたとおり Frobenius ノルムに関して/""""%,,and %)%),$,"$' %")（７）)$!

（１４）

ここで，!#%#&，$$&!%"!である．この Lagged vec-が成り立つ．したがって，%)は各 Elementary matrix %)tor を用いて Trajectory matrix %を次のように定義する．
!
"
111/ 1!
"#$%&%11/ 1#$$"!&%" 1&%&
11/ 1!""
!$ .
%$!!"
$%1/"
%
#$$""&（８）%&%/
00 10&#$
1
1% 1%"! 1%"" /&が %にどれだけ寄与するかを示す指標となる．主要な %)に対応する %)だけ選び，そこから新たな時系列データを構成することで，主要項と非主要項を分離することができる．本稿では主要項をトレンド成分とし非主要項を揺動成分として時系列データを分解する．
以上の手続きの中で，"& に対する Lagged vector の長さ%には任意性があるのがわかる．%によって Trajectory

Trajectory matrix は Hankel 行列となっていることに注意

matrix が変化し，したがって，特異値分解で得られる Ele-す る．Hankel 行 列 と は，行 列 要 素 ()"*が 別 の 行 列 要 素

mentary matrix が変化する．%は分離したい成分に合わせ(+",と次の関係にある任意の行列 !のことを言う．て選択する必要がある．ある周波数を持つ揺動成分につい()"*$(+", （)"*$+",）ては，分離したい揺動成分を含むデータ長以上を設定する（９）

ことが文献［６］では推奨されており，筆者の経験からもそある %が定まると，"& から Hankel 行列である Trajectoryのように設定するのが適切だと考える．

matrix %への変換は 一 対 一 線 形 写 像 と な る．し たがっ

また，時系列データに，ある程度大きい振幅を持つ揺動て，もし %がいくつかの Hankel 行列に分解できるなら，そ成分が二つ以上含まれている場合は，これらの揺動成分をれは"& がいくつかの時系列データに分解できることを意

分離できるかどうかが実際の解析では問題になりうるだろ

１ -#.行列 ! の各成分を（-#.）次元ベクトル 'に展開した時，,!," と 'のユークリッドノルムは等しく，この意味で「自然な」ノルムと考えられる．

２６３

Journal of Plasma and Fusion Research Vol.96, No.5 May 2020う．こういった分離の問題は，特異スペクトル解析法ではある．図２に（a）生成したテスト信号，（b）テスト信号を構「Separability」として議論される．いま興味があるのは，成する各項の成分を示した．このような時系列データに対二つの正弦波 !",$"+$#+()'()'"""(および !#,""#(が分しては，Fourier 変換によって揺動成分のみ取り出すこと離できるかどうかである．ここで，変数 +は時刻でありそが難しい．図３に（a）計算された特異値，（b）１番目の大き$$#かつれ以外の変数は定数である．文献［６］では，$"&さを持つ成分を用いて再構成した信号，（c）
２〜５番目の!"&$!#で あ れ ば，時 系 列 デ ー タ の デ ー タ 長 # お よ び

大きさを持つ成分を用いて再構成した信号を示した．図３

Lagged vector の長さ "を十分に長くとることで，これら（b）ではトレンド成分が，図３（c）
では揺動成分が得られて二つの正弦波を分離できることが示唆されている．Sepa-いることがわかる．この再構成結果からわかる通り，特異rability について詳しくは，文献［６］の２．
３節，文献［７］の１．
５スペクトル解析法はノイズ除去にも適用できる．

節および６．
１節にて解説されている．

３．
２ 実験データへの適用例
LHD のプラズマ実験で得られた軟 X 線信号に特異スペ

３．テスト信号および実験データへの適用例

クトル解析法を適用し，揺動成分を分離した例を示す．

３．
１ テスト信号への適用例

図４（a）にプラズマに水素ペレットを入射した直後の軟 X線信号を示した．この軟 X 線信号には，電子温度の回復にテスト信号に対する特異スペクトル解析法の適用例を２よる緩やかに増加するトレンド成分と振幅が時間的に変わ例示す．

る揺動成分の２成分がある．この時間的に変化する振幅は

例１．,$,##)
()'
##+
(

Snake-like 構造［２］をなす成分である．Snake-like 構造の空!
#%とする．この$!
,$,##)()'##+(の正弦波に対して )間的な構造を研究するためには，トレンド成分と揺動成分関数の値を!%+%$!の区間に*$"!!の点数で等間隔に与え時系列データ!# を作る．この時系列データに対して"$%!として，特異スペクトル解析法を適用する．図１に（a）作成したテスト信号，（b）１番目・２番目の大きさを持つ特異値を用いて再構成した信号，（c）全特異値を用いて再構成した信号を示す．対応する特異値は図１（d）に示した．この例では，１番目と２番目の主要な特異値で元の時系列データが再現されていることがわかる．１番目と２番目の特異値が同等の大きさになっているが，これは特異スペクトル解析法の特性である．すなわち，ある正弦波は同じ程度の寄与を持つ２つの Elementary matrix に分解される［６］．
%(例２．,$$,##)
()'
#,'#+
(#'-+'!&#+("'#+")*(%(,$$,##)()'#,'という関#+(#'-+'!&#+("'#+")*(!
!&，'$!
数 に 対 し て， )!
!
"，&，&$!
$"，$$#!，%$!
($"!
&と す る．こ の 関 数 の 値 を !%+%$!の 区 間 に*$"%!の点数で等間隔に与え時系列データ!# を作る．この時系列データに対して "$#&として，特異スペクトル解

図２ （a）
生成したテスト信号，（b）
テスト信号を構成する各項の成分．

析 法 を 適 用 す る．こ の テ ス ト 信 号 の う ち，%$,##)()'##+(#'-+'!&#+(は振幅・周波数が時間的に変(化する正弦波であり，'#+は時間的に増加するトレンドで

図３ （a）
L ＝ 28 とした時のテスト信号（例２）に対する特異値を大きい順に８つ示した．（b）
１番目の大きさを持つ特異値を用いて再構成した信号，（c）
２〜５番目の大きさを持つ特異値を用いて再構成した信号．（b），（c）には，比較のためにテスト信号を構成する成分（trend，periodic）
を示した．

図１ （a）
テスト信号，（b）
１番目・２番目の大きさを持つ特異値を用いて再構成した信号，（c）
全特異値を用いて再構成した信号．（d）
L ＝ 50 とした時のテスト信号に対する特異値を大きい順に８つ示した．

２６４

Research and Technology Note

Separation of Fluctuation Component with Singular Spectrum Analysis

T. Bando and S. Ohdachi

最後に関連する文献を紹介する．まず，Fourier 変換，を分離する必要がある．
図４（b），（c）に特異スペクトル解析法または Fourier 変Wavelet 変 換 に 対 し て は い く つ も 教 科 書 が あ る が，換によるバンドパスフィルタを用いて揺動成分の取り出しWavelet変換については［８］が網羅的で詳しい．本学会誌にを試みた例を示す．図４（c）では，観測された揺動成分の周おいても過去に［９，１０］などの解説がある．特異スペクトル波数は600 Hz〜700 Hzであるので，400 Hz〜900 Hzのバン解析法を用いたトレンド成分と揺動成分の分離に関し，日ドパスフィルタを軟 X 線信号に適用した．バンドパスフィ

本語で読める教科書的な説明は［１１］にある．より詳しい文ルタを行う際には，時間方向の窓関数として，始端付近の献としては［６，７，１２］があげられる．特異値分解のコンパク揺動成分の変化を抑えるため方形波を用いている．特異ストで教科書的な説明は［５，１３］にある．特異値分解は多次元ペクトル解析法を適用した図４（b）の場合では上手く揺動データに対しても利用でき，画像データに対する特異値分成分のみを取り出せているが，バンドパスフィルタを適用解を用いた揺動成分の分離手法が本学会誌にて解説されてした図４（c）の場合では，Gibbs 現象と思われる振動が発生

い る［１４］．特 異 ス ぺ ク ト ル 解 析 法 は，１次 元 デ ー タ をしており時系列データの始端と終端における揺動成分の取Trajectory matrix に変換して特異値分解を適用するとい得に失敗していることがわかる．

う点で，多次元データに対する特異値分解を１次元データに適用した例とみることもできよう．

４．まとめ
謝

本稿では，特異スペクトル解析法を用いたトレンド成分

辞と揺動成分の分離方法およびその適用例について述べた．

本稿は筆者（坂東 隆宏）の博士論文［１５］の第５章の一部

特異スペクトル解析法を用いることで，周波数・振幅が時を修正，加筆したものである．博士論文を審査いただいた

間的に変動している揺動と，緩やかに変動しているトレン

稲垣 滋博士，榊原 悟博士，武智 学博士，藤堂 泰博士にド成分とを分離することができる．Fourier変換やWavelet感謝いたします．また，本稿の内容の一部は核融合科学研

変換と異なり，特異スペクトル解析法は決まった基底を用究所の一般共同研究 NIFS07KLPH004 による支援を受けたいないので様々な信号に対して柔軟な分解を行うことがでものです．

きる．一方で，特異値分解を行うため，長大なデータ列を処理すると膨大な計算コストがかかる点には注意が必要で参考文献ある．テスト信号および実験的に得られた軟 X 線信号に特［１］L. Marrelli et al., Rev. Scientic Instrum. 72, 499 (2001).
［２］T. Bando et al., Phys. Plasmas 25, 012507 (2018).
［３］T. Bando et al., Nucl. Fusion 58, 082025 (2018).
［４］A. Zhigljavsky, Statistics and Its Interface Volume 3, 255(2010).
［５］金谷健一：線形代数セミナー：射影，特異値分解，一般逆行列（共立出版，２０１８）．
［６］N. Golyandina and A. Zhigljavsky, Singular Spectrum
Analysis for Time Series (Springer, 2013).
［７］N. Golyandina et al., Analysis of Time Series Structure: SSAand Related Techniques (Chapman and Hall/CRC, 2001).
［８］新 誠一，中野和司：図説ウェーブレット変換ハンド

異スペクトル解析法を用いることで，緩やかに変動するトレンド成分と振幅・周波数が時間的に変動する揺動成分とを分離可能であることを示した．
このように特異スペクトル解析法は柔軟に分離を行うことができるが，実際のところ，Fourier 変換・Wavelet 変換で事足りてしまうことが多い．しかしながら，Fourier 変換・Wavelet 変換で対応できないような例（たとえば，利用できるデータ点が少なく，トレンドが緩やかに変化し，揺動成分の振幅・周波数が時間的に変動する場合［２，３］）では特異スペクトル解析法は有効に働く．

図４

軟 X 線 信 号 へ の 特 異 ス ペ ク ト ル 解 析 法・Fourier 変 換（FFT）の適用例．（a）
水素ペレット入射直後の軟 X 線信号．
（b）
軟 X 線信号に特異スペクトル解析法を適用し揺動成分の取得を試みたもの．（c）
軟 X 線信号に Fourier 変換によるバンドパスフィルタを適用し揺動成分の取得を試みたもの．

図５ （a）
Trajectory matrixの第２行（row2）および第２列（clm2）の成分．（b）
Trajectory matrix の第２１行（row21）および第２１列（clm21）の成分．

２６５

Journal of Plasma and Fusion Research Vol.96, No.5 May 2020

ブック（朝倉書店，２００５）．
［９］辻 義之 他：プラズマ・核融合学会誌 85, 620 (2009).
［１０］犬塚 博：プラズマ・核融合学会誌 74, 1023 (1998).
［１１］馬杉正男：信号解析‐信号処理とデータ分析の 基 礎‐
（森北出版，２０１３）．
［１２］H. Hassani, J. Data Science 5, 239 (2007).
［１３］池辺八洲彦 他：現代線形代数‐分解定理を中心として‐（共立出版，２００９）．
［１４］西野信博，大舘 暁：プラズマ・核融合学会誌 86, 648(2010).
［１５］坂東隆宏：「LHDにおける重水素プラズマのMHD不安定性に関する実験研究」総研大物理科学研究科博士論
文，２０１
８年３月．
（総合研究大学院大学学術情報リポジ
トリの Web サイトよりダウンロード可能）．

図５では，（a）Trajectory matrix の第２行および第２列の成分，（b）Trajectory matrix の第２１行および第２１列の成分を示し た．&#% 行 列 'および そ の 行 列 要 素 1*"+に つ い て，，
では '"$'&&"%(とした時，図５（a）
1
11!
!
!
"1('#"
""
#"
#"
#"
$"
#"
"(を 順 に 示 し，図５（b）で
1
11!
!
!
"
1'""
#"
#"
#"
$"
#"
#(""は，'
1
111!
!
!
"
111!
!
!
"
1'
#""
""
#""
#"
#""
$"
""
#""
#"
#""
$"
#""
#"(を順に示#"""""(している．図５（a），（b）それぞれで，行と列は同じデータ列である こ と が わ か る が，こ れ は Trajectory matrix がHankel 行列であり 1*"+$,"(が満たされている,"-'*"+$1ためである．

行 列 'に 対 し ,$!',#を最小とする
Hankel 行列 $ の導出

付録．３'を .#/実行列とする．ここでは Hankel 行列 $ の中

付録．１!
( #!
(#" #"! の証明

で,$!',$!',# を最小とする $ を求める．,# は常に正#なので，この問題は,$!',#を最小とする $を求める問題"を .#/実行列とする．"*の定義と Frobenius ノルムと同値である．以降，$ と 'の成分をそれぞれ )*"+，1*"+の定義から
.

/

.

/

#
#
#
#!*"'
"*,
# #'
#*!*"'((,#$# # '"-$#"* ,** ,,$"-$",$"-$"と書く．ここで，

（１５）

.

/

##$!',$$,),"-!1,"-(#$# # 'である．さて，（２１）,$"-$"と す る．求 め る $ は，$ の 最 小 値 を 与 え る)*"+'"%*%.""%+%/(を求めることで得られる．
さて，同じ 0$*"+をとる )*"+を含む項をまとめると，$$# $0（２２）0#$0')*"+(),"-!1)*"+!1$ # '３）,"-($ # ',"-( （２（１６）"',(&!0/.と書ける２．各々の $0')*"+(を最小とする )*"+'*"+$0(を求/###!*"'!*("*(# #'("-$# # ','* ,,$"-$",$"-$"（１７）

めることで，$ を最小とする $ が得られる．とできる．ここで，% と & の各列は正規化されたベクトル

そこで，なので，.

#!*(#',$" and,$"($0(#) !1 ($!# '"$!
()*"+ ()*"+ ',"-(&!0 *"+ ,"-/#"*(#'-$"-$"（１８）であるので，式（２４）で得られる )*"+により $0の最小値が得/##!*("*(# #','-$",$"-$"（１９）られる．計算を進めると，となる．この結果と式（１５）から，求める関係式##"*,,#$#*(#) !1 ()*"+! # 1$!
!# '"$#+!!0+,"-"()*"+ ',"-(&!0 *"+ ,""',(&!0（２５）（２０）となる．ここで +!0+は集合 !0の元の個数である．ゆえに，が示された．

付録．２)*"+$ " # 1,"!0++"',(&!0Hankel 行列である Trajectory matrix の具体的な様子（２６）が一意に得られる．式（２６）で求めた )*"+を成分として持つHankel 行列である Trajectory matrix の具体的な様子を本$ により $ は最小となる．

稿３．
１節の例２で作成する Trajectory matrix の場合で示す．

２（２４）を満たす )*"+'*)*"+(は下に凸の関数"+$0(を求める．$0'であることを用いると式（１７）は，.

"',(&!0!0$),"%,"-$*'("+$0*であるので，.

#$ は Hankel 行列であるので，)*"+$),""+$,"(であることに注意されたい．
*-'

２６６