修士論文
次世代超大型望遠鏡の広視野多天体補償光学のためのトモグラフィック波面再構成の精度および計算速度に関する研究

The accuracy and computing speed of tomographic
wavefront reconstruction for Multi-Object Adaptive Optics
of Extreme Large Telescope

東北大学大学院 理学研究科 天文学専攻 修士２年
大野 良人
平成 24 年度

要旨
TMT をはじめとする次世代超大型地上望遠鏡にとって、大気揺らぎの影響を補正する補償光学
は必須の技術である。また、近年では次世代望遠鏡に向けた次世代補償光学の検討、開発も活発に行われている。特に遠方の銀河の形成や進化を解明するためには暗い、遠い天体を観測するための集光力、内部を詳細に調べるための空間分解能、統計的な議論を行うのためのサンプル数を稼ぐ効率性が要求され、多天体補償光学と TMT の組み合わせはこれらの要素をすべてクリアすることができる。
広視野多天体補償光学では複数のガイド星を用いて複数の方向の大気揺らぎを測定し、その情報からトモグラフィーの手法を用いて大気揺らぎを３次元的に推定する。この推定値を用いて各天体方向に最適化な補正を行うことで、広視野内の複数の天体に対して高分解能な同時観測を達成し、観測効率を向上することができる。現状でこの多天体補償光学を実現するためにネックとなっているのが、主にオープンループ制御に伴う高精度化と膨大な情報をを用いてリアルタイム補正を行う
ためのシステムの高速化の２点である。さらに、現在考えられているサイエンスケースに対応する
ためにはさらなる広視野化が必要である。本研究ではこれらの課題を克服するための新しい推定手法を検討した。
まず、各天体ごとに補正する前に全方向共通の DM で地表層を先に補正する手法を検討した。

TMT が建設されるマウナケア山頂のサイト調査の結果から、低高度の地表層の揺らぎのパワーが全体の揺らぎのパワーの半分以上を占めていることがわかっている。その地表層を先に補正することで、それぞれの WFS が測定する位相差や DM が補正する位相差を小さくすることができる。また、今回は地表層の推定方法としてトモグラフィー計算とは別に、複数のガイド星の測定値の平均値から求める手法を検討した。こうのようにトモグラフィーをは別で地表層を推定することでトモグラフィー計算で求める大気揺らぎの層の数を減らし、推定精度を向上させることができると考えた。検証の結果、理想的な場合ではこの手法により、すべての層をトモグラフィー計算で求める従来の手法よりも推定精度を向上させることができた。しかし、実際に複数のガイド星の平均値を用いて地表層を推定した場合には推定精度は大きく変わらなかった。これは平均値による地表層の推定精度が従来のトモグラフィーを用いた手法と同程度であることが原因であると考えられる。この手法を用いて推定精度を向上させるためには地表層をより高精度に推定する手法が必要である。
次に風速の情報を用いてある時間の推定値を次の計算に反映させる手法を開発した。過去の測定値を反映させ、トモグラフィー計算に用いる情報量を増やすことで、推定精度の向上が期待される。今回は風速の情報が既知であり、かつ大気ゆらぎ自身は時間変化しない ﬂozen frow の仮定が成り立っている理想的な状態で検証を行った。その結果、理想的な場合ならばこの手法によって推定精度が大きく向上がすることがわかった。さらにこの手法の現実性についても議論した。風速の推定手法や ﬂozen frow の仮定の妥当性に限界により多少の誤差が生じる可能性はあるが、それを考慮しても従来の手法より推定精度の向上が期待されることがわかった。
どちらの手法でも目標とする精度、視野を達成することはできなかったが、推定精度を従来の手法よりも大きく向上させることができた。また、本研究の計算には高速化を狙って 「GPGPU」による並列計算を取り入れている。

目次
第 1 章 広視野多天体補償光学

1.1
1.2

補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
広視野補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2.1
1.2.2
1.3

超大型望遠鏡と広視野補償光学 . . . . . . . . . . . . . . . . . . . . . . . .

23

1.2.3 広視野補償光学の種類 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
多天体補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 多天体補償光学実現の意義 . . . . . . . . . . . . . . . . . . . . . . . . . . .

366

サイエンスケース . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
要求される性能 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7

1.3.4 多天体補償光学の課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.5 多天体補償光学の現状 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
本論文で検討した新しい手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11
11
12

1.4.1
1.4.2

12
13

1.3.2
1.3.3

1.4

従来の SCAO の限界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
2

GLAO+MOAO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
風速の情報を用いる手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 2 章 波面再構成

2.1

仮定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2
2.3
2.4

大気揺らぎと WFS で測定する波面の関係 . . . . . . . . . . . . . . . . . . . . . .

2.5
2.6

最小分散推定 (Minimum Variance Reconstruction) . . . . . . . . . . . . . . . . .

2.7
2.8

開口の影響

WFS の測定値と波面の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
問題の定式化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
大気ゆらぎの共分散行列 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
風速の情報を用いる手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 3 章 シミュレーションの詳細

3.1

シミュレーションの設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2
3.3
3.4

シミュレーションの流れ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
大気揺らぎモデルの再現性の確認 . . . . . . . . . . . . . . . . . . . . . . . . . . .
計算手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4.1
3.4.2
3.4.3

共役勾配法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

GPGPU を用いた並列計算 . . . . . . . . . . . . . . . . . . . . . . . . . . .
共役勾配法の収束条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14
14
14
16
17
18
20
21
21
24
24
26
27
29
29
30
30

第 4 章 結果

4.1

最小分散推定によるトモグラフィック波面再構成の結果 . . . . . . . . . . . . . . .

4.1.1
4.1.2
4.1.3
4.1.4
4.1.5
4.2

4.3

4.4

推定結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

IRMOS と NFIRAOS. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ガイド星の数との関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3132

ガイド星の距離との関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1 の結果のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38
39
39

MOAO+GLAO の結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 理想的な場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41

4.2.2 平均値を用いた結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.3 4.2 の結果のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
風速を用いた結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
42
43

4.3.1
4.3.2

すべての大気ゆらぎの層に同じ風速を与えた場合 . . . . . . . . . . . . . .

43
44

4.3.3
4.3.4
4.3.5

風速のモデルで計算した場合

ある高さの大気ゆらぎの層にだけ風速を与えた場合 . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . .
風速の推定に誤差がある場合 . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 の結果のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45
46
47

計算速度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.4.1

収束条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49
49

4.4.2

各過程の計算に要する時間 . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

第 5 章 議論

5.1
5.2

31
31

52

誤差評価 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ガイド星とガイド星の間で推定が悪くなることについて . . . . . . . . . . . . . . .
5.2.1 各高さでの推定精度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52
54
54

5.3

5.2.2 上層と下層での推定精度 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
GLAO+MOAO の結果について . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57
60

5.4

5.3.1 結果が Normal Method と変わらないことについて . . . . . . . . . . . . .
5.3.2 GLAO の可能性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
風速を用いた結果について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60
61
64

ノイズの影響を抑える手法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

64
65

5.4.1
5.4.2

手法の現実性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

第 6 章 結論

69

謝辞

70

Appendix A AdaptiveOptics

71

A.1 望遠鏡による観測 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.1.1 PSF と回折限界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.1.2 大気揺らぎの影響の定量的な評価 . . . . . . . . . . . . . . . . . . . . . . .

71
71
72

A.2 大気揺らぎ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2.1 大気ゆらぎの空間スペクトル . . . . . . . . . . . . . . . . . . . . . . . . .

74
74

A.2.2 大気ゆらぎの強度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

A.2.3 大気揺らぎの波長依存性 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2.4 風速のモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75
75

A.2.5 大気ゆらぎのパラメーター . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3.1 ガイド星 (Guide Star) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

76
78
78

A.3.2 波面センサー (WsveFront Sensor) . . . . . . . . . . . . . . . . . . . . . . .
A.3.3 波面再構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81
83

A.3.4 可変形鏡 (Deformable Mirror) . . . . . . . . . . . . . . . . . . . . . . . . .
A.3.5 制御 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3.6 補償光学の現状 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

83
84
85

Appendix B Mathematics
B.1 行列微分を用いた計算の証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

86
86

B.2 共役勾配法 (Conjugate Gradient Method) . . . . . . . . . . . . . . . . . . . . . .
B.2.1 共役勾配法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87
87

B.2.2 共役勾配法の実装 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3 前処理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.3.1 対角スケーリング . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

88
89
89

第 1 章 広視野多天体補償光学

第1章

広視野多天体補償光学

本章ではまずはじめに、導入として補償光学および広視野補償光学について説明する。次に本研究で検討した多天体補償光学について、その必要性と実現における課題を述べる。最後にそれらの課題を打破するために本研究で検討した新しい手法について触れる。補償光学に関する一般的な知識については Appendix. A に記述した。

1.1

補償光学

地球大気の揺らぎの影響は地上観測の質を制限する大きな要因の一つである。補償光学 (Adaptive

Optics, AO) はこの大気揺らぎの影響をリアルタイムで補正し、望遠鏡本来の性能での観測を実現するシステムである。特に口径が直径数 m を超えるような大型望遠鏡を使用する場合、その高い分解能を生かすためには AO による補正が必要不可欠である。実際、すばる望遠鏡をはじめとする
稼働中の 8∼10m 望遠鏡には AO システムが装備されており、様々な科学的成果を出している。さらに 2020 年代に運用が予定されている TMT(Thirty Meter Telescope) をはじめとする次世代超大型望遠鏡 (Extreme Large Telescope, ELT) でも、ほぼすべての装置に AO が関係している。近年ではこの ELT に向けた新しい AO システムが世界中で検討、開発されている。
補償光学は大まかに大気揺らぎの測定、大気揺らぎの推定、大気揺らぎの補償の 3 つの要素で構成されている。
大気揺らぎの測定
まず最終的に補正を適用するターゲット天体の近くにある明るい星 (ガイド星、Guide Star、

GS) の光の位相 (波面) を波面センサー (WaveFront Sensor、WFS) で観測する。GS の光は
地球大気中を通過する際に大気の (屈折率の) 揺らぎの影響を受け、その光の位相が乱れる。
その位相波面の乱れを WFS で測定することで、GS の光が通ってきた経路の大気揺らぎの情報を得ることができる。
大気揺らぎの推定

WFS で測定した情報から、補償すべき大気揺らぎを推定する。一般的にこの推定は連立方
程式を行列の形式で解く手法が用いられる。
大気揺らぎの補償
補償は表面の形を変形させることのできる可変形鏡 (Deformable Mirror, DM) で行う。反射
したときに推定した大気揺らぎをキャンセルするような形に DM の表面を変形させ、補正したい天体をこの DM を通して観測することで大気揺らぎの影響を補正する。
現在使用されている AO のほとんどは、1 つの GS と WFS から大気揺らぎを測定するシステムである。本論文ではこれを単層共役補償光学 (Single Conjugate Adaptive Optics, SCAO) と呼ぶことにする。図 1.1 は SCAO のシステムの模式図である。

1

第 1 章 広視野多天体補償光学

1.2. 広視野補償光学

図 1.1: 1 つの GS と WFS を用いる SCAO の模式図。測定波面全体に広がる傾き (Tip/Tilt) の成分は、その他の高次の揺らぎとは別の鏡で補正される。

1.2

広視野補償光学

本章の冒頭でも述べたように、現在様々な次世代の AO システムが検討、開発されている。そのうちの 1 つとして広視野補償光学 (Wide Field AO, WFAO) がある。ここではこの WFAO について説明する。

1.2.1

従来の SCAO の限界

従来の SCAO では 1 つの GS 方向の情報を用いる。そのため実際に補正する天体と GS が離れている場合、推定する大気揺らぎと補正する大気揺らぎが異なるので補償精度が落ちてく。図 1.1でも見られるように特に上層でターゲットの光と GS の光が通過する経路が異なるため、それぞれの光が受ける大気揺らぎの影響も異なる。その結果、補償精度が悪化するのである。このように、

SCAO では補正が適用できるのは GS の周辺だけであり、一度に補正できる視野が非常に狭いという制限がある (isoplanatic angle, A.2)。この制限を打破する新しいシステムが WFAO である。
先に述べたように上層での光の経路の違いが一度に補正できる視野の広さを制限している。そこ
で WFAO では一度に複数のガイド星を用いて複数の方向の大気揺らぎを測定し、その情報から大気揺らぎを高さごとに分けて推定する。大気揺らぎを高さごとに分解して扱うことで、方向による経路の違いを考慮し一度に広い視野を補償することが可能となる。これは CT スキャンなどで用いられるトモグラフィーの手法を計算に取り入れることで達成される。

WFAO には補正可能な視野の広さや補正精度が異なる複数の種類が検討されている。どの種類を採用するかは達成すべきサイエンスを考慮した上で判断される。どの WFAO でも従来の SCAO

2

第 1 章 広視野多天体補償光学

1.2. 広視野補償光学の補正範囲 (K バンドで数十秒角) よりは十分広視野の補償が行える。WFAO は ELT 時代の新しいサイエンスを切り開く装置の 1 つとして期待されており、さまざまな開発、検討が行われている。

1.2.2

超大型望遠鏡と広視野補償光学

