修士論文
すばる望遠鏡 RAVEN 観測データによる大気揺らぎトモ
グラフィック測定の精度限界の評価

Evaluation of accuracy limit in atmospheric turbulence
tomographic measurements based on Subaru RAVEN dataset
東北大学大学院

理学研究科
井上欣彦
令和元年

天文学専攻

修士 2 年

概要
次世代の高精度レーザー補償光学において複数のガイド星と波面センサの測定データからトモグラ
フィーの手法を用いて、ターゲット天体方向の波面を推定する手法の開発・研究は重要となっている。波面推定法としては大気揺らぎ高度プロファイルに基づくモデルを用いたものをはじめとして、
オンスカイ観測データから直接推定するものなどが提案されている。その中で観測データに基づいた直接的な波面推定法はモデルで考慮されていないパラメータの有無を検証するうえで重要である。
一方で、データに基づく波面推定に関する基礎研究は無い。
本研究では、すばる望遠鏡 RAVEN のオンスカイ観測で得られたデータを用いて、最小二乗法によって複数のガイド星の波面センサの測定データからターゲット方向の波面センサの測定データを結びつけるトモグラフィック推定行列 R の推定を行った。まず、R の推定の際に推定に寄与する成分とそれ以外の成分を分けるために特異値分解を行い低特異値の除去を行なった。その中でこれまでの特異値除去の基準に問題があったことを明らかにした。次に、データ数が少ない時に R が推定した時刻の特徴を過学習してしまう問題を L2 正則化の手法を用いて改善することができた。
さらにトモグラフィーの手法を用いた波面推定の推定精度の限界を決める要因として推定に寄与する大気揺らぎの情報量があることがわかった。その中で GS の数と配置が重要であると考え、議論を行った。その結果、情報量を増やすために GS の数を増やすこと、GS の配置をターゲットから近い距離に方向をずらしておくことが重要であることが本手法によっても確認された。また、タイムス
テップデータを用いた推定法の提案を行った。その結果、過去の測定値を反映させ実効的に推定に寄与する情報量を増やすことで推定精度を向上させることができた。最後に最終的に得られる波面推定精度に関する評価を行った。

目次
第 1 章 イントロダクション

1.1

1.2

1.3

3

補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.1.1

原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.1.2

単層共役補償光学とその課題. . . . . . . . . . . . . . . . . . . . . . . . . .

5レーザートモグラフィー補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2.1

原理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.2.2

サイエンスケースと要求性能. . . . . . . . . . . . . . . . . . . . . . . . . .

6

波面再構成

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3.1

シャックハルトマン型波面センサの測定値 . . . . . . . . . . . . . . . . . . .

7

1.3.2

波面センサの測定値と波面位相の関係. . . . . . . . . . . . . . . . . . . . .

8

1.3.3

トモグラフィーモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.3.4

課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

1.4

トモグラフィック測定データに基づく波面推定法 . . . . . . . . . . . . . . . . . . . .

13

1.5

本研究の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

第 2 章 多天体補償光学実証装置 RAVEN

14

2.1

装置仕様 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.2

RAVEN のデータセット . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.2.1

WFS の測定値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.2.2

データ行列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2.3

観測コンディション . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

第 3 章 データに基づく波面推定法

22

3.1

推定法の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

3.2

計算手法の検討 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.2.1
3.2.2

最小二乗法を用いた推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
低特異値の除去 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
25

3.2.3

L2 正則化を用いた推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

タイムステップデータを用いた推定法の検討 . . . . . . . . . . . . . . . . . . . . . .

26

3.3

第 4 章 結果

4.1

4.2

4.3

27

最小二乗法を用いた波面推定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

4.1.1

最小二乗法を用いた波面推定の結果 . . . . . . . . . . . . . . . . . . . . . . .

27

4.1.2

波面センサの測定値について. . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.1.3

静的成分を除いた波面推定の結果 . . . . . . . . . . . . . . . . . . . . . . . .

31

特異値除去を行った時の結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.2.1

最大特異値での基準 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.2.2

累積寄与率での基準 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

L2 正則化による波面推定. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

34

4.4

4.3.1

正則化の項 λ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

4.3.2

L2 正則化による結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

タイムステップデータを用いた波面推定の結果 . . . . . . . . . . . . . . . . . . . . .

36

第 5 章 議論

37

5.1

平均値に現れるパターンについて . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

5.2

特異値除去の基準 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.1

最大特異値の基準 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41

5.2.2

累積寄与率の基準 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

5.2.3

基準の比較 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.3

正則化と過学習 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44
44

5.4

推定精度を決める要因 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.4.1

光学系パラメータ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.4.2

ガイド星の数とガイド星の配置 . . . . . . . . . . . . . . . . . . . . . . . . .

50

5.4.3

タイムステップデータを用いた推定法との関係 . . . . . . . . . . . . . . . . .

56

5.5

推定誤差の評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.6

トモグラフィック推定行列 R の時間依存性 . . . . . . . . . . . . . . . . . . . . . . .

62

第 6 章 結論

69

AppendixA71

A.1 行列の最小二乗解の導出 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

A.2 L2 正則化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

A.3 ゼルニケ多項式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

A.3.1 ゼルニケ多項式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

A.3.2 波面収差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

74

2

第 1 章 イントロダクション

1.1

補償光学

天体からの光は天体が十分遠くにあるため平面波となって地球に到達する。しかし、地上での観測では大気揺らぎによってこの波面が歪んでしまい星像もぼやけてしまう。その結果、感度や空間分解能を損ねてしまう。この大気揺らぎによって歪んだ波面をリアルタイムで補正して、望遠鏡本来の性能 (回折限界) を実現するためのシステムが補償光学 (Adaptive Optics: AO) である。すばる望遠
鏡では 2000 年 12 月に AO がファーストライトをむかえ、それ以降 AO を用いた観測により多くの成果を得ることができている [19]。さらに現在マウナケア山頂に建設が予定されている 30m 級望遠鏡 (Thirty Meter Telescope: TMT) をはじめとする口径が数十メートルを超える次世代超大型望遠鏡 (Extreme Large Telescope: ELT) ではその高い空間分解能を活かすために、この AO が必要不可
欠となってくる。この ELT に向けた AO が世界中で研究、開発されている。以下にすばる望遠鏡で観測した２枚の星像をのせた。左が AO を使用していない時、右が AO を使用したとき。AO によって大きく星像が改善していることがわかる。

図 1.1: AO による星像の改善例。左：AO を使用していないとき。右：AO を使用したとき。観測波
長は K バンド (2.2µm)。(Subaru Telescope HP[19])

1.1.1

原理

補償光学のシステムは波面センサー（Wave Front Sensor: WFS）、リアルタイムコンピューター、
可変形鏡 (Deformable Mirror: DM) の３つの基本的な構成に分けることができる。

3

第 1. イントロダクション

1.1. 補償光学

波面センサー
波面センサーでは望遠鏡に入ってくる光の位相（波面）を測定する。観測対象天体が暗い場合には観測対象天体自体を参照光源として用いることができない。そのため通常は波面を測定するための参照星が必要となってくる。この参照星をガイド星 (Guide Star: GS) と呼ぶ。GS には自然にある天体を GS として用いる自然ガイドスター (Natural Guide Star: NGS) と、レーザーを打ち上げて人工的に仮想の星をつくりそれを GS として用いるレーザーガイドスター(Laser Guide Star: LGS) がある。LGS としては高度 90km 付近のナトリウム層に波長 589nmのレーザーを照射してナトリウムを励起させて発光させるものと高度 20km 付近までの大気のレイリー散乱を利用するものがある。
リアルタイムコンピュータ
波面センサーでの測定値を基にして、観測ターゲット方向の大気揺らぎの 2 次元マップの推定を行う。このプロセスを波面再構成と呼ぶ。また推定した波面の歪みを打ち消すように可変形鏡を制御する。
可変形鏡
薄い鏡面からなる可変形鏡によって乱れた波面の補償を行う。推定された波面が反射された時に、平面波になるように薄い鏡面を変形（各 element の傾きが推定された波面の傾きの 1/2 になるように）させる。鏡の変形は鏡の裏に取り付けられたアクチュエーターと呼ばれる素子に電圧をかけて伸縮させることで行う。
大気揺らぎの変動のタイムスケールは数百〜数千 Hz なので、これらの補正を同程度のタイムスケー
ルでリアルタイムに絶え間なく行う必要がある。このプロセスによって大気揺らぎによって乱れた波面を補償する。以下の図 1.2 に補償光学システムの概念図を示す。

図 1.2: AO 装置の概念図

4

第 1. イントロダクション

1.1.2

1.2. レーザートモグラフィー補償光学

単層共役補償光学とその課題

単層共役補償光学 (Single Conjugate Adaptive Optics: SCAO) は図 1.2 のように１つの GS と 1
つの WFS によってサイエンスターゲット方向の光の波面の乱れを測定、補償するシステムである。

2000 年以降 SCAO はすばる望遠鏡の AO188 と呼ばれる補償光学系をはじめとした 8〜10m 級の望遠鏡に用いられ、多くの成果をあげてきた。その一方で SCAO では、観測ターゲットの近くに十分
明るい NGS がなくてはいけないといった制約や、回折限界に近い補正が達成できる波長が近赤外領域に限られてしまうといった制約がある。そのため SCAO による補償が行えるのは天球面上の数％の領域に限られてしまう。これらの課題は具体的には以下のような要因によるものである。
コーン効果
一つの LGS では LGS 自体が数十ｋｍの高度しかないため、図 1.2 で示したようにターゲット方向からの光で GS からの光が通らない大気の領域ができてしまう。この部分の大気の情報は大気測定にのらない。そのため、波面の補償精度の低下が起こってしまう。この問題はコーン効果と呼ばれ、望遠鏡の口径が大きいほど影響も強くなる。
アンアイソプラナティズム

NGS がターゲットと非常に近い領域にあるときには、この NGS によってターゲット方向の位相の乱れを直接測定することができる。しかし、ターゲットから離れてしまうと異なる大気揺
らぎを通過してしまうため補償性能が低下してしまう。この問題はアンアイソプラナティズムと呼ばれ、SCAO ではこの大気揺らぎの補償が効く視野は近赤外線で 30 秒角程度であると言われている。これは大気揺らぎが同じとみなせる領域がこの 30 秒角程度であるためである。

1.2

レーザートモグラフィー補償光学

1.2.1

原理

レーザートモグラフィー補償光学 (Laser Tomography Adaptive Optics: LTAO) は SCAO における課題であったコーン効果やアンナイソプラナティズムの問題を、トモグラフィーの手法によって取り除き高空間分解能での観測を可能にするシステムである。

SCAO では１つのガイド星と波面センサを用いて補償を行なっていたが、LTAO では図 1.3 のように複数のレーザーガイド星と波面センサを用いる。従来の SCAO では得ることのできなかった大気の情報を含んだ複数の測定データをもとに大気の立体構造を推定し、観測ターゲット方向の波面推定および可変形鏡を用いた補償を行う。この複数の測定値から大気揺らぎの３次元構造を回復するプ
ロセスをトモグラフィック波面再構成と呼ぶ。

5

第 1. イントロダクション

1.2. レーザートモグラフィー補償光学

図 1.3: レーザーガイド星を用いた AO 装置の概念図

1.2.2

サイエンスケースと要求性能

LTAO によって期待されるサイエンスとしては銀河形態進化がある。
銀河は、見かけの形状や明るさに基づいて多くの形態に分類されており、ハッブル分類では巨大銀河は、楕円銀河、レンズ状銀河、渦巻銀河、不規則銀河、特異銀河に分類される。銀河の進化を説明するものとして、階層的集団化モデルが提案されている。これは質量の小さい銀河が初めに生まれ、それらが集まることで大質量の銀河が作られていくというモデルである。しかし、この銀河進化のモデルには観測との矛盾や物理的に不明確な部分が残っている。これを明らかにするためにはより多くの銀河の観測によって、銀河中のガスの内部運動、ガスの金属量や星の分布などを調べることが必要不可欠である [21]。
近年の近赤外線観測によると、星形成率密度は z=0〜1 の時代には増加をし、それ以降の z=1〜3の時代にはその増加が止まっていることがわかっている（図 1.4）[9]。このことから、銀河における星形成は z=0〜1 の時代に活発に行われており、その後銀河の形態が確立していったと考えられる。
したがって銀河の形態進化の解明のためには、銀河が最も激しく進化していると考えられるこの z=1〜3 の時代にわたる銀河内の星の速度構造などを知る必要がある。そのために、銀河内の星からの光の輝線や吸収線を空間分解して観測する面分光装置の導入が望まれる。
しかし、これらの銀河は非常に暗いため観測するためには大きな口径の望遠鏡での観測が必要となる。また、遠方にあるため高空間分解能で観測する必要がある。そこで、LTAO による高空間分解能が期待されている。
上記の銀河の形態進化の歴史を探るためのプロジェクトとして TMT で近赤外広視野多天体分光装置 (InfraRed Multi-Object Spectragraph: IRMOS) の導入が予定されている。
この TMT-IRMOS のサイエンスからの要求では z ＝ 2 ー 4 の銀河の構造を分解して観測するためには、H バンドでエンサークル・エンスクエアエネルギーで 50%以上（0.05”×0.05”）が必要となってくる [4]。これは SR 比で 0.5、波面誤差で 211nm(H バンド) に対応する [10]。

6第 1. イントロダクション

1.3. 波面再構成

図 1.4: 星形成密度 (Hopkins et al,2006[9])

1.3

波面再構成

本節では WFS の測定値から大気揺らぎの構造を推定する波面再構成の手法について説明する。

1.3.1

シャックハルトマン型波面センサの測定値

波面センサには様々なものがあるが、ここでは今回の波面測定で使われたシャックハルトマン型波面センサについて説明する。
シャックハルトマン型波面センサ (Shack Haltman Wave-Front Sensor: SH-WFS) では入射瞳に置いたマイクロレンズアレイによって入射光を分割する。理想的な平面波が入射した際にそれぞれの小開口を通って結像する位置（スポット）を基準として、乱れた波面が入ってきたときのスポットの移動量を測定する。この波面の位相の乱れから、GS からの光が通ってきた大気の揺らぎの情報を得ることができる。このスポットのずれは波面の傾きに対応しており、全小開口での傾きを測定すること
で波面の乱れを再構成することができる。

7

第 1. イントロダクション

1.3. 波面再構成

図 1.5: SH-WFS による波面センシングの概念図。平面波が入ってきたときに各小開口によって光が CCD に結像するスポットを基準点として、歪んだ波面の各小開口上でのスポットの移動量を測定
する。

1.3.2

波面センサの測定値と波面位相の関係

図 1.6: SH-WFS の測定値と波面位相点の幾何学的配置 (Fried Geomtry)

図 1.7: SH-WFS の測定値と波面位相点における位相の関係。幾何学的に関係を求める。

8

第 1. イントロダクション

1.3. 波面再構成

x 方向のスポットの移動量 sxi,j と望遠鏡開口上での波面位相 φ の関係を求める。位相の変化 ϕ[rad]と光学的距離の変化 d[m] の間には、ϕ = 2πd/λ という関係がある。そこで、波面位相 φ を基準となる平面波からの高さ方向のずれと考えると、ｘ方向についてスポットの移動量 sxi,j とそれに関係する
位相点における波面位相 φ との関係は図 1.7 から次のように表せる。

d tanθ =φi+1,j + φi+1,j+1
φi,j + φi,j+1
−
2
2

(1.1)

次にスポットの移動量 sxi,j と θ の関係を求める。ピクセルスケールが 0.38”/pix なので図 1.7 より

θ = 0.38” · sxi,j = 0.38 ·1π·· sx [rad]3600 180 i,j(1.2)という関係が求まる。今、θ は微小量なので tanθ ∼ θ の近似を用いるとtanθ ∼ θ = F · sxi,j(F =0.38 · π)
3600 · 180(1.3)

式 1.2 と式 1.3 より

1(φi+1,j + φi+1,j+1 − φi,j − φi,j+1 )2dFが求まる。y 方向についても同様にして
sxi,j =