WFAO は ELT で有効であると考えられている。ここではその理由について簡単に触れる。まず、広視野 AO を実現するためには複数の GS が必要となる。しかし、GS となるような明るい星が密集している領域は限られている。そのため自然にあるガイド星 (Natural Guide Star, NGS) を用いることは困難である。そこで、基本的にはレーザーで人工的に作るレーザーガイド星 (Laser

Guide Star, LGS) を使用することが検討されている。しかし、LGS を用いた場合は揺らぎの傾き成分を測定する別の星 (Tip/Tilt NGS, TT-NGS) が必要となる。つまり、LGS を用いても結局のNGS が必要となる。しかし、TT-NGS の明るさに対する要求は口径が大きくなるほど緩和される。
そのため現在の 8m 級望遠鏡よりも口径が直径 30m を超えるような超大型望遠鏡の方が暗い星をTT-NGS として選ぶことができるので sky coverage が広くなる。さらに口径が大きいほど１つのGS で一度に多くの大気揺らぎの情報が得られる。そのため特に上空での各ガイド星の光の経路の重なりが増え、広い視野の大気揺らぎを十分に測定できる。図 1.2 は LGS を 6 個配置したときに、これらの LGS でどのくらいの領域の大気揺らぎを測定できるかを示した図である。横軸は視野中心からの距離、縦軸はその方向の大気揺らぎのうち何%が LGS によって測定されているかを表している。LGS で測定されているどうかは本論文中で用いた大気揺らぎモデルの高度で計算している (3 章)。図中の赤が直径 8m、青が直径 30m の口径の望遠鏡を用いた場合の結果である。明らかに青の 30m を用いたほうが視野中心から離れた領域の大気揺らぎまで LGS によって測定されていることがわかる。
Coverage Ratio ( 6GS, 240")

Coverage Ratio [%]

100
80
60
40
20

8m
30m

0

0

50

100

150

Distance [arcsec]

200

250

300

図 1.2: 6 個の LGS を視野中心に１つ、その周りに 240”離して 5 つ配置したときに、各方向の大気揺らぎのうち何%が LGS によって測定されているかを示した図である。横軸は視野中心からの角度
距離である。赤が口径が直径 8m の場合、青が直径 30m の場合の結果である。

1.2.3

広視野補償光学の種類

さらに WFAO には大きくわけて３つの種類が検討されている。これらの違いは主に補正方法の違いである。広視野の大気揺らぎを複数の GS の情報からトモグラフィーの手法を用いて推定する

3

第 1 章 広視野多天体補償光学

1.2. 広視野補償光学点はすべて同じである。それぞれについて簡単にまとめる。
多層共役補償光学 (Multi Conjugate Adaptive Optics, MCAO)
複数の DM で高さごとに補正を行う。各 DM は補正したい高さに共役な位置に設置する。推定した視野全体に対して高精度な補正をすることができるので、広がった天体に対しての補償が可能となる。しかし、推定した広い視野の大気揺らぎを完璧に補正するためには、その広い視野に対応した素子数が非常に多い DM が複数必要となる。これは技術的に大きな制限となるため、他の二つに比べると補正視野は狭くなる。TMT の初期補償光学として搭載される NFIRAOS はこの MCAO のシステムを採用している。
地表層補償光学 (Ground Layer Adaptive Optics, MCAO)
地表に近い層は理想的にはすべての方向で共通な揺らぎの成分である。GLAO ではこの地表層のみを補正する。上層の大気揺らぎは補正しないので MCAO や SCAO に比べて補正精度は落ちる。しかし、全方向で共通な地表しか補正しないので 1 枚の DM のみで広い視野に対応できる。さらに観測から地表層ゆらぎのパワーは大気揺らぎ全体の半分以上を占めることがわかっており (Els et al. [9])、支配的な地表層の成分を取り除くだけでも AO 無しの状態に比べて大幅な改善が期待される。
多天体補償光学 (Multi Object Adaptive Optics, MOAO)

MOAO で視野内の複数の天体方向に対して補正を行う。つまり上の２つの WFAO とは違い、補正を行うのは推定した視野全体ではなく目標天体の方向のみである。このように天体の方向だけ選択的に補正することのメリットは、まず各天体に対して DM が１つで高精度な補正が行えることである。また、MCAO のように広い視野ではなく天体の方向だけの補正なので大規模な DM は必要ない。大規模な DM を用意する技術的な制限が緩和されるため、

MCAO よりも広視野から補正する天体を選ぶことができる。その代わりに補正する天体の数だけ DM を用意しなければならないこと、測定方向と観測方向が異なるので Open Loopでの制御が必要となる。
図 1.3∼ 図 1.5 は各 WFAO システムの模式図である。

4

第 1 章 広視野多天体補償光学

1.2. 広視野補償光学

図 1.3: MCAO。複数の DM で視野全体を補償する。

図 1.4: GLAO。１つの DM で地表層のみ補償する。

5

第 1 章 広視野多天体補償光学

1.3. 多天体補償光学

図 1.5: MOAO。天体ごとに補償する。Open Loop 制御である。

1.3

多天体補償光学

ここでは WFAO の中でも本研究で検討している MOAO について説明していく。

1.3.1

多天体補償光学実現の意義

遠方銀河の性質、構造、進化などを議論するためには、暗い銀河まで観測するための集光力、より細かい構造まで調べるための高分解能、さらに統計的な議論を行うためのサンプル数が重要となってくる。ELT と AO を用いれば集光力、高分解能は達成される。しかし、国際協力のもとで成り立っている ELT では各研究者が観測できる時間は非常に限られている。一度に補正できる視野が狭い従来の SCAO ではサンプル数を稼ぐのに膨大な時間がかかってしまう。これらすべての要求を満たすことができるのが ELT と MOAO の組み合わせである。
では、前章で説明した他の WFAO ではなくなぜ MOAO なのか。まず、1 つに補正する天体を選ぶことのできる視野の広さである。MCAO では視野全体を補正するための技術的な制限により、

MOAO よりも狭い視野の補償しかできない。次に推定精度である。GLAO は非常に広い視野を確保できるが、大気揺らぎのうち低層のゆらぎのみの補償なので精度として不十分である。一方、

MCAO は低層のゆらぎも上層のゆらぎも補正するが、視野全体にたいして補償を行うため、各天体にしぼって補正を行う MOAO と比べて高精度を達成するのが困難である。また、遠方に行くほど銀河は小さくなるので各天体に対する補償範囲は狭くてもよく、重要なのは多天体性となってくる。このように遠方銀河を観測するにあたって MOAO は非常に強力なツールであるといえる。

6

第 1 章 広視野多天体補償光学

1.3.2

1.3. 多天体補償光学

サイエンスケース

MOAO で期待されるサイエンスとして、ここでは Thirty Meter Telescope Detailed Science
Case([27],2007) でまとめられているものを簡単に紹介する。TMT では現在第２期装置の１つとして MOAO と多天体面分光装置 (Near-InfraRed, Multi-Object Spectrometer, IRMOS) の組み合わせを検討している、面分光装置 (Integral Field Unit, IFU) は空間情報と波長情報を同時に取得でき、銀河の空間的な構造を調べるのに非常に適した装置である。さらに AO を用いることでより暗い、細かい構造を観測することができるようになる。MOAO+IRMOS はこの IFU 観測を多天体化したものであり、大幅な観測効率の向上が期待できる装置である。具体的なサイエンスケースとしては以下の２つがまとめられている。

z=1∼5 における形成中の銀河の空間構造の調査
宇宙の星形成史の研究から、銀河の成長は赤方偏移が 1 から 3 の時代に最も活発であったことが知られている。また、この時代に銀河中心のブラックホールの活動生がピークを迎えることも知られており、この時代の銀河の力学的構造や形状を詳細に調べること銀河形成を理解する上で非常に重要であるといえる。TMT と IRMOS を用いることでこれまでの数 kpc
スケールの観測から、数 100pc スケールまで分解能が向上し、さらに多天体観測により銀河の空間的な情報を効率よく得ることができる。この情報から銀河形成を引き起こす要因として Major Merger が支配的なのか、それとも別のガス降着過程によるものかという謎に対してさらに詳細な議論をすることが可能となる。また、今までは分解能や観測の深さの関係で議論が困難であった Minor Merger の影響についても迫ることができると期待される。

z=7 を超える宇宙での初代の銀河の調査、宇宙再電離の解明
これまで観測することが困難であった z=7 を超えるような宇宙での初代銀河の調査も TMT+IRMOSで期待されるサイエンスの１つである。この時代は宇宙再電離が起こったとされる時代でもある。宇宙再電離を引き起こす原因の１つとして銀河内の星形成による紫外線放射が考えられており、宇宙再電離のメカニズムを解明するためにもこの時代の銀河の統計的な調査が必要である。この時代の銀河は現在の宇宙の銀河に比べて非常にコンパクトであることが示唆されているため、z が 7 を超えるような宇宙に存在する銀河に対しては、「面分光観測」というよりは「分光観測」という役割が強いかもしれない。

1.3.3

要求される性能2009 年にまとめられた TMT の Science-Based Requirements Document([27]) で述べられている AO に関する要求性能は以下の表の通りである。
これらの要求が妥当であるかを簡単に確認する。まず IFU の視野について考える。図 1.6 はある
スケールが各時代でどのくらいの角度サイズに対応するかを表した図である。赤線は 10kpc に対応するの見かけの大きさを表してる。また、青線は Ono et al. [17] で得られた銀河のサイズ進化を角度で表したものである。Ono et al. [17] では HST のデータから z が 2-8 までの銀河の有効半径と赤方偏移の関係を議論しており、今回はそこから得られた各時代の銀河の典型的な有効半径の２倍に対応する角度サイズを図 1.6 にプロットしている。この図から現在考えられている IFU の視
野 (約 3”) で z が 10 までの 10kpc のスケールの天体はカバーできることがわかる。また、銀河の
サイズは宇宙初期になるにつれて小さくなっていく傾向にあるので 10kpc よりも小さくなる。つまり、z>1 の典型的な銀河は視野 3”ほどで観測できることがわかる。

7

第 1 章 広視野多天体補償光学

1.3. 多天体補償光学

NFIRAOS (MCAO)
波長域0.6-2.5µm

補正視野2.3’ 以上

補正精度 (波面残差 RMS)

133nm 以下 (視野 30”以内)
190nm 以下 (視野 2.3’ 以内)
160-250nm (視野 2.3’ より外側)

補正精度 (160mas スリット内に入る Energy)

視野全体で 30% (J-band)
視野全体で 50% (K-band)

表 1.1: NFIRAOS に関するサイエンスからの要求

IRMOS(MOAO)
波長域

0.6-2.5µm
5’ 以上 (天体を選ぶことができる視野)

推定視野

1”-5” (実際に補正をする視野、IFU の視野)補正できる天体数

一度に 10 から 20 天体補正精度 (0.05”の正方形内に入る Energy)

50% (J-band)

表 1.2: IRMOS(MOAO) に関するサイエンスからの要求次に推定する視野、つまり天体を選ぶことのできる視野の広さについて考える。多天体性を生かすためには大気揺らぎを推定する視野内に十分な数の対象天体があることが求められる。ここでは特に観測対象が少ないと予想される z > 7 の天体の数を見ていく。まず IRMOS の IFU で達成されるであろう限界等級を求め、その値と観測的に得られている光度関数 (Luminosity Function,LF) を比べてある視野内で期待される天体の数を見積もる。今回は Bouwens et al. [4] で得られたz∼7 と z∼8 の Liman Break 銀河 (LBG) の LF を使用する。図 1.7 の赤線が z∼ 7(上)、z∼8(下)の LF を明るいものから積分した結果をプロットしたものである。限界等級の計算では積分時間を10 時間、スリットに対応するサイズは IFU の 1 ピクセルのサイズである 0.05”×0.05”の正方形とした。図中の黒の折れ曲がり線は、ある視野内に要求した数の天体が見つかるかどうかを示している。つまり黒のプロットが赤線の LF よりも右下にあれば、視野内に要求した数の天体が見つかることが期待される。視野は直径 10’(約 78arcmin2 ) と直径 5’(約 19arcmin2 ) の２パターン、要求する天体数は 10 個を考えている。それぞれの黒のプロットがどの視野、天体数に対応するかは図中に載せてある。 今回は観測する天体が点光源の場合に視野内に 10 個以上見つかるかどうかを考える。表 1.2 に書かれている 0.05”×0.05”の Ensquared Energy が 50%という精度目標が達成されるものと考えて計算すると図 1.7 の Point と書かれたプロットの結果になる。点光源の場合はSN=10、R=2000 で計算を行った。図より 10-20 個の複数天体の同時観測というた多天体性を生かすためには、直径 5’(19arcmin2 ) の視野では狭いことがわかる。そこで本研究では直径 10’ の視野で高精度な推定を行うことを目標とする。

8

第 1 章 広視野多天体補償光学

1.3. 多天体補償光学

3
10kpc

Angular Size [arcsec]

2.5

2Re from Ono et al (2012)

2
1.5
1
0.5
0

1

2

3

4

5

z

6

7

8

9

10

図 1.6: 各時代での見かけの大きさを表した図。赤線は 10kpc、青線は Ono et al. [17] で求められた各時代の銀河の有効半径の 2 倍に対応する角度サイズをプロットしたもの。

9

第 1 章 広視野多天体補償光学

-19

-18.5

1 Kimlative LF From Bouwens et al 2011

-2.5

0.5
0

Point FoV=19arcmin2
N=10 SN=10 R=2000

-0.5

2

Point FoV=78arcmin
N=10 SN=10 R=2000

-1

-3

-3.5
-4
-4.5

-1.5

-5

texp=10h D=30m-5.5

z = 6.4 - 7.2
25

25.5

26

-22

-21.5

-21

26.5

27

27.5

28

28.5

MUV [mag]-20.5 -20 -19.5

-19

-18.5

mAB [mag]

1 Kimlative LF From Bouwens et al 2011

-18
-2
-2.5

0.5
0

Point FoV=19arcmin2
N=10 SN=10 R=2000-0.5

Point FoV=78arcmin2
N=10 SN=10 R=2000

-1

-6

-3

3

-2

-2.5

log10 (Number / sq.arcmin / z=7.4-8.6)

-18
-2

3

-21

-3.5
-4
-4.5

-1.5

-5

texp=10h D=30m

-2

-5.5

z = 7.4 - 8.6

-2.5
25

25.5

26

26.5

27

mAB [mag]

27.5

log10 (Number / Mpc )-21.5

MUV [mag]-20.5 -20 -19.5

28

log10 (Number / Mpc )log10 (Number / sq.arcmin / z=6.4-7.2)-22

1.3. 多天体補償光学

28.5

-6

図 1.7: 光度関数と視野内の天体数の関係を表した図。上の図が z∼7、下の図が z∼8 の結果である。赤線は光度関数を明るい方から積分したもの。

10

第 1 章 広視野多天体補償光学

1.3.4

1.3. 多天体補償光学

多天体補償光学の課題

MOAO を実現するにあたりいくつか大きな課題が存在する。まず推定の計算速度である。AO
はリアルタイムで補正を行うため大気の変動よりも十分早い速度で測定、推定、補正のループを回す必要がある。当然、そのループに含まれている推定計算も高速で行わなければいけない。しかし、これまでの SCAO とは違い複数の GS の情報を用いて、さらに広視野の大気揺らぎを推定するため莫大な計算量が必要となる。例えば、現在 Keck 望遠鏡で使用されているシャックハルトマン
WFS(Shack Hartmann WFS, SH-WFS) のサブアパーチャーのサイズは 0.5m であるので、口径が10m の Keck 望遠鏡では 1WFS につき 400 個の subaperture が装備されている。これを口径 30mの TMT に適用すると 3600 個の subaperture をもつ SH-WFS が必要となる。つまり、WFS １つにつき TMT は Keck 望遠鏡の 9 倍の情報量を扱う計算になる。また、WFAO では複数の WFS を使用するため、WFS の数だけその差は大きくなる。さらに推定する大気揺らぎの領域も広がるので、そこでも計算量が増大する。これを打破するためには計算するコンピュータの進歩はもちろんのこと、高速な新しい推定の手法を考える必要がある。
次に推定精度である。MOAO では測定方向と補償方向が異なるため、Open Loop での制御が必要となる。Open Loop 制御では Closed Loop 制御とは異なり、１回の測定や推定や補償の精度が最終的な精度に大きく影響する。そのため、WFS や DM の高精度化はもちろんのこと、高精度な推定計算が要求される。世界的には MCAO のための推定計算の研究がよく行われている (Ellerbroek[8], Gilles et al. [14])。MCAO は Closed Loop 制御なので、計算アルゴリズムを簡易化して精度よりも速度を向上させる手法が最近の話題となっている。しかし、MOAO では計算精度も重要になってくるため、計算速度と精度を両方達成する計算手法が必要である。
最後の課題として、WFS で測定すべき位相差の範囲、DM で補正すべき位相差の範囲をが大きくなることである。口径が大きくなるほど大きな位相の乱れを測定・補償しなければいけないので、それに対応できる WFS、DM が必要である。

1.3.5

多天体補償光学の現状

世界での MOAO の開発の現状を簡単に紹介する。MOAO はヨーロッパ南天天文台 (European
Southern Observatory, ESO) で検討されている Very Large Telescope(VLT) の次世代多天体面分
光装置 FALCON(Fiber optics spectrograph with Adaptive optics on Large ﬁelds to Correct at

Optical and Near-infrared) のプロジェクトの中で初めて提案されたと言われている (Francois Assemat [10], Hammer et al. [23])。さらに現在では前述した TMT の IRMOS、ESO で計画されている直径 42m の口径の超大型望遠鏡 E-ELT の多天体面分光装置 EAGLE(Cuby et al. [6]) が MOAOを用いた装置として検討されている。

MOAO は実際の装置としてはまだ実現していないが、ESO で計画されている EAGLE プロジェ
クトの一環として MOAO の試験装置 CANARY(Gendron et al. [11]) が 2010 年からスペインのCanary Islands にある直径 4.2m の口径を持つ Wiiliam Hershel Telescope で試験されている。この CANARY は Open Loop 制御とトモグラフィック波面再構成の実現可能性を試験する目的で作成され、比較のために MOAO の他にも GLAO と SCAO の機能も兼ね備えている。2011 年に３つの NGS を用いた場合での試験ではストレル比 0.2 を達成しており、誤差の見積もりも行われている。今後は LGS を用いての試験が行われる予定である。

CANARY は比較的口径が小さい望遠鏡で試験されているので、主に技術的な試験にとどまる。
さらにサイエンスを含めた MOAO の試験装置として、ビクトリア大学を中心に試験装置 RAVEN

11

第 1 章 広視野多天体補償光学

1.4. 本論文で検討した新しい手法が開発されている。RAVEN は 2014 年にすばる望遠鏡での on-sky 試験が予定されている。RAVEN
では 3 つの NGS と１つの LGS を用いてトモグラフィック波面再構成を行う。
このように複数の MOAO の試験装置の計画が進んでいるが、最終的な 30m 級の望遠鏡で運営する装置の実現のためには上記の課題を克服するためのさらなる開発、検討の必要がある。

1.4

本論文で検討した新しい手法

上記の課題を克服するために本研究では以下の二つの手法について検討を行った。

1.4.1

GLAO+MOAOWFS の測定範囲や DM の補正範囲への要求を緩和するために図 1.8 のように、GLAO 的に Closed
Loop 制御で地表層を先に補正する DM を組み込む。地表層は全方向で共通の揺らぎなので、1 枚の DM ですべての方向に対して補正が行われる。TMT のサイト調査の結果 (Els et al. [9]) より、
マウナケア山頂の大気揺らぎは地表層の寄与が全体の 50%以上であることがわかっている。この
地表層成分を最初に補正することで、その後の WFS や各方向用の DM で測定、補償する位相のゆらぎを小さくできる。さらに各方向の DM を取り外す、もしくは表面をフラットにすることで

GLAO としての観測も行うことも可能かもしれない。その場合、GLAO を用いたサーベイ観測と
MOAO を用いた高分解能な観測を１つの装置でカバーすることができるかもしれない。本研究では地表層を求める方法として、すべての GS の測定値の平均値を用いる手法を検討した。地表層は共通成分であるので、理想的にはすべての方向の平均値として扱えるはずである。このようにトモ
グラフィーの計算とは独立に地表層を求め、トモグラフィー計算で求める大気揺らぎの層を減らすことで推定結果の向上につながるかどうかを検証する。また、平均値から求めた推定値を用いた場
合で GLAO 的な補正を行った場合に、どれだけの性能が期待されるかについても検討を行った。

12

第 1 章 広視野多天体補償光学

1.4. 本論文で検討した新しい手法

図 1.8: GLAO+MOAO。全方向で共通成分である地表層を 1 枚の DM で Closed Loop 的に補正する。

1.4.2

風速の情報を用いる手法

次に風速の情報を取り入れた手法を説明する。MOAO の推定精度を決める要因の 1 つとして計算に用いる情報量が大きく影響している。この情報量を増やす最も単純な方法は、使用する GS の数を増やすことである。参照する GS の数が多いほど、様々な方向の大気揺らぎの情報が得られる。そのため、推定に使える情報量が笛て推定精度も向上する。しかし、一度に打ち上げられる

LGS の数は限られている。それに伴い推定精度も制限されてしまう。そこで、本研究では風速の情報を取り入れて推定に用いる情報を増やし、推定精度を向上させる新しい方法を検討した。
多くの観測から大気揺らぎは複数の層として考えられることがわかっている (Azouit and Vernin[1], Hardy [24])。さらにいくつかの研究ではある程度の時間内なら大気揺らぎは時間変化しない”frozen ﬂow”とみなすことができ、この時間内であれば大気の時間変化は大気揺らぎの層の風速による移動だけで表すことができることが示されている (Schock and Spillar [20])。このような事
実から少し前の測定の時に GS の光が通過した大気中の経路の時間変化を風速の情報から推測し、少し前の測定値を次の計算に含めることで計算に用いる情報を増やし、推定の精度を向上させることができると考え、検証を行った。

13

第 2 章 波面再構成

第2章

波面再構成

本章では本研究の核である大気揺らぎを推定する波面再構成について説明する。

2.1

仮定

まず推定を定式化するために以下のことを仮定する。

• 大気揺らぎの影響は「光が大気中のある点を通過した時に光の位相がどれだけ変化するか」で表す。

• 大気揺らぎは複数の高さの無限に薄い層として考える。これは観測的に示されている。
• 各層の大気揺らぎは連続的である。
• 大気ゆらぎの時間変化は各層が風によって平行移動することのみで表現する。大気揺らぎ自身の時間変化はここでは考えない。この過程を frozen ﬂow と呼ぶ。
• 単純化するために天体の光は大気中で回折など起こさずに位相のみが変化すると考える。よって、シンチレーションの効果もないと考える。
• 大気の屈折率の波長依存性は非常に小さく、波面の形状はすべての波長で共通であるとする。
(Hardy [24])2.2

大気揺らぎと WFS で測定する波面の関係

大気揺らぎが WFS の測定値にどのような影響を与えるかを考えていくために、図 2.1 のような状況を設定する。このときの座標は望遠鏡の口径の中心を原点とする。この図では大気揺らぎが 2層、WFS が 2 個にしてあるが、一般的に考えて Nlayer 枚の大気揺らぎの層と Nwf s 個の WFS がある状況を考える。
[i]i 番目の大気揺らぎの層の位相点の数を Nφ 、高さを h[i] とする。このとき i 番目の大気揺らぎ[i][i] [i][i]の位相点は φk = φ[i] (xk , yk , h[i] ) {k = 1, . . . , Nφ } と表す。
j 番目の WFS の波面の位相点の数を Nw[j] 、対応するガイド星の方向を (θx[j] , θy[j] ) とする。
WFS は口径と conjugate させているので、口径上にあると考えてよい。WFS の k 番目の位相点を w[j]k = w[j] (xk[j] , yk[j] , 0) {k = 1, . . . , Nw[j] } とする。また、ガイド星が NGS である場合、この WFS の i 番目の大気揺らぎの層上での位相点は[i]wk[j] = w[j] (xk[j] + h[i] θx[j] , yk[j] + h[i] θy[j] , h[i] ) {k = 1, . . . , Nw[j] }

14

(2.1)

第 2 章 波面再構成

2.2. 大気揺らぎと WFS で測定する波面の関係

図 2.1: ガイド星 (または WFS)[赤丸] の位相点と大気揺らぎの位相点 [黒丸] の関係。各高さで必ずしも一致するわけではないので大気揺らぎの情報を内挿しなければいけない。

ここでガイド星は視野中心から非常に近い (θx[j] , θy[j]  1) なので、tan(θ) ∼ θ の近似を用いた。
一方、LGS の場合は高度が高くなるにつれて WFS の位相点の間隔が狭くなっていくので、[i]wk[j] = w[j] (xk[j]Hlgs − h[i]Hlgs − h[i]+ h[i] θx[j] , yk[j]+ h[i] θy[j] , h[i] ) {k = 1, . . . , Nw[j] } (2.2)
Hlgs
Hlgs

となる。この後は特に LGS、NGS を区別せず
[i]

[i]

[i]

wk[j] = w[j] (xk[j] , yk[j] , h[i] ) {k = 1, . . . , Nw[j] }(2.3)

として進める。
[i]

図 2.1 で見られるように、ガイド星の受ける影響 wk[j] を求めるためには大気揺らぎの情報を内挿しなければいけない。今回は線形補完で内挿を行う。線形補完では補完したい点の周囲の 4 点を用いて、各点から補完したい点までの距離で重み付けをして計算する。大気の揺らぎ層の位相点の間隔を d0 、補完に使用する大気揺らぎの４点の位相点を k1 ∼ k4、補完したいガイド星の位相点を k としたとき、次のように補完を計算する。

 

[i]
[i]
[i]
[i]
xt − xk[j]yt − yk[j]
∑
[i]
1 −
 × 1 − φ[i]
wk[j] =
t
d0
d0

(2.4)

t=k1,k2,k3,k4

これをガイド星のすべての位相点に対して計算を行う。このとき以下のように行列の形式にまとめることができる。








[i]
φ1


 [i] [i]
 w2[j] φ2 



[i] 
=
A


.. 
..
[j] 




.


 . 
[i]
[i]
φNφ[i]
wN [j]
[i]

w1[j]

w[j]

[i]

[i]

w[j] = A[j] φ[i]
[i]

(2.5)

φ[i] は i 番目の大気揺らぎのベクトル、w[j] は j 番目のガイド星の i 番目の大気揺らぎを通過した[i]

後の波面ベクトル、A[j] はそれらの関係をつなぐ補完行列である。

15

第 2 章 波面再構成

2.3. WFS の測定値と波面の関係

最終的に WFS で測定される波面 w[j] はすべての大気揺らぎの影響の合計であるから、∑

Nlayerw[j] =∑

Nlayer[i]w[j] =[i]A[j] φ[i]

(2.6)

i=1

i=1

となる。

2.3

WFS の測定値と波面の関係

次に WFS の測定値と波面の関係を考える。補償光学では測定値と波面の離散的な位相点のコン
フィギュレーションとして主に図 2.2 の３つが考えられている。本研究で用いる Shack Haltmann

図 2.2: 補償光学で用いられる WFS の測定値と波面の位相点のコンフィギュレーション。×印が波面の位相点、矢印が WFS の測定値を表している。

WFS(SH-WFS) を用いた測定方法では Fried Geometry との対応がよい。
次に Fied Geometry の場合に各波面の位相点と各測定値の関係を考えていく。Fried Geometryでは一つの測定値に対して、その周りの 4 つの位相点が関係する。実際に図 2.3 のような 2×2 のsubaperture を持った SH-WFS の場合で説明をしていく。 このとき、波面の位相点は 3 × 3 = 9

図 2.3: 2×2 の subaperture を持つ SH-WFS。各 subaperture は quad cell。

個、測定値は (2 × 2) × 2(x と y)= 8 個である。このとき各測定値は波面の位相点の値を用いて次16

第 2 章 波面再構成

2.4. 問題の定式化のように表される。

)
1 ( [1]sy = −w11 − w12 + w21 + w22
2)
1 ( [2]sy = −w12 − w13 + w22 + w23
2)
1 ( [3]sy = −w21 − w22 + w31 + w32
2)
1 ( [4]sy = −w22 − w23 + w32 + w33
2

1(−w11 + w12 − w21 + w22 )2
1(−w12 + w13 − w22 + w23 )
s[2]
x =
2
1
(−w21 + w22 − w31 + w32 )
s[3]
x =
2
1
s[4]
(−w22 + w23 − w32 + w33 )
x =
2
s[1]
x =

(2.7)

この連立方程式を行列とベクトルを用いて書くと以下のようになる。

 [1] sx
−1 1 [1] sy 
−1 −1 s[2]  0 −1 x 
 [2] sy  1  0 −1
 = 
s[3]  2  0
0
 x 

 [3] 

sy 
0
0
 

s[4] 
0
0
 x 

[4]
sy
0
0

0

−1

1

0

0
1
−1

1
0
0

1
−1
1

0
1
1

0
0

−1 1
−1 −1

0
0

0
0

0
0

−1 1
−1 −1


w11
0
0 0 

 w12 
00 0 
w13 



0
0 0  w21 
0
0 0 w22 
−1 1 0  w23 
1
1 0 
w31 
0 −1 1 w32 
0
1 1
w33
s = Gw



(2.8)

(2.9)

ここで w は波面の位相点の値のベクトル、s は測定値のベクトル、G はその二つを関係を表す行列である。なお、本論文では数式中の太字の大文字は行列を、太字の小文字はベクトルを表しているとする。測定値の数を ns × ns × 2 = Ns 、波面の位相点の数を nw × nw = Nw とすると、Fried Geometry では ns = nw − 1 なので Ns = (nw − 1)2 × 2 という関係が成り立つ。これよりNw < Ns (nw > 4) となり、一般的に式 (2.9) は求める解よりも測定値 (または方程式) の数が多い
過決定方程式であることがわかる。しかし、ピストン成分やワッフルのようにギザギザした高周波数成分など、Freid Geometry では測定できないモードが存在する。これが計算時に悪影響を及ぼすことがあるため注意する必要がある (Poyneer et al. [18])。G は Ns × Nw の行列である。測定ノ
イズを区別する場合は

s = Gw + η(2.10)のように定式化する。η は測定ノイズのベクトルである。ガイド星と WFS が一つの場合の従来の補償光学の場合はこれらの式を解くことで再構成を行う。

2.4問題の定式化式 (2.6) と式 (2.9) を用いて、WFS の測定値から大気揺らぎを求めるために解くべき関係式を構築していく。j 番目の WFS で測定される測定値ベクトルを s[j] とすると、この２式より i 番目の大気揺らぎと j 番目の WFS の測定値の関係は∑Nlayers[j] = G[j] w[j] = G[j]

i=1

17

[i]

A[j] φ[i](2.11)

第 2 章 波面再構成

2.5. 最小分散推定 (MINIMUM VARIANCE RECONSTRUCTION)と書ける。すべてのガイド星の情報から大気揺らぎを求めるために、上の式をすべてのガイド星で連立させる。




G[1]



∑Nlayer[i]A[1] φ[i]





s[1]

∑Nlayer [i] [i] G[2] i=1 A[2] φ s[2]   .  = .  .
 .  
..





s[NGS ]
∑Nlayer [i]G[NGS ] i=1 A[NGS ] φ[i]


i=1

[1]

G[1] A[1] G A[1][2] [2]


=
..


.


[1]
G[NGS ] A[NGS ]

[2]

[N

G[1] A[1]layer[2]G[2] A[2]

···

[N
]
G[2] A[2]layer

..
.

..

.

..
.

···

G[NGS ] A[Nlayer
GS ]

[2]

G[NGS ] A[NGS ]

[N



]

···

G[1] A[1]

]



φ[1]


  φ[2] 




..



.


[N
]
layer
 φ

ここで [ ] の表記は行列の拡張を表している。この式を次のように簡単に書く。

s = G0 φ(2.12)

これで大気ゆらぎと WFS 測定値の関係式が導けた。この連立方程式を φ について解けば再構成が行える。しかし、ここで注意しなければいけないことがある。図 2.1 でも見られるようにどのガイド星でも測定されない大気揺らぎの領域が存在する。特に上層のガイド星よりも外側では大気揺らぎの情報を得ることはできない。式 (2.12) を解く場合に、このような領域を扱うための工夫が必要となる。また、測定ノイズを別に扱う場合は

s = G0 φ + η(2.13)となる。

2.5

最小分散推定 (Minimum Variance Reconstruction)

式 (2.12) や式 (2.13) で測定できないモードや領域を統計的に制限して最適な解を求める手法として Ellerbroek [8] で最小分散推定が提案された。最小分散推定は推定値と真の値の残差の期待値が最小になるような推定値を求めていく手法であり、解やノイズの統計的性質を計算に含めることができる。本研究ではこの方法を採用する。
最小分散法推定の推定行列を R、最小分散推定による大気揺らぎの推定値を φmv とする。このとき

φmv = Rs

(2.14)と書くことができる。最小分散推定の定義より、R は次の条件をみたすような行列となる。

min ||φmv − φ||2R

18

(2.15)

第 2 章 波面再構成

2.5. 最小分散推定 (MINIMUM VARIANCE RECONSTRUCTION)

ここで h i は期待値を表している。では、具体的な R を求めていく。

||φmv − φ||2 = trace= trace= trace{{{(φmv − φ)(φmv − φ)T}

φmv φTmv − φmv φT − φφTmv + φφT}

φmv φTmv − φmv φT − φφTmv + φφT}{}= trace R ssT RT − R sφT − φsT RT + φφT
R が式 (2.15) を満たす場合は ∂( ||φmv − φ||2 )/∂R = 0 が成り立つので
0=

)
)
)
∂ (
∂ (
∂ (
R ssT RT −R sφT −
φsT RT∂R∂R∂R= 2R ssT − 2 φsT
上の式の変形は Appenix B を参照。よって(

R = φsT

ssT)−1(2.16)

ここで式 (2.13) の式より、

ssT = (G0 φ + η)(G0 φ + η)T= G0 φφT G0T + G0 φη T + ηφT G0T + ηη T
φsT = φ(G0 φ + η)T= φφT G0T + φη Tとなる。大気揺らぎのピストン成分は再構成に特に影響しないので φ の平均値を 0、つまり hφi = 0とみなすことができる。さらに大気揺らぎの共分散行列を φφT = Cφφ とする。また、測定ノイズを平均が 0 で分散 σ 2 のガウシアンノイズと仮定すると hηi = 0、 ηη T = σ 2 I となる。ここでI は単位行列を表す。さらに大気揺らぎと測定ノイズは独立であるとする。これらを考慮すると、
ssT = G0 Cφφ G0T + σ 2 I
φsT = Cφφ G0Tとなり、式 (2.16) に代入すると()−1
R = Cφφ G0T G0 Cφφ G0T + σ 2 I()−1−1= G0T G0 + σ 2 Cφφ
G0T(2.17)−1となる。1 段目から 2 段目の式変形は 1 段目と 2 段目の式に左から G0T G0 + σ 2 Cφφ、右から

G0 Cφφ G0T + σ 2 I を作用させたときに等しくなることから求めることができる。これで、)−1(−1
φmv = G0T G0 + σ 2 Cφφ
G0T s
よりトモグラフィック波面再構成が計算できる。

19(2.18)

第 2 章 波面再構成

2.6

2.6. 大気ゆらぎの共分散行列

大気ゆらぎの共分散行列

ここでは大気の共分散行列 Cφφ について考える。まず各層の大気揺らぎは独立なので、この共
分散行列は

 [1]
Cφφ

Cφφ = 

0



0
..





.

(2.19)[N]
Cφφlayer[i]のようなブロック対角行列である。ここで Cφφ は i 番目の大気揺らぎの層の共分散行列を表す。

Gilles et al. [13]、Gilles et al. [14] より推定する大気揺らぎが周期的であるという境界条件を仮定することで各大気の共分散行列を以下のように計算することが可能になる。
Cφφ = F −1 Λ[i] F[i](2.20)

ここで κ は空間周波数であり、揺らぎのスケールを ` とすると κ = 2π/` と表される。この仮定により低周波の誤差が発生するが、推定する大気揺らぎのサイズが口径に比べて十分大きければ特に問題はないことがわかっている (Gilles et al. [13])。ここで F は離散的フーリエ変換、F −1 = F ∗ /Nは逆フーリエ変換を表している。Λ は各要素が大気の空間的なパワースペクトル (Appendx A) から決まる対角行列である。アウタースケール を考慮した場合、次のように計算される。

(−5/3
Λ[i] (κ) = fi (0.1517)2 r0|κ|2 +

1)−11/6L20(2.21)

fi は各大気揺らぎの強度の割合であり、∑Nlayer0 < fi < 1

fi = 1(2.22)

i=1を満たす。
式 (2.19)、式 (2.20)、式 (2.21) から共分散行列を計算できるが、この共分散行列はほぼすべての要素が 0 以外である密行列となる。行列演算において密行列によるかけ算は膨大な計算量を必要とするため、Ellerbroek [8] でこの共分散行列をほとんどの要素の値が 0 である疎行列の形式に近似的に変形させる手法が提案された。本研究でもこの手法を用いる。

Appendix A で述べているように Kolmogorov の大気モデルを仮定した場合、位相ゆらぎのパ
ワースペクトルは ∝ κ−11/3 となるが、この近似では ≈ κ−12/3 = κ−4 と近似する。これにより共[i]
分散行列 Cφφ は離散的なラプラシアン L[i] を用いて次のように近似することができる。
(

[i]

Cφφ

)−1

≈ c2i L[i]T L[i](2.23)

ここで離散的なラプラシアン L は次のような形になる。

(i = j)
 −1

Lij =




1/ni
0

(i 6= j 、i と j 番目の位相点が隣り合う場合。)(i 6= j 、i と j 番目の位相点が隣り合わない場合)(2.24)

ここで ni は i 番目の位相点と隣合う位相点の数である。普通は ni = 4 であり、端の方では 2 や 3になる。c2i は各大気揺らぎの重みを表す定数である。本研究では c2i = 1/fi とする。Ellerbroek [8]20第 2 章 波面再構成

2.7. 開口の影響で、この近似による誤差は大きくても 1%であると考えられている。最終的な共分散行列は以下のような形になる。

 −1
f1 L[1]T L[1]
LT L = 

0

0
..

.
−1
fN
L[Nlayer ]T L[Nlayer ]
layer






(2.25)

この近似により解くべき式は

(
)−1
σ2 T
0T 0
φmv = G G +
G0T s
L L
fi

(2.26)

となる。

2.7

開口の影響

これまで WFS と測定した波面の関係は正方形のグリッドで考えてきた。しかし、実際は望遠鏡の口径は円、もしくは副鏡の影によるドーナッツ型である。また、副鏡を支えるスパイダーの影も移るのでそれらの影響を考えなければいけない。そのために開口行列 M を定義する。
{
1 (測定値が範囲内の場合
M=
0 (測定値が範囲外の場合(2.27)

開口の影響を考える場合は M を用いて、式 (2.9) は

M s = M Gw

(2.28)となる。この式を用いて計算すると、最終的に最小分散推定で解くべき式式 (2.26) は

φmv =

=

(
)−1
σ2 T
G0T M 0T M 0 G0 +L L
G0T M 0T M 0 s
fi
(
)−1
σ2 T
G0T M 0 G +L L
G0T M 0 s
fi

(2.29)となる。ここで M 0 は M を GS の数だけ対角に並べたような行列であり、M 0T = M 0 、M 0 M 0 =M 0 である。本研究では円口径の場合で計算を行っていく。

2.8

風速の情報を用いる手法

本研究で新たに開発した風速の情報を用いる手法について説明していく。
まず図 2.4 のような状況で、ある時間 t1 で各 GS が通過した経路の時間変化について考える。時間変化は風による大気揺らぎの層の移動によるものと考えているので、i 番目の大気揺らぎの風ベ
[i]

[i]

クトルを v [i] = (vx , vy ) とすると、時刻 t1 から t2 の間 ∆t = t2 − t1 に大気の層が移動する距離はv [i] ∆t となる。図 2.4 の左図のように時刻 t1 で GS を測定した場合、この GS が通過する各大気の層上での座標は式 (2.3) で与えられる。時刻 t2 では図 2.4 の右図のように大気の層の移動によってt1 での位置からずれる。このときの座標は
[i]

[i]

[i]

wk[j] (∆t) = w[j] (xk[j] + ∆tvx[i] , yk[j] + ∆tvy[i] , h[i] ) {k = 1, . . . , Nw[j] }21(2.30)

第 2 章 波面再構成

2.8. 風速の情報を用いる手法

図 2.4: 大気揺らぎの層の時間変化による、GS の通過した経路のずれ。

と書くことができる。
式 (2.30) から ∆t 時間経ったときの座標を計算できれば、図 2.4 のように時刻 t2 では t2 での測定値 (青) と ∆t 前の測定値 (赤) を用いて計算することができる。つまり、実質 GS が 2 倍になったと考えられる。このとき ∆t 前での測定値と t2 での大気揺らぎの関係式は式 (2.3) の代わりに式(2.30) を用いて同様に計算することができる。
∑

Nlayerw[j] (−∆t) =∑Nlayer[i]w[j] (−∆t) =

i=1

[i]

A[j] (−∆t)φ[i]i=1s[j] (−∆t) = G[j] w[j] (−∆t)
よってこの ∆t 前の測定値を考慮すると



[1]

···

G[1] A[1] (0)



..

.

 

 

 
[1]

 

  G[NGS ] A[NGS ] (0)

 

 =
 s (−∆t)   [1]  G[1] A[1][1] (−∆t)

 
.

 
.
.

 
..

.

s[NGS ] (−∆t)
[1]
G[NGS ] A[NGS ] (−∆t)
s[1] (0)..
.
s[NGS ] (0)



[

s(0)

..

···
···
..

.

···

]

s(−∆t)

.

[
=

[N

]

G[1] A[1]layer (0)









φ[1]

[Nlayer ]
G[NGS ] A[NGS ] (0)   φ[2] 




.

..

[Nlayer ]G[1] A[1](−∆t) 

[N
]
layer
 φ

..

.


[N
]
G[NGS ] A[Nlayer
(−∆t)
GS ]

G0 (0)..
.

]

G0 (−∆t)φとなる。さらに最小分散推定で解くべき式は式 (2.29) を変形させて[][]M 0 G0 (0)s(0)
00G (0, −∆t) =s (0, −∆t) =M 0 G0 (−∆t)
s(−∆t)

22(2.31)

第 2 章 波面再構成

2.8. 風速の情報を用いる手法

{
}−1
σ2 T
T
T
0
0
φmv = (G (0, −∆t)) G (0, −∆t) +L L(G0 (0, −∆t)) s0 (0, −∆t)fi

(2.32)と書くことができる。また、さらに前の情報も使用したい場合は G0 (0, ∆t, −2∆t) というように拡張していけばいい。

23

第 3 章 シミュレーションの詳細

第3章

シミュレーションの詳細

シミュレーションの設定

3.1

シミュレーションの設定について説明していく。
望遠鏡
望遠鏡の口径は直径 30m とする。望遠鏡の口径は円形とする。副鏡の影は考えない。TMTの主鏡は分割鏡であるが、今回は単純化のため理想的な 1 枚鏡とする。また、光学系によって発生する収差は無いとする。
観測条件
観測地はマウナケア山頂 (4km) とする。また、観測方向は天頂方向のみとする。
推定範囲
直径 10’ の正方形の視野に対して推定を行う。
ガイド星

GS は Na LGS を用いる。高さは 92km であるが、マウナケア山頂 (4km) を仮定しているので計算上は 88km となる。実際は LGS から波面全体に広がる傾き成分の揺らぎは測定できないので、TT-NGS を考える必要がある。しかし、今回は単純化のため LGS から傾き成分が測定できると仮定する。つまり、TT-NGS の存在は考えない。また、LGS の測定波長は589nm とする。GS と GS の距離や GS の数の違いを見るために GS の配置は視野中心に 1
つ、残りをその周りに円上に等間隔に配置する (図 3.1 の左図、中心の図)。さらに IRMOSで検討されている配置についても計算を行っていく (図 3.1 の右図)。

図 3.1: シミュレーションで用いる LGS のコンフィギュレーション。左は LGS が 5 つの場合、中心はLGS が 7 つの場合の例である。右は IRMOS で検討されている配置である。推定精度を評価す
るときは LGS と重なる方向 (青線上)、LGS と LGS の間の方向 (赤線上) で評価を行う。

24

第 3 章 シミュレーションの詳細

3.1. シミュレーションの設定

WFS
WFS は SH-WFS を仮定する。各 subaperture のサイズは Keck 望遠鏡で用いられているサイズと同じ 0.5m とする。また今回は WFS による測定と DM による補正の間に時間差はな
いものとする。

DM
実際は DM の各アクチュエーターの影響関数などを考慮しなければいけないが、今回はまず推定のみの精度を検討するため考慮しない。よって DM は推定値を完璧に再現する理想的なものを考える。
大気揺らぎ
大気揺らぎのモデルは TMT のサイト調査 (Els et al. [9]) で得られたマウナケア山頂のモデルをもとに、ドームで発生するシーイングを考慮したモデルを用いる。r0 = 0.156m(@500nm)
であり、大気揺らぎの層の高さと強度は表 3.1 の通りである。一番の特徴は 0m の地表層が高度 [m]

0

500

1000

2000

4000

8000

16000

強度の割合

0.59600.09630.03250.03720.08690.06840.0826

表 3.1: 大気モデル

全大気揺らぎの 60%ほどを占めていることである。また、大気揺らぎはカリフォルニア大学の Matthew C. Britton によって開発された補償光学シミュレーションコードである Arroyo
を用いて作成する。用いる大気揺らぎのパワースペクトルは von Karman のパワースペクト
ルとする。また、outer scale の値は 30m とする。inner scale については影響が少ないので考慮しない。図 3.2 は作成したモデル大気ゆらぎの図である。推定する大気揺らぎの位相点の図 3.2: Arroyo で作成したモデル大気揺らぎ。スケールを合わせてあるので、地表層が非常に強いことがわかる。

サンプリングは 0m∼2000m で 0.25m、それ以上で 0.5m という値を用いた。下層では各 GSの光の経路の重なりが大きいのでサンプリングを細かくした。これ以上細かくしても、WFSのサンプリングで制限されるため結果に影響はなかった。計算時間が増えるだけである。
最初に用意するモデル大気のサンプリングは 0.05m とする。これは今回用いる r0 の値よりも十分小さく、0.05m よりも細かいスケールの揺らぎの影響は非常に小さいので無視することができる。また、各計算ごとに大気揺らぎは作り直す。
風速のモデル 一般的に風速は wind shear により、9km∼10km くらいの高さで最も強くなるガウ
シアン的なプロファイルを持っていることが知られている (Appendix A)。Hardy [24] では

25

第 3 章 シミュレーションの詳細

3.2. シミュレーションの流れ

次のようなモデルが提唱されている。
[ (
)2 ]
[
]
z cos(ζ) − HTv(z) = vG + vT exp −× sin2 φ + cos2 φ cos2 ζ
LT

(3.1)

ここで、vG は低い層での風速、vT は対流圏界面付近での風速、ζ は観測している高度角、HT は対流圏界面の高さ、LT は対流圏界面の厚さ、φ は望遠鏡が向いている方向と風向の角度である。今回は vG =5m/s、vT =20m/s、ζ = 0、HT =6km、LT =5km、φ = π/2 とした。
この場合の風速モデルは図 3.3 のようになる。風向きはモデルに従って計算ごとにランダム

Wind Model
25
20

Telescope at 4km

Wind Velocity [m/s]

22.5
17.5
15
12.5
10
7.5
5
2.5
0

0

2.5

5

7.5

10

12.5

Altitude [km]

15

17.5

20

図 3.3: 風速のモデルに与える。
表 3.2 にシミュレーションで用いた各設定、値についてまとめる。

3.2

シミュレーションの流れ

本研究のシミュレーションの流れを図 3.4 に示す。まず、各値を設定する。その後、モデル大気
揺らぎと最小分散推定の行列を作成する。次に各 LGS の光の位相が大気揺らぎを通過することで受ける影響を計算し、WFS で測定される位相波面を求める。このとき、有限の高さにある LGS の
光は大気中をコーン上に通過してくることを考慮する。さらに求めた波面を測定値に変換し、そこから最小分散推定を用いて大気揺らぎの推定を行う。推定した大気揺らぎを最初に準備したモデル
大気揺らぎを比較することで推定精度の評価を行う。推定精度の評価は、WFS で測定できない小さなスケールの揺らぎの影響も考慮するため WFS よりも細かいサンプリング (0.05m) で行う。まず、評価したい方向にターゲット天体があると仮定し、そのターゲットの光がモデル大気揺らぎから受ける影響を計算し、 WFS で測定される波面を求める。次に推定した大気揺らぎを用いて同様に WFS で観測される波面を求める。推定値から求めた波面を補償量とし、大気モデルから求めた波面との差をとる。この残差が補正しきれなかった推定誤差に対応する。この残差波面から RMSと PSF を計算する。本来はこの評価を視野全体で行うのが望ましいが、ここでは計算時間の関係

26

第 3 章 シミュレーションの詳細

3.3. 大気揺らぎモデルの再現性の確認

口径

30m

観測地

Mauna Kea 山頂 4km

観測方向

天頂方向

ガイド星

Na LGSLGS の高度92km

LGS の観測波長

589nm

推定視野

直径 10’WFS

SH-WFSWFS の subaperture のサイズ

0.5m

r0 (@500nm)

0.156mOuter Scale30m

推定する大気揺らぎの位相点のサンプリング

0.5m (4000m 以上) 、0.25m (2000m 以下)

モデル大気の位相点のサンプリング

0.05m

表 3.2: シミュレーションの設定

から 2 方向に対して評価を行う。1 つは視野内で最もいい精度が出ることが期待される GS と重なる方向 (図 3.1 の青線)、もう 1 つは視野内で最も精度が悪くなると期待される GS と GS の間の方
向 (図 3.1 の赤線) で評価する。
さらに大気揺らぎの時間変化を再現するためにこれらの操作を繰り返し行う。その際、最初に与えた風速に従い、各高さの大気揺らぎを動かしていく。今回は 1 ステップを 0.1s のサンプリングとし、それを 50 回繰り返し、5 秒分の計算を行った。最終的に RMS の値は全ステップの平均値、PSF は全ステップの合計を用いる。またその PSF からストレル比、FWHM、Ensquared Eneregyを計算する。

3.3

大気揺らぎモデルの再現性の確認

今回のシミュレーションが実際の大気揺らぎを再現しているかを確認するためにサイト調査で得られているシーイングとシミュレーションで計算したシーイングの値を比較する。シーイングはある程度の時間露出した PSF から計算される。計算を行った波長は 500nm、800nm、1250nm、1650nm、2200nm である。図はそれぞれの 5s 露出の PSF の画像である。計算は各波長で 10 回ずつ行い、その平均値を用いる。シーイングの測定は IRAF の imexam を用いて行う。このとき
moﬀat 関数でフィッティングを行う。図 3.5 がその結果である。横軸は波長で、縦軸はシーイング
サイズである、黒四角がシミュレーションの結果で、緑の５角形がサイト調査で得られた観測値である。また、赤線は Kolmogorov のパワースペクトル仮定して観測値から他の波長でのシーイングを計算した値を示している。Kolmogorov の理論に従えばシーイングは長波長になるほど λ−0.2 で小さくなっていく (Appendix A)。この図からシミュレーションは観測値、また観測値から計算した解析的な値と比べてシーイングを小さく見積もっているように見える。
この結果はアウタースケールの影響で説明できることが Tokovinin [22] でわかっている。Tokovinin[22] では口径が r0 よりもずっと大きいような望遠鏡での観測の場合に、アウタースケール L0 を考慮した von Karman スペクトルから計算したシーイング vK がアウタースケールを考慮しないコ
ロモゴロフのパワースペクトルから求めたシーイング 0 よりも系統的に小さくなることが数値シ

27

第 3 章 シミュレーションの詳細

3.3. 大気揺らぎモデルの再現性の確認

図 3.4: シミュレーションの流れ

ミュレーションで求めらている。さらにこのとき vK と 0 の関係は L0 /r0 > 20 であれば式 (3.2)
で近似的に精度良く表すことができることも求められている。

(

vK
0

)2

(
≈ 1 − 2.183

r0
L0)0.336(3.2)

この式から L0 =30m、r0 =0.156m のときではアウタースケールの影響でシーイングが λ =500nm
で 20%、λ =2200 で 40%も小さくなることがわかる。実際にアウタースケールの影響を修正した
シーイングと波長の関係は図の青線になる。今回のシミュレーションの結果はこの青線によく一致しており、アウタースケールの影響でコロモゴロフのパワースペクトルを仮定して解析的に求めた
シーイング (赤線) よりもシミュレーションの結果が小さくなったと考えることができる。
また、サイト調査で得られるシーイングは主に DIMM などの口径が r0 に比べてそんなに大きくない装置で測定されている。そのためアウタースケールの影響は特に受けないので、今回のシミュ
レーションの結果とは一致していない。実際に David J. E. Floyd and Prieto [7] では DIMM で測定されたシーイングよりも 6.5m の望遠鏡で観測した星のサイズが小さくなるという結果が得られており、アウタースケールを 25m と仮定して式 (3.2) の関係を考慮するとその違いが説明できると結論づけている。よって、今回のシミュレーションの結果が測定値に比べて小さくなったことも現実的であると考えられるので、特に補正は行わない。

28

第 3 章 シミュレーションの詳細

3.4. 計算手法

1
Kolmogorov
Tokovinin(2002)
Simulation
Observation

0.9

FWHM [arcsec]

0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

500

800

1250

1650

Wavelength [nm]2200

図 3.5: 波長とシーイングの関係を示した図。黒四角がシミュレーションから計算した結果。緑の５角形
は TMT のサイト調査で得られた観測値。赤線はコロモゴロフのパワースペクトルを仮定した解析的な値。青線は Tokovinin [22] のアウタースケールを考慮した場合の結果。青線とシミュ
レーションの結果が良く一致している。シミュレーションではシーイングを観測値よりも小さく見積もっているが、それはアウタースケールの影響で説明できることがわかる。

3.4

計算手法

3.4.1

共役勾配法

最小分散の計算は推定行列 R を s に掛けることで行われる。R の要素数は (推定する大気揺らぎの位相点の総数)×(測定値の総数) となる。推定する視野や GS の数にもよるがだいたい推定する大気揺らぎの位相点の総数は 105 以上、測定値の総数は 104 以上となる。つまり、少なくとも約109 の要素数を持つ行列となる。例えばこの行列を double 型で確保する場合、8GB ものメモリが必要となる。さらに R を計算するためには式 (2.29) を見てもわかるように逆行列を計算する必要がある。このような大規模な行列の逆行列の計算は非常に時間がかかるため現実的ではない。
このような大規模問題に対しては、一般的に直接逆行列を求めない方法が用いられる。特に頻繁に採用される手法に共役勾配法 (Conjugate Gradient Methed, CGM) がある。詳細は Appendix

B で述べるが、CGM では Ax = b という問題を x について解く場合に逆行列 A−1 を計算することなく、行列 A とベクトル b のかけ算を繰り返すことで反復的に求めることができる。
さらに、R を構成する行列は要素のほとんどが 0 であるような疎行列である。そこで必要なメ
モリ量を減らすために 0 以外の値のみを確保する形式を用いる。本研究では Compressed Sparse

Row(CSR)1 という形式を用いて各行列を格納する。疎行列とベクトルの演算は高速なライブラリが多く存在する。
また、CGM は係数行列に対してある前処理を行うことで収束回数を減らすことができる。これを前処理付き共役勾配法 (Preconditioned Conjugate Gradient Method, PCGM) と呼ぶ。いくつかの研究では AO の大規模な波面再構成の計算にこの PCGM を用いており、例えば Bardsley [3]では前処理として正則化をかけた最小自乗法のための行列を作用させる Least Square PCGM が用
いられ、Gilles et al. [14] では多重格子法を前処理として用いた Multi Grid PCGM が検討されて1 Compressed Row Storage(CRS) と呼ばれる場合もある

29

第 3 章 シミュレーションの詳細

3.4. 計算手法いる。今回の計算では前処理として最も簡単な対角スケーリングを用いる。対角スケーリングの詳
細は Appendix B で述べる。
実際の計算ではまず式 (2.29) を以下のように変形する。
[][ ]M 0 G0s0Ḡ =s̄ =ḠT Ḡφmv = ḠT s̄
σL
0

(3.3)

まず ḠT s̄ = z を計算し、次に ḠT Ḡφmv = z に対して対角スケーリングを前処理とした PCGMを適用し、φmv を求める。

3.4.2

GPGPU を用いた並列計算

行列とベクトルのかけ算は一般的に並列計算と相性がいい。本研究では計算の高速化を達成する
ために、General-Purpose computing on Graphics Processing Units(GPGPU) を用いた並列計算を取り入れる。GPU は一般的に画像処理を専門とする演算装置であり、単純な計算に対しては強力な威力を発揮する。GPGPU はこの GPU を様々な計算に生かす手法であり、近年急速に発展してきた。特に NVIDIA 社の提供する Compute Uniﬁed Device Architecture(CUDA) は GPU 向けの C 言語の統合開発環境であり、C 言語ライクな言語を採用しているため非常に汎用性の高いものとなっている。また、CUDA 用の行列演算や高速フーリエ変換のライブラリも無料で用意されている。
今回はこの CUDA とそのライブラリを用いて共役勾配法の並列化を行った。用いたライブラリ
は密行列演算ライブラリ CUBLAS と疎行列演算ライブラリ CUSPARSE である。用いた GPU や

CUDA の情報を表 3.3 にまとめる。
GPGPU
種類

Tesla C2070

GPU の数

1

コア数

448

プロセッサ周波数

1.15GHz

倍精度浮動小数点性能（ピーク時）

515Gﬂops

単精度浮動小数点性能（ピーク時）

1.03Tﬂops

専用メモリ合計

6GB

メモリ帯域

144GB/sec

CUDA
Version3.2

表 3.3: シミュレーションの設定

3.4.3

共役勾配法の収束条件

今回の CGM での収束条件としては相対誤差を用いている。k 回目の反復の推定値を φkmv とす
ると、残差 rk は rk = z − ḠT Ḡφkmv で定義される。収束条件は相対誤差 |rk |/|z| がある値 r よりも小さくなった場合、収束したと判断される。検証の結果、今回の計算では r = 10−3 程度で十分であった。詳細は 4 章で述べるが、この場合の計算の反復回数は 250-300 回である。

30

第 4 章 結果

第4章

結果

本章ではシミュレーションの結果を説明する。ここで、それぞれの方法を区別するために普通のトモグラフィーの手法を Normal Method、平均値を用いた手法を Average Method、風の情報を用いた手法を Wind Method と呼ぶことにする。

4.1

最小分散推定によるトモグラフィック波面再構成の結果

まず Normal Method で計算を行い、現状での推定精度、問題点について見ていく。

4.1.1

推定結果

図 4.1 は GS が 6 個で周辺の GS の中心からの距離が 240離れた場合で、推定された大気揺らぎの図である。図 4.1 の 1 段目は与えたモデル大気ゆらぎ、2 段目は推定された大気揺らぎ、３段
目はモデルと推定値の差である。この図から LGS がカバーしている中心付近の領域では測定値に

図 4.1: LGS が 6 個、中心からの距離が 240” のときの各高さでのモデル大気ゆらぎ、推定した大気揺
らぎ、モデルと推定値の残差の図。

基づいた推定がなされているのに対し、LGS がカバーしていない外側では大気揺らぎの情報が無いので推定が行われていないことがはっきりとわかる。今回の結果の場合、GS からの情報が無い

31

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

中心から 240”よりも外側では正しい推定が行われない。つまり、重要なのは GS で測定されている 240”よりも内側の領域の推定精度である。
次に 50 ステップ繰り返したときの最終的な視野内の各方向での推定精度を見ていく。図 4.2 は視野内で推定残差の RMS(左) とストレル比 (左) の 50 ステップの平均値を示した図である。波面
残差は [nm] で表している。明らかに GS の位置での推定精度が良いことがわかる。これは推定が

図 4.2: 視野内の補正後の残差の RMS とストレル比のマップ。図の直径は 10’ である。

GS からの測定値に合うように行われているため、GS 方向に最適化されていることが原因である。
また外側の情報が無い領域では推定が正しく行われていないので推定残差が大きくなっている。このような GS 方向で特に良くなる傾向は Neichel et al. [16] でのフーリエ変換を用いた場合の結果でも見られている。

4.1.2

IRMOS と NFIRAOS

ここから推定精度についてさらに定量的に考えていく。まずは IRMOS と NFIRAOS で検討されている GS の配置で計算を行い、表 1.1 と表 1.2 で要求される精度と比較していく。IRMOS のGS の配置は前章の図 3.1 の右図で示した通りである。NFIRAOS では 6GS を視野中心に１つ、周りに 35”だけ離して５つ配置することが検討されている。
図 4.3 はそれぞれの場合での推定残差の RMS の平均値を表した図である。赤線は IRMOS、青線は NFIRAOS の結果である。また、Open Square は GS と重なる方向 (図 3.1 の青線上)、Filled

Square は GS と重ならない方向 (図 3.1 の赤線上) を表している。NFIRAOS は GS が視野の中心に集まっているので、視野中心では波面残差が非常に小さいが中心から距離が離れるにつれて急激に推定精度が悪くなる。図 4.3 の黒線は表 1.1 の NFIRAOS に対する要求精度をプロットしたものであり、この黒線より推定残差の RMS を小さくすることが望まれる。今回は理想的な場合での計算結果にもかかわらず、この要求を完全には満たすことができていない。赤線の IRMOS の場合は

NFIRAOS に比べて GS の位置が広がるため中心での残差は大きくなるが、NFIRAOS よりも広い範囲で一様に補正ができていることがわかる。

32

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

Residual WFE (IRMOS and NFIRAOS)

500

WFE [nm]

400

300

200

IRMOS
NFIRAOS
Requirement for NFIRAOS

100

0

0

50

100

150

200

Distance [arcsec]

250

300

図 4.3: 視野中心からの距離と推定残差の RMS を示した図。赤線は IRMOS、青線は NFIRAOS の結果である。Open Square は GS と重なる方向 (図 3.1 の青線上)、Filled Square は GS と重ならない方向 (図 3.1 の赤線上) を表している。黒線は NFIRAOS に対する要求精度 (表 1.1) を示している。

次に補正後の PSF を見ていく。次のページの図 4.4 は IRMOS と NFIRAOS の J バンド (λ =1250nm)、

H バンド (λ =1650nm)、K バンド (λ =2200nm) での 5 秒分を足し合わせた PSF である。縦軸は視野中心からの距離を示している。図 4.3 の波面残差と同じように視野中心から離れるほど PSF
広がりぼやけていくのがわかる。また、短波長ほど回折限界が小さいので視野中心での PSF の広
がりは小さい。しかし、同じ波面残差でも波長に対する割合が短波長ほど大きい (例えば波面残差が 200nm のとき K バンドでは約 λ/10、J バンドでは約 λ/6) ので、視野中心からはなれた場所で得られる PSF は短波長ほど広がり淡くなっている。
この PSF から IRAF の imexam で測定した FWHM を図 4.5 に、PSF から計算したストレル比を図 4.6 に、ある領域にどれだけの光が含まれているかを表す Ensquared Energy (EE) を図 4.7 に示す。
まず図 4.5 では赤が IRMOS、青が NFIRAOS の FWHM の測定結果を示している。図 4.3 の波
面残差の RMS のプロファイルでは NFIRAOS(青) は中心から 35”以上、IRMOS は 150”以降で急速に残差が大きくなるのに比べて、FWHM のプロファイルはどの波長でもそれより広い範囲で回折限界の FWHM に近い性能が出ていることがわかる。つまり、FWHM を回折限界像に近づけるために必要な推定精度は比較的緩いことがわかる。NFIRAOS(青) の場合どの波長でも中心から90”くらいまで、IRMOS(赤) の場合だいたい 210”くらいまではほぼ回折限界像の FWHM と一致しているので、図 4.3 と合わせて考えると波面残差が約 350nm 以下であれば回折限界の FWHMが得られることがわかる。
次にストレル比 (図 4.6) と EE(図 4.7) について見ていく。EE は 0.05”×0.05”の開口で計算している。どちらの図も図 4.3 の波面残差の RMS と似たようなプロファイルをしている。この図と図4.5 を比べると、FWHM は小さくても残差が大きいほど中心部より外側に光が広がっていき、中心部の光量は大きく減っていくことがわかる。これは図 4.4 の図からでも見ることができる。また、上でも触れたように FWHM はある程度の範囲で回折限界像と等しくなっており特に波長による大きな違いはない。しかしストレル比や EE では同じ波面残差でも波長による違いが非常に大きく出ており、短波長ほど33

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

図 4.4: 各波長での IRMOS と NFIRAOS で得られる PSF。縦軸は視野中心からの距離を示している。
画像のスケールは回折限界の PSF の最大値と最小値に合わしている。また対数スケールで表し
ている。

34

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

FWHM of IRMOS and NFIRAOS (Jband)
30m Difflaction Limit (FWHM 8.8mas)
8.2m Difflaction Limit (38.4mas)
IRMOS
NFIRAOS

0.14

FWHM [arcsec]

0.12
0.1
0.08
0.06
0.04
0.02
0

0

50

100

150

Distance [arcsec]

200

250

300

FWHM of IRMOS and NFIRAOS (Hband)
30m Difflaction Limit (FWHM 11.5mas)
8.2m Difflaction Limit (50.6mas)
IRMOS
NFIRAOS

0.14

FWHM [arcsec]

0.12
0.1
0.08
0.06
0.04
0.02
0

0

50

100

150

Distance [arcsec]

200

250

300

FWHM of IRMOS and NFIRAOS (Kband)
30m Difflaction Limit (FWHM 15.4mas)
8.2m Difflaction Limit (67.5mas)
IRMOS
NFIRAOS

0.14

FWHM [arcsec]

0.12
0.1
0.08
0.06
0.04
0.02
0

0

50

100

150

Distance [arcsec]

200

250

300

図 4.5: 上から J バンド (1250nm)、H バンド (1650nm)、K バンド (2200nm) での IRMOS(赤) と NFIRAOS(青) の FWHM の視野中心からの変化を表している。黒線は 30m の回折限界像の FWHMと 8.2m の回折限界である。ここで、30m の回折限界像の FWHM は Appendix A の式から計算される回折限界のサイズよりも小さくなっている。これは回折限界の定義が FWHM ではな
く最初に 0 になる暗環までの距離なので、FWHM で考えるとそれよりも小さくなる。四角の種類は評価している方向の違い (白四角が GS と重なる方向) を表している。

35

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

Strehl Ratio (NFIRAOS)

100

Kband
Hband
Jband

Strehl Ratio [%]

80

60

40

20

0

0

50

100

150

200

Distance [arcsec]

300

Strehl Ratio (IRMOS)

100

Kband
Hband
Jband

80

Strehl Ratio [%]

250

60

40

20

0

0

50

100

150

200

Distance [arcsec]

250

300

図 4.6: 上から NFIRAOS、IRMOS の各波長での PSF から計算したストレル比と視野中心からの距離の関係を表す図である。色は波長の違いを、四角の種類は評価している方向の違い (白四角が GSと重なる方向) を表している。

36

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

Ensquared Energy 0.05"x0.05" [%]Ensquared Energy 0.05"x0.05" [%]第 4 章 結果

Ensquared Energy 0.05"x0.05" (NFIRAOS)

100

Kband
Hband
Jband

80

60

40

20

0

0

50

100

150

200

Distance [arcsec]

250

300

Ensquared Energy 0.05"x0.05" (IRMOS)

100

Requirement for IRMOS
Kband
Hband
Jband

80

60

40

20

0

0

50

100

150

200

Distance [arcsec]

250

300

図 4.7: 上から NFIRAOS、IRMOS の各波長での 0.05”×0.05”の四角内に全体の何%の光が入かを表す
Ensquared Enegy(EE) と視野中心の関係を示した図である。色は波長の違いを、四角の種類は評価している方向の違い (白四角が GS と重なる方向) を表している。IRMOS の図中の黒線は
表 1.2 の直径 5’ 以内で J バンドの EE が 50%以上という要求精度を示している。

37

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

波面残差に対して敏感であることがわかる。図 4.7 の下図から、IRMOS は要求される精度をぎり
ぎりではあるがクリアしていることがわかる。しかし、本研究で目標としている視野 10’ を達成するためには外側での推定精度が足りず、改良が必要なことがわかる。
以上の結果より IRMOS に対する表 1.2 の要求精度を達成するためには、波面誤差の RMS を

200nm 以下にする必要があることがわかる。この後はこの波面残差の RMS200nm を達成すべき目標としていく。

4.1.3

ガイド星の数との関係

IRMOS と NFIRAOS の結果から、推定精度を決める要因として GS の数と配置が大きく影響していることがわかる。最初に述べたように推定視野を直径 10’ に拡張していくことを考えていくた
めに、まず は GS の数の影響を見ていく。図 4.8 の図は周辺の GS の距離を 240 で固定して、GSの数を変化させたときの推定精度を示した図である。横軸は中心からの視野中心からの角度距離、

Residual WFE (GS=6,8,10, R=240")
550

6GS, 240"
8GS, 240"
10GS, 240"
Requirement 200[nm]500

WFE [nm]

450
400
350
300
250
200
150

0

50

100

150

200

250

300

Distance [arcsec]
図 4.8: 縦軸は推定残差の RMS、横軸は視野中心からの角度距離。周辺の GS の距離を 240”に固定してガイドの数を変化させた場合の推定精度を示している。色の違いは GS の数の違いである。四角の種類は評価している方向の違い (白四角が GS と重なる方向) を表している。

縦軸は推定残差の RMS である。まず IRMOS の場合は視野 150”内に隙間を埋めるように GS を配置していたので評価する方向の違いは見られなかった。図 4.8 では視野中心の GS 以外の GS を同じ距離に配置しているため、GS と重なる方向 (白四角) と重ならない方向 (黒四角) の差がはっきりと出ている。また中心の GS と外側の GS の間での推定精度も GS の方向に比べて悪くなる。

GS の数による違いはまず視野全体での波面残差が GS が多いほど小さくなることがあげられる。
これは LGS が有限の高さにあることによる影響である。今回の場合 Na LGS を仮定しているので
LGS は望遠鏡から 88km の高さにある。天体はほぼ無限遠に存在するため天体からの光は平行光で届くが、有限の高さにある LGS の光は円錐状に望遠鏡に届く、そのため LGS を用いた場合は特に上層で天体の光が通過する経路を LGS の光ではカバーできない領域が生じることによる誤差が生まれる。この影響をコーン効果、もしくは Focal Anisoplanatism と呼ぶ。この影響は口径が大きくなるほど強くなり、口径 30m、K バンドではコーン効果によって生じる誤差は 400nm ほどにもなる (Appendix A)。しかし、今回の場合は複数の LGS を用いているので、それぞれの LGS が各 LGS が測定できない領域を補い合いこの影響が小さくなる。実際に図 4.8 では波面残差の RMS

38

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

が 400nm よりもずっと小さくなっている。これも ELT で複数の LGS を用いた補償光学が必要な理由の１つである。しかし、GS の数が減るとその効果が薄れ GS の方向でも残差が大きくなっていく。また、視野中心の GS 方向より外側の GS 方向の精度が悪いのは視野中心の GS は他の GSに囲まれているのでコーン効果の影響をより抑えられているためである。さらに外側の GS 方向よりもその少し内側で推定残差が小さくなるのは、GS よりも内側でのみ他の GS の情報によりコーン効果が抑えられているからである。
この結果より外側の GS を 240”離したときには、GS を 10 個置けば GS と重なる方向では半径250”くらいまでで推定残差 200nm 以下を達成できる。しかし、GS と重ならない方向や 250”よりも外側では目標値を達成することができていない。

4.1.4

ガイド星の距離との関係

次に GS の配置する距離の影響を見ていく。図 4.9 は GS の数を 8 つに固定して、周辺の GS を配置する距離を変化させた結果を示している。まず、GS の配置する距離が近いほど視野全体でのResidual WFE (GS=8, R=180",240",300")
550

8GS, 180"
8GS, 240"
8GS, 300"
Requirement 200[nm]500

WFE [nm]

450
400
350
300
250
200
150

0

50

100

150

200

250

300

Distance [arcsec]
図 4.9: 縦軸は推定残差、横軸は視野中心からの距離。GS の数を 8 個に固定し、周辺の GS の距離を変化させた場合の推定精度を示している。色の違いは GS の数の違いである。四角の種類は評価している方向の違い (白四角が GS と重なる方向) を表している。

推定精度が良くなる。これは先ほどと同じくコーン効果が原因であり、GS を密に配置した方がそれぞれの GS がカバーし合う面積が大きくなるためコーン効果の影響が抑えられる。また密に配置するほど中心の GS と外側の GS の間での推定精度の悪化も弱くなる。

4.1.5

4.1 の結果のまとめ

まず IRMOS での GS の配置で計算を行い、J バンドでの 0.05”×0.05”の EE が 50%以上という要求精度が波面残差の RMS で 200nm 以下に対応することを確認した。また、GS の数や配置する
距離の推定精度に対する影響を調べ、GS を密に配置するほど視野全体の推定精度が向上することがわかった。高度の低い層では各 GS の光は同じような経路を通過するため、密に配置することによるメリットは高度の高い位置で様々な方向の大気揺らぎの情報が得られることである。つまり、高高度の大気揺らぎの情報をいかに増やすかが推定精度の向上につながることが考えられる。しか39

第 4 章 結果

4.1. 最小分散推定によるトモグラフィック波面再構成の結果

し中心の GS と外側の GS の間は WFS が測定している領域には含まれているため、推定精度はもう少し良くできると考えられる。次章ではこの GS 方向以外の方向で推定精度が悪くなる原因についてもう少し詳細まで議論する。

40

第 4 章 結果

4.2

4.2. MOAO+GLAO の結果

MOAO+GLAO の結果

本節では MOAO+GLAO システムのための Average Method の結果を示していく。

4.2.1

理想的な場合

まず、仮に地表層を 100%取り除けた場合、どのくらいの精度向上が期待されるかを見積もった、
ここでは地表層は完全に省いて 6 層でトモグラフィー計算を行った。図 4.10 がその結果で、上図は GS と重なる方向、下図は重ならない方向のプロファイルである。図の丸印が 6 層で計算した理想的な Averege Method の結果である。全体の 60%近くを占める地表層を除いてトモグラフィー
Ideal Average Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Ideal Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150
100

0

50

100

150

200

250

300

Distance [arcsec]

Ideal Average Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Ideal Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150
100

0

50

100

150

200

250

300

Distance [arcsec]
図 4.10: 縦軸は推定残差、横軸は視野中心からの距離。周辺のガイド星の中心からの距離を 240 に固定
して、ガイド星の数を変化させた。四角印は Normal Method の結果、丸印は理想的な Average
Method(6 層でトモグラフィー計算) の結果である。上の図は GS と重なる方向、下図は重ならない方向のプロファイルである。

計算を行ったので、理想的な Average Method の方が推定精度が良くなっている。特に GS が多くなるほどその傾向が強い。この結果から精度良く地表層を除くことができれば、Average Methodによるトモグラフィー計算の推定精度の向上が期待されると考えられる。

41

第 4 章 結果

4.2.2

4.2. MOAO+GLAO の結果

平均値を用いた結果

次に実際にすべての GS の平均値から地表層の推定値を求める Average Method 行った結果について見ていく。図 4.11 がその結果である。ここでは四角印は Normal Method の結果、三角印はAverage Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150

0

50

100

150

200

250

300

Distance [arcsec]

Average Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150

0

50

100

150

200

250

300

Distance [arcsec]
図 4.11: 縦軸は推定残差、横軸は視野中心からの距離。周辺のガイド星の中心からの距離を 240 に固定して、ガイド星の数を変化させた。四角印は Normal Method の結果、三角印は Average
Method の結果である。上の図は GS と重なる方向、下図は重ならない方向のプロファイルである。

Average Methoid の結果である。図からわかるようにどちらの方法でも特に優位な差はなかった。
また、GS の数も変えても特に変化は見られなかった。

4.2.3

4.2 の結果のまとめ

理想的な場合では Average Method を用いることによって精度が向上した。しかし、実際に平均値を用いて計算を行うと Normal Method と比べて大きな違いは見られなかった。次の章ではこの理由について考察し、また平均値による地表層の推定値を用いて GLAO としてどのくらいの性能が期待されるかを考えていく。

42

第 4 章 結果

4.3

4.3. 風速を用いた結果

風速を用いた結果

ここでは Wind Method の結果を見ていく。ここではそのとき測定した情報を”現在の測定値 、前に測定した情報を”過去の測定値4.3.1と呼ぶことにする。

すべての大気ゆらぎの層に同じ風速を与えた場合

まずシンプルに考えるためにすべての大気揺らぎの層に同じ風速を与え、風速の強さを変化させて風速の強さの影響を見る。今回は 1m/s、3m/s、5m/s、7m/s、10m/s の風速を与えた場合で、1s 前の測定値を用いて計算する。図 4.12 は GS を 8 個で中心から 270 の位置に配置した場合のWind Method (8GS, 270")

400

Normal Method
Wind 1m/s
Wind 3m/s
Wind 5m/s
Wind 7m/s
Wind 10m/s

350

WFE [nm]

300
250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

Wind Method (8GS, 270")

400

Normal Method
Wind 1m/s
Wind 3m/s
Wind 5m/s
Wind 7m/s
Wind 10m/s

350
300

WFE [nm]

300

250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

300

図 4.12: すべての大気揺らぎの層に同じ風速を与えた場合の Wind Method の結果。ここでは GS の数は８個で、外側の GS を配置した距離は 270”とする。赤は Normal Method の結果、他は風速の情報を用いた場合の結果である。上図は GS と重なる方向、下図は重ならない方向のプロファイルである。

結果を表した図である。赤が Normal Method、他が Wind Method の結果である。この計算では

1 ステップのみ計算を行った。明らかに風速の情報を用いて計算を行ったほうが推定精度が向上している。また、風速が強いほど推定精度が向上していることがわかる。風速が強いほど過去の測定

43

第 4 章 結果

4.3. 風速を用いた結果

値は現在の測定値からずれた位置の情報となる。よりずれた位置の情報を用いることで、現時刻での測定だけでは情報が薄い箇所を補うことができるので風速が強い方が精度が向上したと考えられる。この計算では外側の GS を中心から 270”離した位置に配置したので、最も上層の大気揺らぎがある 16000m では外側の GS の光は中心から約 21m 離れた経路を通過する。風速が 10m/s
で 1s 前の情報を使うということは 10m だけずれた情報を用いることになるので、ちょうど GS とGS の間の経路を過去の測定値で埋めることができる。その結果、GS と重なる方向 (図 4.12 の上図) では半径 5’ の視野内で目標である「波面残差の RMS が 200nm 以下」を達成できている。GSと重ならない方向 (図 4.13 の下図) でも大きな推定精度の向上が見られる。

4.3.2

ある高さの大気ゆらぎの層にだけ風速を与えた場合
Wind Method (8GS, 270")

400

Normal Method
Wind 5m/s (All layers)
Wind 5m/s (low layers)
Wind 5m/s (Middle layers)
Wind 5m/s (High layers)

350

WFE [nm]

300
250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

300

Wind Method (8GS, 270")

400

Normal Method
Wind 5m/s (All layers)
Wind 5m/s (low layers)
Wind 5m/s (Middle layers)
Wind 5m/s (High layers)

350

WFE [nm]

300
250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

300

図 4.13: ある高さの大気揺らぎのみに風速を与えて Wind Method を用いた場合の推定残差。赤は Normal
Method、青はすべての大気揺らぎに 5m/s の風速を与えた結果、それ以外が下層、中層、上層それぞれにのみに風速を与えて Wind Method を用いた結果である。上図は GS と重なる方
向、下図は重ならない方向のプロファイルである。

次にある高さの大気揺らぎだけに風速を与えた場合を考える。7 層ある大気揺らぎのうち 0m と500m を下層、1000m と 2000m と 4000m を中層、8000m と 16000m を上層としてそれぞれに風速44

第 4 章 結果

4.3. 風速を用いた結果を与え、Wind Method の手法を用いて計算を行う。つまり風速を与えてない高さの層は過去の測定値と現在の測定値が等しくなる。ここでも計算は 1 ステップのみ行う。図 4.13 の赤は Normal

Method、青はすべての大気揺らぎに 5m/s の風速を与えた結果、それ以外が下層、中層、上層それぞれに風速を与えて Wind Method を用いた結果である。下層、中層、上層のどの層に風速を与えて Wind Merhod を用いた場合でも Normal Method に比べて残差が小さくなっているが、明らかに上層のみに風速を与えた場合の推定残差が他に比べて非常に小さくなっていることがわかる。
下層や中層ではどの角度から来た光もほぼ同じ経路を通過するのに対し、上層では各 GS の光は別の経路を通過するので得られる情報が薄くなりがちである。そのため今回の結果では Wind

Method で下層よりも上層の情報を増やしたほうが推定精度が向上したと考えられる。つまり Wind
Method で上層の情報を増やすことは 4.1 で見た GS を密に配置することに対応するといえる。また、すべての大気揺らぎを風速でずらして Wind Method を使用した結果 (青線) よりも、上層のみ
を風速でずらして Wind Method を用いた結果 (水色) のほうが推定精度が向上している。これは下層では過去の測定値が風速によりずれるほど最終的な補正に関係する領域からはみだしていくので、逆に計算に悪影響を及ぼしていると考えられる。以上より Wind Method を用いる場合は上層で風速が強く、下層では風速が弱いような状況で大きな推定精度向上が期待できることがわかる。

4.3.3

風速のモデルで計算した場合

次に図 3.3 の風速のモデルを仮定した場合で Wind Method を用いて計算を行う。今回は 50 ス
テップ計算し、5 秒分の推定残差を求めた。図 4.14 がその結果である。ここでは 0.01s 前、0.05s
前、0.1s 前、0.2s 前、0.5s 前の情報を用いたときの結果を載せている。より過去の情報を用いた方が風速によるずれ量が大きいので推定精度が向上していくが、0.5s 前の情報を用いた時は 0.2s 前の情報も用いた場合より推定残差が大きくなっている。これは図 4.13 で見たように下層の風速によるずれが大きくなったためであると考えられる。ここで用いた風速のモデルは図 3.3 の図からもわかるように上層で風速が強く、下層で風速が弱くなるので Wind Method に適した状況であるので、そこまで昔の情報を用いなくても推定精度の向上が期待される。

45

第 4 章 結果

4.3. 風速を用いた結果

Wind Method (8GS, 270")

400

Normal Method
0.01s
0.05s
0.1s
0.2s
0.5s

350

WFE [nm]

300
250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

Wind Method (8GS, 270")

400

Normal Method
0.01s
0.05s
0.1s
0.2s
0.5s

350
300

WFE [nm]

300

250
200
150
100

0

50

100

150

Distance [arcsec]

200

250

300

図 4.14: 図 3.3 のモデルから風速を与えた場合の結果。赤が風速の情報を用いない Normal Method の結果、他が風速の情報を用いる Wind Method の結果である。Wind Method は 0.01s 前、0.05s
前、0.1s 前、0.2s 前、0.3s 前の情報を用いたときの５つの結果を載せている。

4.3.4

風速の推定に誤差がある場合

さらに風速の推定に誤差がある場合を考える。風速は前節のモデルの値を用いて、0.1s 前の情報で Wind Method を計算する。風速の推定誤差による影響はどのくらい過去の測定値を使用する
かで変わってくる。そこで、ここでは風速の推定誤差によってずれる距離で風速の推定誤差の大きさを表す。まずすべての大気揺らぎの層が風速の推定誤差によって考慮していた位置より 0.1m、0.2m、0.5m ずれたときの推定精度を計算する。今回は 0.1s 前の測定値を使っているので、風速の推定誤差にすると 1m/s、2m/s、5m/s となる。図 4.15 がその結果である。青線は風速の推定誤差のない Wind Method の結果、赤線は Normal Method の結果、その他が風の推定誤差を入れた
Wind Method の結果である。風速の推定誤差によるずれが 0.1m 以下のときはその影響はほとん
ど見られない。しかし、風速の推定誤差がそれ以上大きくなると徐々に推定精度も落ちていき、風速の推定誤差によるずれが 0.5m のときは Normal Method よりも推定精度が悪くなる。
次に大気揺らぎを上層 (16000m、8000m)、中層 (4000m、2000m、1000m)、下層 (500m、0m)

46

第 4 章 結果

4.3. 風速を用いた結果

Wind Method

400

Normal Method
350

No Noise

WFE [nm]300

Noise = 0.1m

250Noise = 0.2m

200

Noise = 0.5m

150
100

0

50

100

150

200

Distance [arcsec]

250

300

図 4.15: ノイズを与えた場合の結果。青線はノイズが無い場合の Wind Method の結果、赤色は Normal
Method の結果、その他が風速の推定にノイズを入れた Wind Method の結果である。

の３つに分けてそれぞれに風速の推定誤差をのせ、高さごとの風速の推定後さの影響を見ていく。
図 4.16 がその結果であり、上から下層、中層、上層の結果である。この図から同じ誤差をのせても下層の推定残差が中層や上層と比べて大きくなることがわかる。これは今回用いている大気揺らぎモデルが 0m で 60%近いパワーを持つためであると考えられる。他の高さでは風速の推定誤差によるずれが 0.5m になっても Normal Method からは十分精度が向上している。この結果から下層でどれだけ風速を正確に推定できるかが重要となることがわかる。

4.3.5

4.3 の結果のまとめ

風速の推定誤差がない理想的な場合では Normal Method に比べて Wind Method を使用した方がに大きく推定精度が向上しこの手法が有効であることがわかった。とくに今回用いてる風速のモ
デルは上層で強い風速を持つので、より Wind Method に適しているといえる。しかし、風速の推定誤差が大きい場合は Normal Method よりも精度が悪くなる。次の章では風速の推定誤差による影響を抑える手法と、この Wind Method が現実的に使えるかどうかを議論する。

47

第 4 章 結果

4.3. 風速を用いた結果

Wind Method (Low)400

Normal Method
350

No Noise

WFE [nm]300

Noise = 0.1m

250Noise = 0.2m

200

Noise = 0.5m

150
100

0

50

100

150

200

Distance [arcsec]

250

300

Wind Method (Middle)

400

Normal Method
350

No Noise

WFE [nm]300

Noise = 0.1m

250Noise = 0.2m

200

Noise = 0.5m

150
100

0

50

100

150

200

Distance [arcsec]

250

300

Wind Method (High)400

Normal Method
350

No Noise

WFE [nm]300

Noise = 0.1m

250Noise = 0.2m

200

Noise = 0.5m

150
100

0

50

100

150

200

Distance [arcsec]

250

300

図 4.16: 上から下層、中層、上層にノイズを与えた場合の結果。青線はノイズが無い場合の Wind Methodの結果、赤色は Normal Method の結果、その他が風速の推定にノイズを入れた Wind Methodの結果である。

48

第 4 章 結果

4.4. 計算速度

計算速度

4.4

この章の最後に計算速度について説明する。

4.4.1

収束条件

計算速度に大きく関係するものとして CGM での計算の反復回数がある。この反復回数を決めるのは主に収束の早さと収束条件のきつさである。ここでは GS が７個、中心からの距離を 180”
離した場合で収束条件を考えていく。今回の計算で収束条件は相対誤差を用いている。k 回目の推定値を φkmv とすると、式 (3.3) から相対誤差は rk = ||s̄ − ḠT Ḡφkmv ||2 で定義される。この相対誤差がある値 r 以下になったとき収束したと考え計算を終了する。ここでは r の値を変えて実際に計算を行い、r の値と精度の関係について見ていく。図 4.17 は各 r の値での 16000km の大気揺

図 4.17: 収束条件 r を変えたときの 16000km の推定値の図。n は反復回数。赤線は GS がカバーしているおおよその領域。

らぎの推定値の画像である。当然だが r の値を小さくする、つまり収束条件を厳しくすると反復回数が増加する。また推定された画像を見てみると GS がカバーしている領域では推定がすぐに行
われ、GS がカバーしていない領域では徐々に内側の値が外側にのびていくような推定となっている。今、赤線の外側で測定値は無いので推定は共分散行列のみで行われているためどんなに計算を行っても正しい推定は行われない。重要となってくるのは赤線の内側の推定であり、図 4.17 では

r = 10−3 ではすでに赤線の内側はほぼ推定されているように見える。
次に各方向での推定残差を見ていく。図 4.18 は横軸は視野中心からの角度距離、縦軸は各方向
Iteration Criteria (7GS , R=180")600-2550 r=10r=10-3500 r=10-4-5
450 r=10-6
r=10

400
350
300
250
200

r=10-2r=10-3500 r=10-4-5
450 r=10-6
r=10
550

Residual WFE [nm]Residual WFE [nm]

600

400
350
300
250

0

25

50

75 100 125 150 175 200 225 250 275 300

Distance from the center [arcsec]

200

0

25

50

75 100 125 150 175 200 225 250 275 300

Distance from the center [arcsec]

図 4.18: 収束条件 r を変えたときの各方向の推定値残差。横軸は評価した方向、縦軸は残差の RMS である。色の違いは r の値の違いである。左側の図は GS と重なる方向、右の図は GS と GS の間の方向での結果である。

49

第 4 章 結果

4.4. 計算速度での推定残差の値であり、各方向での推定精度を表した図である。この図から、 r = 10−2 では推定しきれずに精度が他と比べてやや悪くなるが、r ≤ 10−3 の場合では特に GS より内側でほとん
ど推定結果が変化していないことがわかる。図 4.17 と図 4.18 より、まず GS より情報がある赤線より内側で推定が行われ、その後外側の推定が行われていると考えることができる。しかし赤線より外側では GS からの情報はまったくないので計算を繰り返しても正しい推定がされることはない。つまり、内側の推定が終わった頃に計算を打ち切るのが最も効率がいい。今回はこの結果から

r = 10−3 を採用する。このときの計算の反復回数は GS の数や配置によらず 300 回前後であり、計算時間にすると約 1.5 秒であった。Wind Method の場合、行列の大きさ、測定値の量が増えるた
め 3s ほどかかった。

4.4.2

各過程の計算に要する時間

さて、ここではさらに詳しく計算時間について見ていく。今回の計算はA.

測定値 s̄ を CPU メモリから GPU メモリに移す

B.

GPU メモリ上で (ḠP )T s̄ = z

C.

GPU メモリ上で (ḠP )T (ḠP )x = z に対して共役勾配法を行う。

D.

GPU メモリ上で φmv = P x により推定値を求める。

E.

最後に GPU メモリ上の推定値のデータを CPU メモリ上に移す。

共役勾配法の詳細は Appendix B を参照。

という過程で行われる。ここで P は対角スケーリングの前処理行列である (Appendix B)。また、ḠP や (ḠP )T は事前に計算してある。実際に各過程でかかった時間を表 4.1 にまとめる。表 4.1
から計算時間の大部分は行列とベクトルのかけ算が占めていることがわかる。特に共役勾配法の反復回数が増加するほど行列とベクトルのかけ算の回数が増えるため、全体の計算時間も増加する。
今回の計算では反復回数は 250 回〜300 回であるので、共役勾配法で約 1.5〜1.8s の計算時間が必要となっている。前節で見たように今回計算にかかった時間は約 1.5s であるので共役勾配法を解くのにかかる時間がほぼ全体の計算時間となっていることがわかる。
この結果から計算時間を減らすためにはまず共役勾配法の反復回数を減らすことが必要である。
今回用いた対角スケーリングによって、反復回数は約半分になった。この対角スケーリングは前処理の中でもっともシンプルな手法であり、そこまで大きな効果は期待できない。共役勾配法の前処理は盛んに研究されている分野であり、様々な手法があるのでいろいろ試してみる必要がある。また、より高性能なパソコンや GPU を用いて基礎的な計算時間の向上も重要である。

50

第 4 章 結果

4.4. 計算速度

過程

計算時間 [ms]A. CPU から GPU へ s̄ のデータ転送
T

B. z = (ḠP ) s̄

0.64
2.58

C. 共役勾配法 (Appendix B 参照)
for 文より前

for 文

C.0. B = |z|0.3

C.1. t = Ax0

5.5

C.2. r0 = z − t

0.54

C.3. p0 = r0

0.005

C.4. R0 = |r0 |0.35

C.5. 収束の判定

0.001 以下

C.6. qk−1 = Apk−1

5.4

C.7. ek = pk−1 · qk−1

0.32

C.8. αk−1 = Rk−1 /ek

0.001 以下

C.9. xk = xk−1 + αk−1 pk−1

0.004

C.10. rk = rk−1 − αk−1 qk−1

0.004

C.11. Rk = |rk |0.3

C.12. 収束の判定

0.001 以下

C.13. ck = rk · qk−1

0.3

C.14. βk = −ck /ek

0.001 以下

C.15. pk = rk + βk pk−1

0.004

for 文内の合計

6.332

D. φmv = P x

0.91

E. GPU から CPU へ φmv のデータ転送

0.71

表 4.1: 各過程に要する計算時間。ここで Ax = (ḠP )T (ḠP )x を表している。測定できる計算時間は
1ms 単位であるので、ここでは同じ操作を 100 回繰り返したときにかかった計算時間を 100 で割った数を載せている。また、100 回の計算でも 1ms 以下の計算時間である場合は 1000 回の計算を行った。

51

第 5 章 議論

第5章

5.1

議論

誤差評価

まずトモグラフィー計算による誤差の寄与を評価する。今回の計算で影響していると考えられる
誤差はトモグラフィー計算の推定誤差、DM の分割鏡のサイズよりも小さなスケールのゆらぎによる影響 (フィッティングエラー) がある。コーン効果による影響は GS の配置や数に依存するのでト
モグラフィーによる推定誤差に含めるとする。よって今回の計算での誤差は以下のように書くことができる。

σtotal =√2
σtomo+ σf2 it

(5.1)

Hardy [24] によればフィッティングエラーによって発生する位相誤差 σf it [rad] は式 (5.2) で記述される。

(
σf2 it = aF

d
r0

)5/3
(5.2)

ここで d は DM で補正できる揺らぎの最小スケール (分割鏡のサイズ、アクチュエイターの間隔) である。今回の場合は DM は推定結果を完璧に再現するものを考えているので、WFS の subapertureのサイズが d となる。aF は DM の構造によって決まる定数である。今回の DM は波面の局所的な傾きを推定値に従って補正するだけであり、隣り合うアクチュエイター同士の影響はないと考えている。このような場合では Hardy [24] より aF =0.18 となることがわかっている。r0 の値は観測値 0.156m(@500nm) からコロモゴロフのパワースペクトルを仮定して計算すれば、図 3.5 より

0.9231m となる。この場合、式 (5.2) を用いて計算すると、フィッティングエラーは
√
(
)5/3
0.5
σf it = 0.18
= 0.2545[rad] = 89.1210[nm]0.9231(5.3)となる。３章で見たようにアウタースケールの影響を考慮した場合に得られるシーイングは図 3.5のように小さくなり、このシーイングから r0 を計算すると 1.4143m となる。当然、この値からフィッティングエラーを計算すれば式 (5.3) の値より小さくなる。しかし、アウタースケールの影
響を考慮して変化するのは大きなスケールの揺らぎであり、フィッティングエラーに影響するような小さなスケールの揺らぎは特に変化しないと考えられるため、フィッティングエラーに関しては
コロモゴロフを仮定して計算した式 (5.3) の値を採用する。シミュレーションでは、完璧に推定ができた場合の残差がフィッティングエラーとなるはずである。実際に計算を行い、視野全体で平均すると 80.3711[nm] となり、式 (5.3) の値とコンシステントである。この結果から、今回の目標値である波面残差の RMS が 200nm 以下を達成するために満たすべきトモグラフィー計算での誤差√は σtomo = 2002 − 80.37112 = 183.1406[nm] であることがわかる。つまり 200nm 以下という目
標値の場合、フィッティングエラーを小さくすることよりもトモグラフィー計算をいかに精度良く行うかが重要であることがわかる。フィッティングエラーを除いて図 4.3 の IRMOS と NFIRAOSの推定残差の図を書き直すと図 5.1 のようになる。NFIRAOS の場合の最小の推定残差は 100nm、

IRMOS の場合 150nm くらいである。IRMOS の場合は推定残差が支配的であることがわかる。
52

第 5 章 議論

5.1. 誤差評価

Residual WFE (subtracted Fitting Error)

500

WFE [nm]

400

300

200

IRMOS
NFIRAOSFitting Error

100

0

0

50

100

150

200

Distance [arcsec]

250

300

図 5.1: Fitting Error を分離して、IRMOS と NFIRAOS で生じる推定誤差をプロットしたもの。

この他に今回考慮しなかったものに WFS の測定と DM による補正の間の時間差で生じる Temporal

Error がある。時間差が生じる原因としては計算にかかる時間や WFS の露出時間である。Hardy[24] によれば遅延時間 τ が発生した時、風速 v で動く 1 層の大気揺らぎによって生じる波面誤差
σT D は√()5/3
0.427vτ
σT D = 28.4
[rad]
(5.4)
r0
と書ける。複数の層がある場合は各層から受ける影響の自乗和の平方根で最終的な影響を計算できる。また、各大気揺らぎの層の r0 は Appendix A の r0 の定義より r0 = r0 × (w[i] )−3/5 となる。
[i]

ここで w[i] は i 番目の大気揺らぎの層の揺らぎパワーの割合である。各大気揺らぎの風速を図 3.3で、揺らぎのパワーの割合を表 3.1 から与えて計算すると図 5.2 のような図が書ける。横軸は遅延Error Estimation
Requirement

200
175

WFE [nm]

150
125
100
75
50
25
0

0

1

2

3

4

Temporal Error
Limit for Tomographic Error
Fitting Error
5
6
7
8
9

Time Delay [ms]10

図 5.2: 横軸は遅延時間 [ms]、縦軸は発生する波面誤差の RMS[nm] である。赤線は遅延時間に対して生じる Temporal Error、緑線は先に計算した Fitting Error 80.3711[nm] を示している。青線
は個 Temporal Error と Fitting Error が決まった時、今回の目標値である 200nm を達成する
ために満たすべきトモグラフィー計算の推定残差の値である。

時間 [ms]、縦軸は波面残差の RMS[nm] である。赤線は生じる Temporal Error の値を示しており、

53

第 5 章 議論

5.2. ガイド星とガイド星の間で推定が悪くなることについて

遅延時間が大きくなるほど誤差が大きくなることがわかる。緑の線は先に計算したフィッティング
エラーの値である。Temporal Error とフィッティングエラーを考慮した時、目標値
200nm を達成√

2002 − σT2 D − σf2 it ) は青線のようになる。遅延時間が増えるほど Temporal Error 増加するのでそれに伴いトモグラフィー
計算で許される推定残差も小さくなる。図 5.1 より、今回の計算で IRMOS が達成している推定残するために満たすべきトモグラフィー計算の推定残差 (つまり σtomo =差は約 150nm であるので 200nm を達成するためには最低でも遅延時間を 5ms 以下にする必要がある。さらに視野の外側まで考えるとさらに遅延時間を小さくする必要がある。遅延時間を小さくする手段としてはまず WFS の露出時間を短くすることが考えられるが、WFS の露出時間を短くするほど WFS の subaperture で得られる星像の SN が下がり、今度は WFS の測定誤差が増大する危険がある。WFS の測定と推定計算は平行して行うシステムもあるため、はっきりとした値は確定できないが、Temporal Error と測定誤差を抑えるためにはトモグラフィー計算にかかる時間を最高でも 5ms 以下にすることが必要であるといえる。

5.2ガイド星とガイド星の間で推定が悪くなることについて

Normal Method で中心の GS と周辺の GS の間の領域で推定精度が悪くなることについてさらに掘り下げて考えていく。

5.2.1

各高さでの推定精度

次のページの図 5.3 は GS と重なる方向での高さごとの残差を表した図である。例えば一番左上の図は 0m での推定残差であるが、0m はすべての方向で共通なので方向 (横軸) で違いは無い。
すべての高さで推定残差が補正前の波面誤差 (灰色) よりも小さくなっているが、補正量は各高さでばらつきがある。特に 500m、1000m の高さの補正量が他の層に比べて小さい。今回用いている大気モデル (表 3.1) では 1000m の高さの揺らぎのパワーが最も小さいので、推定値もそれに伴い小さくなりすぎている可能性がある。500m の層は２番目に強いパワー (約 10%) を持っているので、補正量が小さいのは揺らぎのパワーが原因とは考えにくい。図 4.1 の２段目の 500m の推定値を見てみると、明らかに 500m のモデルよりも 0m のモデルを再現している。この２層は非常に高さが近いので、推定が分離できていないことが予想される。また 0m は非常に大きい揺らぎのパ
ワー (60%) をもつため、500m の推定値もそれに引きずられて 0m の大気揺らぎを再現している可能性が高い。
図 5.3 は GS と重なる方向の推定残差であり、すべての層を合計して考えると図 4.3 の青いプロッ
トと一致する。しかし、図 5.3 ではいくつかの層で図 4.3 の全体の推定残差よりも大きい残差が残っている。また、図 4.3 で見られるような GS と GS の間で推定が悪くなるような傾向は図 5.3 では見られない。

√
σ12 + σ22で与えられる。同様にそれぞれの高さの推定誤差が独立である場合、最終的な推定誤差 σALL は式
あるノイズ σ1 と σ2 が独立な場合、これらの２つの影響の合計のノイズは２乗和の平方根(5.5) のように計算される。

σALL

v
u 7
u∑
σ2
=t
i

(各高さの RMS が独立である場合)

i=1

54

(5.5)

第 5 章 議論

5.2. ガイド星とガイド星の間で推定が悪くなることについて

ここですべての大気揺らぎの層を合わせて考えた図 4.3 の最終的な推定誤差を σreal と表記する。
実際に式 (5.5) の結果から σALL を計算してみると図 5.4 のようになる。明らかに σALL 6= σreal であり、各高さの推定値が関係していることがわかる。

55

第 5 章 議論

430
Residual WFE [nm]

800
700
600
500
400
300

Residual WFE [nm]Layer 0m

270
265
260
255
250
245
240
235
230
225
220

460

Layer 1000m

Layer 4000m
Residual WFE [nm]Residual WFE [nm]

400
380
360
340
320

Residual WFE [nm]0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

380
370
0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

Layer 2000m

360
340
320
300
280
0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

Layer 8000m

400
350
300
250
200
150
100

0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

Layer 16000m

400

4GS

350

6GS

300

8GS

250

10GS

200

No Correction

150
100

390

450

420

450

400

260

440

300

410

380

0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

Layer 500m

420

360

0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

Residual WFE [nm]Residual WFE [nm]9005.2. ガイド星とガイド星の間で推定が悪くなることについて

0 25 50 75 100 125 150 175 200 225 250 275 300
Distance [arcsec]

図 5.3: GS と重なる方向での各高さでの推定残差を示した図。灰色のプロットは補正前の大気揺らぎのRMS である。色の違いは GS の違いである。各図で横軸のスケールが異なるので注意。

56

第 5 章 議論

5.2. ガイド星とガイド星の間で推定が悪くなることについて

1000
900
800

σALL

WFE [nm]

700

4GS
6GS
8GS
10GS

600

σreal

500
400
300
200
100

0

50

100

150

Distance [arcsec]

200

250

300

図 5.4: 丸印はすべて層の合計の RMS。四角印は式 (5.5) を用いて計算した値である。色の違いは GSの数の違いである。明らかに σALL 6= σreal であるので、推定が各高さ方向に分解できていない
ことがわかる。

図 5.3 から見られる GS の数による違いは、16000m と 8000m の 2 層で GS の数が多いほど推定残差が小さくなっている。図 4.3 で見られた GS の数が増えるほど視野全体の精度が向上する理由の１つは、この上層の推定精度の向上であることが考えられる。また、GS の数を増やしても低高度の層を分離して推定することができないことがわかる。

上層と下層での推定精度

5.2.2

前節から各高さの推定が強く関係しており、推定が分解できていないことがわかった。本節では上層と下層に分けてさらに推定の高さ方向のプロファイルについて見ていく。下から５層を下層、残りの上２層を上層として考えていく。まず下層の σreal (5 層で計算) と σALL (5 層で計算) を比べてみる。その結果が図 5.5 である。
σreal

Low Layer (0m, 500m, 1000m, 2000m, 4000m)

4GS
6GS
700 8GS
10GS
600

800
700

WFE [nm]WFE [nm]

800

σALL

500
400

600
500
400

4GS
6GS
8GS
200 10GS

300

300

200
0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

図 5.5: 左図は σreal (5 層で計算)、右図は σALL (5 層で計算) の結果。明らかに σALL 6= σreal である。

左図は σreal (5 層で計算)、右図は σALL (5 層で計算) の結果であるが、σALL 6= σreal であること

57

第 5 章 議論

5.2. ガイド星とガイド星の間で推定が悪くなることについてが一目瞭然である。これは先ほども述べたように推定が高さごとに分解できていないことを示している。また、GS の数が増えるほど推定残差の値は下がっていくので、下 5 層をまとめて考えると精度良く推定が行えていることがわかる。
次に上層についても同様に見ていく。図 5.6 は図 5.5 と同様に上層に対して作成した図である。や
σreal

450

450

4GS
6GS
8GS
10GS

400
350
300
250
200
150

σALL

4GS
6GS
8GS
10GS

400

WFE [nm]WFE [nm]High Layer (8000m, 16000m)

350
300
250
200

0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

150

0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

図 5.6: 図 5.5 と同様に上層対して作成した図。

や σreal のほうが小さいが、σALL ≈ σreal であるようにも見える。つまり、8000m の層と 16000mの層の推定誤差は互いに独立であり、推定が分離できていると考えることができる。また、下層と同様に GS が数が多いほど精度が向上している。しかし、GS と GS の間の精度が特別悪くなる傾向は見られない。
最後に下層と上層が独立かどうかを見ていく。GS が６個と 10 個の場合について、σreal (7 層で計算) と σALL (下層と上層から計算) を求めていく。図 5.7 が計算結果である。σreal (紫) と σALL (緑)は一致せず、下層と上層でも推定が混ざっていることがわかる。下層と上層を合わせて考えるとGS と GS の間が悪くなるが、GS 方向は特に変化がない。このことから GS 方向ではすべての層の影響を足すと誤差がキャンセルされるように推定が行われていることがわかる。しかし、GS 以外の方向では GS の方向に最適化している影響を受けて精度が悪くなってしまうということ起きている。GS を増やした場合全体的な推定残差が小さくなる。GS を増やして特に変わることは上層での情報量であるので、GS を増やすと図 5.6 で見たように上層の推定精度が向上し、それに伴って
図 5.5 で見られるように下層の推定精度も向上するというこが考えられる。また、前章の風速をある層だけに与えて Wind Method を用いた場合 (図 4.13) でも上層の大気揺らぎの情報を増やした
ときが最も推定残差が小さくなることから、やはり上層でどれだけ密に情報を得られるかが推定精度を決める大きな要因となっている。

58

第 5 章 議論

5.2. ガイド星とガイド星の間で推定が悪くなることについてUpper Layers and Lower Layers 6GS R=240"

550

Upper Layers RMS(8000m-16000m)
Lower Layers RMS (0m-4000m)
σALL
450
σreal

Residual WFE [nm]

500

400
350
300
250
200
150

0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

Upper Layers and Lower Layers 10GS R=240"

550

Upper Layers RMS(8000m-16000m)
Lower Layers RMS (0m-4000m)
σALL
450
σreal

Residual WFE [nm]

500

400
350
300
250
200
150

0

25

50

75

100 125 150 175 200 225 250 275 300

Distance [arcsec]

図 5.7: 上層と下層の結果から σreal と σALL を計算した。上の図が GS が 6 個の場合、下の図が GS が10 個の場合である。

59

第 5 章 議論

5.3. GLAO+MOAO の結果について

5.3

GLAO+MOAO の結果について

5.3.1

結果が Normal Method と変わらないことについて

４章では理想的な場合として 0m の地表層完全に取り除いて計算を行った。しかし、平均値で求める場合は WFS の subaperture より小さなスケールの揺らぎは測定できないので、フィッティン
グエラーは必ず生じる。0m の地表層で生じるフィッティングエラーは式 (5.3) と地表層でのフリードパラメータ r0 (0.5960)−3/5 より 69.8nm と見積もることができる。図 4.10 の図にこの 0m でのフィッティングエラーを追加すると図 5.8 のようになる。丸印が理想的に地表層を抜いて計算した
Ideal Average Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Ideal Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150
100

0

50

100

150

200

250

300

Distance [arcsec]

Ideal Average Method

400

6GS, 240"
8GS, 240"
10GS, 240"
Normal Method
Ideal Average Method
Requirement 200nm

WFE [nm]

350
300
250
200
150
100

0

50

100

150

200

250

300

Distance [arcsec]
図 5.8: 図 4.10 に 0m で予想されるフィッティングエラーを追加した図。四角が Nomarl Method、丸印が理想的な Average Method である。上の図は GS と重なる方向、下図は重ならない方向のプ
ロファイルである。

結果に 0m のフィッティングエラーを追加したプロファイルであるが、0m のフィッティングエラーを考慮しても Normal Method よりも十分に推定残差が小さくなっている。実際に平均値を用いた
とき (図 4.11) に Normal Method と Averge Method に違いが無かったのは、平均値を用いた地表層の推定の精度が影響していると考えられる。
実際に GS が８個の場合に Normal Method と Average Method の地表層成分の推定残差と、図60

第 5 章 議論

5.3. GLAO+MOAO の結果について

5.7 と同様に下層 (0m 4000m) と上層 (8000m,16000m) に分けた推定残差を図 5.9 に示す。図の四Average Method and Normal Method (8GS, 240")

600

Average Method

550

Normal Method

500

ALL

WFE [nm]4500m Layer

400

High Layers

350

Low Layers

300
250
200
150

0

50

100

150

200

Distance [arcsec]

250

300

図 5.9: 四角印が Average Method、丸印が Normal Method の結果であり、どちらとも GS が 8 個で外側の GS の距離が 240”で計算を行っている。赤は７層すべてでの推定残差、青は 0m の地表層の推定残差、緑は上層の推定残差、紫は下層の推定残差を示している。

角が Average Method、丸印が Normal Method の結果である。まず青色の地表層の推定残差を見てみると Normal Method の方がやや推定残差が小さいが、どちらの手法でも大きな残差が残っている。Normal Method の場合は下層の推定値は高さごとに分離ができていないのでそれぞれの高さでの推定残差は大きく見えることを Section 5.2 で見た。Average Method の場合も平均値をとったときに０ m の大気揺らぎのみではなく 500m、1000m の大気揺らぎも含まれており、0m での推定残差としては大きくなったと考えられる。次に上層 (緑)、下層 (紫) の推定残差を見るとどちらの手法もほぼ同じような結果になっている。これは Average Method で平均値をとって別で地表層を推定しても、上でも述べたように平均値には地表層以外の大気揺らぎも含まれていて、結局推定が分離できていないので Normal Method と同じような結果になったと考えられる。その結果最終的な残差 (赤) を見ても特に違いは見られなくなる。
今回の結果から地表層を別で推定してトモグラフィー計算から外す手法を用いることで、普通にすべての大気揺らぎでトモグラフィー計算をするよりも推定残差が小さくなることがわかったが、これは地表層をどれだけ精度良く推定できるかに依存し、平均値を用いた場合は地表層だけを別で推定することができず精度の向上は見られなかった。

5.3.2

GLAO の可能性

次に Normal Method と Average Method で GLAO 的な補正を行った場合の補正精度を検証する。これまでの結果から下層の推定値は高さごとに分解ができてないことがわかっているので、下３層 (0m、500m、1000m) を用いて補正を行った。500m と 1000m に関しては視野中心方向での推定値を用いる。図 5.10 がその結果である。今回は半径 10’ の視野内で評価した。赤が Nomal

Method、青が Average Method の結果であるが、２つの手法で最も異なる点は傾き成分の推定残差 (×印) である。Normal Method のほうが傾き成分の推定残差が大きくなっているが、これは７層のトモグラフィー計算の結果の下３層だけを用いていることが原因である。複数の層のトモグラ
フィー計算では全体に広がる傾き成分は高さごとで分解することはできず、７層を合わせて考慮する必要がある。そのため下３層のみの結果では、傾き成分に関しては正しい推定値とはいえないの61

第 5 章 議論

5.3. GLAO+MOAO の結果について

GLAO WFE (8GS, 240")

1100

Normal Method

1000
900

Average Method

WFE [nm]

800
700
600
500
400
300
200

0

100

200

300

400

Distance [arcsec]

500

600

図 5.10: 下 3 層の視野中心方向の推定値で補正を行った場合の結果。赤が Normal Method、青が Average
Method の結果である。×印が低次成分 (傾き) のみの推定残差、四角印が高次成分の推定残差、丸印が全成分の推定残差である。

で残差が大きくなっている。Average Method では 0m の推定値は全 GS の測定値の平均値を用いているため Normal Method よりは残差が小さくなっている。この影響で全体の残差 (丸印) も２つの手法で大きく差がでてしまっている。しかし、高次の成分の推定残差 (四角印) を見るとそこまで差は無いので、２つの手法による GLAO 的な補正精度に差はないと考えられる。また、傾きの成分を正しく推定するように手法を検討しなければない。どちらの手法も視野中心方向が最も精度が良くなる。これは中心方向のプロファイルを用いているためであると考えられる。

GLAO 的な補正で計算した PSF から IRAF の imexam を用いて測定した FWHM を図 5.11 に示す。赤は Normal Method、青は Averag Method を示している。誤差棒は視野内の最大値と最小
GLAO FWHM (8GS, 240")

0.45

Normal Method
Average Method
Seeing

0.4

FWHM [arcsec]

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

Jband

Hband

Kband

図 5.11: GLAO 的な補正を行ったときの FWHM。横軸は波長で赤は Normal Method、青は Average
Method を示している。三角印、または丸印は視野内での平均値、誤差棒は視野内での最大値と最小値を示している。黒四角は Seeing 時の PSF から測定した FWHM である。

値を表している。全体的に Average Method の方がいいのは傾き成分の残差が小さいためである。
どの波長でも Seeing よりも FWHM が小さくなっており、GLAO 的な補正が効果的であることが
わかる。

62

第 5 章 議論

5.3. GLAO+MOAO の結果について

Ensquared Energy 0.2"x0.2" [%]Ensquared Energy 0.2"x0.2" (8GS, 240", Jband)
60

Normal Method

50

Average Method

40

Seeing

30
20
10
0

0

100

200

300

Distance [arcsec]

400

500

600

Ensquared Energy 0.2"x0.2" [%]Ensquared Energy 0.2"x0.2" (8GS, 240", Hband)
60

Normal Method

50

Average Method

40

Seeing

30
20
10
0

0

100

200

300

Distance [arcsec]

400

500

600

Ensquared Energy 0.2"x0.2" [%]Ensquared Energy 0.2"x0.2" (8GS, 240", Kband)
60

Normal Method

50

Average Method

40

Seeing

30
20
10
0

0

100

200

300

Distance [arcsec]

400

500

600

図 5.12: 0,2”×0.2”の正方形で測定した Esquared Energy。上から J バンド、H バンド、K バンドの結果である。黒線は Seeing の PSF から測定した結果である。

63

第 5 章 議論

5.4. 風速を用いた結果について

最後に 0.2”×0.2”の正方形内の Ensquared Energy を図 5.12 に示す。全波長で中心方向が最も良くなっており、Seeing に比べて視野中心で約 2.5-3 倍、中心から 10”離れたところで約 2 倍のゲインがあることがわかる。

5.4

風速を用いた結果について

5.4.1

ノイズの影響を抑える手法

まず風速の推定誤差の影響を抑えるために、過去の測定値に対して重み付けを行った。式に両辺から次のような対角行列を左から掛ける。

1

..

.



Γ=






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

1
γ
..

.

(5.6)

γ
ここで γ は重みである。すると次にような形になる。
[
]
[
]
s(0)
G0 (0)=ΓφΓs(−∆t)
G0 (−∆t)
この式を式 (3.3) と同様に以下のような形にして解く。
[][]
ΓG0 (0, −∆t)Γs0 (0, −∆t)
Ḡ =
s̄ =
σL
0

(5.7)

ḠT Ḡφmv = ḠT s̄

(5.8)

最初に風速の推定誤差が無い場合でこの式を用いて計算を行い、重みによる振る舞いを調べた。
図 5.13 は 1s 前の情報を用いた Wind Method の場合に、重み付けを行った結果である。赤線がWeight Wind Method

400

Normal Method
γ2 = 1.00

350

γ2 = 0.75

WFE [nm]300

γ2 = 0.25250

γ2 = 0.10
γ2 = 0.01

200
150
100

0

50

100

150

200

Distance [arcsec]

250

300

図 5.13: 過去の測定値に対して重み付けした結果。ここでは風速の推定誤差は無い。

64

第 5 章 議論

5.4. 風速を用いた結果について

Normal Method で計算した場合で、その他が γ の値を変えて Wind Method を計算した結果である。γ を小さくしていくと過去の情報の重みが小さくなっていくため Normal Method の結果に近づいていく。
次に実際に風速の推定誤差がある場合に適用する。風速の推定誤差によって大気揺らぎの層が考えている位置よりも 0.5m ずれた場合に、γ の値を変化させて重み付けの影響を検証する。図 5.14

Weight Wind Method

400

Normal Method

350

No Noise2

γ = 1.00

WFE [nm]3002

γ = 0.502502

γ = 0.25
2

γ = 0.1200

γ2 = 0.01
150
100

0

50

100

150

200

Distance [arcsec]

250

300

図 5.14: 風速の推定誤差によるずれ 0.5m であるときの重み付けの結果。赤が Normal Method の結果、
青が風速の推定誤差がないときの W い nd Method の結果である。

がその結果の図である。赤が Normal Method の結果、青がノイズが無い場合の Wind Method の結果、他が風速の推定誤差がある場合に重み付け行った結果である。重み付けを行わないとき (緑、γ = 1) わないときでは、Normal Method よりも推定残差が大きかったが、重み付けを行うことで特に GS と GS の間 (中心から 150”) ほどで Normal Method に比べて推定残差が小さくなっている。γ の値を小さくしすぎると、Normal Method の結果に近づいていくが、γ の値をうまく調整
すれば風速の推定誤差が大きい場合でも Wind Method である程度の推定精度向上が期待できると
いえる。

5.4.2

手法の現実性最後にこの Wind Method が現実性について考えていく。今回の計算では大気揺らぎ自体は変化
せず、さらに風速の情報がわかっている理想的な場合での結果である。しかし実際にこの手法を用いるためには以下の２つのことを考慮しなければいけない。

1. 大気揺らぎの”frozen ﬂow”の仮定がどこまで正しいか。
2. 風速の推定手法
今回用いたような風速の情報を用いる手法は、具体的な手法は少し異なるが系外惑星の直接撮像をターゲットとした超高精度な補償光学や可視の補償光学でも検討されている (Johnson et al.

[15], Poyneer et al. [19])。このような補償光学では先に述べた Temporal Error を限界まで減らす必要がある。そこで風速の情報から WFS の測定と DM の補正の間の遅延時間中の大気揺らぎの65

第 5 章 議論

5.4. 風速を用いた結果について

時間変化を予測し、Temporal Error を小さくすることが考えられている。ここではこれらの先行研究の結果から、本節の冒頭であげた２点について考えていく。
まず frozen ﬂow の仮定がどれだけの時間保たれるかについて考えていく。frozen ﬂow の仮定が保たれるタイムスケールは、Wind Method においてどれだけ過去の測定値を用いることができる
かに大きく影響する。4 章で見たように Wind Method で推定精度を向上させるためには特に上層で現在の測定からずれた位置の測定結果を用いる必要がある。それを達成するためには風速が強くなるか、なるべく過去の測定値を用いるかどちらかである。前章の結果から今回の風速モデルの場
合、0.1s 前の測定値を用いることができれば十分な推定精度の向上が見られた。ここではこの 0.1s
を 1 つの基準として考えていく。まずはこのタイムスケールを理論的に概算していく。この概算は国立天文台ハワイ観測所の高遠氏の概算を参照した。大気揺らぎの時間変化は乱流運動により大きなスケールの揺らぎが小さなスケールの揺らぎに崩れていくことによって生じる。ある任意のス
ケール ` での典型的な乱流運動の速度を v` とすると、このスケールの大気揺らぎが完全に崩されるタイムスケールは

τ` ∼

`
v`

(5.9)

と記述されるであろう。より小さなスケールに崩されていく過程で、大きなスケールから小さなス
ケールへ伝達される単位時間当たり、単位質量当たりの乱流の運動エネルギー  は

 ∼ v`2 /τ` = v`3 /`(5.10)と表される。この値はスケールによらず一定であると考え、乱流のエネルギー入力が起こるスケー
ル L での典型的な速度を V とすると

V 3 /L = v`3 /`(5.11)( )1/3`V
L(5.12)となる。これよりv` =したがって、スケール ` の大気揺らぎが乱流によって崩され変化するタイムスケールは

τ` ∼`2/3 L1/3
V(5.13)と見積もることができる。ここで L はアウタースケール L0 に対応するので今回の場合 L = 30mとなる。さて、frozen ﬂow の定義を”ある時間の大気ゆらぎの構造の 50%以上が保たれている状態”とすれば、そのタイムスケールは τf (`) = τ` /2 と考えることができる。図 5.15 は横軸を揺らぎの空
間スケールとし、いくつかの風速の場合の frozen ﬂow が成り立つタイムスケールをプロットした図である。色の違いは考えている風速の違いである。今回考えている風速のモデルでは下層で 10m/s
前後、上層で約 25m/s の風速を持つ。また測定できる揺らぎのスケールは WFS の subaperture のサイズである 0.5m より大きいスケールである。まず下層では風速が 10m/s なので、図 5.15 から

subaperture サイズの揺らぎが frozen ﬂow として見なせるタイムスケールは 0.1s である。大気揺
らぎのパワースペクトルより、揺らぎのスケールが小さいほど揺らぎのパワーも小さくなるので、
subaperture サイズの大気揺らぎは測定できる揺らぎの中で最も揺らぎのパワーが小さい。より大きなパワーを持つ揺らぎの frozen ﬂow であるタイムスケールはさらに長いため、下層では Wind
Method に必要な 0.1s 内で十分 frozen ﬂow の仮定が成り立っていると言える。一方、上層では風速は 25m/s と大きいため subaperture サイズの揺らぎの frozen ﬂow であるタイムスケールは約40ms と非常に小さい。この場合 0.1s で frozen ﬂow として見なせるのは 2m 以上のスケールを持66

第 5 章 議論

5.4. 風速を用いた結果についてTime scale of Frozen Flow

WFS subaperturetime scale of frozen flow [s]

1
0.5
0.3
0.2
0.1

V = 5m/s
V = 10m/s
V = 15m/s
V = 20m/s
V = 25m/s
V = 30m/s

0.05

0.02
0.1

0.25

0.5

1

2.5

5

spatial scale of turbulance [m]1030

図 5.15: 揺らぎの空間スケールと対応する frozen flow が成り立つタイムスケール。色の違いは風速の違いである。

つ揺らぎであり、上層では Wind Method で 0.1s 前の測定値を用いた場合 2m 以下のスケールを持
つ揺らぎによる誤差が発生する可能性がある。風速がもっと小さければタイムスケールは大きくなるが、その代わりに Wind Method で推定精度を向上させるためには 0.1s より過去の測定値を用いる必要があるため、この 2m 以下のスケールの揺らぎによる影響は風速によらず生じることがと考えられる。
次に frozen ﬂow の仮定のタイムスケールを観測的な結果から検証していく。”frozen ﬂow”の存
在は Gendron and Lena [12] や Schock and Spillar [20] によって観測的に見つかっている。Schockand Spillar [20] では 3.5m と 1.5m の望遠鏡で WFS を用いて frozen ﬂow であるタイムスケールを見積もった。彼らはある時間に測定した波面とそれから少し時間が経った波面の相関をしらべることでどのくらいの波面の構造がどれだけの時間保たれているかを調べた。その結果、90%の構造が保たれている時間は約 20ms、50%の構造が保たれている時間は 60ms-100ms であると結論している。また、彼らは同時に大気揺らぎの層の数と風速も推定している。推定した大気揺らぎの層の数は約 5 層で、風速のプロファイルは 1 層だけ約 40m/s の強い風速を持ち、残りの層は 5-10m/sの風速を持つという結果であった。 Johnson et al. [15] では口径が 1m の望遠鏡と 5m の望遠鏡でWFS を用いて測定されたデータを用いたシミュレーションで大気揺らぎの 90%の構造保たれている時間は 20ms-130ms という値を求めた。これらの観測的な結果はすべての大気揺らぎの層を合計して考えた結果であり、さらに口径より小さなスケールの大気揺らぎをまとめた結果である。そのため先に見積もった理論的な概算と比較するのは難しいが、やはり 0.1s 前後が frozen ﬂow の仮定が成り立つ限界であると考えられる。また、frozen ﬂow のタイムスケールは観測する日によって大きくばらついていることもわかる。
次に風速の推定に関して考えていく。風速の推定は Johnson et al. [15] や Poyneer et al. [19] で行われている。Poyneer et al. [19] では Keck 望遠鏡や Gemini 望遠鏡の AO で取得された WFS のデータを用いて複数の層の風速の推定を行い、5m/s-30m/s の風速を検出し、median の値は 10m/sであった。彼らは風速の時間変化についても言及しており、彼らの結果では風速の時間変化は 10sの間に RMS で 0.5m/s であった。彼らはシミュレーションで彼らの風速の推定手法の精度を求めており、この風速の変化が RMS で 0.5m/s という結果が推定誤差ではないことを確認している。

67

第 5 章 議論

5.4. 風速を用いた結果について

この結果から風速の情報を使用する場合、その情報の更新頻度の要求はそれほど厳しくないことがわかる。また、彼らの手法の精度は風速が 10m/s の場合で推定誤差が 0.05m/s、風速が 20m/sの時は 0.1m/s という非常に高精度な推定が行えることが確認された。しかし、この手法では高さ方向の情報はわからないため、今回試した Wind Method 用の手法としては使用できない。その他に [2] では光の intensity の変化から風速を高さごとに推定する手法を数値シミュレーションによって検討しており、推定している高さは 1000m 付近と浅いが非常に精度よく高さごとに推定できている。４章の結果では風速の推定による位置のずれが 0.1s ほどであれば推定精度にそこまで影響しないという結果が出ている。0.1s 前の測定値を用いる場合は、風速の推定誤差にして 1m/s である。つまり風速の推定に必要な精度は Wind Method の場合はそこまで厳しくないので、精度に関してはここで紹介した先行研究の手法でも十分達成されている。問題は高さごとに推定する手法が必要なことである。その手法としてまずはトモグラフィー計算の推定結果から見積もる手法を検討していきたい。
以上のことをふまえて、まず frozen ﬂow のタイムスケールに関しては 0.1s の時間内では特に上層の小さなスケールの大気揺らぎが大きく変化する可能性がある。今回は frozen ﬂow の定義として”ある時間の大気揺らぎの構造の 50%が保たれているタイムスケール”としたが、大気揺らぎがどのくらい変化したら最終的な推定残差にどのくらい影響を与えるかを見積もってからこの定義を決めなければいけない。また、どのくらいのスケールの大気揺らぎまでが frozen ﬂow として保たれているべきかもしっかりと考慮していかなければいけない。風速の推定に関してはこれまでの手法で精度としては十分であるため、あとは高さごとに推定する手法を検討する必要がある。しかし、今回の計算結果から多少のノイズがのったとしても普通にトモグラフィー計算を行うよりはWind Method を用いることで推定精度の向上が期待される。

68

第 6 章 結論

第6章

結論

本論文では次世代超大型望遠鏡に向けた多天体補償光学のためのトモグラフィック波面再構成について検討を行った。多天体補償光学を実現するにあたり現状で大きな課題となっている高精度化、高速化を達成するために今回新たに開発した手法と検証結果は以下の通りである。

1. 各天体ごとに補正を行う前に、全方向で共通成分である地表層を GLAO 的に補正するシス
テムを検討した。補正を２段にすることで DM が補正する位相差や WFS の測定する位相差の大きさを小さくすることができると考えている。また、地表層の推定方法としてトモグラ
フィー計算とは別に、複数のガイド星の測定値の平均を用いる方法をシミュレーションを用いて検証した。理想的に地表層を取り除いた場合にはすべての大気揺らぎをトモグラフィー計
算から求めるよりも推定精度の向上が確認されたが、実際に平均値を用いて地表層を求めた
場合には推定精度に違いはなかった。これは平均値による地表層の推定精度がトモグラフィーを用いた場合とほとんど変わらないことが原因であり、この手法によって推定精度を向上させるためには高精度に地表層を取り除く必要がある。また、GLAO 的な補正を行った場合の補正精度についても検証を行い、直径 20’ の視野内で 0.2”×0.2”の正方形内の Ensquared

Energy が Seeing で制限される値より 2 3 倍向上することがわかった。
2. 風速の情報を用いて過去の測定値と現在の測定値を合わせて情報を増やし、計算の精度を向上させる手法を開発、検証した。その結果、frozen ﬂow の仮定と風速の情報がある程度の精度でわかっていればこの手法により大きな精度向上につながることがわかった。また、手法の現実性についても議論した。今回の風速のモデルで、この手法に必要な 0.1s の時間内では特に上層の小さなスケールの揺らぎに対して frozen ﬂow の仮定が成り立たない可能性があり、多少の誤差が生じる可能性がある。風速の推定に関しては高さごとに推定する手法が必要である。
いずれの手法でも理想的な場合では大きな推定精度の向上が見られ、その結果推定視野を直径 5’以上に広げることができた。しかし、目標である直径 10’ の視野内で波面残差 200nm 以下を達成
することはできなかった。
今後の展望としては以下のことを行っていきたい。

• TT-NGS の存在、WFS の測定誤差、DM の形状、光学系の収差などを考慮してシミュレー
ションを発展させ、より現実的な状況での推定の精度や推定視野について議論する。
• 目標とするサイエンスからの要求精度や、予想される成果をさらに正確に見積もる。
• 実際に複数の WFS で観測したデータを用いて、frozen ﬂow の過程の妥当性の検証や風速の推定手法の開発を行う。
• 今回開発した手法を実際に望遠鏡を用いて試験を行う。そのために現在始動している MOAOの試験装置開発のプロジェクトなどに積極的に参加する。

69

謝辞
まず本修士論文をはじめ、この３年間の勉学や研究に関して丁寧なご指導を頂いた担当教員の秋山正幸准教授に感謝いたします。また、たくさんの海外出張や発表の機会を与えてくださったことにも心より感謝いたします。来年度からは博士後期課程の学生としてお世話になります。まだまだ未熟な部分もあると思いますが、早く一人前になれるよう精進していきますので今後ともよろしくお願いいたします。
国立天文台ハワイ観測所の大屋真氏にはゼミや研究会を通じて多くのアドバイスや幅広い知識をいただきました。ここに感謝に意を表します。
また、この 1 年間頑張ってこれたのは天文学専攻の先輩、同期、後輩の皆様のおかげです。本当にありがとうございました。
最後に修士、博士への進学を快く受け入れ、応援してくれた両親に感謝いたします。

Appendix A AdaptiveOptics

Appendix A
AdaptiveOptics
Appendi A では補償光学の一般的な知識について説明する。

A.1

望遠鏡による観測

我々は望遠鏡を用いて、天体を観測する。まずはこの望遠鏡を用いることでどのように天体が観測されるのか、また大気揺らぎがどのような影響を与えるのかを簡単に説明する。

A.1.1

PSF と回折限界

望遠鏡によってどのように天体が観測されるかを説明する。我々が観測する星や銀河などの天体
は地球から離れているため、外部の収差の影響が無い場合、天体の光は平面波として望遠鏡に入射すると考えることができる。望遠鏡に入射した光は望遠鏡の光学系を通過し、焦点面上に結像される。非常に遠くの点光源を観測した場合、得られる像は無限に小さい点となるべきだが、実際は回折の影響で有限の広がりを持った像として観測される。このとき像は点広がり関数 (Point Spread

Function, PSF) で表すことができる。PSF は観測波長や光が受けた収差の影響や望遠鏡の口径で決まり、開口面での位相面をフーリエ変換することで得られる。収差の影響が無い理想的な場合、PSF は望遠鏡の口径と観測波長のみで決まる。図 A.1 に回折限界の PSF を載せる。このときのPSF の大きさ θD (中心から第一暗環までの距離) が望遠鏡の解像力、分解能の理論的な限界値となる。これを回折限界 (Diﬀraction Limit) と呼ぶ。回折限界は望遠鏡の性能を特徴づける重要な値であり、観測波長が λ[m]、望遠鏡の口径が D[m] のとき回折限界の大きさ θD [rad] は次のように近似的に求めることができる。

λ
[rad]
(A.1)
D
式 (A.1) から口径 D が大きいほど、また観測波長 λ が短いほど回折限界 θD が小さくなり、空間分解能が向上することがわかる。図 A.1 には口径が 30m と 8m の場合の回折限界像が載せてあるが、θD ∼ 1.22

30m の方がより小さい領域の光が結像されている。これがより大きな口径を持つ望遠鏡を作るメ
リットの一つである。図 A.2 は回折限界と観測波長、口径の直径の関係を示した図である。この
図より K バンドの場合、すばる望遠鏡 (口径 8.2m) の回折限界は ∼67[mas]、TMT(口径 30m) は∼18[mas] となり、TMT の方が約 3.5 倍の空間分解能を持つことがわかる。
しかし、地上観測の場合、PSF は光学系の収差や大気揺らぎの影響を受ける。図 A.1 の右側の
大気揺らぎの影響を受けたときの PSF を見てみると、短い露出時間の場合では焦点面で光は一点集まらず回折限界くらいの大きさに分散してしまう。これをスペックルノイズと呼ぶ。また、長時
間露出の場合、口径の大きさに関わらず大気揺らぎの強さで決まるサイズに星像が広がっているこ

71Appendix A AdaptiveOptics

図 A.1: 観測波長 λ = 2200[nm] の場合の PSFとがわかる。つまり、地上観測の場合どんなに口径を大きくしても、大気揺らぎの影響を補正しな
い限り大気揺らぎでリミットされる分解能しか得られない。大気揺らぎの影響を補正する補償光学
は地上の大型望遠鏡にとって必須のシステムであるといえる。

A.1.2

大気揺らぎの影響の定量的な評価

収差や大気揺らぎの影響を定量的に表すものとして光の位相の二乗平均平方根 (Root Mean Square,

RMS) を考える。これは平均値からのばらつきを表している。位相の平均値 σφ [rad] を持ったある連続的な波面 f (x, y) {x, y ∈ Ω} の RMS は√∫2{f (x, y) − σφ } dΩ
Ω
∫
RMS =
[rad]
(A.2)
dΩ
Ω
と書ける。また n 個の点で表せる離散的な波面を考えた場合は√∑n2i=1 (fi − σφ )
[rad]
RMS =
n

(A.3)

と書くことができる。位相の変化 φ[rad] は光学的な距離の変化 d[m] と同じであり、φ = 2πd/λ の関係でつながっている。同様に位相の RMS σφ と光学的な距離の RMSσd の関係は σφ = 2πσd /λとなる。
さらに収差の最終的な像への影響の評価方法としてストレル比 (Strehl Ratio, SR) を考える。SR
の定義は「回折限界の PSF のピークの値と位相が歪んだ時の PSF のピークの値の比」である。つまり、収差によってどれだけ PSF が分散してしまったかを表しており、回折限界像の PSF のピー
クの値を I0 、位相が乱れたときの PSF のピークの値を I とすると SR = I/I0 と書くことができる。SR は 0 ∼ 1 の値をとり、1 のときが理想的な場合である。光学的には SR≥ 0.8 のときはほぼ

72

Appendix A AdaptiveOptics

0.12

K-band (2200nm)
H-band (1650nm)
J-band (1200nm)
I-band (800nm)
V-band (550nm)

0.11
Diffraction Limit [arcsec]

0.1

Subaru 8.2m

0.09
0.08
0.07
0.06
0.05

TMT 30m

0.04
0.03
0.02
0.01
0

5

10

15

20
25
30
Diameter of Aperture [m]

35

40

45

図 A.2: 回折限界と口径、観測波長無収差であると言われている。また、SR は歪んだ光の波面の RMS 平均から以下のように近似的に計算することができる。

{ (
)2 }
( 2)
2πσd
SR ≈ exp −σφ = exp −
λ

(A.4)

は σφ ≤ 2[rad] の範囲で正しいと考えられている。図 A.3 は各波長でのストレル比と RMS の関係を示した図である。短波長ほど同じ RMS[nm] でも SR の値が小さいので、短波長の方が大気揺ら
1K-band (2200nm)
H-band (1650nm)
J-band (1200nm)
I-band (800nm)
V-band (550nm)

0.9
0.8

Strehl Ratio

0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

0

100

200

300

400

500

600

700

800

900

1000

RMS [nm]

図 A.3: 各波長でのストレル比

ぎの影響を受けやすいことがわかる。現在稼働している補償光学は赤外域での補償を行っており、条件が良い場合で K-バンド (2200nm) で SR∼0.6 ほどの精度が達成されている。

73

Appendix A AdaptiveOptics

大気揺らぎ

A.2

次に補償対象である大気揺らぎについて説明する。ここで我々が大気揺らぎと読んでいるのは大
気中の温度ムラ (つまり、屈折率ムラ) のことである。この後の補償光学について理解するために、さらに補償光学のデザインを決める上でも大気揺らぎについて理解することは重要である。

A.2.1

大気ゆらぎの空間スペクトル

大気揺らぎの空間周波数のパワースペクトルはコロモゴロフの乱流理論 (Kolmogorov’s model)でよく表せることがわかっている。コロモゴロフの理論では揺らぎのエネルギーは大きいスケール
から小さいスケールへと分配されていく。地球の大気の場合、日中の間に比較的大きなスケールの揺らぎが太陽光からエネルギーを受け取る。大きなスケールの大気揺らぎが持つエネルギーは、風によって混ぜられ徐々に小さなスケールの揺らぎに分割されていく。このような過程で地球大気には温度のムラが生じる。温度のムラは屈折率のムラと等価であり、この屈折率が不均一な大気中を天体の光が通過することで位相が変化し、天体の光は収差を持った状態で観測される。
コロモゴロフの理論に従うと、大気中の屈折率の１次元のパワースペクトル ΦN (κ) は ΦN (κ) ∝−5/3

κを満たす。3 次元の場合は Tatarski によって解析的にΦN (κ) ∝ κ−11/3(A.5)となることが示された Hardy [24]。
しかし実際には式 (A.5) で表すことのできる大気揺らぎの空間スケールには限界があり、その空
間スケールの上限値を outer scale L0 、下限値を inner scale `0 と言う。outer scale より大きなス
ケールの揺らぎは等しいパワースペクトルを持ち、innner scale より小さなスケールでは指数関数
的に値が小さくなるような振る舞いをすると言われている。一般的に inner scale 付近では大気揺
らぎのパワーは非常に小さく、inner scale の影響はほとんど無視できる。しかし、outer scale 付
近では揺らぎは大きなパワーも持つため、outer scale の効果を考慮したスペクトルが提案された(Strohbehn [21])。
ΦN (κ) ∝ (κ2 + κ20 )−11/6(A.6)

である。これを Karman Spectrum と呼ぶ。ここで κ0 = 2π/L0 は outer scale によって決まる空間周波数である。図 A.4 は修正されたパワースペクトルと Kolmogorov のパワースペクトルを比較した図である。観測的には inner scale は数 mm と見積もられている。一方、outer scale は数 mから 100m 以上と大きく変化することがわかっている。outer scale の大きさは大気揺らぎの強さに大きく影響する。
また Hardy [24] より、位相の変化のパワースペクトルが求められている。コロモゴロフスペク
トルの場合と Karman スペクトルの場合でそれぞれ以下のように書くことができる。

Φφ (κ) =Φφ (κ) =0.023 −11/3
κ
5/3
r0
0.023
5/3

r0

(κ2 + κ20 )−11/6

74

(A.7)

(A.8)

Log of normalized power of reflactive index

Appendix A AdaptiveOptics

10

5

1

spatial scale [m]
0.1

0.01

0.001

Kolmogorov Power Spectrum
Modified Power Spectrum

0-5

Inner Scale = 2mm-10Outer Scale = 10m

-15
-20

-1

0

1
2
Log of spatial wave number [1/m]34

図 A.4: 修正された Spectrum と Kolmogrov Spectrum。縦軸は規格化した屈折率のパワー、横軸は空間周波数。この図では L0 = 10[m]、`0 = 2[mm] で計算している。

A.2.2

大気ゆらぎの強度

2
大気揺らぎの強さは大気構造定数 CN(h) で表される。大気構造定数は任意の２点の屈折率の違2いを空間平均することで計算することができる。また、CNは海抜からの高度の関数である。通常、
2
地球大気は地上から 1km 付近に最も強い揺らぎの層を持つ。ここでの CNの平均値は地上付近で10−15 ∼ 10−14 m−2/3 、1km 付近で 10−16 m−2/3 であると観測されている。現在稼働している大型
望遠鏡はこの層よりも高い山の山頂などに建設されているため、この層の影響は無視できる。こ2の層よりも上では 8km∼10km くらいまで CN= 10−17 m−2/3 の強さの大気揺らぎが存在する。この高度では大気ゆらぎは複数の薄い層上に発生することが観測されている。Azouit and Vernin [1]での観測結果を以下に示す。さらに 10km 付近には対流圏界面 (対流圏と成層圏の界面) が存在し、ここでは急速な風向、風速の変動である wind shear の影響でしばしば強い大気揺らぎが発生することがある。また、10km 以上では大気揺らぎの強さは急激に減衰する。そして 20km∼25km よりも上では大気揺らぎの影響はほとんどないと考えることができる。

A.2.3

大気揺らぎの波長依存性

大気の屈折率のゆらぎは可視域から赤外域あたりまでは変化しないと言われている。そのためこの波長域では測定される波面は波長によらず共通であり、可視域で測定を行い、赤外域で補正を行うことが可能となる。

A.2.4

風速のモデル

風は大気揺らぎの時間変化に大きく寄与する。また、対流圏界面では wind shear より強い大気
揺らぎが発生するなど大気揺らぎの強さとも関係している。気球による測定により、風の強さは対流圏界面付近で最も大きな値を持ち、その高度を中心としたガウシアン的な高さプロファイルを持つことが示唆されている。Azouit and Vernin [1] での観測結果を以下に示す。また、一般的なモデ

75

Appendix A AdaptiveOptics2
図 A.5: 気球を用いた大気構造定数 CN(h) の観測結果 (Azouit and Vernin [1])。薄い層状になっている。

ルとして次のような式で表される (Hardy [24])。
[ (
)2 ]
[
]
z cos(ζ) − HT× sin2 φ + cos2 φ cos2 ζv(z) = vG + vT exp −
LT

(A.9)

ここで、vG は低い層での風速、vT は対流圏界面付近での風速、ζ は観測している高度角、HT は
対流圏界面の高さ、LT は対流圏界面の厚さ、φ は望遠鏡が向いている方向と相対的な風の角度で
ある。

A.2.5

大気ゆらぎのパラメーター

次に大気揺らぎの状態を表すいくつかのパラメーターを紹介する。これらのパラメータは補償光学のデザインを決める上でも非常に重要なものである。
フリード パラメーター (Fried Parameter) r0
フリードパラメーター r0 は大気揺らぎの強さを表す指標として用いられる。r0 の次元は長さであり、波面の RMS が 1[rad] となるような口径の直径を表している。つまり、r0 の値が2大きいほど大気揺らぎの強さは弱い。r0 は次のような式で CNから計算できる。

[
]−3/5
∫
2
2
r0 = 0.423k (sec ζ) dh CN (h)(A.10)

ここで、k は波数 k = 2π/λ、ζ は観測している高度角である。この式から r0 は波長に対して λ6/5 の依存性を持つことがわかる。一般的に r0 の値は λ = 500[nm] での値が用いられ、海抜付近での値は数 cm 以下、マウナケア山頂などの条件のいい場所では約 15cm くらいである。また、λ =500nm で r0 =15cm のとき K バンド (2200[nm]) では約 90cm となる。

76

Appendix A AdaptiveOptics

図 A.6: 気球を用いた風速の観測結果 (Azouit and Vernin [1])。ガウシアン的なプロファイルになっている。

シーイング (Seeing) θseeing
大気揺らぎの天体観測への影響の大きさを表すものとしてシーイングが使用される。シーイ
ングは大気揺らぎの影響による星像の広がりの大きさを表している。図??の右側の PSF のサイズに対応する。シーイングは θseeing ∼ λ/r0 で計算される。望遠鏡の口径が r0 よりも小さい場合は、PSF のサイズは望遠鏡の回折限界で決まる。しかし、望遠鏡の口径が r0 よりも大きい場合、口径のサイズに関わらずシーイングほどの分解能しか得られないことがわかる。マウナケア山頂でのシーイングは r0 = 0.156m(@500nm) から計算すると約 0.6” である。

Isoplanatic Angle θ0
isoplanatic angle θ0 は大気揺らぎが方向によってどれだけ異なるかを表すパラメーターである。θ0 は 2 方向の大気揺らぎの差が RMS で 1[rad] となるような角度を表しており、次の式
で書ける。

]−3/5
[
∫
2
5/3
2
8/3
dh CN (h)hθ0 = 2.914k (sec ζ)(A.11)

式からわかるように θ0 に波長依存性があり、短波長ほど θ0 が小さくなる。マウナケア山頂
では可視域でだいたい θ0 ∼ 数秒角、赤外域で θ0 ∼ 数十秒角であり可視域での大気揺らぎの方向による違いが大きいことがわかる。

Tilt Isoplanatic Angle θT A
これは大気揺らぎのうち傾き成分についての isoplanatic angle である。大気揺らぎの傾き成分は口径に比べて大きいスケールの揺らぎの影響であるので、θ0 に比べて θT A は大きい。2方向の角度が小さい時 (口径 D に対して D/40, 000 以下の角度)、θT A は以下のような式で

77

Appendix A AdaptiveOptics

書ける。

[
2

θT A = 0.668k D

−1/3

]−1/2

∫
2
dhCN
(h)h2

(A.12)

Greenwood Frequency fG
これは大気の時間変化の周波数を特徴づけるパラメーターであり、以下の式で表される。

[∫
2fG = 0.102k sec(ζ)

]3/5
2
dhCN
(h)v 5/3 (h)(A.13)

補償光学

A.3

補償光学は大気揺らぎの影響を補正するために、大気揺らぎの測定、推定、補償の３つの大きな過程で構成されている。
大気揺らぎの測定
まず大気揺らぎを補正するために最終的に補正を適用するターゲットの近くにあるガイド星(Guide Star, GS) と呼ばれる明るい星を観測する。ガイド星の光は望遠鏡に届くまでに大気
揺らぎの影響を受けるので、その位相を波面センサー (WaveFront Sensor, WFS) で測定することで、ガイド星の光が通ってきた経路の大気揺らぎの情報を得ることができる。
大気揺らぎの推定
波面センサーで測定した情報から、補償すべき大気揺らぎを推定する。一般的に推定は連立方程式を行列の形式で解くことに帰着する。
大気揺らぎの補償
補償は表面の形を変形させることのできる可変形鏡 (Deformable Mirror, DM) で行う。反射
したときに推定した大気揺らぎをキャンセルするような形に DM の表面を変形させ、補正したい天体をこの DM を通して観測することで大気揺らぎの影響を補正する。傾き成分は変化が大きいので別の鏡で補償する場合が多い。
これらの行程を大気揺らぎの時間変化よりも十分早く繰り返すことで常に大気揺らぎの影響をリア
ルタイムで補正し続ける。この節では補償光学を構成する各要素について簡単に紹介する。

A.3.1

ガイド星 (Guide Star)

精度良く大気揺らぎを測定するためには、ガイド星はできるだけ明るい必要がある。また、ガイド星と補正したい天体がなるべく近くにあることが望ましい。ガイド星としては自然の天体を用いる Natural Gaude Star (NGS) とレーザーで人工的に作る Laser Guide Star (LGS) がある。
最も簡単な場合は補正したい天体が明るい星である場合である。この場合ターゲット天体自身を

NGS として扱え、測定方向と補正方向が一致するため最も最適な補正が行える。ターゲット天体が暗かったり、広がっている場合はそれ自身をガイド星としては使用できないため、別に近くの明
るい星をガイド星として扱う必要がある。このときターゲット天体とガイド星の距離が離れるほ
ど測定方向と推定方向が異なっていくので、最終的な補償精度は悪くなっていく。これを Angular78

Appendix A AdaptiveOptics

Anisoplanatism と呼び、この影響によって発生する平均的な誤差は解析的に以下のように表すことができる。
( )5/3
θ(A.14)σθ2 =θ0θ はガイド星とターゲットの角度距離、θ0 は isoplanatic angle である。この式からガイド星とター
ゲット天体の距離が θ0 だけ離れた時に σθ =1[rad] の誤差が発生する。これは SR=0.37 に相当する。A.2 述べたように観測条件の良いマウナケア山頂でも θ0 は可視域で数秒角、近赤外域で数十秒角であるので、一つのガイド星の情報で精度良く補正できる範囲も可視域で数秒角、近赤外域で数十秒角程度であることがわかる。図 A.7 と図 A.8 はすばる望遠鏡の AO188 での NGS の距離、
明るさと性能の関係を表した図である。図 A.8 からもガイド星の距離が離れるにつれて性能が悪く

図 A.8: すばる望遠鏡の AO188 での
NGS の距離と SR の関係。
([25])

図 A.7: すばる望遠鏡の AO188 での
NGS の明ると SR の関係。
([25])

なっていることがわかる。また isoplanatic angle はおよそ 10[arcsec]∼25[arcsec] であることがわかる。この図では測定誤差などの他の影響も含まれているので実際はもう少し良いかもしれない。
さらに図 A.7 では NGS が暗くなるにつれて補正精度が悪くなっている。このように精度良く補正を行うためには近くに明るい星が必要だが、そのような天体は限られているため NGS を用いた補償光学を適用できる領域は非常に限られてくる。図 A.9 は SR と銀緯と sky coverage の関係を示した図である。明るい星が多い銀河中心でも SR=0.3 を達成できるような領域は 1%以下となっており、非常に限られた領域でしか NGS による補償光学が使えないことがわかる。
この問題を解決するために明るい星が無い領域では LGS を用いて補償光学を行う。LGS としては主にレイリー散乱による LGS とナトリウムの励起を利用したものが使用されている。

Rayleigh LGS
レイリー散乱の散乱断面積は光の波長が長くなるほど小さくなっていくので、一般的に紫外
線や可視光のレーザーが用いられる。またレイリー散乱の強さは大気密度に依存する。高度が高くなるほど大気が薄くなっていくので、レイリー散乱による LGS の高さは 20km∼30km
くらいである。

Sodium LGS (Na LGS)
地球大気には約 92km 付近にナトリウムの層があることがわかっている。この層に 689nm のレーザーを照射し、ナトリウムを励起させることで人工的にガイド星を作ることができる。

79

Appendix A AdaptiveOptics

図 A.9: 各 SR を達成するガイド星の sky coverage の割合。横軸は銀緯。([28])ナトリウムの励起による光の放出は、レイリー散乱に比べて効率がいい。また、LGS の高度
もレイリー散乱による LGS に比べて 92km と非常に高い。
このように人工的にガイド星を作ることで補償光学を適用できる領域を大幅に増やすことが可能となった。現在稼働している大型望遠鏡では高度が高い Na LGS が採用されている。図 A.10 はマウ
ナケア山頂で LGS のためのレーザーが打ち上げられている様子である。また図 A.11 は Na LGS

図 A.10: ハワイ、マウナケア山頂でレーザーガイド星を使用する様子。すばる望遠鏡 (左端)、Keck 望
遠鏡 (左から２番目)、Gemini North 望遠鏡 (右端) からレーザーが出ているのが見える。ここでレーザー光が見えるのはレイリー散乱の影響である。([25])を用いた場合の補償光学の模式図である。
しかし、LGS を使用する場合に注意しなければいけないことがある。ほぼ無限遠にある NGS と違い LGS は望遠鏡から比較的に近い高さにある。そのため図 A.11 でも見られるように LGS からの光はコーン上に望遠鏡に届いてしまう。そのため LGS からの光では測定できない大気揺らぎがあるので補償精度を悪化させる原因となる。この効果は Focal Anisoplanatism と呼ばれ、その影響によって生じる誤差は解析的には以下のように記述される。

(
σF2 A =

Dd0

)5/3
(A.15)

ここで D は望遠鏡の口径、d0 は大気の状態、観測波長、観測高度、LGS の高さなどのに依存する値である。マウナケア山頂で Na LGS を用いた場合、d0 の値は K バンドで 20∼25m くらいである。
この式から現在マウナケア山頂にある 8m 級の望遠鏡では、K バンドでの Focal Anisoplanatism

80Appendix A AdaptiveOptics

図 A.11: Na LGS+TT-NGS を用いた補償光学。

の影響は σF2 A = 0.15 ∼ 0.2[rad2 ] となる。これは SR がだいたい 0.8 くらいに相当するので、あまり影響が大きくないと言える。しかし口径が 30m である TMT の場合は Focal Anisoplanatism の影響のみで SR が 0.3 以下まで落ちてしまうため、大きな誤差の原因となる。
もう一つの注意点として、LGS はそれ自身が大気揺らぎの影響で動いてしまうことである。望遠
鏡から打ち上げられたレーザーは大気中を通るため、LGS の高さに到達する前に大気揺らぎの影響
で進行方向が曲がってしまう。さらに上空で散乱したあと再び曲がったパスで望遠鏡に戻ってくる
ため、LGS 自身で口径全体に広がるような傾き成分 (Tip/Tilt) を測定できない。そのため、LGSを用いる場合は別で傾き成分を測定する NGS が必要となる。この NGS のことを Tip Tilt NGS、

TT-NGS と呼ぶ。TT-NGS は大気の傾き成分を測定するだけなので NGS よりも暗い星を使用することができる。もちろん TT-NGS は LGS やターゲット天体と近いことが要求されるが、上でも述べたように大気揺らぎの傾き成分は高次の成分よりも比較的広い範囲で一定であるので、NGSと比べて距離が離れることによって生じる誤差は小さい。距離が離れることによる傾き成分が異なることを Tilt Anisoplanatism と呼び、その影響は次のような式で書くことができる。
()2
θ2
σT A =
θT A(A.16)

ここで θ は TT-NGS とターゲット天体の角度距離、 θT A は Tilt isoplanatic angle である。注意しな
ければいけないのは、この式は TT-NGS とターゲット天体の距離が比較的近い場合 (θ ≤ D/40, 000)
のみ成り立つ。図 A.12 と図 A.12 は TT-NGS の距離、明るさと性能の関係を示した図である。この二つの図を図 A.7 と図 A.8 と比較した時、TT-NGS の方が明るさ、距離共に NGS と比べて制限が緩いことがわかる。そのため TT-NGS の存在を考慮しても、LGS を用いた方が sky coverageが大きくなる。

A.3.2

波面センサー (WsveFront Sensor)

ガイド星の光の波面を測定する波面センサーにはいくつかの種類があるが、ここでは最も使用されている波面センサーの一つであるシャックハルトマン波面センサー (Shack-Hartmann WFS、

81

Appendix A AdaptiveOptics

図 A.12: TT-NGS の明るさとストレル比の関係 ([25])図 A.13: TT-NGS の距離とストレル比の関係 ([25])SH-WFS) について説明する。
SH-WFS は大気揺らぎの測定だけでなく、光学系の歪みを測定する装置としてもよく用いられる。検出器の前にマイクロレンズアレイ (Micro Lens Array, MLA) を置くことによって、入射光を分割し局所的な波面の傾きを測定することができる。図 A.14 は SH-WFS の原理を簡単の示した図である。まず左側の図のように、SH-WFS に理想的な平面波が入射した場合は MLA によっ

図 A.14: Shack Hartmann WFS の原理て分割され、それぞれのレンズ (ここでは subaperture と呼ぶ) で検出器上に結像される。このときの像の位置は各 subaperture の位置と対応する。次に右図のように位相が乱れた波面が入射した場合を考える。入射した波面は subaperture で分割され像を作るが、像の位置は分割された波面の傾きに対応して平面波の場合の位置から移動する。この各 subaperture で作られた像の位置の移動距離を測定することで、局所的な波面の傾きの情報を得る。この情報から波面を計算することで大気揺らぎの再構成が完了する。
現在使用されている WFS は可視域用のものがほとんどである。検出器としては CCD が最も用いられている。理由としては LGS が可視域で光ること、さらに、一般的に星は可視光で明るいた
め可視域の方が NGS として利用できる星が多いことがあげられる。大気の屈折率のゆらぎの波長依存性は小さいため測定を可視で行い、補正を赤外で行ってもよい。実際、現在使われている補償
光学はそのようなシステムである。

82

Appendix A AdaptiveOptics

SH-WFS で測定できる空間周波数成分の限界値は subaperture のサイズ、つまり MLA のピッチによって決まる。subaperture よりも小さな揺らぎは測定することができないので、より細かな大気揺らぎを測定する場合はより細かいピッチを持った MLA を使用しなければ行けない。しかし、

subaperture を小さくするほど各 subaperture に入ってくる光の量が減るため、作られる像も暗くなり像の位置の測定精度が悪くなる。特に補償光学の場合はリアルタイムで補償を行うため、波面センサーの露出時間は非常に短く、subaperture のサイズが光の量を決める。subaperture のサイズはこの二つのバランスを考えて決める必要がある。また光量を稼ぐために各 sabaperture は 2×2のピクセルが割り当てられる quad cell であることが多い。

TT-NGS 用の波面センサーは傾き成分のみを測定できればいいので subaperture は一つで十分である。この場合、subaperture の大きさは口径と等しくなるため高次の成分を測定する WFS に比べて光量が大幅に増える。TT-NGS は NGS に比べて暗い星を利用できる。
次に各 subaperture で生じる測定ノイズを考える。測定誤差としてはフォトンノイズによるラン
ダムな誤差と光学系のアライメントや検出器の特性による系統的な誤差がある。ここでは前者について考える。quad cell を用いた点光源の位置測定時の誤差については先行研究で調べられている。

aperture のサイズが d[m]、観測波長が λ[m]、ガイド星の角度サイズが θ[arcsec] であるとき、生じる一方向 (x,y 方向のどちらか) の位置測定の誤差は


[( )
]1/2
2

3λ
π

2


+θ

 8SN R
2d
σα =
[(
]1/2
)2



π
3λ

2

+θ

8SN R
2r0

(r0 > d)[rad]

(A.17)(r0 < d)

この式はそのまま各 subaperture で生じる測定誤差としてみなすことができる。ここで SNR(Signal-

to-Noise Ratio) は各 subaperture での天体からの光とノイズの比であり以下のように記述できる。
SN R = √Nobj[]
2 + e2
Nobj + Npix Nbg(A.18)

Nobj はガイド星からの光子数、Nbg は subaperture のピクセル数 (quad cell の場合 Npix =4)、Nbg
は天体以外の背景光からの光子数、e は各ピクセルの Readout ノイズである。これらの式より、ガイド星が明るいほど測定ノイズは小さく、また、ガイド星の像が小さいほど測定ノイズが小さいことがわかる。

A.3.3

波面再構成

波面センサーのデータから DM の変形量を決めるために、大気揺らぎを再構成する必要がある。
一般的には行列とベクトルの演算やフーリエ変換によって計算される。詳細は本論文中で説明するが、リアルタイム補正のために高速で計算する必要がある。

A.3.4

可変形鏡 (Deformable Mirror)

補償光学では波面の位相の乱れを補正するために可変形鏡 (Deformable Mirror, DM) と呼ばれる表面の形を変形させることのできる鏡を用いる。DM の表面を反射したときに光の乱れた位相を

83Appendix A AdaptiveOptics

平行光に戻すような形に変形させ、この DM を通して観測を行うことで補正を行う。補償光学で用いられる DM としては傾きを補正する Tip-Tilt 鏡とそれよりも高次の揺らぎを補正する DM(一般的にはこちらを DM と呼ぶ) がある。ここでは後者について説明する。DM の稼働方式はいくつかあるが主に天文学で用いられているのは以下の２つである。
ピエゾ (Piezo)
最もよく用いられているものである。この方式では電圧をかけると変形する圧電素子 (ピエゾ素子) の特徴を用いて鏡を変形させる。この方式のメリットとしては得られる変形量が他の方式に比べて大きいということである。逆でメリットとしては装置自体が比較的大きくなってしまうことがあげられる。

Micro Electro Mechanical Systems(MEMS)
MEMS 技術を用いた方式で主に静電力によって鏡を変形する。特徴としてはピエゾに比べて用意に小型化が可能である。また、ピエゾの比べて安価であることもメリットである。デメリットとしては得られる変形量が小さいことである。
本研究で検討する多天体補償光学では補正する天体の数だけ DM を用意しなければ、安価で小型な MEMS-DM を用いることが検討されている。しかし、現在の MEMS-DM で達成されている変形量では不十分であるため、大きなストロークを持つ MEMS-DM の開発が必要である。また、近年では望遠鏡の副鏡を可変形鏡にする動きもある。こうすることで光学系内の鏡の枚数を減らし、光学系からの熱放射や光学系による光の損失を抑えることができる。
可変形鏡で考えなければいけないことは素子数、変形量、応答速度、用いることができる波長域、鏡表面の滑らかさである。ここでは素子数について考えていく。DM の 1 素子の大きさ以下のスケールの揺らぎは補正できないため、高精度な補償を達成するためには素子のサイズは小さい(つまり素子数が大きい) ほうが望まれる。しかし、実際には WFS の subaperture の大きさや要求精度などとのバランスを考えて決めていかなければならない。この補正できない小さな揺らぎによって生じる誤差はフィッティングエラーと呼ばれ、Hardy [24] よりコロモゴロフのパワースペク
トルを仮定すれば以下のような式で表すことができる。

(
σf2 it = aF

d
r0

)5/3
(A.19)

ここで d は可変形鏡の素子のサイズである。この式より素子のサイズを固定した場合、短波長ほど

r0 が小さいので生じる誤差が大きいことがわかる。また、aF は DM のタイプによって決まる定数
であり、鏡が連続的か分割鏡か、隣り合う素子同士の影響などから決まる。
実際にすばる望遠鏡の補償光学の観測波長の違いによる性能を表したのが図 A.15 である。この
図から短波長ほど SR が下がっていることがわかる。これは短波長ほど DM で補正できない小さな揺らぎの影響が大きくなることが１つの原因である。

A.3.5

制御

大気揺らぎは時間とともに変化していくので、補正を繰り返すことが必要となる。その制御方法には大きく分けて Closed Loop 制御と Open Loop 制御の２種類がある。

Closed Loop 制御

84

Appendix A AdaptiveOptics

図 A.15: すばる望遠鏡の補償光学 AO188 の性能。観測波長による SR の違い。([25])。

図 1.1 と図 A.11 は Closed Loop 制御になる。Closed Loop 制御では WFS の前に DM を置
き、補正後の波面の残差を測定する。その残差をキャンセルするように、一つ前の DM の変形量にフィードバックをかけて補正をアップデートする。この制御方法のメリットとしては、
まず常に補正後の残差を測定するので波面センサーで測定可能な位相差の範囲が小さくてよ
いことである。また一般的にフィードバックをかける際にゲイン (0 ∼ 1 の値) を掛けて少しずつ補正をしていくため、Open Loop に比べて測定・推定精度が悪くても、ループを高速に繰り返すことで十分な精度を保つことができる。デメリットとしてはゲインを掛けて少しず
つ補正を行っていくので、Open Loop に比べて高速にループを回す必要がある。また、ゲインの値を間違えると値が発散してしまい間違った補正を行うことがある。

Open Loop 制御
Open Loop 制御では DM で補正前の生データを測定し、大気揺らぎの絶対値を推定する。その推定値の逆位相がほぼ DM の変形量と等しくなる。1 回の測定・推定・補正の結果が最終的な精度の直接反映されるため、Closed Loop に比べて高精度な測定、推定、補正が必要となる。
現在稼働しているシステムでは主に Closed Loop が採用されている。理由としては測定や推定の精度に対する要求が低いことがあげられる。

A.3.6

補償光学の現状

現在世界中のほぼすべての 8m∼10m 級の大型望遠鏡には一つの NGS、もしくは LGS ＋ TT-NGSを用いた補償光学が搭載されている。達成されている精度としては K バンド観測で NGS 用いた場
合 SR∼0.6、LGS を用いた場合 SR∼0.4 ほどである。特に補償光学を用いた高分解能と面分光装置を組み合わせた銀河の内部構造の探査や、コロナグロフを用いた系外惑星の直接撮像などで大きな成果をあげている。
近年では、現在建設が予定されている TMT をはじめとした次世代超大型地上望遠鏡に向けた新しい補償光学の検討や開発が活発に行われている。

85

Appendix B Mathematics

Appendix B
Mathematics
B.1

行列微分を用いた計算の証明

ここでは次の式を証明していく。

∂ trace(XAX T )
= 2XA
∂X

(B.1)

∂ trace(BX T )
=B
∂X

(B.2)

∂ trace(XB T )
=B
(B.3)
∂X
ここで X を m × n の行列、A を n × n の対称行列、B は m × n の行列である。各行列の (i, j) 成
分を xij 、aij 、bij とする。
まず、行列微分について説明する。φ(X) を X を含むスカラーとすると φ(X) の X での微分は以下のように定義される。



∂φ

∂φ

···


 ∂x11∂φ(X) 
.
=
 ..
∂X

 ∂φ

..


∂x1n 
.. 

. 

∂φ 

.

···

∂xm1



∂xmn

次に XAX T の (i, j) 成分は

(

n ∑
n
∑
)
XAX T ij =xik akt xjt
k=1 t=1Tと書けるので、trace(XAX ) はtrace(XAX T ) =

n ∑
n ∑
m
∑

xpk akt xpt

k=1 t=1 p=1

これより、

[

( n n m)]∑∑∑∂∂ trace(XAX T )=xpk akt xpt
∂X
∂xij
ij
t=1 p=1
k=1

=

n
∑

(xik akj + xik ajk )

k=1
n
∑

=2

k=1

86

xik akj(B.4)

Appendix B Mathematicsとなる。最後の式変形では A が対称行列であることから、ajk = akj を用いた。この結果は (2XAT )ijに等しいので、すべての要素について考えると式 (B.1) が証明された。
同様に trace(BX T ) はtrace(BX T ) =

m ∑
m
∑

bkt xkt

k=1 t=1となるので、

[

(m n
)
]
∑∑
∂
∂ trace(BX T )=bkt xkt
∂X
∂xij
ij
t=1
k=1

= bij
となる。よって式 (B.2) が証明された。また trace(BX T ) = trace(XB T ) であるから式 (B.3) も
成り立つ。

B.2

共役勾配法 (Conjugate Gradient Method)

本節では共役勾配法について説明する。

B.2.1

共役勾配法

次のような連立一次方程式を考える。

Ax = b(B.5)

A は n × n の正定値対称行列、x と b は n の要素をもつベクトルである。この式の解を求めることは、変数 x の行列関数 φ(x)φ(x) =1 T
x Ax − xT b
2

(B.6)

を最小化する x がを求めることと同義である (Kawakami [29])。共役勾配法はこの式 (B.6) の最小値を効率的に求め、それによって式 (B.5) を解く手法である。初期値 x0 から出発して反復的に最小値を求めることを考える。k + 1 回目の修正における x を

xk+1 = xk + αk pk(B.7)とする。ここで pk は修正ベクトルであり、αk は修正の大きさを決める値である。αk は φ(xk+1 )を最小にするように、つまり∂φ(xk+1 )=0∂αkを満たすように決める。この式を解くと α は
αk =pTk (b − Axk )
pT rk
= Tk
T
pk Apkpk Apk(B.8)と書ける。ここで rk = (b − Axk ) は k 回目の修正の残差とする。共役勾配法では探索方向 pk が、

k − 1 回目までのすべての残差 ri (i = 0, . . . , k − 1) の一次結合で表せるように選ぶ。Kawakami [29]より、{
p0 = r0 = b − Ax0(B.9)
pk = rk + βk−1 pk−1
87Appendix B Mathematicsとなる。また、

βk−1 = −rkT Apk−1
pTk−1 Apk−1(B.10)である。k + 1 回目の残差 rk+1 は式 (B.8) より

rk+1 = b − Axk+1 = b − A(xk + αk pk ) = rk − αk Apk(B.11)

さらに式 (B.9) より、

αk =pTk rk(rk + βk−1 pk−1 ) rk
|rk |2==
pTk ApkpTk ApkpTk Apk(B.12)となる。式 (B.7)、式 (B.9)、式 (B.10)、式 (B.11)、式 (B.12) を用いて解を反復的に求めていく。
収束条件としては相対残差ノルム |rk |/|b| がある値よりも小さくなったときに収束したとする。

B.2.2

共役勾配法の実装

さて、ここからは実際に共役勾配法を実装していくことを考える。まず初期値 x0 はなるべく真の解に近いものが望ましい。しかし、真の解がわからない場合の方が多い。そのときは x0 はすべて 0 のベクトルを用いる。プログラムとしては以下のようになる。

0. B = |b|
1. t = Ax0
2. r0 = b − t
3. p0 = r0
4. R0 = |r0 |5. 収束条件の確認 R0 /B
条件を満たしていれば x0 が解で計算終了

for k = 1, · · ·6.

qk−1 = Apk−1

7.

ek = pk−1 · qk−18.

αk−1 = Rk−1 /ek

9.

xk = xk−1 + αk−1 pk−1

10.

rk = rk−1 − αk−1 qk−1

11.

Rk = |rk |12.

収束条件の確認

Rk /B

条件を満たしていれば xk が解で計算終了

13.

ck = rk · qk−1

14.

βk = −ck /ek15.

pk = rk + βk pk−1
end for88

Appendix B Mathematics

ここで t と qk は作業用のベクトル、Rk は k 回目の残差ノルム、B は b のノルム、ek と ck は作業用の定数である。

B.3

前処理

ここでは共役勾配法に用いる前処理について説明する。共役勾配法では係数行列 (式 (B.5) の A)に対して事前に処理をすることで収束速度を向上させることができることがしられている。この操作を前処理といい、前処理を用いた共役勾配法を一般的に前処理付き共役勾配法と言う。

B.3.1

対角スケーリング

これは本研究でも用いている手法であり、もっとも簡単な前処理である。次のような前処理行列

P を作用させ、係数行列 A の対角成分をすべて 1 にすることで収束速度が向上する。
 √1/ a11


P =




√
1/ a22..

.

√
1/ ann(B.13)

Ax = bP AP P −1 x = P b
Ãx̃ = b̃

(B.14)

ここで Ã = P AP 、x̃ = P −1 x、b̃ = P b である。また、P −1 は P の対角成分の逆数をとった行列である。これより、式 (B.14) から共役勾配法によって x̃ を計算し、x = P x̃ から最終的な解を求めることで、対角スケーリングによる前処理付き共役勾配法を計算できる。

89

参考文献

[1] Azouit, M. and Vernin, J., 2005. “Optical Turbulence Proﬁling with Balloons Relevant
to Astronomy and Atmospheric Physics”, Publications of the Astronomical Society of the
Pacific, 117, 536-543.
[2] Banakh, V. A. and Marakasov, D. A., 2007. “Wind proﬁling based on the optical beam
intensity statistics in a turbulent atmosphere”, J. Opt. Soc. Am. A, 24, 3245-3254.
[3] Bardsley, J. M., 2008. “Wavefront Reconstruction Methods for Adaptive Optics Systems
on Ground-Based Telescopes”, SIAM J. Matrix Analysis Applications, 30, 67-83.
[4] Bouwens, R. J., Illingworth, G. D., Oesch, P. A., Labbe, I., Trenti, M., van Dokkum,
P., Franx, M., Stiavelli, M., Carollo, C. M., Magee, D., and Gonzalez, V., 2011. “ULTRAVIOLET LUMINOSITY FUNCTIONS FROM 132 z∼7 AND z∼8 LYMAN-BREAKGALAXIES IN THE ULTRA-DEEP HUDF09 AND WIDE-AREA EARLY RELEASESCIENCE WFC3/IR OBSERVATIONS”, ApJ, 737, 90.
[5] Conan, R., Bradley, C., Lardire, O., Blain, C., Venn, K., Andersen, D., Simard, L., Vran,J.-P., Herriot, G., Loop, D., Usuda, T., Oya, S., Hayano, Y., Terada, H., and Akiyama,M., 2010. “Raven: a harbinger of multi-object adaptive optics-based instruments at theSubaru Telescope”, Proc. SPIE, 7736.
[6] Cuby, J.-G., Morris, S., Parr-Burman, P., Lehnert, M., Evans, C., Fusco, T., Jagourel, P.,
Mignant, D. L., Myers, R., Rousset, G., Schnetler, H., Amans, J.-P., Assemat, F., Beard,
S., Cohen, M., Dipper, N., Ferrari, M., Gendron, E., Gimenez, J.-L., Hastings, P., Hubert,
Z., Hugot, E., Laporte, P., Leroux, B., Madec, F., McGregor, H., Morris, T., Neichel, B.,
Puech, M., Robert, C., Rolt, S., Swinbank, M., Talbot, G., Taylor, W., Vidal, F., Vivs,
S., Vola, P., and Wells, M., 2010. “EAGLE: a MOAO fed multi-IFU NIR workhorse for
E-ELT”, SPIE Proc, 7735.
[7] David J. E. Floyd, J. T.-O. and Prieto, G., 2010. “Seeing, Wind, and Outer Scale Eﬀects
on Image Quality at the Magellan Telescopes”, Publications of the Astronomical Societyof the Pacific, 122, 731-742.
[8] Ellerbroek, B. L., 2002. “Eﬃcient computation of minimum-variance wave-front reconstructors with sparse matrix techniques”, J. Opt. Soc. Am. A, 19, 1803-1816.
[9] Els, S., Travouillon, T., Schck, M., Riddle, R., Skidmore, W., Seguel, J., Bustos, E.,
and Walker, D., 2009. “Thirty Meter Telescope Site Testing VI: Turbulence Proﬁles”,
Publications of the Astronomical Society of the Pacific, 121, 527-543.

90

Appendix B Mathematics[10] Francois Assemat, F. H., Eric Gendron,, 2006. “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., 376, 287-312.
[11] Gendron, E., Vidal, F., Brangier, M., Morris, T., Hubert, Z., Basden, A., Rousset, G.,
Myers, R., Chemla, F., Longmore, A., Butterley, T., Dipper, N., Dunlop, C., Geng, D.,
Gratadour, D., Henry, D., Laporte, P., Looker, N., Perret, D., Sevin, A., Talbot, G., andYounger, E., 2011. “MOAO ﬁrst on-sky demonstration with CANARY”, A&A, 529, L2.
[12] Gendron, E. and Lena, P., 1996. “Single layer atmospheric turbulence demonstrated byadaptive optics observations”, Astrophys. Space Sci, 239, 221-228.
[13] Gilles, L., Vogel, C. R., and Ellerbroek, B. L., 2002. “Multigrid preconditioned conjugategradient method for large-scale wave-front reconstruction”, J. Opt. Soc. Am. A, 19, 18171822.
[14] Gilles, L., Ellerbroek, B. L., and Vogel, C. R., 2003. “Preconditioned conjugate gradientwave-front reconstructors for multiconjugate adaptive optics”, Applied Optics, 42, 52335250.
[15] Johnson, L. C., Gavel, D. T., and Wiberg, D. M., 2011. “Bulk wind estimation and prediction for adaptive optics control systems”, J. Opt. Soc. Am. A, 28, 1566-1577.
[16] Neichel, B. B., Fusco, T. T., and Conan, J.-M. J., 2009. “Tomographic reconstruction for
wide-ﬁeld adaptive optics systems: Fourier domain analysis and fundamental limitations”,
J. Opt. Soc. Am. A, 26, 219-35.
[17] Ono, Y., Ouchi, M., Curtis-Lake, E., Schenker, M. A., Ellis, R. S., McLure, R. J., Dunlop,
J. S., Robertson, B. E., Koekemoer, A. M., Bowler, R. A. A., Rogers, A. B., Schneider,
E., Charlot, S., Stark, D. P., Shimasaku, K., Furlanettoa, S. R., and Cirasuolo, M., 2012.
“Evolution of the Sizes of Galaxies over 7<z<12 Revealed by the 2012 Hubble Ultra Deep
Field Campaign”, ApJ, submitted , arXiv:1212.3869.
[18] Poyneer, L. A., Gavel, D. T., and Brase, J. M., 2002. “Fast wave-front reconstruction inlarge adaptive optics systems with use of the Fourier transform”, J. Opt. Soc. Am. A, 19,
2100-2111.
[19] Poyneer, L., van Dam, M., and Vran, J. P., 2009. “Experimental veriﬁcation of the frozen
ﬂow atmospheric turbulence assumption with use of astronomical adaptive optics telemetry”, J. Opt. Soc. Am. A, 26, 833-846.
[20] Schock, M. and Spillar, E. J., 2000. “Method for a quantitative investigation of the frozen
ﬂow hypothesis”, J. Opt. Soc. Am. A, 17, 1650-1658.
[21] Strohbehn, W. J., 1968. “Line-oS-sight wave propagation through the turbulent atmosphere”, Proc. IEEE, 56, 1301-1318.
[22] Tokovinin, A., 2002. “From Diﬀerential Image Motion to Seeing”, Publications of the Astronomical Society of the Pacific, 114, 1156-1166.
91

Appendix B Mathematics[23] Hammer, F., Say‘ede, F., Gendron, E., Fusco, T., Bur-garella, D., Cayatte, V., and Conan, J., 2002. The Scientific Requirements for Extremely Large Telescopes, Cape Town:
Cambridge University Press.
[24] Hardy, J. W., 1998. Adaptive Optics for Astronomical Telescope, New York: Oxford University Press.
[25] “Subaru Telescope AO188”, [http://www.naoj.org/Observing/Instruments/AO/performance.html].
[26] “NVIDIA HP”, [http://www.nvidia.co.jp/page/home.html].
[27] “TMT HP”, [http://www.tmt.org/].
[28] “W.M. Keck Observatory Adaptive Optics”, [http://www2.keck.hawaii.edu/optics/ao/].
[29] Kawakami,

I.,

2009.

“Fundamentalof[http://www7.ocn.ne.jp/ kawa1/].

92

Numerical

Computation(inJapanese)”,