syi,j =

1
(φi,j+1 + φi+1,j+1 − φi,j − φi+1,j )
2dF

(1.4)

(1.5)

が求まる。
これを拡張して、以下の図 1.8 ような N ×N 個のスポット測定点とそれに対応する (N +1)×(N +1)
個の位相点がある場合を考える。

図 1.8: N × N のサブアパーチャーを持つ SH-WFS。

このときスポットの移動量 si,j と波面位相 φ を結び付ける関係行列は以下のような形になる。

9

第 1. イントロダクション

1.3. 波面再構成







 x 
−1 − 1 0 . . . 0 1 1 0
...
0 
s1,1



..
 . 


.
 .. 






..
 x



 sN,N .

= 1 
 sy  2dF ..
 1,1 


.
 . 


 . 


.


 . 
..



y

sN,N
0
...
0 − 1 1 0...0 − 1 1





φ1,1
φ1,2
..
.
φ2,1
φ2,2
..
.
φN,N
φN,N +1..
.
φN +1,N



























(1.6)

φN +1,N +1
測定値ベクトルを s = [sx1,1 , · · · , sxN,N , sy1,1 , · · · , syN,N ]T と定義して、s と φ を結びつける (1.6) 式の右辺の左の行列を G と表すことにすると (1.6) 式は(1.7)

s = Gφと書ける。今スポットの移動量 S は測定値として与えられる。また、行列 G は上で求めたように

WFS の測定値と波面位相点における幾何学的な関係によって与えられる。これらを用いて波面位相点における波面 φ を最小二乗法によって以下のように求めることができる。
（行列の最小二乗解については Appendix 参照）

φ = (GT G)−1 GT s

1.3.3

(1.8)

トモグラフィーモデル

ここまで、従来の SCAO のように一つの GS と WFS の測定値から波面位相を再構成する手法について述べてきた。以下では大気揺らぎと複数の WFS の測定値を結びつけるトモグラフィーモデル
について Martin, 2016[12] を参考に記述する。
まず、トモグラフィーモデルでは大気揺らぎを無限に薄い独立な層の足し合わせで表現できると考える。そこで、Nlayer 個の大気揺らぎの層が異なる高度に存在していると仮定して、NGS 個の GSの測定値から各層の大気揺らぎ Φl (1 ≤ l ≤ Nlayer ) を計算する。
今 WFS は望遠鏡の主鏡に共役にさせているため、WFS は開口上にあると考えてよい。p 番目のWFS に着目して、この WFS の i 番目の副開口における波面位相を φip と表すと、φip は GS 方向の各大気層における i 番目の副開口の位置に対応する大気揺らぎの足し合わせとして表現できる。l 番
目の大気層における i 番目の副開口の位置に対応する座標を δipl 、i 番目の WFS の開口上での座標をri とすると δipl は{δipl =hl
αp hl + ri (1 − hLGS
)

(LGS)

αp hl + ri

(N GS)p(1.9)と書ける。ここで LGS の高度を hLGS
、大気揺らぎの l 番目の層を高度 hl 、GS が中心軸から αp = (αp,xp、αp,y ) だけずれた方向を向いているとした。また |αp | ≪ 1 として tanαp ≈ αp の近似を用いた。LGSは有限な高さにあるため位相点の間隔が狭くなっている。

10

第 1. イントロダクション

1.3. 波面再構成

図 1.9: l 番目の大気揺らぎ層と p 番目の WFS の測定値の関係。

変数

説明

Φl

l 番目の大気層における大気揺らぎ

φpp 番目の WFS の望遠鏡開口上での大気揺らぎ

ri

i 番目の副開口の望遠鏡開口上における座標αpp 番目の GS の中心軸からの角度のずれ

hll 番目の大気層の高度

δipl

p 番目の WFS の i 番目の副開口の l 番目の大気層に対応する座標

hLGS
p

LGS の高度表 1.1: トモグラフィーモデルにおける変数

この δipl を用いると φip は以下のように記述できる。

∑

Nlayerφip =

Φl (δipl )(1.10)

l=1

次に、l 層目の大気揺らぎ Φl から αp 方向の揺らぎを切り出す行列を Apl と定義すると p 番目のWFS の開口上での揺らぎ φp は∑Nlayerφp =Apl Φll=111(1.11)

第 1. イントロダクション

1.3. 波面再構成と書ける。ここで、関係式 (1.7) を用いて p 個の WFS の測定値 sp と波面位相 φp を結びつけると、∑Nlayer

sp = Gp φp = Gp

Apl Φl(1.12)

l=1と書ける。式 (1.12) より WFS の測定値と大気揺らぎの関係式は以下のように導出される。








s1
s2
..
.





 
 
=
 
 

sNGS

G1 A11G1 A12···G1 A1Nlayer

G2 A21G2 A22···G2 A2Nlayer···

GNGS ANGS Nlayer..
.
GNGS ANGS 1 GNGS ANGS 2








Φ1
Φ2
..
.








(1.13)

ΦNlayer

式 (1.13) 右辺第１項を G′ と定義し、NGS 個の WFS の測定値を並べた行列を S = [s1 , s2 , · · · , sNGS ]および大気揺らぎを表す行列を Φ = [Φ1 , Φ2 , · · · , ΦNGS ] と定義すると、WFS の測定値と大気揺らぎの関係は最終的に以下の式で記述される。

S = G′ Φ(1.14)

ここから大気揺らぎ Φ は以下の式を用いて求めることができる。

Φ = (G′ G′ )T G′T S

1.3.4(1.15)

課題

LTAO における課題としてこのトモグラフィーアルゴリズムの開発・研究がある。望遠鏡の大型化や MOAO の導入に伴い、鏡やセンサーの大きさや数は増えている。そうした中で、望遠鏡のデザインの問題はもちろんのこと、大量のデータに対してトモグラフィー計算を行う新しいアルゴリズムの開発・研究が必要不可欠になっている [7]。
その中でも実際の観測下における状態を反映したモデルを立てられるかが課題となってくる。Martin,

2017[13] で述べられているトモグラフィーモデルで用いられるパラメータのリストを以下にのせた。
• Nlayer : 独立な大気層の数• hl : 大気層の高度• r0 (hl ) : 大気層の強さ• L0 (hl ) : アウタースケール• ガイド星配置などの観測パラメータ• DM や WFS などの光学系パラメータ• 各 WFS の測定ノイズ
トモグラフィーモデルではこれらのパラメータを用いた計算が行われるが、Frabrice.2010[2] の中で述べらているように実際のシステムにおける DM や WFS の光学的な関係や装置由来の光学収差のモデルを立てることは非常に難しい。トモグラフィーモデルでは、このモデルと実際のシステムとのズレをいかにしてなくすかが大きな問題となってくる。また、トモグラフィーモデルにおける大気
揺らぎに関連するパラメータが実際のオンスカイ観測下での状態を反映できているのかを検証する必要がある。

12

第 1. イントロダクション

1.4

1.4. トモグラフィック測定データに基づく波面推定法

トモグラフィック測定データに基づく波面推定法

本論文の中でデータに基づく波面推定とは、実際のオンスカイ観測で得られた複数の GS 方向の測定データとターゲット方向の測定データを線形関係で結びつけて、GS 方向の測定値から直接的にターゲット方向の波面を推定する手法を指すこととする。
この測定データに基づいた手法の利点としては、光学系由来の収差や装置内部の位置関係などのキャリブレーションが可能であるということが挙げられる。また、大気揺らぎのプロファイルなどのパラメータはデータの中に内包されているため明示的にこれらのインプットを与える必要がない。
以上は Osvorn et al, 2012[8] でも述べられている。

Yamazaki, 2015[6] の中でデータに基づく波面推定が従来のモデルに基づいた手法よりも良い推定
結果を残すことがあることがわかった。また、オンスカイ観測データに基づいた手法では従来のトモグラフィーの手法では考慮されていなかったオンスカイの特性を反映したパラメーターの存在が示唆された。
しかし、実際の観測下ではターゲット方向からの光は弱いためターゲット方向の測定値を得ることはできない。そのため、データに基づく波面推定の目的としては、上記のトモグラフィーモデルにおいて装置のキャリブレーションができていない問題や、実際の観測下での大気のパラメータなどを抽出できていない問題を明らかにするための位置付けと考えている。また、観測データに基づいた波面
推定でどれほどまで推定精度を上げられるかを検証することによってモデルに基づいた推定における精度限界の一つに指標になると考えている。
波面推定法の詳しい内容については３章で述べる。

1.5

本研究の目的

本研究では上記の研究を踏まえ、すばる望遠鏡 RAVEN による大気揺らぎ測定データを用いて、波面推定法の検討を行い、推定精度を向上させる方法を提案し推定精度限界の評価を行うことを目的とする。これによって実際のトモグラフィー推定における精度限界の一つの指標を提示するとともに、従来の手法と比較することでこれまでの推定で考慮されていなかったパラメータについての議論
を行う。

13

第 2 章 多天体補償光学実証装置 RAVEN
RAVEN はビクトリア大学の AO チームによって設計された多天体補償光学 (Multi-Object Adaptive
Optics: MOAO) の実証装置で、2014 年 5 月、2014 年 8 月、2015 年 6 月にすばる望遠鏡のナスミス
焦点に取り付けられてオンスカイ観測が行われた。RAVEN では３つの NGS と１つの LGS を用いて、２つの独立なサイエンスチャンネル１および２に対してトモグラフィック波面再構成の試験を行い、大気揺らぎの推定が問題なく行えているかを実証することが目的とされた。

2.1

装置仕様

RAVEN の光学設計および光学ブロック図を以下に示した。

図 2.1: RAVEN の光学設計 (Lardière et al., 2012[11])14

第 2. 多天体補償光学実証装置 RAVEN

2.1. 装置仕様

図 2.2: RAVEN の光学ブロック図 (Lardière et al., 2012[11])これらのシステムを Lardière et al., 2012[11] を参考に以下で詳しく紹介する。

• キャリブレーションユニット
キャリブレーションユニット (Calibration Unit: CU) では大気揺らぎを生成して AO システムのテストおよびキャリブレーションを行う。2.4′ の視野に渡って 7×7 の配置の星および軸上のLGS を作り出す。
• 開ループ波面センサー
RAVEN では３つのピックオフミラーが 3.5′ の視野を自由に可動できるようになっており、３方向の NGS からの光を導入して開ループ波面センサー（Open-Loop WFS: OL-WFS）を用いた測定を行う。本論文ではこれら３つの OL-WFS を OL-WFS1、OL-WFS2、OL-WFS3 と呼ぶことにする。

OL-WFS の仕様を以下の表にまとめた。
有効ピクセル

128×128

ピクセルサイズ24×24µm

副開口

10×10

副開口あたりの有効ピクセル

12×12

ピクセルスケール

0.38′′ /pix

視野

4.8′′

表 2.1: OL-WFS の仕様• 閉ループ波面センサー
閉ループ波面センサー (Closed-Loop WFS: CL-WFS) は２つのサイエンスチャンネルからの光を測定する。これら CL-WFS を用いて AO システムのキャリブレーションや観測中に DM

15

第 2. 多天体補償光学実証装置 RAVEN

2.2. RAVEN のデータセットの形状の測定を行う。本論文ではこれら 2 つの CL-WFS を CL-WFS1、CL-WFS2 と呼ぶことにする。

• サイエンスリレー
RAVEN は２つのサイエンスチャンネルを持っており、ρーθ ステージに取り付けられた２つのピックオフミラーによってそれぞれのサイエンスチャンネルに光が導入される。それぞれのサ
イエンスチャンネルは望遠鏡の焦点に共役な DM、ビーム合体装置、イメージローテーターから構成されている。DM は 11×11 のアクチュエーターで制御されており、ビーム合体装置により２つのサイエンスチャンネルからくる補正された光が赤外線撮像分光装置 (Infrared Camera

and Spectrograph: IRCS) に送られる。

2.2

RAVEN のデータセット

2.2.1

WFS の測定値

RAVEN では３つの自然ガイド星と１つのレーザーガイド星のデータを用いたトモグラフィーを行った。その補償精度の実証のために２つの天体に対してこのトモグラフィーを行なった。今回は、３つの自然ガイド星方向の OL-WFS1〜3 および２つのサイエンスターゲット方向の CL-WFS1〜2のデータを用いた。

RAVEN では SH-WFS が採用されており、SH-WFS は 10×10 の小開口によって分割されている。
SH-WFS は望遠鏡の開口および副鏡によって図 2.3 のようにマスクされている。そのため、実際には図 2.3 に示した 72 個の小開口上での測定値が記録されている。図 2.3 における番号は各 WFS のデータ記録順を表している。この際、１つの小開口に対してそれぞれ x 方向、y 方向のスポットの基準点からの移動量が記録される。

図 2.3: RAVEN のマスク形状とデータ記録順

16

第 2. 多天体補償光学実証装置 RAVEN

2.2. RAVEN のデータセット

p 番目の WFS における測定値ベクトルは以下のような形になっている。



sx1
 . 
 .. 

 sx72  − srefsp = p sy1 
 . 
 . 
 . 
sy72

(2.1)

ここで、srefは基準ベクトルを表しおり、あらかじめ CU で平面波を作り出した時のスポットのp
位置でキャリブレーションが行われている。

2.2.2

データ行列

RAVEN のオンスカイ観測で得られたデータセットを図 2.4 に示す。RAVEN のデータの 1 行目から 72 行目までには OL-WFS1 の x 方向のスポットの移動量が pix 単位で図 2.3 の番号の順番で記録されている。その次の 73 行目から 144 行目には OL-WFS1 の y 方向のスポットの移動量が同様の順
番で記録されている。その次からは OL-WFS2、OL-WFS3、CL-WFS1、CL-WFS2 の測定値が同様に記録されている。１回あたりの測定に対してこれら 720 個のデータが記録され、各時刻に対する
これら 720 個の測定値が時系列順に並んだ形となっている。

図 2.4: RAVEN のデータセット (data size:720×42000)

本研究では、RAVEN のオンスカイ観測で得られた各データセットに対してまず Rotation Matrix
をかけている。今各 WFS の測定値は各サイエンスパスに対して回転してしまっているため、この

Rotation Matrix をかけることで同じ座標系へ変換している。この行列は各 WFS の座標系を CL-WFSの座標系に合わせる 720×720 の行列で CL-WFS1 の座標系に合わせる Rot1 と CL-WFS2 の座標系に合わせる Rot2 の二つがある。
また、AO では tip-tilt 成分の補正については高次の成分とは別に独立で補償することが多い。tip-tiltとは大気揺らぎによる天体の重心の移動成分である。この大気揺らぎによる重心の移動成分を補正す

17

第 2. 多天体補償光学実証装置 RAVEN

2.2. RAVEN のデータセット

ることで、最も大きい空間スケールの大気揺らぎの影響を大きく抑えることが可能である。AO では
この重心の移動成分は観測視野内の明るい星を使った測定を行い、tip-tilt 鏡による補正を行なっている [21]。そこで今回は tip-tilt 成分について、あらかじめデータから取り除いた上で計算を行った。

18

第 2. 多天体補償光学実証装置 RAVEN

2.2.3

2.2. RAVEN のデータセット

観測コンディション

今回解析に用いたデータセットについて、観測時のコンディションを以下にまとめた。
大気揺らぎ高度および大気揺らぎ強度分布については SLOpe Detection and Ranging(SLODAR)と呼ばれる技術を用いて推定が行われている。
データ名20140514043709

観測日時2014 年 5 月 14 日

観測開始時刻4 時 37 分 9 秒

観測時間172.8 秒

フレームレート

250/秒

データサイズ (CL-WFS,OL-WFS の測定値)

720×43199

大気揺らぎ高度[0 4.5 6 7.5 9 12]km

大気揺らぎ強度分布[0.507 0.134 0.063 0.040 0.103 0.150]

アウタースケール

40m

フリード長

0.224m

表 2.2: データ 20140514043709 の観測コンディション。

図 2.5: データ 20140514043709 の観測における各ガイド星とサイエンスターゲットの配置

19

第 2. 多天体補償光学実証装置 RAVEN

2.2. RAVEN のデータセット

データ名20150626063402

観測日時2015 年 6 月 26 日

観測開始時刻6 時 34 分 2 秒

観測時間259.8 秒

フレームレート

150/秒

データサイズ (CL-WFS,OL-WFS の測定値)

720×36565

大気揺らぎ高度[0 9 10.5 12]km

大気揺らぎ強度分布[0.725 0.014 0.039 0.022]

アウタースケール

40m

フリード長

0.185m

表 2.3: データ 20150626063402 の観測コンディション。

図 2.6: データ 20150626063402 の観測における各ガイド星とサイエンスターゲットの配置

20

第 2. 多天体補償光学実証装置 RAVEN

2.2. RAVEN のデータセット

データ名20150627105819

観測日時2015 年 6 月 27 日

観測開始時刻10 時 58 分 19 秒

観測時間

494.7 秒

フレームレート

150/秒

データサイズ (CL-WFS,OL-WFS の測定値)

720×69042

大気揺らぎ高度[0 2.5 7.5 10 12.5]km

大気揺らぎ強度分布[0.570 0.208 0.065 0.018 0.134]

アウタースケール

40m

フリード長

0.175m

表 2.4: データ 20150627105819 の観測コンディション。

図 2.7: データ 20150627105819 の観測における各ガイド星とサイエンスターゲットの配置

21

第 3 章 データに基づく波面推定法
以下、時刻 t における GS 方向の WFS1〜3 の測定値ベクトルを at 、サイエンスターゲット方向の

WFS の測定値ベクトルを bt と表す。経験的トモグラフィック波面推定では GS 方向の測定値 at からサイエンスターゲット方向の測定値 bt を推定する行列を R を実際のオンスカイの観測データを用いて求める。

3.1

推定法の概要トモグラフィック波面推定の流れを以下の図 3.1 に簡単にまとめた。

図 3.1: 解析の流れ

トモグラフィック推定行列 R の決定
図 3.1 のように GS 方向の WFS の測定値およびサイエンスターゲット方向の測定値に対して、時刻 t1 から時刻 t2 の期間 (duration) のデータセットを取り出し、それぞれ行列 A、B とする。この22

第 3. データに基づく波面推定法

3.1. 推定法の概要

期間のデータを R を訓練したという意味で、訓練データと呼ぶことにする。

sOL−W F S1 (t = t1 ) · · · sOL−W F S1 (t = t2 )
A = [at1 · · · at2 ] =  sOL−W F S2 (t = t1 ) · · · sOL−W F S2 (t = t2 ) 

sOL−W F S3 (t = t1 ) · · · sOL−W F S3 (t = t2 )B

= [bt1 · · · bt2 ] = [sCL−W F S (t = t1 ) · · · sCL−W F S (t = t2 )]データに基づいた波面推定では、このある期間におけるデータセット A、B を推定行列 R を用いて以下の式でモデル化する。

B = RA
この関係式からトモグラフィック推定行列 R の推定解 R̂ を求める。本研究では、このトモグラ
フィック推定行列 R の決定方法の提案を行い、推定精度の向上を図る。

トモグラフィック推定行列 R による推定
次に上で求めたトモグラフィック推定行列 R̂ を訓練データ以降の時刻 t における GS 方向の測定
値 at に適用することで、時刻 t におけるターゲット方向の測定値 bt を推定する。この bt の推定値を b̂t と表す。

b̂t = R̂at
訓練データ以外の時刻でターゲット方向の測定値を予想するとういう意味で、これらの時刻を予測データと呼ぶことにする。

推定精度の評価
最後に推定誤差を評価するために slope 測定値を望遠鏡開口上での波面位相に復元する波面再構成を行う。(1.3.2 を参照)

φ = (GT G)−1 GT s

(3.1)

この波面再構成のプロセスによってオリジナルの測定値 bt および推定誤差 b̂t − bt から波面位相を回復する。
波面の平均値からのずれを評価するために本論文では二乗平均平方根 (Root Mean Square: RMS)を用いる。位相 φi の平均値を φ とすると RMS は√∑N
2
i=1 (φi − φ)
RM S =
N

[m]

(3.2)

と書ける。今回のデータでは波面の位相点は 96 個の点で表しているので N=96 である。
また、元の波面と GS 方向の測定値から推定した波面の誤差を推定誤差 (Wave Front Error: WFE)と表して以下のように定義することにする。
√∑N2i=1 (φi − φ̂i )
WFE =
N

[m]

本論文では以上の二つの指標を用いて、推定精度を評価していく。

23

(3.3)

第 3. データに基づく波面推定法

3.2. 計算手法の検討

3.2

計算手法の検討

3.2.1

最小二乗法を用いた推定

まず初めに、最小二乗法によって推定行列 R の推定解 R̂ を求める。最小二乗法による推定ではコ
ストファンクション F を以下のように定義し、これを最小化するような R の解を求める。

F = ∥B − RA∥2F(3.4)

この推定解は数学的に求めることができ、解は以下のようになる。（行列の最小二乗解の導出については Appendix を参照。）

R̂ = BAT (AAT )−1
また、本研究の中でこの最小二乗解を計算する際には以下のような行列 A の特異値分解を用いた。


A




= [w1 w2 · · · wm ] 





0···0
v1t

..
..   t 
v 
.
. 
 .2 ..
..   . .
.  . 

γ1
γ2
..

.

vnt

γm 0 · · · 0= W ΓV T(3.5)

ここで w および v はそれぞれ m 次元空間、n 次元空間を張る正規直交基底であるのでこれを用ると

R̂ は
R̂

= BAT (AAT )−1= BV ΓT W T (W ΓV T V ΓT W T )−1= BV ΓT W T (W ΓΓT W T )−1= BV ΓT W T W (ΓΓT )−1 W T= BV ΓT (ΓΓT )−1 W T


γ1


..


.






γm 

= BV 

 0
···
0 

 .
.. 

 .
. 
 .
0
···
0

 1
γ

 1




= BV 
 0

 .
 .
 .
0

..

.

···
···



1
γ12

..

 T
W


.
1
2
γm





1 
γm 
WT
0 

.. 

. 
0

によって求めることができる。特異値分解は matlab の関数 svd を用いて行なった。

24

(3.6)

第 3. データに基づく波面推定法

3.2.2

3.2. 計算手法の検討

低特異値の除去

最小二乗法による推定では、非正方行列 A を特異値分解して計算を行った。しかし、ここで注意
しなくてはならないこととして、行列 A によっては特異値が 0 になる、あるいは他の特異値と比較
して絶対値が極端に小さくなるものが存在するということである。

(3.6) 式で表される R の最小二乗解は、行列 A の特異値分解で以下の形で表すことができる。
∑

min(m,n)R̂ = Bk=1

1
vk wkTγk(3.7)

ここで、m は GS 方向の測定点に対応するので m=432、n は訓練データの長さに対応する。今、訓練データとしては測定点の値を決めるために十分長くとる場合を考えて、m < n とする。また、Bを l × n 行列とおいた時 l=144 である。この時 R̂ は l × m 行列となり、R̂ の i 行目 (1 ≤ i ≤ 144) の行ベクトル r̂i は B の i 行目の行ベクトル bi を用いて以下のように表せる。
m∑1r̂i = biγk

k=1

vk wkT(3.8)

ここでベクトル bi を n 次元空間をはる正規直交基底 v を用いて
n
∑

b̂i =

αk vkT(3.9)

k=1と表すと、式 (3.8) は

r̂i =

m
∑
αk
k=1

γk

wkT

(3.10)と表せる。ここで、B に計測誤差の影響が含まれている場合、特異値 γk の中に 0 に近いものがある場合、 αγkk によって、小さい特異値に対応する成分に含まれる誤差が増幅されてしまうことがわかる。今、特異値は絶対値から大きいものから順番に γ1 ≤ γ2 ≤ · · · ≤ γm のように並んでいる。そこ
で p 番目 (1 ≤ p ≤ m) 以降の項を無視して、

r̂ =

p
∑
αk
k=1

γk

wkT

(3.11)とすることで低特異値に対応するモードを切り捨て、計測誤差に対する解の不安定性を取り除く。

3.2.3

L2 正則化を用いた推定

最小二乗法は多くの機械学習の基礎となるアルゴリズムである一方で、ノイズに対して過適合してしまうというデメリットがある。これはデータ数に対してモデルが複雑すぎることが原因である。
つまり、計算に用いるデータ数が十分ではない時にトモグラフィック推定行列 R の推定解が訓練デー
タに合わせすぎてしまい、予測データに対する波面推定の際に推定精度が劣化してしまう。このモデルの複雑さを制御するのが正則化を用いた手法である。そこで、正則化の項を計算に取り入れた波面推定を行い、従来の最小二乗法を用いた時と比較して推定精度を向上させられるか検証を行った。

25

第 3. データに基づく波面推定法

3.3. タイムステップデータを用いた推定法の検討

最小二乗法による推定ではF = ∥B − RA∥2F(3.12)

を最小とするような R を求めていた。
これに対し、今回用いた L2 正則化では上記の式 3.12 に R のノルムの二乗の項を加えることで Rの解に制限をつける。データ数が少ない時に第一項を小さくするために R のパラメータを合わせす
ぎてしまう問題を新しく加えた第二項で R のパラメータの値にペナルティをかけることで防ぐ。

F = ∥B − RA∥2F + λ∥R∥2F(3.13)

この解は解析的に求まり次式で与えられる。(Appendix 参照)

R̂ = BAT (AAT + λI)−1

3.3(3.14)

タイムステップデータを用いた推定法の検討

ここまでは計算手法に関する検討を行なってきたが、最後に推定に用いるデータに関する検討を行なう。トモグラフィーの計算では複数のガイド星の測定値から、大気の高さ方向の構造に分解して３次元構造を推定するので、推定に寄与する特徴量を増やすことによって推定精度を向上させることができると考えられる。そこで短い時間内であれば測定値には高い相関があると考え、タイムステッ
プ ∆t だけずらしたデータセット A−∆t をA−∆t = [at1 −∆t · · · at2 −∆t ](3.15)と定義し、以下のように A に測定値 A−∆t を追加した A′ を用いて推定を行うことで推定精度を上げられないかを検討を行なう。

[
′

A =

A
A−∆t

]
(3.16)

つまり、R を上記の A′ を用いて次のように推定する。

R̂ = BA T (A A T )−126(3.17)

第 4 章 結果

4.1

最小二乗法を用いた波面推定

4.1.1

最小二乗法を用いた波面推定の結果

まず、最小二乗法を用いた波面推定を行なった。
データ 20140514043709 に対して ch1 方向と ch2 方向に対して時刻 90 秒から 130 秒を取り出し訓練データとして波面推定を行った時の波面推定誤差及びオリジナルの波面の波面誤差を以下の図 4.1に示した。ここからオリジナルの波面に対して、推定波面誤差が小さくなっていることが確認できる。
次に訓練データを 70 秒、100 秒、130 秒と増やして波面推定を行った際に、波面推定誤差の変化を比較した結果を図 4.2 に示す。図 4.2 から訓練データにおける推定波面誤差は予測データにおける推定波面誤差よりも小さいことがわかる。また、訓練データを長く取った時に訓練データと予測データの波面推定誤差の差分が小さくなっていく傾向も確認できる。これはトモグラフィック推定行列 R が訓練データを用いて決定されているため、訓練データ内の特徴を反映しているためだと考えられる。
次に、訓練データの長さと波面推定誤差の関係を調べるために 5 秒間隔で訓練データを増やしていったときの予測データに対する波面推定誤差の平均および分散の結果を図 4.3 に示した。本論文における波面推定の目的としては訓練データを用いてトレーニングした推定行列 R から、その後の GS
方向の測定に対して適用することでサイエンスターゲット方向の波面を推定することであった。したがって、今回関心があるのは予測データに対する波面推定精度である。本論文では Yamazaki, 2015との比較を行うために訓練データの終点を 130 秒に固定して、135 秒から 140 秒における波面推定誤差の平均値を評価している。

図 4.1: データ 20140514043709 に対して波面推定を行った時の波面推定誤差のプロット（左図：サイエンスパス１。右図：サイエンスパス２）緑色プロットはオリジナルの波面の RMS、赤色プロットは同時刻のデータを用いて推定を行った際の波面推定誤差を表している。

27

第 4. 結果

4.1. 最小二乗法を用いた波面推定

図 4.2: データ 20140514043709 に対して波面推定を行った時の波面推定誤差のプロット（左：サイエンスパス１、右：サイエンスパス２）上から順に訓練データを 70 秒、100 秒、130 秒と増加させていった時の波面推定誤差を表している。

28

第 4. 結果

4.1. 最小二乗法を用いた波面推定

図 4.3: データ 20140514043709 に対して訓練データを増加させていったときの波面推定誤差の推移の結果を表している。（左：サイエンスパパス、右：サイエンスパス２）29

第 4. 結果

4.1.2

4.1. 最小二乗法を用いた波面推定

波面センサの測定値について

WFS ごとの測定値の性質を調べるために、各時刻の WFS の測定値の平均の分布をプロットしたものを図 4.4 に示した。次に、各 WFS の各測定点における測定値の時間平均をとった時の分布を図4.5 に示した。

図 4.4: データ 20140514043709 の各 WFS の各測定点における測定平均の分布。単位は pix。

図 4.5: データ 20140514043709 の各時刻の WFS の測定値の平均の分布。単位は pix

30

第 4. 結果

4.1. 最小二乗法を用いた波面推定

各時刻における測定値の平均は正規分布に近い分布になっており、これを RMS に直して考えると標準偏差 177nm のばらつきを持っている。各時刻における測定値の平均は重心の移動量に対応している。つまりこれは全体的な波面の持っている傾き、すなわち tip-tilt を表していると考えられる。
さらに、図 4.5 から各時刻の WFS の測定値の平均が 0 ではないことは各 WFS の測定値に固定的な収差が乗っていることを示唆している。この収差は装置由来の収差であると考えられる。
また図 4.4 の中に小さい飛び値を持っていることが確認できるが、これは何かしらのノイズによる影響であると考えられる。図 4.2 を見ても、ch1 方向と ch2 方向共に 60 秒あたりで、オリジナルの波
面の RMS の値にこのノイズに起因すると見られる大きな飛びが見られる箇所がある。これらがチャンネル方向によらずに同時刻に見えたことから、これは WFS の測定に起因するノイズではなく外部の要因によるものだったと推測できる。さらにはこの影響が波面推定にも残っていることから飛び値
はあらかじめ取り除いて推定する必要があると言える。

4.1.3

静的成分を除いた波面推定の結果

前節で、各 WFS の測定値の平均をとることによって tip-tilt 成分と測定値に固定的にのっている
収差の存在が明らかになった。tip-tilt 成分についてはあらかじめ取り除いて計算を行っていたが、固定的にのっている収差に関してもあらかじめデータから取り除くことで大気揺らぎに由来する成分を取り出して推定が行えると考え、オリジナルのデータを用いた場合と各 WFS の測定値の時間平均を引いたデータを用いた比較を行った。
その結果を図 4.6 に示した。青いプロットがオリジナルのデータに対して波面推定を行ったときの、赤のプロットがオリジナルのデータから各 WFS ごとの時間平均値を引いたデータを用いて波面推定を行ったときの結果を表している。この結果から時間平均を引いたデータを用いた時が推定誤差
としては 10nm ほど小さくなる結果となった。
これは推定に寄与する次元を落とすことで推定精度をあげることができたと考えられる。この結果
から装置由来の収差成分が推定精度を劣化させる原因になりうることがわかった。そのため、装置由来の収差成分に関してはあらかじめなくなるようにキャリブレーションをする必要があると言える。

図 4.6: 青 : オリジナル . 赤: オリジナルデータから全時間平均を引いたデータを用いて波面推定を行った時の波面推定誤差。（左：サイエンスチャンネル１。右：サイエンスチャンネル２）。

31

第 4. 結果

4.2. 特異値除去を行った時の結果

4.2

特異値除去を行った時の結果

4.2.1

最大特異値での基準

まず、特異値除去の基準として以下の yamazaki,2015 の中で使われていた最大特異値で基準化した値を用いた。

γp
γ1

(4.1)

この基準で 10−9 、0.01、0.05、0.1 以下で特異値を除去した時の結果を図 4.7 に示した。この結果
から訓練データが短い時には特異値除去の基準を大きくした時に推定結果が良くなり、訓練データが長くなるにつれて特異値除去の基準を小さくした方が推定結果が良くなる傾向が確認できる。
これは訓練データが短い時にはノイズの影響が大きく、これを取り除いた方が推定精度として良くなるが十分なデータ数が確保できた場合にはより多くの情報を推定に取り入れることが望ましいことを表している。

図 4.7: 特異値除去の基準を 10−9 、0.01、0.05、0.1 にしてトモグラフィー推定を行った時の波面推定
誤差とデータ記録時間の関係。(左：サイエンスチャンネル１。右：サイエンスチャンネル２)。

4.2.2

累積寄与率での基準

次に特異値除去の基準として累積寄与率を用いた場合の推定結果を示す。累積寄与率は以下の式
で与えられるように、第 p 個までの特異値がどれだけデータの情報を説明しているかを示す尺度で
ある。
p
∑
k=1
m
∑

γi
(4.2)
γi

k=1

以下の図 4.8 にデータ 20140514043709 に対して累積寄与率基準で 40%、80%、100% を基準に、それ以降の特異値を除去した場合の波面推定誤差の結果を示した。この結果から累積寄与率基準を大きくして特異値除去をすると、訓練データが長い場合に対しては推定結果はよくなる一方で、訓練データが短い時に対しては悪くなるという傾向が確認できる。また、累積寄与率 100% で基準をと32

第 4. 結果

4.2. 特異値除去を行った時の結果

る、つまり全ての特異値を計算に含めてしまうと訓練データが短い時に推定精度が非常に悪くなる
ことがわかる。
また、累積寄与率 80% および 100% の基準で推定したトモグラフィック推定行列 R の構造の比較を図 4.9 に示した。累積寄与率 80% で特異値を切った場合には R の構造が見えてくるが、全ての情報を使って推定した R（右図）を見るとノイズの影響で構造が見えなくなっていることがわかる。

図 4.8: データ 20140514043709 に対して累積寄与率 40%、80%、100% で特異値を切った時の波面
推定誤差の比較。（左：サイエンスチャンネル１。右：サイエンスチャンネル２。）

図 4.9: 推定期間 80 秒〜130 秒から推定した R の構造。左は累積寄与率 80%、右は 100%で低特異値除去を行って計算した結果を表している。

33

第 4. 結果

4.3

L2 正則化による波面推定

4.3.1

正則化の項 λ4.3. L2 正則化による波面推定

正則化の項をつけて波面推定を行った時の正則パラメータ λ の大きさと波面推定誤差の関係を図4.10 に示した。この図 4.10 の上段は訓練データ 100 秒、下段は訓練データ 10 秒に対して推定を行った時の結果を示しており、ここから訓練データによって λ の値とその際の波面推定誤差の振る舞いも変化することがわかる。そこで、次に訓練データの長さと予測データに対する波面推定誤差を最小にする λ の関係を図 4.11 に示した。

図 4.10: 正則化の項をつけてトモグラフィーを行った時の正則パラメータ λ と波面推定誤差の関係。
上段は訓練データ 10 秒、下段は訓練データ 100 秒に対して推定を行った時の結果（左：サイエンスチャンネル１、右：サイエンスチャンネル２）

34

第 4. 結果

4.3. L2 正則化による波面推定

図 4.11: 訓練データの長さと予測データに対する波面推定誤差を最小にする λ の関係。（左：サイエ
ンスチャンネル１、右：サイエンスチャンネル２）

4.3.2

L2 正則化による結果

上で求めた各訓練データで波面推定誤差を最小にする正則化パラメータ λ を用いて波面推定を行い、この結果を従来の最小二乗法を用いた推定と比較したものを以下の図 4.12 に示した。ここからλ の値が 0.1 よりも小さいデータ記録時間 25 秒以上に対しては正則化の項を入れて推定しても、従来の最小二乗法を用いた推定と同程度の結果しか得られないことがわかる。一方で、データ記録時間が短い場合には、正則化の項を入れた推定が上回る結果となった。これは訓練データが短い時には、訓練データに対する過学習の問題が生じてくるためだと考えられる。

図 4.12: 正則化の手法を用いて波面推定を行った場合と従来の最小二乗法を用いて波面推定を行った
場合の波面推定誤差の比較（左：サイエンスチャンネル１、右：サイエンスチャンネル２）

35

第 4. 結果

4.4

4.4. タイムステップデータを用いた波面推定の結果

タイムステップデータを用いた波面推定の結果

タイムステップデータを追加して波面推定を行った時の結果を図 4.13 に示した。横軸は何 frame 前のデータを追加したかを表しており、30〜130 の 100 秒間を訓練時刻として 130〜135 秒における波面推定誤差の平均値を示している。この結果から、タイムステップデータを追加すると、より過去のデータを追加していくごとに波面推定誤差が小さくなっていく傾向が確認できる。そして、15frame（60 ミリ秒前）のデータを追加して推定した時に波面推定誤差が最小になる。この時の波面推定誤差
としてはタイムステップデータを用いない時と比較して 15nm 小さくなっている。そして、これよりもさらに過去の測定データを追加した場合には波面推定誤差は徐々に大きくなっていく。

図 4.13: タイムステップデータを追加して推定した波面推定誤差の結果。横軸は何 frame 前のデー
タを追加したかを表している。30〜130 秒の時刻を訓練時刻として用いた (左：サイエンスチャンネ
ル１、右：サイエンスチャンネル２)

36

第 5 章 議論

5.1

平均値に現れるパターンについて

4.1.2 で見つかった固定的なパターンについて、各データセットに対してその形を開口上での波面の形に直したものと各開口上で Zernike 多項式展開した時の各モードの強度を図 5.1 に示した。座標
系はサイエンスチャンネル１に合わせた。Zernike 多項式については Appendix を参照。
図 5.1 から各データセットの各 WFS の測定値の平均値には同様のパターンが見られることがわかった。この各 WFS の測定値の平均値に現れる固定的なパターンは各サイエンスパス内の光学系で生じるシステマティックな収差に由来すると考えられる。各 WFS においてこのパターンに類似性が見られるのは各サイエンスパスにおける共通部分で生じた収差の影響を反映しているからだと考えられる。CL-WFS ごとの差異は向いている方向の違いによるものだと考えられる。固定的なパター
ンの RMS を以下の表 5.1 に示した。データ 20140514043709 のサイエンスチャンネル２における大きな収差に関しては、サイエンスターゲットをピックアップする光学系で大きな収差が生じてしまっている可能性がある。
データセットの違いによる固定パターンの比較をすると、いずれのデータセットに対しても 6 番目のモードに対応するゼルニケ係数で大きな値を持っていることがわかった。これは非点収差に対応するモードを表している (Appendix 参照)。今すばる望遠鏡の主鏡は 261 個のアクチュエータで支持され、主鏡セルに３点の固定機構によってつながれている。主鏡の鏡面は鏡面測定装置によって計測され、鏡面を補正するための支持力がアクチュエータにフィードバックされる [18]。今回、全てのデー
タセットに現れた非点収差は 3 点の固定機構につながれたことによる主鏡の変形に起因している可能性があるが、さらなる検証が必要である。
この固定パターンを引くことで波面推定精度が上がったのは各サイエンスパスで異なる望遠鏡由来の収差の成分 σsystematic を推定から取り除くことができたからであると考えられる。図 4.3 の緑のプロットで示したようにオリジナルデータの RMS はサイエンスチャンネル１方向では 554nm、
サイエンスチャンネル２方向では 831nm であった。方向の違いはあるものの、ほとんど同じ波面を測定しているにも関わらず RMS に大きな違いが見られたのも望遠鏡由来の収差が波面測定にのっていたからだと考えられる。したがって、オリジナルデータから表 5.1 の固定パターンを引いたものが大気揺らぎ由来の収差 σatmosphere であると考えられる。これらの収差は独立であるので
222から計算すると、データセット 20140514043709 に対して大気揺
σoriginal= σatmosphere+ σsystematic

らぎ由来の収差はサイエンスチャンネル１方向では 399nm、サイエンスチャンネル２方向では 435nmと求まる。サイエンスチャンネルごとの差異は WFS の測定ノイズや方向の違いによるものであると考えられる。

OL-WFS1

OL-WFS2

OL-WFS3

CL-WFS1

CL-WFS2

データ 20140514043709

419nm

440nm

425nm

384nm

708nm

データ 20150626063402

375nm

419nm

337nm

370nm

410nm

データ 20150627105819

437nm

482nm

381nm

428nm

451nm

表 5.1: 各データセットの各 WFS の測定値の固定パターンの収差37

第 5. 議論

5.1. 平均値に現れるパターンについて

図 5.1: データ 20140514043709 における CL-WFS1〜3 および OL-WFS1,2 の各 WFS の測定値平均から計算した開口上での固定パターンおよび開口上でのパターンを Zernike 多項式展開した時の各モードの振幅38

第 5. 議論

5.1. 平均値に現れるパターンについて

図 5.2: データ 20150626063402 における CL-WFS1〜3 および OL-WFS1,2 の各 WFS の測定値平均から計算した開口上での固定パターンおよび開口上でのパターンを Zernike 多項式展開した時の各モードの振幅39

第 5. 議論

5.1. 平均値に現れるパターンについて

図 5.3: データ 20150627105819 における CL-WFS1〜3 および OL-WFS1,2 の各 WFS の測定値平均から計算した開口上での固定パターンおよび開口上でのパターンを Zernike 多項式展開した時の各モードの振幅40

第 5. 議論

5.2

5.2. 特異値除去の基準

特異値除去の基準

本節では特異値除去の基準の妥当性について考えるために、最大特異値で基準化した場合と寄与率分布で考えた場合を比較しながら議論を行っていく。
この基準を考えるにあたってそれぞれオリジナルのデータを用いた時と、時間平均を引いたデータを用いた時の推定結果を比較した。

5.2.1

最大特異値の基準

Yamazaki, 2015[6] の中で、『特異値を小さい側から除去した時の経験的トモグラフィック再構成行
列による波面誤差が増加に転じる限界を特異値除去の基準とする』と述べられている。そこで、まず特異値除去の基準として最大特異値で基準化した時にこの中で使われた 0.05 という値を用いた推定を行った。その結果を以下に示す。

図 5.4: 最大特異値で基準化した時に 0.05 を基準として特異値を除去した時の推定結果。青：オリジ
ナルのデータを用いた推定結果。赤：時間平均を引いたデータを用いた推定結果。(左：サイエンスチャンネル１。右：サイエンスチャンネル２）

ここで平均値を引いて推定を行った時に、結果の改善が見られたがこれを特異値分布の観点から議論していく。図 5.5 に訓練データ 100 秒に対するオリジナルのデータおよび時間平均を引いたデー
タの A の特異値分解を行なった際の特異値を大きさ順に並べた。他と比べて明らかに小さい特異値
が 6 個存在することが確認できる。
また、表 5.2 に、オリジナルのデータおよび時間平均を引いたデータの A に対して特異値分解を行なった際の特異値分布を示した。左はオリジナルのデータに対する特異値の分布、右は時間平均を引いたデータに対する特異値の分布を示している。なお各特異値は最大の特異値で規格化している。この表から平均値を引いた場合には最大の特異値が小さくなって相対的に特異値分布が前半によっていることがわかる。

41

第 5. 議論

5.2. 特異値除去の基準

図 5.5: 特異値を大きいものから順番に並べたプロット。(左：オリジナルのデータ。右：時間平均を引いたデータに対してそれぞれ特異値を求めた)。

1∼0.5

5

10

30

70

100

130

5

10

30

70

100

130

1

1

1

1

1

1

1∼0.5

2

2

4

4

5

5

−1

45

78

140

169

191

202

−1

12

15

20

21

22

23

∼ 10

∼ 0.05

42

53

67

75

81

84

∼ 0.05

78

102

123

117

106

101

241

∼ 10

−2

215

187

123

104

95

90

∼ 10

−3

74

50

32

30

27

26

−4

∼ 10
∼ 10

−2

∼ 10

−3

∼ 10

−4

19

13

11

10

9

9

∼ 10

11

7

4

2

2

2

∼ 10−5

1

1

1

1

1

1

∼ 10−5

1

0

0

0

0

0

0

∼ 10

−6

0

0

0

0

0

0

∼ 10

−7

0

0

0

0

0

0

∼ 10

−8

0

0

0

0

0

0

∼ 10

−9

0

0

0

0

0

0

−10

∼ 10

−6

∼ 10

−7

∼ 10

−8

∼ 10

−9

215
136

0
0
0
0

239
104

0
0
0
0

247
79

0
0
0
0

246
72

0
0
0
0

243
69

0
0
0
0

67

0
0
0

−10

0

0

0

0

0

0

∼ 10

0

0

0

0

0

0

≤ 10−10

6

6

6

6

6

6

≤ 10−10

6

6

6

6

6

6

∼ 10

表 5.2: データ記録時間、10、30、70、100、130 秒に対して行列 A を特異値分解した時の特異値の分布。(左：オリジナルのデータ。右：時間平均を引いたデータに対してそれぞれ特異値を求めた)。
また特異値は最も大きいもので割ることで基準化している。

42

第 5. 議論

5.2.2

5.2. 特異値除去の基準

累積寄与率の基準

主成分分析などでは全体の情報の 8 割程度が含まれていれば良いという考えのもと、8 割を超えた
ところで切ってしまうのが一般的である。本論文においてもこの累積寄与率 80%を用いた推定を行った結果を以下の図 5.6 に示した。この結果を見ると、時間平均を引いたデータを用いた推定誤差が小さくなっているものの最大特異値の基準を用いた時ほど大きな差異は見られない。また、各訓練デー
タの長さに対応する行列 A の特異値の累積寄与率分布を表 5.3 にまとめた。

図 5.6: 累積寄与率基準 80%で特異値除去をした時の推定結果。青：オリジナルのデータを用いた推定結果。赤：時間平均を引いたデータを用いた推定結果。(左：サイエンスチャンネル１。右：サイ
エンスチャンネル２）

5

10

30

70

100

130

5

10

30

70

100

130

∼ 10%

2

3

3

4

4

4

∼ 10%

4

5

7

8

8

9

∼ 20%

9

11

14

15

16

17

∼ 20%

12

15

19

21

22

23

∼ 30%

20

25

29

31

33

34

∼ 30%

24

30

36

39

40

41

∼ 40%

36

43

50

53

55

56

∼ 40%

41

48

57

60

62

63

∼ 50%

57

66

75

79

81

83

∼ 50%

62

72

82

87

89

90

∼ 60%

83

95

107

111

113

115

∼ 60%

88

100

114

119

121

122

∼ 70%

118

131

145

150

153

154

∼ 70%

123

137

152

157

159

161

∼ 80%

165

179

194

199

201

202

∼ 80%

170

184

199

204

206

207

∼ 90%

236

249

259

263

264

265

∼ 90%

240

252

263

267

268

269

∼ 100%

432

432

432

432

432

432

∼ 100%

432

432

432

432

432

432

表 5.3: データ記録時間、10、30、70、100、130 秒に対して行列 A を特異値分解した時の特異値の累積寄与率分布。(左：オリジナルのデータ。右：時間平均を引いたデータに対してそれぞれ特異値を求めた)。

43

第 5. 議論

5.2.3

5.3. 正則化と過学習基準の比較

表 5.3 をみると同数の特異値に対してはオリジナルのデータよりも平均値を引いたデータの累積寄与率が小さくなっていることがわかる。しかし、これらの間に多少の違いはあるものの、累積寄与率が 80%を超えてくる特異値の個数についてはほとんど相違がない。つまり、平均値を引くことによって累積寄与率分布が大きく変動するという傾向は見られない。一方で、表 5.2 からわかるように平均値を引くことによって、最大特異値に対する各特異値の相対的な大きさの分布は前半にずれてしまうことがわかる。これは最大特異値の寄与率が平均値という特徴を引くことによって小さくなってしまうことに起因すると考えられる。今、Yamazaki.2015[6] で用いられていた特異値除去の基準を見てみるとオリジナルのデータでは 0.05 以上の特異値は、データ記録時間 100 秒については 104 個であるのに対して、平均値を引いたデータに対しては 302 個になっている。この値と表 5.3 における値を比較してみると、特異値除去の基準 0.05 という値はオリジナルのデータでは累積寄与率 40%程
度でそれ以外の 60%の情報を切ってしまっていることになる。一方で、平均値を引いたデータに対してはほぼ累積寄与率 80%を達成しており、十分な水準に達していると言える。まとめると、特異値除去の基準として最大特異値で基準化した値を用いることは推定に寄与している情報まで落としてしまっている可能性があるため、特異値除去の基準としては累積寄与率で判断するのが妥当である
ことがわかった。

5.3

正則化と過学習

以下の図 5.7 は訓練データの長さと波面推定誤差の関係を表している。青いバーが 130〜135 秒の予測データにおける推定誤差の結果を、赤いバーが 125〜130 秒の訓練データに対する波面推定誤差を、また緑の線が元の波面の平均の RMS を表している。ここから訓練データが短いときに、予測
データに対する推定誤差が非常に大きいことがわかる。これは R のパラメータに対して訓練データが少なすぎるために R の自由度が大きく、訓練データに合わせすぎてしまうために起こる過学習の問題であると考えられる。つまり予測データに対して正しく推定できないということになる。これを解決するためには一つとしてデータ数を増やすという方法が考えられる。これは手元にある訓練
データの量が少ないことが理由でデータの傾向を十分に表現しきれていない場合があるため、学習データの数を増やすことが過学習を避け、汎化性能を高めることへ効果的である。実際この結果からもわかるように、訓練データを増やすと予測データに対する推定誤差が小さくなっていく。
次に各訓練データで波面推定誤差を最小にする正則化パラメータ λ を用いてトモグラフィック推定
行列 R を推定し、訓練データに対する推定結果を最小二乗法を用いた時の結果と比較したものを以下の図 5.8 に示した。訓練データが 20 秒以下の短い時には正則化の手法を用いることで、図 4.12 に示したように予測データに対する波面推定誤差は改善していた。一方で、図 5.8 から訓練データに対する波面推定誤差は悪くなっていることがわかる。この結果からも正則化の項を計算に含めたことに
よって訓練データに対する過学習が改善されたことがわかる。
これをトモグラフィック推定行列 R の形から見ていく。図 5.9 はそれぞれトモグラフィック推定行
列 R における 49 番目の CL-WFS の x 方向の測定値に寄与する OL-WFS1〜OL-WFS3 の x 方向の各測定値の重みを表している。上段は最小二乗法を用いて訓練時間５秒で推定を行った結果、中段は正則化の項を加えて訓練時間５秒で推定を行った結果、下段は訓練時間 100 秒で累積寄与率 80%で低特異値除去を行った時の結果を示しているの結果を示している。ここで、トモグラフィック推定行
列 R の構造において、地上層の影響は CL-WFS と同じ副開口に対応する位置に現れる。また、高層
大気の影響はサイエンスターゲット方向と各 GS の方向の違いだけずれて現れる。Yamazaki, 2015[6]では RAVEN で用いた大気のモデルとガイド星の配置に基づいて、トモグラフィック推定行列 R に現れる高層大気の影響を推定した。これを図 5.10 に示した。ここでは 49 番目の副開口 (図の中で赤44

第 5. 議論

5.3. 正則化と過学習でプロット) に着目し、地上層からはじまり各高層大気層の影響がトモグラフィック推定行列 R の中で、３つの OL-WFS の測定値のどの箇所の寄与として現れるかをマップしている。したがって、ト
モグラフィック推定行列 R の推定が正しく行えているかの指標として R の構造に図 5.10 で示したよ
うなパターンが現れるかを確認すれば良い。
図 5.9 上段の結果から訓練時間が短い時には訓練時間の特徴を過学習してしまい、R 本来の構造が見えなくなってしまっていることがわかる。一方、図 5.9 中段の正則化の項を含んだ推定の結果では、上の結果で見られた R 全体にわたる大きな寄与がなくなっている。特に、CL-WFS1 に対する

OL-WFS1 の寄与と CL-WFS2 に対する OL-WFS2 の寄与には大気のモデルから推定されるパターンが確認できる。このことからも正則化の手法によって R が汎化ができていることがわかる。しかし、正則化の項を入れた推定を行った場合であっても、サイエンスターゲット方向の推定に無関係と考えられる場所の測定値に強い寄与が見られるなど十分に R の構造を推定しきれているとは言えない。
以上をまとめると、トモグラフィック推定行列 R の構造を推定する際に訓練時間が短い場合には正則化の項を含めることにより、訓練時間に対する過学習の影響を抑えることができた。一方で、図5.9 下段の訓練時間 100 秒で推定を行った時の結果で、大気のモデルから推定されるパターンをより鮮明に確認することができるように R の推定の際に十分な訓練時間が必要となることも改めて確認
された。

45

第 5. 議論

5.3. 正則化と過学習

図 5.7: 訓練データの長さと波面推定誤差の関係。青いバーは 130〜135 秒の予測データにおける波面推定誤差の結果、赤いバーが 125〜130 秒の訓練データに対する波面推定誤差の結果を示している（左：サイエンスチャンネル１、右：サイエンスチャンネル２）

図 5.8: 正則化の手法を用いてトモグラフィーを行った場合と従来の最小二乗法を用いてトモグラ
フィーを行った場合の波面推定誤差の比較。訓練データに対する推定誤差を表している（左：サイエ
ンスチャンネル１、右：サイエンスチャンネル２）

46

第 5. 議論

5.3. 正則化と過学習

図 5.9: 49 番目の WFS の測定値に寄与する OL-WFS1〜OL-WFS3 の各測定値の重み。上段は最小二乗法を用いて訓練時間５秒で推定を行った結果、中段は正則化の項を加えて訓練時間５秒で推定を行った結果、下段は訓練時間 100 秒で累積寄与率 80%で低特異値除去を行った時の結果を示している（左：サイエンスチャンネル１、右：サイエンスチャンネル２）

47

第 5. 議論

5.3. 正則化と過学習

図 5.10: 大気モデルおよびガイド星の配置から考えられるトモグラフィック推定行列 R に現れるパターン (From Yamazaki, 2015[6])48第 5. 議論

5.4

5.4. 推定精度を決める要因

推定精度を決める要因

従来のトモグラフィーモデルでは様々なパラメータを仮定してモデル計算が行われているが、本研究で用いたデータに基づいた波面推定において各パラメータが推定精度にどれほど影響してくるのかについて考察していく。

5.4.1

光学系パラメータ

以下の図 5.11 はデータに基づいた波面推定において回転補正をかけて補正を行なったときとかけずに補正を行なった結果を示している。この図を見ると、回転補正をかけて座標系をサイエンスチャ
ンネル１に合わせたときと、回転補正をかけなかったときの結果に差異は見られない。つまりデータに基づいた波面推定では系の回転に対するキャリブレーションが自動的に実行されていることがわかる。また、この他にも 5.1 で見られた装置内部の光学系に由来する収差についても、時間平均を引くことによってあらかじめ取り除いて推定を行っている。この光学系由来の収差を時間平均を引くことでキャリブレーションした時としなかった時の波面推定誤差を比較したものが図 5.11 であった。その結果、推定精度としては 15nm の違いが見られた。以上のような装置由来のパラメータについてもデータに基づいた波面推定ではモデルを立てることなく自動的に計算の中で処理を行うことが可能である。一方で、トモグラフィーモデルにおいては DM や WFS の光学的な関係や装置由来の収差のモデルを立てる必要がある。モデルと実際のシステムにズレが生じている場合、この部分でデータに基づく波面推定の結果がモデルに基づいた推定結果よりも良くなってくる一つの要因になっていると考えられる。

図 5.11: 回転補正の有無。データ 20140514043709 に対して紫：回転補正をかけずに波面推定を行った時の結果。緑：ch1 の座標系に合わせて波面推定を行った時の結果を示している（サイエンスチャ
ンネル１方向）49

第 5. 議論

5.4.2

5.4. 推定精度を決める要因

ガイド星の数とガイド星の配置

本節では GS の数および配置と推定精度の関係について考えていく。そこでまず、これまで CL-

WFS1〜CL-WFS3 の測定値を用いてターゲット方向の波面推定を行ってきたところに対して、１つの WFS の測定値のみを用いて波面推定を行った場合の結果を図 5.12 に示した。また、図 5.13 に２つの WFS の測定値を用いた場合の波面推定結果を示した。
図 5.12 のデータ 20140514043709 の結果では、１つの GS から推定した場合には、サイエンスター
ゲット１方向に対しては WFS １の測定値を用いた時がもっとも推定誤差としては小さくなった。これは図 2.5 の GS の配置からもわかるように、CL-WFS １に対応する GS の配置がサイエンスチャン
ネル１方向にもっとも近くなっていることに起因すると考えられる。逆にサイエンスターゲット２
方向に対しては WFS １の測定値を用いて推定した時の推定誤差が最も大きくなっている。これはサイエンスターゲット方向２と CL-WFS1 の方向が大きくずれてしまっているために起こるアンナイ
ソプラナティズムの問題を表している。次に２つの WFS の測定値を用いた場合の図 5.13 のデータ

20140514043709 の結果では、サイエンスチャンネル１方向では WFS1、WFS2 を用いた場合、サイエンスチャンネル２方向では WFS2、WFS3 を用いた場合に波面推定誤差が最も小さくなった。これは、サイエンスチャンネル１に対しては WFS3 方向の向きのずれが最も大きいため WFS3 以外の
２つの WFS の測定値を用いた場合に、サイエンスターゲット方向の大気の情報をより多く計算に取り入れることができているからだと考えられる。サイエンスチャンネル２方向に対しても同様に考えることができる。

図 5.12:波面推定を行う際に１つの WFS の測定値のみから推定した結果。上段はデータ

20140514043709 に対する結果、下段はデータ 20150626063402 を示している（左：サイエンスチャ
ンネル１、右：サイエンスチャンネル２）50

第 5. 議論

5.4. 推定精度を決める要因

図 5.13: 波面推定を行う際に２つの WFS の測定値から推定した結果。上段はデータ 20140514043709
に対する結果、下段はデータ 20150626063402 を示している (左：サイエンスチャンネル１、右：サ
イエンスチャンネル２）

51

第 5. 議論

5.4. 推定精度を決める要因

ここからは GS の数および配置と推定精度の関係をデータ 20140514043709 およびデータ 20150626063402
を用いてより詳細に調べていく。そのために、実際の想定とは違うが実験的に OL-WFS1、OL-WFS2、

OL-WFS3、CL-WFS1 の測定値からサイエンスチャンネル２方向の波面推定および波面推定誤差の評価を行った。なお、今回の測定では DM の制御は行なっていないため CL-WFS1 を擬似的に NGSとみなして推定を行っても問題ない。
データ 20140514043709 に対して訓練時刻 0〜130 秒の 130 秒間を用いた場合と、推定データ

20150626063402 に対して訓練時刻 0〜200 秒の 200 秒間を用いた場合の波面推定誤差と用いたガイド星の関係を以下の表 5.4 および表 5.5 にまとめた。また、この結果について以下で簡単に述べる。
データ 20140514043709
まず GS を１つから４つと増やしていった時の波面推定誤差について、それぞれで最も推定精度が高くなった時の波面推定誤差および用いた GS についてまとめる。一つの GS を用いた推定では CL-WFS2 を用いた場合に 334nm(結果 2)、2 つの GS 用いた推定では OL-WFS2、

OL-WFS3 を用いた場合と OL-WFS2、CL-WFS1 を用いた場合に 312nm(結果 8、結果 9)、3
つの GS 用いた推定では、OL-WFS2、OL-WFS3、CL-WFS1 を用いた場合に 296nm(結果 14)、

4 つの GS 用いた推定では 295nm(結果 15) という結果となった。
1 つの GS を用いた推定では角度距離が 1.328’ 離れた GS １つを用いた場合と 0.670’ 離れたGS1 つを用いた場合 (結果 1 と結果 2) では波面推定誤差に 51nm もの差が見られた。次に、1
つの GS を用いた場合と 2 つの GS を用いた場合の結果を比較すると、1 つの GS を用いた推定に対して、OL-WFS1 を選んで追加した時には波面推定誤差としては 2〜8nm の改善しか見られな買った（結果 2 と結果 5、結果 3 と結果 6、結果 4 と結果 7）。一方で、残りの 3 つの GSの中から 2 つの GS を選んだ場合には GS を 2 つにすることで波面推定誤差に 22〜43nm の改善が確認できた（結果 2 と結果 8、結果 2 と結果 9、結果 3 と結果 8、結果 3 と結果 10、結果 4と結果 9、結果 4 と結果 10）。また、OL-WFS1 以外の 3 つの GS について GS を２つから 3 つに増やした時には 16〜22nm の改善が見られた (結果 8 と結果 14、結果 9 と結果 14、結果 10と結果 14)。最後に OL-WFS2、OL-WFS3、CL-WFS1 に対して OL-WFS1 を追加して GS を4 つに増やして推定を行った場合には波面推定誤差として 1nm の改善しか見られなかった (結果 14 と結果 15)。
データ 20150626063402
上と同様に最も推定精度が高くなった時の波面推定誤差および用いた GS についてまとめる。
一つの GS を用いた推定では CL-WFS1 を用いた場合に 510nm(結果 4)、2 つの GS 用いた推定では OL-WFS3、CL-WFS1 を用いた場合に 454nm(結果 10)、3 つの GS 用いた推定では、

OL-WFS2、OL-WFS3、CL-WFS1 を用いた場合に 425nm(結果 14)、4 つの GS 用いた推定では 413nm(結果 15) という結果となった。

1 つの GS を用いた推定では角度距離が 1.478’ 離れた GS を用いた場合と 0.658’ 離れた GS を用いた場合 (結果 1 と結果 4) で波面推定誤差に 53nm の差が見られた。GS の数と推定精度の関係では結果 10 と結果 11 を除いて GS が多いほど波面推定誤差が小さくなっていることが確認できる。GS を 1 つから 2 つに増やした時にはおよそ 50nm の改善が見られ、さらに 1 つ増
やして GS を 3 つにするとさらに 40nm ほど改善が見られた。最後に GS を 3 つから 4 つに増やした場合でも 10〜20nm ほどの改善が確認できた。

52

第 5. 議論結果

5.4. 推定精度を決める要因

ガイド星の数

ガイド星および角度距離 (’)

OL-WFS1

OL-WFS2

1

1

2

1

3

1

4

1

5

2

1.328

6

2

1.328

7

2

1.328

8

2

0.670

9

2

0.670

10

2

11

3

1.328

0.670

12

3

1.328

0.670

13

3

1.328

14

3

15

4

OL-WFS3

波面推定誤差 (nm)CL-WFS1

1.328

385
0.670

334
0.610

343
0.888

0.670

329
0.610

335
0.888

353

0.610
0.610

1.328

355

312
0.888

312

0.888

318

0.610

309
0.888

312

0.610

0.888

315

0.670

0.610

0.888

296

0.670

0.610

0.888

295

表 5.4: ガイド星との距離と波面推定誤差の関係。データ 20140514043709 のサイエンスチャンネル
２方向の波面推定に対する結果。なお、訓練時刻は 0〜130 秒の 130 秒間を用いた。

結果

ガイド星の数

ガイド星および角度距離 (’)

OL-WFS1

OL-WFS2

1

1

2

1

3

1

4

1

5

2

0.804

6

2

0.804

7

2

0.804

8

2

1.478

9

2

1.478

10

2

11

3

0.804

1.478

12

3

0.804

1.478

13

3

0.804

14

3

15

4

OL-WFS3

波面推定誤差 (nm)CL-WFS1

0.804

539
1.478

569
0.521

520
0.658

1.478

499
0.521

488
0.658

0.521
0.521

0.804

510

465
482

0.658

463

0.658

454

0.521

456
0.658

434

0.521

0.658

430

1.478

0.521

0.658

425

1.478

0.521

0.658

413

表 5.5: ガイド星との距離と波面推定誤差の関係。データ 20150626063402 のサイエンスチャンネル
２方向の波面推定に対する結果。なお、訓練時刻は 0〜200 秒の 200 秒間を用いた。

53

第 5. 議論

5.4. 推定精度を決める要因

次に以上の結果を踏まえて GS の数および配置と波面推定精度の関係について考察していく。

1 つの GS を用いる場合
表 5.4 および表 5.5 の結果 1〜結果 4 から GS の配置をサイエンスターゲット方向からに遠くにおくほど波面推定誤差が劣化していく傾向があることがわかった。特に 1’ 以上離れた GS を用いた場
合にはそれ以外の GS を用いた場合と比べて推定精度として 30nm 近く悪くなっていた。これは GSの配置を遠くにおくことによって、ターゲット方向からの光が通るパスで GS がカバーする領域が小さくなるためである。また、GS を離しすぎることで推定に寄与しない大気の情報が測定値に多くのってきてしまう影響もあると考えられる。この推定に寄与しない情報がトモグラフィック推定行列

R の計算を行う中で推定精度を劣化させる原因になっている可能性がある。
この結果からトモグラフィック波面推定においてサイエンスターゲットと GS の距離が推定精度に大きく関わってくることがわかった。

複数の GS を用いる場合
表 5.4 および表 5.5 から GS を増やすほど波面推定精度は良くなる。このことからトモグラフィック
波面推定において GS の数が推定精度に大きく関わってくることがわかる。そこで以下では波面推定における推定精度の限界を達成するために必要な GS の数について考えていく。GS の数を 1 つから 4
つまで増やすごとに推定精度としてはデータ 20140514043709 に対しては 22nm、16nm、1nm、データ 20150626063402 に対しては 56nm、29nm、12nm 向上した。この中でデータ 20140514043709 に対して結果 14 と結果 15 を比較した時に、GS の数を３つから 4 つに増やしても推定精度が 1nm しか向上していない。この要因としては、1 つ目にはこの 295nm という値が本手法による精度限界になっていること、2 つ目にガイド星の配置によるものの２つの可能性が考えられる。
仮に 295nm が本手法による推定精度の限界を示しているとすると、このデータ 20140514043709の波面推定に対して GS ３つの情報で大気の構造の十分な推定が行えていることになる。しかし、4.4で見たようにタイムステップデータを用いた推定によって推定精度が 284nm まで向上した。次節で詳細に考察するが、この結果から大気揺らぎの推定を行うための情報として 3 つの GS の測定値では不十分であったと考えられる。また、データ 20150626063402 に対しては GS を 4 つにすることによって推定精度が 12nm 向上していることからも、GS を 3 つによる波面推定の精度が限界に達していないことが示唆される。したがって、データ 20140514043709 に対して結果 14 と結果 15 に違いが見られなかったのは、GS の数による推定精度の限界ではなく GS の配置によるものであると考えられる。CL-WFS2 に対する OL-WFS1 の配置と CL-WFS2 に対する CL-WFS1 の配置を比較するとほぼ同方向で角度距離としても 0.5 しか離れていないことがわかる。結果 4 と結果 7 を比較すると CL-WFS1 に OL-WFS1 の測定値を追加した場合でも波面推定精度としては 2nm の改善しか見られない。一方で、OL-WFS1 以外の３つの GS の結果を比較すると GS を 1 つから 2 つに増やした時に波面推定誤差が 22〜43nm 小さくなった。また結果 7 と結果 2 および結果 3 を比較しても、

OL-WFS1、CL-WFS1 の 2 つを用いた推定 (結果 7) よりも OL-WFS2、OL-WFS3 の 1 つを用いた
推定 (結果 2 および結果 3）の精度が良い結果となった。このことから OL-WFS1 はサイエンスター
ゲット方向から離れすぎていることと、CL-WFS1 と同じ方向に配置されているために追加の情報としてほとんど意味をなさないと考えられる。
この結果から複数の GS を用いる場合にはサイエンスターゲットに対する配置が推定精度に関わってくることがわかった。

54

第 5. 議論

5.4. 推定精度を決める要因

以上をまとめると GS の配置としてサイエンスターゲット方向の波面推定に寄与する情報をより多く含むように、サイエンスターゲットと近い位置に方向をずらしておくのが良いことがわかった。これによって、GS によってカバーされる推定に関係する領域が大きくなり推定精度が良くなる。
それでは、波面推定における推定精度の限界を達成するためには幾つの GS が必要になってくるのだろうか。これは、推定する大気揺らぎが高さ方向にどれだけ構造を持っているかによって決まってくると考えられる。つまり、大気揺らぎの層が多いほどそれらの層の足し合わせとして現れる測定
値を用いて縮退を解くために必要な GS の測定値も増えてくると考えられる。RAVEN の大気モデ
ルのシミュレーションでは、星像に影響を与える大気層を 7 つに分けている [14]。MOAO 実証装置

CANARY のデータを用いた研究ではトモグラフィーの手法によって限界まで精度を高めるためには少なくとも 10 個程度のガイド星が必要であると述べられている [13]。したがって、これらの研究か
らも３つの GS の測定値に基づいた波面推定では推定に必要な情報が不足しており、推定精度として十分な性能を引き出すためには数十個の GS を配置する必要があると考えられる。
また、トモグラフィーモデルでは大気を独立な無限に薄い層の足し合わせで表現できると仮定しているため、モデルにおいてはこの離散的な大気層の数の仮定による誤差が生じてきてしまうと考えられる。モデルにおいて何層の大気揺らぎを仮定した時に十分な推定精度を得られるかについて
は、本論文におけるデータに基づいた手法を用いてさらに GS を増やした時の波面推定誤差の振る舞いから検証していきたい。

55

第 5. 議論

5.4.3

5.4. 推定精度を決める要因

タイムステップデータを用いた推定法との関係

次に、4.4 で見たようにタイムステップデータを用いた推定によって推定精度が上がったことについて考察していく。タイムステップデータを加えて波面推定を行なった時に波面推定誤差が小さくなったことから、このタイムステップデータが大気揺らぎの推定に寄与する情報を持っていると考えられる。以下の図 5.15 はある時刻の CL-WFS2 の 144 個の測定値およびそこから ∆t[frame] ずらした時刻での CL-WFS2 の 144 個の測定値の相関を表している。各行が 72 個のセンサの x 方向の測定
値および y 方向の測定値、各列が ∆t[frame] だけ前の測定値に対応している。

∆t = 0、0 つまり同時刻の測定値同士の相関を見ると、対角成分の値は 1 となっているほか、x 方
向同士の測定値に弦のような成分が確認できる。これは x 方向同士に隣り合う測定値が強い相関を持っていることを表している。また、y 方向同士の測定値では対角成分の周りにも強い成分が確認できる。これは y 方向同士に隣り合う測定値が強い相関を持っていることを表している。また、∆t を増やしていくとある時刻の y 方向の測定値と 10 から 15frame 前の y 方向の測定値との間に新しく弦のような成分（x 方向に隣り合う成分）に強い相関が現れてくることがわかる。以下ではこのタイム
ステップデータとの間に新しく現れた強い相関について考えていく。
今いくつかの観測から 60〜100nm の時間内であれば大気揺らぎは時間変化しない frozen ﬂow”と考えることができ、この時間内であれば大気の時間変化は大気揺らぎの層の風速による移動だけで表せることがわかっている [16]。この短いタイムスケールにおける大気揺らぎの時間変化を表した概念図 5.14 である。短いタイムスケールで大気の各層の変動が風による移動だけであるという”frozen

ﬂow”の仮定をすると、t = t1 における測定値は t = t1 − ∆t の大気揺らぎが風速で移動したぶんだけずれた状態の情報として得られる。そのため大気の構造が保持されたまま風速によって x 方向に移動した場合には t = t1 と t = t1 − ∆t の x 方向にずれた測定値の間に相関が見えてくると考えられる。ある時刻の y 方向の測定値と 10 から 15frame 前の y 方向の測定値との間に新しく現れた弦のような成分はこれを示していると考えられる。

図 5.14: 時間経過による大気揺らぎの層の変化

56

第 5. 議論

5.4. 推定精度を決める要因(a) ∆t = 0(b) ∆t = 1(c) ∆t = 5(d) ∆t = 10(e) ∆t = 15(f) ∆t = 30(g) ∆t = 60(h) ∆t = 250

図 5.15: ある時刻の CL-WFS2 の 144 個の測定値およびそこから ∆t[frame] ずらした時刻での CL-

WFS2 の 144 個の測定値の相関マップ。左上から右下にかけて順番に ∆t = 0、1、5、10、15、30、60、250 だけずらした時刻との相関を表している。

57

第 5. 議論

5.4. 推定精度を決める要因

次に、タイムステップデータを加えた時に推定精度が上がったことについて、図 5.14 から議論していく。NGS 方向の測定値を考えると、図 5.14 の左の時刻 t = t1 においてサイエンスターゲット方
向からの光が通過する領域のうちピンクで示した領域の大気の情報は NGS の測定値にはのらない。
しかし、”frozen ﬂow”の仮定をすると、ピンクで示した領域は時刻 t = t1 − ∆t においては同じ状態
で少しずれた位置に現れる。これは風が図 5.14 の方向に速度を持っているとした場合、現在の時刻(t = t1 ) では測定にのらなかった大気の情報を含んだ測定値が時刻 t = t1 − ∆t のタイムステップデー
タから得られることを表している。そのため、現在の測定値に加えて過去のタイムステップにおける測定値を計算に取り入れることで、測定点を実効的に増やすことができ、これによって波面推定精度が向上したと考えられる。また、15frame 前のデータを用いた場合に最も推定精度が向上したことについて、図 5.15 の相関マップでは 1frame 前の測定値との相関に関してはこの成分が見えないことから、短いタイムスケールでは風速による大気の移動が十分ではなく、追加の情報としても不十分である可能性がある。逆に、30、60frame 前の測定値との相関に関してはこの成分が見えなくなるだけでなく、同じセンサ同士に対応する対角成分の相関も弱くなっていることから長いタイムスケールでは
大気の構造が崩れていき、推定に寄与する追加の情報として意味のないものになっていってしまうと考えられる。これは”frozen ﬂow”の仮定のタイムスケールとも一貫している。この二つの効果によるトレードオフで最も推定誤差が小さくなるのが 15frame 前のデータを用いた時であったと考えられる。
また図 5.16 は各 WFS の測定値の自己相関をプロットしたものである。図 5.16 に示したように

RAVEN のデータからも大気揺らぎ由来の測定値に対する自己相関（赤いプロット）が時間とともに小さくなっていっていき大気の状態が崩れていく様子が確認できる。また自己相関が 0.7 以上を保つ
タイムスケールがおよそ 15frame(60ms) であることがわかる。青のプロットが 1 秒後に対しても大きな相関を示しているのは、オリジナルのデータが 5.1 で議論した装置由来の固定的な収差を含んでいるからである。

図 5.16: WFS 測定値データの自己相関。左上から CL-WFS1、CL-WFS2、CL-WFS3、OL-WFS1、

OL-WFS2 の測定値に対する結果を示している。また青のプロットはオリジナルデータに対して自己相関をとった結果、赤のプロットは時間平均を引いたデータに対して自己相関をとった結果を表している。

58

第 5. 議論

5.5. 推定誤差の評価

推定誤差の評価

5.5

最後に、将来の LTAO のデザインに向けて本研究のデータセット 20140514043709 に基づく波面推定で得られた推定誤差の定量的な評価を行う。まず最終的な MOAO の波面推定誤差を以下の式で与えられるように様々な要因から生じる独立な収差項の二乗和で表して評価する Independent Terms

method(IT method) が Martin, 2017[13] の中で提案された。この中で実際に観測的に求めた波面推定誤差と比較することにより IT mehod による波面推定誤差の評価の妥当性が示された。そこで本論文においても、この IT method を用いて最終的な波面推定誤差 σerr を評価することにする。
2
2
2
2
2
2
σerr
= σAliasing+ σN
oise + σT omography + σF itting + σImplementation(5.1)

各誤差項について Andersen, 2012[14] を参考に考えていく。まずは、各誤差項の概要を以下に記述した。

σAliasing
エリアシング誤差。サンプリング間隔が荒いことでのこの間隔の 1/2 以上の高次の収差を測定できないことによる誤差。

σN oise
読み出しノイズ、背景ノイズに起因する誤差。この項の大きさは GS の明るさとサンプリング
周波数に依存する。

σT omography
トモグラフィック推定行列 R を用いた推定による誤差。ここまで議論してきたようにガイド星の不足などにより十分な推定ができないことや、この後の節で出てくるように大気の構造の変動に伴う R の構造の変動に起因すると考えられる。

σF ittingフィッティング誤差。これは DM が有限な大きさを持っているために、これよりも小さいスケー
ルの揺らぎを再現できないことに起因する。

σImplementation
この他にも AO のループにおいて GS の測定値から補正量を計算して、DM を変形させるまでに時間がかかることで生じる誤差や DM の open-loop 制御に起因する実装上の誤差をまとめて表記した。

closed-loop では DM の後ろに WFS を置いて DM で補正後の波面を測定し、残差を補正する
ようにフィードバックをかけて DM を変形させる。そのため closed-loop では DM の制御に起因する収差は小さい。一方で、トモグラフィック波面推定の場合のように測定方向と補正方向が異なる場合は open-loop 制御を行う。open-loop 制御では残差ではなく波面自体の推定および DM による補正を行うため、DM の変形量が直接補正精度に影響する。この DM の実際の変形量と制御モデルとのずれから生じる収差を表している。

59

第 5. 議論

5.5. 推定誤差の評価

データセット 20140514043709 に対して各誤差項を以下で概算していく。

σF ittingフィッティング誤差は以下のように近似される。

σF2 itting ≈ 0.25(d0 /r0 )5/3

(5.2)

ここで、d0 は DM で補正可能な最小のスケールを表しており、今回の場合は WFS の副開口の大きさ (d0 = 0.8m) がこれに対応する。r0 はフリード長を表しいる。フリード長とは局所的に波面の構造が保たれているとみなせる最小のスケールを表し、今回の場合は観測から波長500nm で r0 = 0.224m という値が得られている。コロモゴロフ乱流を仮定するとフリード長は波長の 6/5 乗に比例するため、H バンド (1650nm) では r0 = 0.654m となる [15]。また、定数項は DM の構造から決まり、今回は [14] に基づいて 0.25 を用いている。したがって、フィッ
ティング誤差は以下のように概算できる。
√
σF itting ≈ 0.25(d0 /r0 )5/3 = 0.591[rad] = 155[nm](5.3)

σAliasing2
エリアシング誤差は σAliasing= 0.1(d0 /r0 )5/3 で見積もることができる [14]。これを用いると、

エリアシング誤差は以下のように概算できる。
√
σAliasing ≈ 0.1(d0 /r0 )5/3 = 0.374[rad] = 98[nm](5.4)

σN oise
測定誤差は天体の明るさ m とサンプリング時間 Fs に依存し天体が明るく、サンプリング時間が長いほど小さくなる。Andersen, 2012[14] の中で、シミュレーションで 3 つの GS を用いた推定において WFS の測定誤差が最終的な波面推定誤差へどれほど影響するかについて見積りが行われた。その結果は m=12、Fs = 500Hz の場合には σN oise = 54nm、m=14、Fs = 180Hzの場合には σN oise = 95nm となった。
今回は測定誤差による影響をこのシミュレーションの結果に基づいて見積もる。データ 20140514043709
では 1 秒間に 250 回の測定が行われている。また天体の明るさを知るために Gaia Archive から本研究で測定された天体の等級 (G バンド) を調べ、その結果を以下にまとめた。ここから、荒い概算ではあるが m=12、Fs = 250Hz としてシミュレーションの結果との比較から

σN oise = 54 ∼ 95 として考えることにする。
WFS

等級 (G バンド)

OL-WFS1

11.943224

OL-WFS2

11.304723

OL-WFS3

11.499274

CL-WFS1

12.402215

CL-WFS2

12.692897

表 5.6: データ 20140514043709 に対して Gaia Archive から調べた各 WFS で測定した天体の等級 (G
バンド)。

60

第 5. 議論

5.5. 推定誤差の評価

σT omography
本研究で評価していた波面推定誤差には、σT omography および σAliasing 、σN oise による誤差項が含まれていたと考えられる。そこで、
22W F E 2 = σT2 omography + σAliasing
+ σN
oise

(5.5)

と考えて、σT omography の概算を行う。
最も推定精度が良かったのはタイムステップデータを用いた推定で、波面推定誤差〜300nm という結果が得られた。この結果および上で概算した σAliasing 、σN oise の値を用いるとトモグラ
フィー誤差は σT omography = 267〜278nm と概算される。
上で概算した誤差項および RAVEN の実験室でのシミュレーションで得られた σImplemetation の
大きさを以下にまとめた。また、最終的にこれらの誤差の二乗和で得られる波面推定誤差 σerr も以下にのせた。
誤差項

大きさ (nmrms)

σT omography

トモグラフィー誤差267〜278

σAliasing

エイリアシング誤差

98

σN oise

測定誤差

54〜95

σF ittingフィッティング誤差

155

σImplemetation

実装上生じる誤差

107

σerr

合計

354

表 5.7: データセットの条件から概算した各誤差項の大きさおよび RAVEN の実験でのシミュレー
ションから得られた各誤差項の大きさ

ここで、概算された LTAO での最終的な収差を評価するためにストレル比 (Strehl Ratio: SR) を考える。SR は点光源を結像した際に、像の中心強度が理想的な回折限界の像に対する比である [21]。

SR は波面収差の RMSφ と波長 λ から以下のように近似することができる。
SR ≈ exp[−(

2πφ
)]
λ

(5.6)

この SR を用いて補正後の波面の収差を評価する。現在の AO は赤外線領域での補償を行なっており、H バンド (1650nm) で考えると上記の波面誤差は SR∼0.26 と求まる。1.2.2 で見た LTAO におけるサイエンスからの要求 SR∼0.5 と比較すると、本手法ではこれを達成することができなかった。

RAVEN のシステムではトモグラフィー誤差を 0 にしたとしても最終的な波面誤差は 218nm となる
ため SR∼0.5 を達成することは難しいと思われるが、今後はさらに GS を増やして推定精度の評価を行いデータに基づいた波面推定の精度限界を検証していきたい。

61

第 5. 議論

5.6

5.6. トモグラフィック推定行列 R の時間依存性

トモグラフィック推定行列 R の時間依存性

トモグラフィック推定行列 R の導出の際にはここまで見てきたように、ある長さの訓練データを用いた推定を行ってきた。この時に GS 方向の測定値と同時刻のサイエンスターゲット方向の測定値を線形関係で結びつけ、大気揺らぎの 2 次元マップの推定を行った。これまでは大気揺らぎの時間変動を考慮していなかったが、現実の観測下では大気揺らぎは時間変動しているためこの変動に応じてR も変化すると考えられる。R は大気揺らぎの構造の特徴量を内包しているため、R の構造の変化を比較することで大気揺らぎの構造の変動を理解することができるはずである。
そこで図 5.17 に同じ観測 (データ 20140514043709) に対して異なる時刻で訓練した R およびそれらの差分をマップした。これを見ると訓練時刻の違いによる R の構造に特徴的な差異は一目では確認できなかった。この R の構造の差分についてサイエンスターゲット 2 方向の 12 番目、17 番目、56
番目、61 番目の x 方向の測定値に対して、異なる訓練時刻 (０〜50 秒および 50 100 秒) から推定した R を用いた場合に OL-WFS1〜OL-WFS3 の x 方向の各測定値の寄与がどれだけ異なっているかを図 5.18 に示した。ここからも訓練時刻の違いによる R の構造の変化のパターンを読み取ることはできなかった。これには、R の推定においてある長さの訓練時間をとってしまっているために、大気の小さな時間変動成分は平均化されて消えてしまっている可能性が考えられる。一方で、予測デー
タに対する訓練データの時刻と波面推定誤差の関係を示した図 5.19 の結果から、予測データに対して直前の時刻におけるデータ (∆ = 0) を訓練データとして用いて R の推定を行った時に波面推定誤差が小さくなることがわかった。このことから R の構造の中に大気揺らぎの時間的に変動している
構造を反映した特徴量が何かしら残って含まれていると推察できる。そこで今後は R の時間に依存した特徴量の抽出を行なっていきたい。

図 5.17: データ 20140514043709 のサイエンスチャンネル 2 方向に対する R の構造。左上：50〜100秒を訓練時刻として推定した R の結果、右上：100〜150 秒を訓練時刻として推定した R の結果、
下：左２つの R について差分をとったマップを示している。

62

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

図 5.18: データ 20140514043709 のサイエンスチャンネル 2 方向の測定値に対して異なる訓練時刻(０〜50 秒および 50 100 秒) から推定した R を用いた場合に OL-WFS1〜OL-WFS3 の各測定値の寄与がどれだけ異なるかを示したグラフ。上から順に 12 番目、17 番目、56 番目、61 番目の測定値に対する結果を示している。

63

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

図 5.19: R の推定時刻と波面推定誤差の関係。横軸は予測データに対して、どれだけ前の訓練デー
タを用いて R の推定を行ったかを表している。今予測データとしては 165〜166 秒を用いているので、例えば ∆ = 0 では 115〜165 秒を訓練時刻として R の推定を行った時の結果を示している。上段はデータ 20140514043709 に対する結果、下段はデータデータ 20150626063402 に対する結果を表している（左：サイエンスパス１、右：サイエンスパス２）∆ = 0 に対する波面推定誤差が最も小さくなっている。

64

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

異なる観測 (データ 20140514043709 とデータ 20150626063402) でサイエンスチャンネル 2 方向に対して推定した R の構造の比較を図 5.20 に示した。また R からサイエンスチャンネル 2 方向の 12
番目、61 番目の x 方向の測定値に対する、OL-WFS1〜OL-WFS3 の x 方向の各測定値の寄与の大きさを 5.18 に示した。
データ 20140514043709 で推定した R の構造では地上層の寄与にあたる同じ位置と比較して高層大気の寄与にあたるずれた位置にも強い成分が見えるのに対して、データ 20150626063402 で推定した R の構造では高層大気にあたる成分が弱くなっている。これはデータ 20150626063402 の観測において大気揺らぎの地上層成分が強いためであると考えられる。そこで、Wilson, 2002[17] で提案された SLODAR という複数の SH-WFS の測定値を用いて大気揺らぎの強度を推定するシステムを使って推定された大気揺らぎの高度プロファイルを以下の図 5.22 に見てみた。なお、これらのデー
タは秋山研究室の先輩大野さんより提供して頂いた。この高度プロファイルの上段および下段はデー
タ 20140514043709 およびデータ 20150626063402 に対する推定結果を表しており、それぞれのデー
タに対して時刻を 60 秒ずらして推定した高度プロファイルの変動についても示した。この高度プロファイルから、高層大気の強度にはデータ間で大きな違いは見られなかった。一方で、地上層成分についてはデータ 20150626063402 で強い強度を持っていることがわかった。R に現れた構造はこの大気揺らぎの構造を反映していると考えられる。
また、データ 20140514043709 から推定した R の構造の中で OL-WFS1 の x 方向の測定値の寄与がほとんど見られない結果となっていた (図 5.20 左)。図 5.6 からもサイエンスターゲット 2 方向の

12 番目と 61 番目の x 方向の測定値に対する OL-WFS1 の寄与がほとんどないことがわかる。これに関しては OL-WFS1 の x 方向の測定値に大きなノイズがあったためではないかと考えている。これを調べるためにデータ 20150626063402 に対して SN 比 7 のホワイトノイズを OL-WFS1 の x 方向の測定値に付加してから推定を行なった。その結果を図 5.6 に示した。データ 201405140437909 と同様に OL-WFS1 の x 方向の測定値の寄与が小さくなっていることがわかる。

図 5.20:

訓練時間 50 秒を用いて推定したトモグラフィック推定行列 R の構造 (左：データ

20140514043709 に対する結果、右：データ 20150626063402 に対する結果)

65

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

図 5.21: データ 20140514043709 とデータ 20150626063402 のサイエンスチャンネル 2 方向に対して推定した R における 12 番目、61 番目の測定値に対する、OL-WFS1〜OL-WFS3 の各測定値の寄与の大きさを示している。

66

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

図 5.22: 大気揺らぎの高度プロファイル。縦軸が大気揺らぎ層の高度、横軸が大気揺らぎの強度を表している。上段はデータ 20140514 に対する推定結果 (左：20140514043659 から 60 秒間に対する結果、右：20140514043759 から 60 秒間に対する結果) 下段はデータ 20150626 に対する推定結果 (左：20150626063402 から 60 秒間に対する結果、右：20150626063502 から 60 秒間に対する結果)

図 5.23: データ 20150626063402 のサイエンスチャンネル 2 方向に対して SN 比 7 のホワイトノイズ
を OL-WFS1 の x 方向の測定値に付加してから推定を行なった結果。

67

第 5. 議論

5.6. トモグラフィック推定行列 R の時間依存性

最後に観測コンディションの推定精度に対する影響について考えていく。まず、データ 20140514043709とデータ 20150626063402 に対して R の推定に用いる訓練時間と波面推定誤差の比較を 5.24 で行なった。訓練時間を 100 秒間とった場合について見ると、データ 20150626063402 の波面推定誤差はデー
タ 20140514043709 で得られる波面推定誤差に対してサイエンスチャンネル 1 方向で 100nm、サイ
エンスチャンネル 2 方向で 150nm 大きくなっていることがわかる。
データ間で推定精度に大きな違いが見られたの原因の一つには上で見たように大気揺らぎの地上層成分の大きさが関係していると考えられる。データ 20150626063402 では地上層の成分が高層の成分に対して大きいため高層大気の影響が埋もれてしまい推定がうまくできていない可能性がある。

MOAO におけるトモグラフィー推定の性能評価を行った Martin, 2017[13] の論文の中でも大気揺らぎの地上層の成分が比較的小さい場合 (Seeing[”]<0.6) には SR ∼ 0.3(@H バンド) を得られたの対して、大気揺らぎの地上層の成分が大きい場合 (Seeing[”]>1.2) には SR ∼ 0.1(H バンド) ほどの精度しか得られないことがわかっている。
二つ目の要因としては大気揺らぎの変動が考えられる。上の図 5.20 で見たように R の構造はガイド星の配置や大気揺らぎの構造などの観測コンディションに応じて変動する。ガイド星の配置は一度の観測においては固定されているが、大気揺らぎは短いタイムスケールで変動している。この変動が大きい場合には R の推定を行う時に用いる訓練時間内でも大気揺らぎの構造が大きく変動してしまう。これによって R の推定精度が劣化がしてしまっている可能性がある。実際に、図 5.22 の 60 秒
間ごとの地上層成分の強度の比較から、データ 20140514043709 では地上層成分を含めて大気の状態が比較的安定しているのに対して、データ 20150626063402 では地上層成分が大きく変動していることがわかる。
以上をまとめるとデータに基づいた波面推定の精度が大気揺らぎの地上層成分の強さや変動に依存している可能性が示唆されたが、観測コンディションと推定精度の関係についてはさらに多くのデータを用いた議論を行なう必要がある。

図 5.24: 訓練データを増加させていったときのデータ 20140514043709 とデータ 20150626063402 に対する波面推定誤差の推移の比較 (左：サイエンスパス１、右：サイエンスパス２)

68

第 6 章 結論
本論文ではすばる望遠鏡 RAVEN による大気揺らぎトモグラフィック測定データを用いて、波面推定法の検討を行った。
本論文では最小二乗法によって複数のガイド星の波面センサの測定データからターゲット方向の波面センサの測定データを結びつけるトモグラフィック推定行列 R の推定を行った。まず、R の推定の際の計算手法に関して訓練データが短い場合に対して以下の二つの手法を導入した。

1. R の推定に際に推定に寄与する成分とそれ以外の成分を分けるために特異値分解を行い低特異値の除去を行い、推定精度を向上させることができた。またその中で Yamazaki, 2015 における特異値除去の基準に問題があったことを明らかにし、累積寄与率に基づいた手法を導入した。

2. 訓練データが短い時に R が推定した時刻の特徴を過学習してしまうことによる推定精度の劣化を確認した。それに対して正則化の手法を用いることで過学習の影響を改善することができた。
次に、訓練データが十分に長い時に、推定精度の限界を決める要因についての議論を行った。その結果、トモグラフィック波面推定の推定精度に関わる大きな指標として推定に寄与する大気揺らぎの情報量があることがわかった。情報量を増やすために

1. GS の数と配置に関する議論を行った。ここでは GS の数と配置が推定精度に大きく関わってくることが本手法によっても確認された。また、GS によってカバーされる推定に関わる領域
を大きくするために、GS の配置をサイエンスターゲットから近い距離におくことと、方向をずらしておくことが重要であることも確認された。

2. タイムステップデータを用いた推定法の提案を行った。その結果、過去の測定値を現在の測定に反映させた推定を行うことで実効的に測定点を増やし、推定精度を向上させることができた。
この際、大気の変動のタイムスケールを考慮するとおよそ 60ms 前の測定値を追加した時に推定精度としては最も大きく向上した。
今後の展望としては本手法を用いてさらに GS を増やした時の波面推定誤差の振る舞いを見ること
で、推定精度限界の評価を行なっていきたいと考えている。ここから必要な GS の数の見積もりなどを行なっていきたい。さらに、データに基づいた波面推定の精度が大気揺らぎの地上層成分の強さや変動に依存し、これらが大きい場合に R による波面推定精度が劣化してしまう可能性が示唆された。そこで、R の時間発展成分についてより詳細な解析を行うことで、よりリアルタイムの大気揺らぎの特徴を反映した R の推定を行なっていきたいと考えている。

69

謝辞
大学院に入った頃の私は、自分の了見の狭さを棚に上げ、基礎的な勉強もせずに、大学から長い間遠ざかっていました。指導教員の秋山正幸教授は、そんな私を見捨てることなく、私が研究の面白さに気づくまで辛抱強く 4 年間も見守ってくださいました。私が研究室を訪ねる度に、勉強不足の私が納得するまで議論に付き合ってくださった秋山教授にはどんなに感謝しても感謝し尽くせません。
最後まで基礎の勉強不足が残り、研究に対する姿勢という面でも未熟な私でしたが、この修士論文への取り組みを通して、研究に対する姿勢や研究の面白さに気づくことができました。人としての生き方という面でも、秋山教授から多くを学び成長できたと思っています。本修士論文をはじめ、学部４年生の時から５年もの間大変お世話になりました。ほんとうにありがとうございました。
また、この修士論文を仕上げる中で、本当に多くの方が力を貸してくださいました。秋山研究室の先輩大野良人さん、山崎公大さんからは本修士論文を書くにあたり必要なデータを提供して頂いた
他、多くのアイデアを頂きました。秋山研究室の先輩方、同期、後輩には発表の機会の度に研究に対する温かい助言や議論を頂きました。天文学専攻以外でも国立情報学研究所の速水謙先生、Ning

Zheng さん、東北大学経済学研究科の松田安昌先生、佐藤宇樹さんには貴重なご助言をいただきました。みなさんの助けがなかったならば、本修士論文を仕上げることはできませんでした。
さらに、研究がなかなか思うように進まず、就活を含め将来の不安で心が折れそうになった時にも多くの方々が支えてくださいました。関わってくださった方々にこの場をお借りして感謝申し上げます。
最後に、４年間私の考えを常に尊重して支え続けてくれた、両親や家族にも感謝いたします。

70

AppendixA

A.1

行列の最小二乗解の導出

今、行列 X から行列 Y を推定する行列 R を考える。

Y = RX(A.1)

この時、行列 R の推定解を最小二乗法によって導出する方法を以下に示す。

行列 R の推定解として、(A.1) 式の両辺の差 Y − RX の大きさが最小になるものを考える。そのために、行列の大きさを表す量として、以下のフロベニウスノルム (Forbinus Norm) を用いる。
√∑ ∑|aij |2(A.2)∥A∥F =ij

また、以下表記を簡単にするために F を次のように定義する。

F = ∥Y − RX∥2F(A.3)

フロベニウスノルムの定義から ∥A∥2F = tr(AAT ) が成り立つことが直ちにわかる。この性質を用
いると

F

=

∥Y − RX∥2F=

tr[(Y − RX)T (Y − RX)]=tr[Y T Y − X T RT Y − Y T RX + X T RT RX](A.4)となることがわかる。この行列のノルム F を最小にする R を行列 R の推定解 R̂ として求める。

R̂ = argmin F

(A.5)

R

この推定解 R̂ は、F を R で微分した時の極値という条件から求めることができる。すなわち、

∂F
∂R

∂
tr[Y T Y − X T RT Y − Y T RX + X T RT RX]
∂R
∂
∂
∂
= −
tr[X T RT Y ] −tr[Y T RX] +tr[X T RT RX]
∂R
∂R
∂R
= 0
=

ここでスカラー F を行列 R で微分することを以下のように定義する。
(
)
∂F
∂F
=
∂R pq
∂Rpq

(A.6)

(A.7)

すなわち、F を R の (p,q) 成分で微分した値を (p,q) 成分にもつ行列ができる。以下、式 (A.6) 式の各項について計算していく。以下、数式はアインシュタインの縮約記法に従って表記するものとする。

71

A.1. 行列の最小二乗解の導出

APPENDIX A.

tr[X T RT Y ]

=

T T
Xij
Rjk Yki = Xji Rkj Yki(A.8)

tr[Y T RX]=

YijT Rjk Xki = Yji Rjk Xki

(A.9)

=

T T
Xij
Rjk Rkl Xli = Xji Rkj Rkl XliT

T

tr[X R RX](A.10)(A.8),(A.9),(A.10) 式より (A.6) 式の各項は∂tr[X T RT Y ]∂Rpq∂=(Xji Rkj Yki )
∂Rpq
∂Rkj
= Xji
Yki
∂Rpq
= Xji δkp δjq Yki(A.6) 式の第一項== Xqi Ypi
T
= Ypi Xiq

= (Y X T )pq(A.11)= tr[Y T RX]∂=(Yji Rjk Xki )
∂Rpq
∂Rjk
= Yji
Xki
∂Rpq
= Yji δjp δkq Xki(A.6) 式の第二項= Ypi Xqi
T
= Ypi Xiq

= (Y X T )pq∂=
tr[X T RT RX]∂Rpq∂(Xji Rkj Rkl Xli )
=
∂Rpq
∂Rkj
∂Rkl
= Xji
Rkl Xli + Xji Rkj
Xli
∂Rpq
∂Rpq
= Xji δkp δjq Rkl Xli + Xji Rkj δkp δlq Xli(A.6) 式の第三項(A.12)= Xqi Rpl Xli + Xji Rpj Xqi
TT= Rpl Xli Xiq+ Rpj Xji Xiq= 2(RXX T )pqと求まる。ただし、

{
δij =

1

(i = j のとき)

0(i ̸= j のとき)

(A.13)

以上の (A.6),(A.11),(A.12),(A.13) 式より∂F= −2Y X T + 2RXX T = 0∂R(A.14)

を満たす R が行列の推定解 R̂ であることから

R̂XX T = Y X T
72(A.15)

A.2. L2 正則化

APPENDIX A.

ここで行列 X は行に対して最大階数を持つので XX T は可逆な正方行列となる。したがって R̂ は次のように求まる。

R̂ = Y X T (XX T )−1

A.2(A.16)

L2 正則化

過学習の問題を抑える手法として L2 正則化を考える。式 (A.3) に対して R の二乗ノルムを加えたF = ∥Y − RX∥2F + λ∥R∥2F(A.17)を考え、このノルム F を最小にする R を行列 R の推定解 R̂ として求める。ここで、λ は正則化パラメータである。つまり、

R̂ = argmin F = argmin[∥Y − RX∥2F + λ∥R∥2F ]R(A.18)R

ここで、最小二乗解の導出と同様に考えると推定解 R̂ は

∂F
∂R

∂
∂
tr[Y T Y − X T RT Y − Y T RX + X T RT RX] −tr[RT R]
∂R
∂R
∂
∂
∂
∂
= −
tr[X T RT Y ] −tr[Y T RX] +tr[X T RT RX] −tr[RT R]
∂R
∂R
∂R
∂R
= 0
(A.19)=

で与えられる R として求まる。式 (A.19) の前半 3 項については最小二乗解の導出のところで考えたので、ここでは 4 項目についてみていく。

tr[RT R] =T
Rij
Rji = Rji Rij(A.20)(A.20) 式より (A.19) 式の 4 項は(A.19) 式の第 4 項=∂
tr[RT R]
∂Rpq∂
λ(Rji Rij )
∂Rpq
∂Rji
∂Rji
+
Rji )
λ(Rji
∂Rpq
∂Rpq
λ(Rji δjp δiq + δjp δiq Rji )

=

2λRpq

=

2λ(R)pq===λ(A.21)となる。(A.8),(A.9),(A.10),(A.21) 式より (A.19) 式は∂F= −2Y X T + 2RXX T + λR = 0∂R(A.22)

したがって、L2 正則化によって求まる推定解 R̂ は

R̂(XX T + λI) = Y X T(A.23)

R̂ = Y X T (XX T + λI)−1(A.24)

を満たす。以上より

73

A.3. ゼルニケ多項式

APPENDIX A.

A.3

ゼルニケ多項式

A.3.1

ゼルニケ多項式

以下では波面収差の関数表現として用いられるゼルニケ多項式 (Zernike Polynomials) について [20]を参考にまとめる。
ゼルニケ多項式の最大の特徴は単位円内部で直交する直交関数系であることで、ゼルニケ多項式を用いると単位円上のすべての関数を表現することができる。光学システムはレンズや望遠鏡の開口に代表されるように円形の素子から構成されているため、波面収差をこのゼルニケ多項式を用いて評価することが多い。
ゼルニケ多項式 Znm (ρ, θ) は極座標 (ρ、θ) で表され、以下の式で定義される。
{
Rnm (θ)cos(mρ)(m ≥ 0)m
Zn (ρ, θ) =|m|
Rn (θ)sin(|m|ρ)(m < 0)(A.25)

ここで n は自然数、m は n ≤ |m| を満たす整数、動径座標 ρ の定義域は 0 から 1 までである。まm
た動径多項式 Rn(θ) は以下の式で定義される。
 n−m

2

(−1)k (n − k)!
 ∑
θn−2k
m
n+m
n−m
Rn (θ) =
k!(
−
k)!(
−
k)!
2
2
k=0



0

(n − m : 偶数)(A.26)(n − m : 奇数)

ゼルニケ多項式は規格化条件によっていくつかの記述法があるが、ここでは本論文で用いた Standard

Zernike 多項式と呼ばれる表記法を用いた。以下の図 A.3.2 および表 A.1 に低次のゼルニケ多項式を
まとめた。

A.3.2

波面収差

ここでは波面収差とゼルニケ多項式の関係について説明する。まず、図 A.3.2 に波面収差の概念図を示した。物体上の１点 P を出た光は球面波として媒質中を伝搬する。Wo は球面波の等位相面を表しており、これを波面と呼ぶ。球面波として伝搬した波面はレンズが理想的に無収差の場合には、レンズを通過後に P の共役点 P’ を中心とする球面波 Wo′ となる。しかし、レンズに収差がある場合にはレンズを通過後に球面波 Wo′ から ∆W だけずれた歪んだ波面 Wi として変換される。この本来の波面からのずれ ∆W を波面収差と呼ぶ。波面収差 ∆W は瞳座標 (ξ, η) の関数として表され、これを W (ξ, η) とおく。
ゼルニケ多項式の座標 (ρ, θ) と瞳座標 (ξ, η) の関係は

ξ = aρcosθ(A.27)

η = aρsinθ

(A.28)と表せるため、波面収差 W(ξ, η) はゼルニケ多項式を用いて以下のように記述される。

W (ξ, η) =∑cj Zj (ρ, θ)j74(A.29)

A.3. ゼルニケ多項式

APPENDIX A.

図 A.1: 波面収差の概念図

各多項式の係数 cj をゼルニケ係数と呼ぶ。各多項式に対応する収差を以下で説明する。また表 A.1にまとめた。
定数項：Z1
収差としての意味はなく無視してよい。

tilt 成分：Z2 、Z3
波面の傾き (tilt) を表しており、物体上の１点の結像を考える時には収差ではなく像の横ずれを示す。
デフォーカス：Z5
物体上の１点の結像を考える時には、像の光軸方向へのずれを示す。
非点収差：Z4 、Z6像面内の直交する 2 方向でピント合う位置がずれる収差を表している。Z4 の項は X 軸方向とY 軸方向でのピントのずれ、Z6 の項は X 軸と ±45 度をなす軸方向でのピントのずれを示す。
コマ収差：Z8 、Z9レンズ系に偏心があった場合に光軸からはずれた視野周辺で発生し、像がボケるだけでなくボケかたが像面内で非対称になる。

75

A.3. ゼルニケ多項式

APPENDIX A.

式番号

n

m

1

0

0

1

ピストン

2

1

1

ρcosθ

傾き (tilt) X 成分

3

1

-1

ρsinθ

傾き (tilt) Y 成分22

ρ cos2θ

非点収差 (0 度と 90 度方向)

02

2ρ − 1

デフォーカス-2

2

ρ sin2θ

非点収差 (±45 度方向)

4
5
6

2
2
2

ゼルニケ多項式

収差

3

7

3

3

ρ cos3θ8

3

1(3ρ2 − 2)ρcosθ

3 次のコマ収差X 成分-1(3ρ − 2)ρsinθ

3 次のコマ収差

Y 成分

9
10

3
3

-3

2

3

ρ sin3θ

表 A.1: Standard Zernike 多項式

図 A.2: ゼルニケ多項式

76

参考文献

[1] Fabrice Vidal, Eric Gendron, and Gérard Rousset, ”Tomography approach for multi-object
adaptive optics,” J. Opt. Soc. Am. A. Vol. 27, pp. A253-A264, 2010
[2] F.Vidal, E. Gendron, M. Brangier, A.Sevin, G.Rousset, and Z.Hubert, ”Tomography reconstruction using the Learn and Apply algorithm,” 1st AO4ELT conference, 07001-p.1-6, 2010
[3] F. Assémat, E. Gendron and F. Hammer, ”The falcon concept: multi-object adaptive optics
and atmospheric tomography for integral ﬁeld spectroscopy. principles and performances on an
8 meter telescope,” Mon.Not.Roy.Astron.Soc. Vol. 376, pp. 287-312, 2006
[4] Masayuki Akiyama, Shin Oya, Yoshito H. Ono, Hideki Takami, Shinobu Ozaki, Yutaka Hayano,
Ikuru Iwata, Kazuhiro Hane, Tong Wu, Tomoyasu Yamamuro and Yuji Ikeda, ”TMT-AGE:
wide ﬁeld of regard multi-object adaptive optics for TMT,” Thirty Meter Telescope Science
Forum, Tucson Arizona,16-19 July, 2014[5] Sergio Luis Suárez Gómez, Jesús Daniel Santos Rodrı́guez, Francisco Javier Iglesias Rodrı́guez
and Francisco Javier de Cos Juez, ”Analysis of the temporal structure evolution of physicalsystems with the self-organising tree algorithm (SOTA): application for validating neural network systems on adaptive optics data before on-sky implementation,” Entropy, Vol. 19, No.
103, pp1-11, 2017[6] K Yamazaki, ”Evaluation of atmospheric tomography method by empirical reconstruction ofSubaru/RAVEN on-sky measurements,” Master thesis, 2015
[7] CARLOS GONZÁLEZ GUTIÉRREZ, MARÍA LUISA SÁNCHEZ RODRÍGUEZ, RAMÓNÁNGEL FERNÁNDEZ DÍAZ, JOSÉ LUIS CALVO ROLLE, NIEVESROQUEÑÍ
GUTIÉRREZ, ”Rapid tomographic reconstruction through GPU-based adaptive optics,”
Oxford University Press, Vol. 27, No. 2, pp214-226, 2018[8] J. Osborn, F. J. D. C. Juez, D. Guzman, T. Butterley, R. Myers, A. Guesalaga and J. Laine,”Using artiﬁcial neural networks for open-loop tomography,” Optics Express, Vol. 20, Issue 3,
pp.2420-2434, 2012
[9] Andrew M. Hopkins and John F. Beacom, ”ON THE NORMALIZATION OF THE COSMIC
STAR FORMATION HISTORY,” The Astrophysical Journal, Vol. 651, pp.142-154, 2006[10] Akiyama, M., Oya, S., Hane, K., Wu, T., and Tokoku, C.,
university,

MOAO activities in Tohoku

Adaptive Optics Systems II., Proc. SPIE, Vol. 7736, 2010[11] Olivier Lardière, Reston Nash, J.-P. Markes, Dave Andersen, Colin Bradley, C
CID3elia Blain, Ronald Desmarais, Darryl Gamroth, Meguru Ito, Kate Jackson, Przemek Lach,
Laurie Pham. ”Final opto-mechanical design of Raven, a MOAO science demonstrator for
Subaru,” Proc. of SPIE Vol. 8447, 844753-p.1-15, 2012
77

参考文献

参考文献

[12] O.A. Martin, C. M. Correia, E. Gendron, G. Rousset, F. Vidal, T.J. Morris, A.G. Basden, R.M.
Myers, Y.H. Ono, B. Neichel, and T. Fusco, ”William Herschel Telescope site characterization
using the MOAO pathﬁnder CANARY on-sky data,” Proc. SPIE. Vol. 9909, Adaptive Optics
Systems V, 99093P, 2016[13] O. A. Martin, É. Gendron, G. Rousset, D. Gratadour, F. Vidal, T. J. Morris, A. G. Basden, R.
M. Myers, C. M. Correia, and D. Henry, ”Wave-front error breakdown in laser guide star multiobject adaptive optics validated on-sky by Canary,” EDP Sciences, A&A 598, A37, pp.1-17,
2017[14] DAVID R. ANDERSEN, KATE J. JACKSON, CÉLIA BLAIN, COLIN BRADLEY, CARLOS
CORREIA, MEGURU ITO, OLIVIER LARDIÉRE, AND JEAN-PIERRE VÉRAN, ”Performance Modeling for the RAVEN Multi-Object Adaptive Optics Demonstrator,” The Astronomical Society of the Paciﬁc. Vol. 124, pp. 469-484, 2012[15] Hardy, J. W., ”Adaptive Optics for Astronomical Telescope,” New York: Oxford Uni- versity
Press., 1998[16] Schock, M. and Spillar, E. J., ”Method for a quantitative investigation of the frozen ﬂow
hypothesis,” J. Opt. Soc. Am. A, 17, pp. 1650-1658. 2000
[17] R.W.Wilson, ”SLODAR: measuring optical turbulence altitude with a Shack Hartmann wavefront sensor,” Monthly Notices of the Royal Astronomical Society, Volume 337, Issue 1, pp.
103-108, 2002[18] 神澤富雄, 伊藤昇, 臼田知史, ”すばる望遠鏡の計測制御技術,” 計測と制御, 第 43 巻, 第 3 号, 2004[19] ”Subaru Telescope HP”,[https://www.subarutelescope.org/][20] 渋谷眞人, 大木裕史, ”回折と結像の光学,” 朝倉書店, 2005[21] 天文学辞典 (日本天文学会)

78