修士論文

補償光学系波面センサーによる
大気ゆらぎの高さ分布のリアルタイム推定Real time estimationof atmospheric turbulence altitude profile
with wavefront sensors in adaptive optics systems
東北大学大学院 理学研究科 天文学専攻 修士２年
大金 原
令和元年

i

要旨
可視光・近赤外線天文観測において、大型望遠鏡の空間分解能を引き上げる補償光学は今や必須技術の一つである。これからの時代において、さらに遠方の宇宙における天体の検出や高い空間分解能での個々の天体の観測、たくさんの天体のサンプルを用いた統計的議論を行ってゆくためには、30m 級の巨大望遠鏡の開発と 8m 級の望遠鏡における高いサーベイ能力が必要となる。それに伴い補償光学の技術にも巨大化、広視野化が要求され、実現の鍵を握るのは複数のレーザーガイド星を用いた補償光学である。
このような補償光学では、複数のレーザーガイド星を用いて複数の方向の大気ゆらぎを測定し、それらの情報からトモグラフィの技術を用いて大気ゆらぎの 3 次元的な構造を推定する。この推定をもとに観測
天体の方向の波面のゆらぎを最適化によって求め、波面補償に用いる。様々な観測方向に対して共通に影響を与える地表層ゆらぎを補償することで広視野補償光学を実現することができる。また、いくつかの天体の方向に対して同時に波面ゆらぎの最適化を行うことで多天体補償光学を実現することができる。しかしながら大気ゆらぎの 3 次元的な構造を推定する過程は悪条件の逆問題を解くことであり、大気ゆらぎの高さ分布 (どの高さにどれくらいの強度のゆらぎがあるか) という事前情報が推定の精度を左右する。本
研究では大気ゆらぎの高さ分布を取得する新しい手法の検証を行った。
新しい手法としては、補償光学で用いられる波面センサーを用いた MASS-DIMM 手法を考えた。波面センサーを用いることによる観測・補償との同一方向性、MASS-DIMM 手法を応用することによるリア
ルタイムな高さ分布取得性が一度に得られると考えたためである。手法の実現性を検証すべく、実際に波面センサー光学系を設計・製作し、大気ゆらぎに起因する星の明るさ・星像位置変動を測定した。
これらを解析することにより、今回の東北大学のサイトにおける大気ゆらぎの高さ分布が典型的な大気構造と矛盾しないことを明らかにし、手法の実現性を実証した。一方でトモグラフィ補償光学の事前情報としての役割を考えた場合、大気ゆらぎ強度の推定誤差、推定の高さ方向の分解能の見積もりが重要である。本手法ではこれらはトレードオフの関係にある。今後トモグラフィ補償光学のシミュレーションに
よって、どれほどの推定誤差の範囲でどれほどの分解能での推定を行う必要があるかについての定量的な見積もりが必要である。

iii

目次

要旨
第1章

1.1

1.2

第2章

2.1

2.2

第3章

i
複数のレーザーガイド星を用いた補償光学

1

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

1

1.1.1

補償光学のシステム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.1.2

補償光学によるゲイン . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1.1.3

補償光学を用いたサイエンス . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.1.4

宇宙望遠鏡との関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

複数のレーザーガイド星を用いた補償光学 . . . . . . . . . . . . . . . . . . . . . . . .

9

1.2.1

LTAO : 可視光領域の補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.2.2

WFAO : 広視野の補償光学 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.2.3

トモグラフィの事前情報としての大気ゆらぎの高さ分布 . . . . . . . . . . . . .

12

大気ゆらぎの高さ分布の推定

15

大気のモデル . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.1

Kolmogorov の乱流理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.2

地球大気の構造. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.1.3

大気ゆらぎに起因する種々の物理量のゆらぎ . . . . . . . . . . . . . . . . . . .

17

2.1.4

大気ゆらぎの強度のパラメータ . . . . . . . . . . . . . . . . . . . . . . . . . .

18

大気ゆらぎの高さ分布の推定手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.2.1

バルーンを用いた直接測定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

2.2.2

2 つの波面センサーの相関を用いる方法 . . . . . . . . . . . . . . . . . . . . . .

21

2.2.3

MASS-DIMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

2.2.4

各手法のまとめ

30. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

波面センサーを用いた MASS-DIMM の提案

31

3.1

Shack-Hartmann 型波面センサーを用いた MASS-DIMM 法の提案 . . . . . . . . . . .

31

3.2

MASS を行うために . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

3.3

DIMM を行うために . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

目次
第4章

4.1

4.2
第5章

iv
東北大学 50cm 望遠鏡と波面センサーによる星の観測

35

光学系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.1.1

波面センサー光学系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.1.2

EMCCD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

東北大学 50cm 望遠鏡を用いた観測 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

新手法の実践および結果

40

5.1

星の変位・明るさの時間変動. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

40

5.2

MASS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

5.2.1

開口パターンの取り方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

5.2.2

Scintillation Index の計算方法 . . . . . . . . . . . . . . . . . . . . . . . . . .

45

5.2.3

明るさゆらぎの相関関数としての Scintillation Index . . . . . . . . . . . . . .

45

5.2.4

Weighting Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

5.2.5

行列の逆解き手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

5.2.6

大気ゆらぎの高さ分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

DIMM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3.1

開口の取り方 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3.2

Diﬀerential Image Motion の計算 . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3.3

構造関数の理論曲線とのフィット . . . . . . . . . . . . . . . . . . . . . . . . .

54

5.3.4

積分時間の補正. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

5.3.5

フリード長とシーイングの推定結果 . . . . . . . . . . . . . . . . . . . . . . . .

55

MASS と DIMM の結果の対応関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

議論

58

MASS の結果に関して . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

58

6.1.1

大気ゆらぎの高さ分布の妥当性 . . . . . . . . . . . . . . . . . . . . . . . . . .

58

6.1.2

大気ゆらぎの強度の推定エラー . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.3

5.4
第6章

6.1

6.2

DIMM の結果に関して. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

6.3

トモグラフィの事前情報取得手法としての機能性 . . . . . . . . . . . . . . . . . . . . .

63

6.4

今後の展望 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

63

結論

64

第7章
謝辞

65

参考文献

66

付録 A

A.1

大気ゆらぎの統計的取り扱い

70

ゆらぎと統計量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

A.1.1 ゆらぎ (Fructuation) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

目次

A.2

v
A.1.2 自己相関関数 (Auto Correlation Function) . . . . . . . . . . . . . . . . . . .

70

A.1.3 相互相関関数 (Cross Corelation Function) . . . . . . . . . . . . . . . . . . . .

71

A.1.4 自己構造関数 (Auto Structure Function) . . . . . . . . . . . . . . . . . . . .

71

A.1.5 相互構造関数 (Cross Structure Function) . . . . . . . . . . . . . . . . . . . .

72

A.1.6 パワースペクトラム (Power Spectrum) . . . . . . . . . . . . . . . . . . . . . .

72

大気ゆらぎによる影響 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

A.2.1 光の状態の記述 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

A.2.2 様々な物理量のゆらぎと統計量 . . . . . . . . . . . . . . . . . . . . . . . . . .

73

1

第1章

複数のレーザーガイド星を用いた補償
光学
1.1 補償光学
1.1.1 補償光学のシステム
地上望遠鏡は地球の大気を通して天体の観測を行う。地球の大気は乱流構造を持っているため、その屈折率は乱流を特徴付ける長さスケールで異なる。この屈折率の空間ゆらぎを大気ゆらぎと呼ぶ。天体から届いた光は大気ゆらぎの影響を受けて平面波から乱れた波面（等位相面）になる。結果として、地上望遠鏡で得られる天体像はぼやけてしまう。補償光学（Adaptive Optics; AO）は、大気ゆらぎの影響による波面の乱れをリアルタイムに補正することで、シャープな天体像を実現するシステムである。補償光学シ
ステムの全体像を図 1.1 に示す。補償光学システムは大きく分けて以下の 4 つのコンポーネントから構成される。

1. ガイド星 (Guide Star; GS)
大気ゆらぎの影響による波面の乱れをモニタするために用いる点光源である。観測天体の付近にある自然の星 (Natural Guide Star; NGS) を用いる場合と、レーザーを打ち上げて上空 90km に存在するナトリウム層を励起発光させることで作る人工の星 (Laser Guide Star; LGS) を用いる
場合とがある。（LGS を用いた場合の補償光学系の概念図は図 1.2 に示した。）GS に要求される
明るさと観測天体への近さは、例えばすばる望遠鏡の補償光学装置 AO188 による典型的な大気状
況下での補償を考えた場合、観測天体から 1 分角以内の mR < 13 の星となる [1]。この要求から、

NGS を用いた場合のスカイカバレッジは全天の 1% 以下にとどまってしまう。一方で LGS を用いた場合のスカイカバレッジは全天の 90% 以上に達する。人工の星を用いているにも関わらずス
カイカバレッジが 100% にならないのは、LGS は波面全体の傾き (Tip/Tilt 成分) に感度を持たないからである。LGS を作るためのレーザー光が大気ゆらぎを空に向かって通過する際に受ける屈折の影響と、LGS からの光が大気ゆらぎを地上に向かって通過する際に受ける屈折の影響とが相殺してしまうことが原因である。したがって LGS を用いる場合でもなお、波面の Tip/Tilt 成分を測定する自然の星 (Tip/Tilt-NGS; TT-NGS) が必要となる。Tip/Tilt の測定には分割数の少な第1章

複数のレーザーガイド星を用いた補償光学

図 1.1

2

補償光学系の概念図。ガイド星が受けている大気ゆらぎの影響を波面センサーで測定し、観測天体が受けている影響の補正に利用する。実際には観測天体もガイド星も無限遠に存在するため、観測天体とガイド星が影響を受ける大気乱流の領域はほとんど同じである。

い波面センサーで十分な光子数を稼ぐことができるため、TT-NGS に要求される明るさはより暗くなり AO188 の場合で mR < 17 である [1]。

2. 波面センサー (Wavefront Sensor; WFS)
大気ゆらぎの影響による波面の乱れを測定するセンサーである。波面センサーには Shack-Hartmann 型やピラミッド型、曲率型などいくつかの種類が存在する。ここでは、代表してShack-Hartmann 型波面センサー (Shack-Hartmann Wavefront Sensor; SHWFS) を紹介する。
SHWFS の概念図を図 1.3 に示す。SHWFS は、マイクロレンズアレイと呼ばれる、小さなレンズが格子状に並んだレンズと検出器から構成される。SHWFS に光が入射すると、マイクロレンズア
レイによって複数の像が格子状に検出器上に結ばれる。この時の各像の位置は、マイクロレンズに入射する光の局所的な波面の傾きで決まる。したがって、平面波が入射した時の像の位置からの変位を測定することで波面の形状を推定することができる。波面の分割数を上げることでより高次の波面形状が測定できる一方で、分割数を上げると 1 つの像あたりの明るさが暗くなり変位の測定が難しくなる。すばる望遠鏡の補償光学装置 AO188 ではシャックハルトマン型ではなく曲率型が用

第1章

複数のレーザーガイド星を用いた補償光学

図 1.2

レーザーガイド星補償光学系の概念図。NGS を用いた補償光学に比べてスカイカバレッジが格段に良くなる。一方で LGS の高さは有限であるため LGS から届く光と観測天体から届く光が通る
経路が異なり、波面測定のエラーとして残る。

いられているが、8.2m の主鏡を 188 個に分割して波面の測定を行なっている。

図 1.3

シャックハルトマンセンサーの概念図。検出器上に結ぶ像の変位の情報から、波面の局所的な傾きを推定することができる。

3

第1章

複数のレーザーガイド星を用いた補償光学

43. 制御系 (Control system)
補償光学における制御系は、波面センサーで取得される画像の解析から波面歪みの情報を取り出し、その歪みを打ち消すように後述の可変形鏡のアクチュエーターに電圧値を送るコンピュータシ
ステムである。大気ゆらぎの構造は時間と共に変化してゆくため、天体を撮像している間高い空間分解能を保つためには常に制御ループを回し続ける必要がある。この時のループの速度は AO188の場合で 1kHz である [2]。

4. 可変形鏡 (Deformable Mirror; DM)
大気ゆらぎの影響による波面の乱れを補正するための鏡である。鏡の裏面に多数のアクチュエー
ターが取り付けられており、それらを独立に動かすことで鏡面形状を変えることができる。反射後に波面の乱れが打ち消されるように鏡面形状が制御される。
これらのサブシステムによって、長い観測時間にわたって高い空間分解能での天体像を取得することができる。

1.1.2 補償光学によるゲイン
ここでは補償光学によって天文観測にどのようなゲインが得られるかについて述べる。まず図 1.4 に示
すのは望遠鏡の角度分解能を波長の関数として表したものであり、望遠鏡が回折限界を達成した時と大気
ゆらぎの影響を受けた場合とを比較した図である。ここで、角度分解能は点像分布関数 (PSF) の半値全幅 (FWHM) としている。

FWHM of PSF [arcsec]

100
10 1
10 2
diffraction limited(D=8m)
diffraction limited(D=30m)
seeing limited(r0=10cm @500nm)
seeing limited(r0=20cm @500nm)

10 3
10 1

100
101
observation wavelength [micron]

図 1.4 角度分解能 (PSF の FWHM) と波長との関係。青系の色の線は回折限界を達成した場合の分解能を、赤系の色の線は大気ゆらぎによって制限された場合の分解能を示している。波長が短くなる
ほど補償光学の持つ潜在的なゲインが大きいことが読み取れる。

第1章

複数のレーザーガイド星を用いた補償光学

5

この図において青系の色の線は回折限界を達成した場合の分解能を、赤系の色の線は大気ゆらぎによって制限された場合の分解能を示している。回折限界像のサイズは観測波長 λ に対して ∝ λ の関係があるのに対して、大気ゆらぎの影響を受けた時の像のサイズは ∝ λ−1/5 の関係がある。したがって両者の差は観測波長が短くなるほど大きくなる。中間赤外線に近い領域における観測ではすばる望遠鏡でも TMT
でも望遠鏡の回折限界が像のサイズを決めるため補償光学は必要とならない一方、可視光・近赤外線における観測では補償光学による潜在的なゲインが非常に大きいことがわかる。
しかし、実際の補償光学の運用において完全な回折限界像が得られるわけではない。波面測定の誤差(photon error) や測定から補正までの時間的遅れ (temporal error)、 波面補償の誤差 (fitting error)、ガ
イド星と観測天体との方向の違い (angular anisoplanatic error) などによってシステム全体のエラーが生じるためである。これを定量的に評価するパラメータとして波面誤差 (wavefront error; WFE) とストレル比 (Strehl ratio; SR) を紹介する。まず、WFE は補償した後の波面 f (x, y) に対して以下の式、√∫WFE[nm] =(f (x, y)− < f (x, y) >)2 dxdy
∫
dxdy

(1.1)

によって定義されている。単位は nm である。f (x, y) は閉ループ制御の場合、波面センサーによって取得される。先に述べた様々な誤差要因に対して、

WFE[nm] =√WFE2photon + WFE2temporal + WFE2f itting + · · ·(1.2)という 2 乗和の関係が成り立つ。ストレル比は PSF のピークの高さの回折限界像のピークの高さに対する比で定義される。光の振幅ゆらぎが無視でき、位相ゆらぎの 3 次以上の項が無視できるとき、ストレル比は近似的に WFE のみの関数でかけることが知れらている [3]。

{ (
)2 }
Peak hight of PSF
2π(WFE[nm])SR =∼ exp −
Peak hight of diﬀraction limited PSF
λ式 1.3 にしたがって、SR と観測波長との関係をプロットしたのが図 1.5 である。

(1.3)

第1章

複数のレーザーガイド星を用いた補償光学

1.0

WFE=100[nm]WFE=200[nm]WFE=300[nm]WFE=400[nm]WFE=500[nm]0.8Strehl ratio

6

0.6
0.4
0.2
0.0

10 1

100
101
observation wavelength [micron]

図 1.5 ストレル比 (SR) と波長との関係。一般的に、SR > 0.2 で補償が効いていると判断され、SR > 0.8 は回折限界像とみなされる。短波長になるほど十分な SR を得るために必要とされる要求波面精度が高くなる。

一般的に大型望遠鏡で補償光学を用いない場合のストレル比は SR < 0.01 である。SR > 0.2 で補償が効いていると判断され、SR > 0.8 は回折限界像とみなされる。このことを踏まえて図 1.5 を見ると、波面誤差の要求精度は近赤外線（2µm）では WFE < 400nm 程度であるのに対して可視光（500nm）ではWFE < 100nm であり、短波長での補償光学の難しさが見て取れる。この理由から、現在の補償光学シ
ステムは基本的に近赤外線領域において運用されている。
次に補償光学が一度に適用できる視野について述べる。図 1.6 は観測波長とアイソプラナティック角（詳細は 2 章にて）との関係を表した図である。アイソプラナティック角とは大気ゆらぎに起因する波面の乱れが同一と見なせる角度領域であり、観測天体とガイド星との離角はこの角度よりも小さくなる必要がある。言い方を変えると 1 つのガイド星によって補償することができる角度領域を表しており、補償光学のゲインが得られる視野を表していると言える。高度分布を含めた大気ゆらぎの状態で決まるため観測
サイトによって異なる。ここでは Skidmore et al. 2009 [4] による TMT のサイト調査の大気ゆらぎの高さ分布測定結果に基づく計算を示している。

複数のレーザーガイド星を用いた補償光学

Isoplanatic angle (arcsec)

第1章

7

Cerro Armazones
Cerro Tolonchar
Mauna Kea

102
101
100
10 1

100
101
observation wavelength [micron]

図 1.6 1 つのガイド星を用いた補償光学の視野（アイソプラナティック角）と観測波長との関係。

Skidmore et al. 2009 [4] によって示された、各サイトの 500nm におけるアイソプラナティック角を元に波長の 6/5 乗に比例する線を描いたものである。1 つのガイド星で達成される補償光学の視野は
近赤外線領域で 10 秒角程度である。

図 1.6 からわかるように、アイソプラナティック角は波長に対して ∝ λ6/5 の関係を持っており、長波長ほど補償光学の視野は広くなる。一方で視野の大きさは近赤外領域 (2µm) でも高々 10 秒角であり、1つのガイド星を用いた補償光学における視野の狭さが伺える。
以上をまとめると、補償光学システムは可視光・近赤外線観測における角度分解能を大きく引き上げる技術であり、特に可視光においてそのゲインは大きい。一方で可視光における補償は近赤外に比べて波面誤差に要求される精度が高く難易度は高い。また、1 つのガイド星で達成される観測視野は 10 秒角と狭く、可視光ではさらに狭くなる。

1.1.3 補償光学を用いたサイエンス
ここでは、前節で述べた補償光学における観測への効果を踏まえて、具体的なサイエンスを紹介する。

1. AGN アウトフローの観測
AGN フィードバックは、AGN から磁場による駆動で吹き出すジェットによって、銀河の周囲のガスを吹き飛ばしたり圧縮加熱したりすることでガスの降着や星の形成を妨げるメカニズムである。星形成銀河の quench やブラックホールと母銀河との共進化などを説明するシナリオとして提唱されている一方で、AGN フィードバックが起こっていることの観測的証拠はまだ乏しいのが現状である。銀河中心部から高速で吹き出すガスアウトフローの観測が補償光学と面分光との組み合
わせによって可能となる。例えば Genzel et al.2014[5] では、VLT+SINFONI（近赤外補償光学つ第1章

複数のレーザーガイド星を用いた補償光学8

き面分光装置）を用いて z ∼ 2 の大質量星形成銀河を観測し、銀河中心の領域において速い（速度の FWHM で ∼ 1500km/s）ガスの流れを検出した。さらにそのようなガスの吹き出しが多くの大質量星形成銀河で起こっていることが示された。

2. 高赤方偏移における銀河の形態観測z ∼ 3 を超えるような遠方銀河に対しては、明るさの問題から面分光のような手法は難しいものの高い空間分解能で撮像観測することによってその形態的特徴を調べることができるようになる。例
えば Akiyama et al.2007[6] では、 Subaru+AO36+IRCS（近赤外補償光学つき撮像装置）を用いて z ∼ 3 のライマンブレイク銀河を観測し、その形態が分厚いガス円盤による高い面密度を持ったものであることが示された。

3. 中間赤方偏移における銀河の面分光観測z ∼ 1 − 2 の中間赤方偏移は銀河の星形成活動が最も盛んになる時期であり、この時代の銀河を調べることは銀河の形成進化史を知る上で非常に重要である。このくらいの光度距離になると銀河面分光により速度場の情報や場所ごとの星形成活動や金属量の違いなどを調べられるようになる。

Schreiber et al.2018[7] では、VLT+SINFONI（近赤外補償光学つき面分光装置）を用いて z ∼ 2の星形成銀河を大規模に面分光サーベイし、星の運動の分布と星形成の分布を調べた。

4. 系外惑星の直接検出
系外惑星の直接撮像に用いられるコロナグラフの技術には、極限補償光学と呼ばれる狭い角度領域(∼ 1”) に対して高いコントラストを実現する補償光学が併用される。Currie et al.2018[8] では、Subaru+SCExAO（近赤外補償光学つきコロナグラフ装置）を用いて系外惑星の直接撮像による検出を行なった。

5. 銀河中心領域の星の観測
高い空間分解能を実現する補償光学は、星々が込み入った銀河中心などの領域における個々の星の観測を可能にする。例えば Mattila et al.2007[9] では、VLT+NAOS+CONICA（近赤外補償光
学つき撮像分光装置）を用いることで、近傍の LIRG のダストに囲まれた銀河中心部において超新星を発見した。また、Schodel et al.2002[10] では、銀河系中心を運動する星々の軌道観測から銀河系中心に SMBH が存在していることを突き止めた。

1.1.4 宇宙望遠鏡との関係
補償光学によって得られる観測のゲインは大気の影響を低減した結果である。したがってそもそも大気の影響のないスペースからの観測を行うことができれば補償光学を使うことなく高い空間分解能での観測ができる。実際にハッブル宇宙望遠鏡はこれまでに可視光・近赤外領域における非常に質の高い観測を可能としてきており、今後ジェームズウェッブ宇宙望遠鏡が打ち上げられることにより宇宙望遠鏡を用いた
天文学の幅がさらに広がることだろう。しかし、宇宙望遠鏡にもデメリットがある。ロケットで打ち上げる際の制限から望遠鏡のサイズや重量、耐衝撃性などに制限がかかること、開発や維持にコストがかかる
こと、メンテナンスのための技術者のアクセスが難しいこと、それに伴い一般に観測期間が短いことなどが挙げられる。したがって地上望遠鏡と宇宙望遠鏡は相補的な関係にあると言え、地上望遠鏡において宇第1章

複数のレーザーガイド星を用いた補償光学

9

宙望遠鏡と同等の性能を実現できる補償光学のシステムは今後も重要な天文観測技術の一つとなるであろう。

1.2 複数のレーザーガイド星を用いた補償光学
前節で述べたように、補償光学による高い空間分解能が期待できる領域は、波長的には近赤外の比較的長い波長、空間的には近赤外線観測で 10 秒角程度に限られる。ここでは、複数のレーザーガイド星を用いることでこれらの問題を解決するシステムについて説明する。

1.2.1 LTAO : 可視光領域の補償光学
短波長になるほど補償光学が難しくなる原因は、図 1.5 に示した通りである。レーザートモグラフィ補
償光学 (Laser Tomography Adaptive Optics; LTAO) は、波面の測定エラーの一つを取り除くことでシステム全体のエラーを小さくし、短波長での補償光学を実現するシステムである。
図 1.2 に示したように、有限の高度にある LGS から届く光の経路は観測天体から届く光の経路と異なるため、大気乱流による影響も異なる。この根本的な波面測定エラーは、LGS から届く光の経路が円錐形になることから「円錐効果」と呼ばれている。LTAO では、複数のレーザーガイド星を観測天体を取り囲むように配置し、観測天体から届く光の経路をカバーすることで円錐効果を克服する。

LTAO の概念図を図 1.7 に示す。図に示すように、4 つの LGS は観測天体が影響を受ける大気乱流の領域を全てカバーしているため、十分な情報取得ができている。一方で、各 LGS がカバーする領域の重なり方が高度により異なるため、測定された LGS の波面の情報から観測天体が受けているであろう波面の歪みの影響を推定する過程が必要となる。この過程では、トモグラフィの原理を用いた推定が行われる。トモグラフィ推定における計算方法の詳細は Appendix に記した。もちろんトモグラフィ推定にも推定のエラーがつくものの、円錐効果によるエラーよりも小さいため LTAO によって測定のエラーを小さくできることがシミュレーションによって試算されている。また、円錐効果は望遠鏡の開口が大きくなるほど顕著になる。したがって次世代の大型望遠鏡では、近赤外における補償光学をする場合にもトモグ
ラフィの技術が必要となることが予想される。LTAO は 8m 級の望遠鏡で可視光領域における補償を実現する技術と、30m 級の望遠鏡の時代における補償光学の基盤技術という 2 つの側面を併せ持った重要な技術の一つである。

第1章

複数のレーザーガイド星を用いた補償光学10

図 1.7 LTAO の概念図。複数 (この図では 4 つ) の LGS が観測天体が影響を受ける大気乱流の領域をカバーするように配置されることによって波面測定のエラーを小さくすることができる。同時に、各 LGS からの波面を測定する波面センサーも複数必要となる。

1.2.2 WFAO : 広視野の補償光学
広視野補償光学 (Wide Field Adaptive Optics; WFAO) は、複数の LGS を LTAO の場合よりも広い間隔で配置することで、測定する大気乱流の領域を広げるものである。LTAO が狭い領域の大気乱流をより高精度に測定し可視光での補償を実現するのに対して、WFAO は測定精度は従来の 1 つの LGS を用いた補償光学と同程度であり近赤外領域の補償に留まるものの広い視野を実現する。WFAO には、大きく分けて以下の 3 つの種類がある。図 1.8 にそれらの概念図を示す。

第1章

複数のレーザーガイド星を用いた補償光学

11

図 1.8 WFAO の概念図。LGS を複数用いる必要上、波面センサーも複数必要になる。GLAO では
視野全域にわたって地表層のみ補償するため可変形鏡は 1 つで十分である。一方で MCAO では上空の層も補償し、MOAO では天体ごとに補償するため可変形鏡も複数必要となる。

1. 地表層補償光学 (Ground Layer Adaptive Optics; GLAO)
地表層（上空数百メートルまで）の大気ゆらぎは様々な観測方向に対して共通に影響を与える大気ゆらぎ成分である。このことを利用し、GLAO では地表層ゆらぎの影響のみを取り除くことで広い視野を実現する。もちろん地表層のみの補償しか行わないため星像の質としては単一のレー
ザーガイド星を用いた従来の補償に比べて劣る。しかし地表層ゆらぎは大気ゆらぎ全体のパワーの50% 以上を持っていることが一般に知られており、地表層ゆらぎの補償だけでも十分な像質のゲインがあると言える。また GLAO は WFAO の中で最も広い視野（近赤外で ∼15 分角）が達成できるシステムである。複数の LGS からの波面を複数の波面センサーで測定しそれらに共通する波面ゆらぎ成分を取り出すことで、地表層成分を推定し補正することができる。実際の運用計画例としては、ULTIMATE-SUBARU において採用されているシステムである。

2. 多層共役補償光学 (Multi Conjugate Adaptive Optics; MCAO)MCAO では、複数の LGS の波面の情報から LTAO と同様にして測定した大気ゆらぎの高さ分布を用いて、各高度ごとに別々の DM によって補償する。各 DM は大気ゆらぎの各高度に光学的に共役な位置に置かれ、ゆらぎの影響を補正する。GLAO と異なり高高度の大気ゆらぎの影響も補正するため像の質は単一のレーザーガイド星を用いた従来の補償と同等のまま、近赤外で ∼1 分角
程度にまで視野を広げることができる。実際の運用計画例としては、TMT の第 1 期補償光学装置

NFIRAOS に採用されているシステムである。
3. 多天体補償光学 (Multi Object Adaptive Optics; MOAO)
MOAO では、複数の LGS の波面の情報から LTAO と同様にして測定した大気ゆらぎの高さ分布を用いて、各天体ごとに別々の DM によって補償する。各 DM は各観測天体の方向の大気ゆらぎの影響を補正する。複数の天体に対して同時に LTAO を行うイメージである。したがって像の質は単一のレーザーガイド星を用いた従来の補償よりも高く、視野は近赤外で ∼10 秒角程度となる。
一方で天体をピックアップすることのできる領域は近赤外で ∼5 分角程度にまで広げることができる。実際の運用計画例としては、TMT の第 2 期装置として提案されている多天体分光器と組み合第1章

複数のレーザーガイド星を用いた補償光学

12

わせた補償光学装置としての提案がある。

1.2.3 トモグラフィの事前情報としての大気ゆらぎの高さ分布
ここまでで紹介した LTAO と 3 つの WFAO に共通するのは、複数の LGS からの光の波面情報から、
大気ゆらぎの 3 次元的な構造を復元し、さらに各観測方向や各観測天体からの光の波面を推定する点である（トモグラフィック波面推定と呼ぶ）。しかしトモグラフィに用いることのできる情報は、数十秒角の離角で打ち上げられた高々数個の LGS からの波面情報であるため、縮退が大きく、数が少ない情報から大気の性質を復元する過程は悪条件の逆問題となる。したがって有効な波面推定を行うためには大気ゆらぎの分布に関する事前情報が不可欠であり、事前情報の質が最終的な推定、ひいては補償光学の精度を左右する。詳しくは 2 章にて述べるが、地球の大気乱流は地表から高度 20km 程度にかけて薄い層状に分布している。事前情報としてはこの「大気ゆらぎ強度の高さ分布」が必要となる。
ここで Fusco & Costille2010[11]、Costille & Fusco2012[12] によって示された事前プロファイルの重要性のシミュレーション結果を紹介する。彼らの計算は以下のようなものである。まずバルーン実験によって取得された詳細な大気ゆらぎの高さ分布から N 個の主要な大気ゆらぎ層を取り出し、これを大気モデルとしてシミュレーション大気を作る。次にその大気を通ってきた LGS や各観測方向からの光が持つ波面を計算する。次に N 個の層の中から Nrec（ただし Nrec ≤ N ）個の層を取り出して事前情報とし、LGS の波面情報を用いてトモグラフィック波面推定を行う。最後にトモグラフィック波面推定によって求められた各観測方向の波面と、シミュレーション上の波面との誤差を計算する。
この計算の中で N と Nrec が変数になっている。N を変えることは大気モデルの質を変えることに相当しており、Nrec を変えることはトモグラフィック波面推定の事前情報の質を変えることに相当している。
図 1.9 は、N を変えながら常に Nrec = N とした状態でトモグラフィック波面推定の推定誤差を調べたものである。横軸に N の値、縦軸に推定波面とシミュレーション上での実際の波面との残差に対する RMS 値を取っている。黒と赤の線は LGS の離角と補償光学のモードがそれぞれ 4.3′ の LTAO,2′の MCAO の場合のシミュレーションになっている。いずれも 30m 級の望遠鏡を仮定している。常に

Nrec = N としているためトモグラフィック波面推定は十分な事前情報のもとで行われている状態であり、この場合の推定誤差はトモグラフィック波面推定が持つ本質的なエラーを意味する。したがってこの
図では 20 以上の大気ゆらぎ層をもつ大気モデルを考えないと推定誤差を過小評価してしまうことが示さ
れている。

第1章

複数のレーザーガイド星を用いた補償光学13

図 1.9 Fusco & Costille2010[11] によるシミュレーションの結果。横軸はシミュレーションの大気モデルに用いるゆらぎ層の数 N 、縦軸は大気モデルが分かっている状態での波面推定エラーを表している。この結果から、20 以上のゆらぎ層を考慮したモデルを用いないと推定エラーを過小評価してしまうことがわかる。

図 1.10 は、N = 250 とした上で Nrec を変えながらトモグラフィック波面推定の推定誤差を調べたものである。横軸に Nrec の値、縦軸に推定波面とシミュレーション上での実際の波面との残差に対するRMS 値を取っている。まず N = 250 であることから大気モデルの質によって推定誤差を過小評価している心配はない。また Nrec が大きくなるにつれて推定誤差が小さくなっていくことは、事前情報の質がトモグラフィック波面推定の誤差に影響を与えていることを意味する。ここでも Nrec > 20 においてト
モグラフィック波面推定の誤差が一定値になっており、Nrec ∼ 20 がエラーを最小限に抑えるために事前
プロファイルに必要とされるゆらぎ層の数になっている。一般的に大気ゆらぎは地表から高度 20km までの間に分布していることを考えるとこの結果は、およそ 1km ごとの分解能で大気ゆらぎの高さ分布を取得する必要があることを示している。

図 1.10 Fusco & Costille2010[11] によるシミュレーションの結果。横軸は事前情報としての大気プ
ロファイルの層数、縦軸は波面推定エラーを表している。この結果から、20 以上のゆらぎ層を考慮した事前プロファイルを用いることで推定エラーを十分小さくすることができることがわかる。

第1章

複数のレーザーガイド星を用いた補償光学

14

また、大気ゆらぎの高さ分布は数十分ほどの時間スケールでの時間変動を伴うことが知られている。図1.11 に示すのは Kornilov et al. 2003[13] で示された大気ゆらぎの高さ分布測定を長時間行なった結果である。横軸に観測時間、縦軸に地表からの高さをとっており各高さでの大気ゆらぎの強度をバーの長さによって示している。ある高度のゆらぎ強度の時間変化を見ると、数十分程度のスケールで大きな変動をしていることが読み取れる。大気ゆらぎ層が風速によって移動することを考えれば観測方向による違いも生まれてくる。

図 1.11 Kornilov et al. 2003[13] で求められた大気ゆらぎの高さ分布の時間変化。ある高度の大気
ゆらぎの強度変動に注目すると、強度が数時間程度のタイムスケールで大きく変化していることがわかる。各高度の振る舞いは独立であり、大気ゆらぎ全体としては数十分程度のタイムスケールでその高さ分布を変えている。

以上をまとめると、トモグラフィック波面推定を用いる補償光学には、観測と同じ方向（視野）における ∼ 1km の分解能での大気ゆらぎの高さ分布を数十分程度のリアルタイム性をもって取得し、事前プロファイルとして更新することが要求されるということになる。

15

第2章

大気ゆらぎの高さ分布の推定
2.1 大気のモデル
2.1.1 Kolmogorov の乱流理論
大気のもつ乱流構造について、以下 Roddier 1981[14] に従って説明する。大気などの流体の流れは、流体のレイノルズ数 Re = V0 L0 /ν0（ただし V0 は流れを特徴付ける速度、L0 は流れを特徴付ける空間ス
ケール、ν0 は流体の動粘性係数）が流れの幾何学的な構造で決まるある閾値を超えると乱流となることが知られている。地球大気の場合には V0 ∼ 1[m/s], L0 ∼ 15[m], ν0 ∼ 1.5 × 10−5 [m2 /s] であり、Re ∼ 106となる。一般的にこのレイノルズ数はよく発達した乱流に相当することから、地球大気は基本的に乱流状態になっていることがわかる。この乱流によって生じる大気屈折率のムラが天体から届く光の位相のムラを作り、集光したときに像の劣化を引き起こす。この大気屈折率のムラのことを「大気ゆらぎ」と呼ぶ。
大気ゆらぎのパワースペクトラムは Kolmogorov の理論 (Kolmogorov 1941; [15]) によってよく表されることが知られている。Kolmogorov のモデルでは、大気は太陽からの輻射によって大きな空間スケー
ルの温度揺らぎを持ち、大気循環によって段々とより小さなスケールの揺らぎへと分割されていき (エネルギーカスケード)、最終的に熱としてエネルギーが散逸する (粘性散逸) と考える。この理論に従うと、
大気ゆらぎの 3 次元パワースペクトラム（大気屈折率のムラの 3 次元パワースペクトラム）は、
2
ΦN (fx , fy , fz ) = 9.7 × 10−3 f −11/3 CN
[m3 ]
(2.1)
√
2
とかける。ここで、fx , fy , fz は各方向の空間周波数 [m−1 ]、f = fx2 + fy2 + fz2 であり、CN[m− 2/3]2は大気屈折率の構造定数と呼ばれる定数である。ここで CNの物理的意味を説明するために構造関数について説明する。
構造関数は相関関数と同様に、ある物理量や物理量同士の相関の度合いを表す関数である。位置の関数としての物理量 f (⃗x) の自己構造関数は次のように定義される。

Df (⃗x) ≡< [f (⃗r) − f (⃗r + ⃗x)]2 >

(2.2)

第2章

大気ゆらぎの高さ分布の推定

16<> はアンサンブル平均を表す。相関関数 Bf (⃗x) =< f (⃗r)f (⃗r + ⃗x) > との関係は、Df (⃗x) ≡< [f (⃗r) − f (⃗r + ⃗x)]2 >(2.3)=< f (⃗r) > + < f (⃗r + ⃗x) > −2 < f (⃗r)f (⃗r + ⃗x) >(2.4)= 2[Bf (⃗0) − Bf (⃗x)]

(2.5)

2

2

となる。一方で物理量 f のゆらぎ ∆f = f − < f > を議論したい場合、D∆f (⃗x) =< [f (⃗r)− < f > −f (⃗r + ⃗x)+ < f >]2 >(2.6)=< [f (⃗r) − f (⃗r + ⃗x)]2 >(2.7)= Df (⃗x)(2.8)となることから、構造関数の値は物理量の平均値に依存しないことがわかる。構造関数は平均値が時間変動するような物理量のゆらぎを議論する上で相関関数よりも適していると言える。大気屈折率の構造関数
は Obukhov 1949[16] によって解析的に求められており、
2
DN (⃗x) ≡< [N (⃗r) − N (⃗r + ⃗x)]2 >= CN|⃗x|2/3(2.9)

2とかける。この式から、CNの値が大きいほど構造関数が急激に大きくなる、すなわち相関長が小さくな
2
はゆらぎ強度の指標となるパラメータであることがわかる。
ることがわかる。CN

2.1.2 地球大気の構造
地球大気のゆらぎは層状に分布していることが知られている。したがって前節で大気ゆらぎの強度を2特徴付けるパラメータとして紹介した CNは地表からの高さ h に依存する関数であると考えることがで2きる。本論文において大気ゆらぎの高さ分布と呼ぶのは、CN(h) のことである。図 2.1 に、Hardy 1998[17] における典型的な大気ゆらぎの高さ分布を示す。

2
図 2.1 典型的な大気ゆらぎの高さ分布。右側のグラフの横軸が大気ゆらぎの強度を表す CNの値、縦軸が地表面からの高さを表している。地表層に最も近い部分でゆらぎ強度は最大であり、高度 20km
付近にかけて弱くなってゆく。しかし高度 10km 付近にはジェット気流に起因する局所的な強度のピークがある。

第2章

大気ゆらぎの高さ分布の推定

17

図 2.1 からわかるように、大気ゆらぎは地上付近で最も強い。これは大気が地表と摩擦を起こすことに起因する。したがって高度が上がるにつれて弱くなってゆき、典型的な大気圧下では上空 20km よりも高い位置にある大気ゆらぎの影響は無視できるほど小さくなる。また、上空 10km 付近にはゆらぎ強度が局所的に高くなるところがある。これは対流圏と成層圏との境界領域 (対流圏界面; tropopose) において吹いているジェット気流に起因するものである。これらの典型的な特徴に加え、緯度や地形、気温などのパラメータが複雑に絡み合うことで観測サイトや季節、昼夜によって異なる大気の特徴が作り出されている。

2.1.3 大気ゆらぎに起因する種々の物理量のゆらぎ
Kolmogorov のパワースペクトラムを仮定することで、大気ゆらぎに起因する様々な物理量のゆらぎのパワースペクトラムを解析的に導くことができる。後述の大気ゆらぎの高さ分布推定に関する予備知識となるものも含まれるためここで紹介する。また導出に関する詳細な計算は Appendix に記した。

1. 大気ゆらぎ層による光路長ゆらぎのパワースペクトラム
前述の通り、大気ゆらぎは薄い層状に分布している。ある 1 つの層（高度 h）に着目し、この層の厚みを ∆h とする。このゆらぎ層を通過することで生じる光路長 l のゆらぎの 2 次元パワースペク
トラム Φl は Kolmogorov パワースペクトラムを層の厚み方向に積分することによって得られ、∫ ∆hΦl (fx , fy ) =2ΦN (fx , fy , fz )dh ∼ ΦN (fx , fy , 0)∆h = 9.7 × 10−3 f −11/3 CN∆h[m4 ] (2.10)0となる。ここで、f =√fx2 + fy2 である。

2. 大気ゆらぎ層による光の位相ゆらぎのパワースペクトラム
光路長のゆらぎによって、天体から届くまっすぐな波面は位相ゆらぎを持つようになる。光の波数を k = 2π/λ とすると位相 ϕ は kl で与えられるため、位相ゆらぎのパワースペクトラム Φϕ は、
2
Φϕ (fx , fy ) = Φl × k 2 = 9.7 × 10−3 k 2 f −11/3 CN∆h[m2 ]

(2.11)となる。

3. 地上へ伝搬した光の位相ゆらぎのパワースペクトラム
我々が観測するのは、位相ゆらぎを持った光が地上へ伝搬した後の光である。伝搬後の光の持つゆらぎのパワースペクトラムは Kolmogorov パワースペクトラムに加えて 2 つの仮定を置くことで計算できる。1 つ目は光の伝搬がフレネル伝搬で記述できること、2 つ目は乱流が比較的弱く各点の位相と平均値との差が 1rad よりも十分小さいこと (弱乱流近似) である。これらの仮定のもとで、地上へ伝搬した後の光の位相 ψ のゆらぎのパワースペクトラム Φψ は以下のようにかける。
2Φψ (fx , fy ) = Φϕ × cos2 (πλhf 2 ) = 9.7 × 10−3 k 2 f −11/3 cos2 (πλhf 2 )CN∆h[m2 ]

(2.12)

4. 地上へ伝搬した光の振幅ゆらぎのパワースペクトラム
位相と同様に、地上へ伝搬した光の振幅ゆらぎについても以下のようなパワースペクトラム Φχ が導かれる。
2Φχ (fx , fy ) = Φϕ × sin2 (πλhf 2 ) = 9.7 × 10−3 k 2 f −11/3 sin2 (πλhf 2 )CN∆h[m2 ]

(2.13)

第2章

大気ゆらぎの高さ分布の推定

185. 星像の位置ゆらぎのパワースペクトラム
星像の位置は、波面の局所的な傾きで決まる。ψ(⃗x) で与えられる波面の光が入射してきた時、検出器から見た光の入射角度は

(

∂ λ
∂ λ
− ∂x
x), − ∂yx)
2π ψ(⃗
2π ψ(⃗)

で与えられる。この関係を用いると、傾きをとる方向 x, y に対する入射角度 α, β のゆらぎのパワースペクトラム Φα , Φβ は次のようになる。
2
Φα (fx , fy ) = λ2 fx2 Φψ = 9.7 × 10−3 (2π)2 fx2 f −11/3 cos2 (πλhf 2 )CN(h)∆h[m2 ](2.14)2Φβ (fx , fy ) = λ2 fy2 Φψ = 9.7 × 10−3 (2π)2 fy2 f −11/3 cos2 (πλhf 2 )CN(h)∆h[m2 ](2.15)

また、ここまでは 1 つの大気ゆらぎ層によって生じる物理量のゆらぎを考えてきたが、前節で見たように実際には複数のゆらぎ層からの寄与を考慮する必要がある。弱乱流近似を仮定するとき、各高度のゆらぎによるパワースペクトラムへの寄与は単純和になることが示せる。これらを踏まえると、星像の位置ゆらぎのパワースペクトラムは、

Φα (fx , fy ) =Φβ (fx , fy ) =∑λ2 fx2 Φψ =

∑

layer

j

∑

∑

λ2 fy2 Φψ =2
9.7 × 10−3 (2π)2 fx2 f −11/3 cos2 (πλhj f 2 )CN(hj )∆hj [m2 ](2.16)
2
9.7 × 10−3 (2π)2 fy2 f −11/3 cos2 (πλhj f 2 )CN(hj )∆hj [m2 ]

j

layer

(2.17)と書ける。

6. 星像の明るさゆらぎのパワースペクトラム
星像の明るさゆらぎは天体によって元の明るさが異なるので明るさの平均値で規格化したゆらぎを考えることとする。この時天体の明るさ I は光の振幅 χ を用いて I(⃗x) ∼ 1 + 2χ(⃗x) と書ける。明るさゆらぎは振幅ゆらぎの 2 倍であるため、パワースペクトラムとしては 4 倍となる。したがって複数の大気ゆらぎ層からの寄与も考慮すると、パワースペクトラム ΦI は、

ΦI (fx , fy ) =∑4 × Φχ =

layer

∑

1.53f

j

−11/3

{

sin(πλhj f 2 )
λ

}2
2
CN
(hj )∆hj [m2 ](2.18)と表せる。

2.1.4 大気ゆらぎの強度のパラメータ2
ここまで、大気ゆらぎの強度を表すパラメータとして大気屈折率の構造定数 CN(h) を用いてきた。しかし、様々な用途に合わせた他の指標も存在する。ここではそれらを紹介する。

1. フリード長 (Fried parameter) r0
大気ゆらぎによって天体から届く光の波面は乱されるが、小さな空間スケールでは局所的に波面の乱れはないものと見なせる。この、波面が維持されていると見なせる最大の空間スケールをフリード長と呼ぶ。したがってフリード長よりも小さな開口サイズを持つ望遠鏡で観測した場合、短い露
出時間においては天体像が大気ゆらぎから受ける影響はないものと見なせる。フリード長は以下の第2章

大気ゆらぎの高さ分布の推定19

式で与えられる。

[

0.423k 2
r0 =
cos γ

]− 35

∫

2
CN
(h)dh

(2.19)

ここで、k = 2π/λ は観測している光の波数、γ は観測方向の天頂角である。式から分かる通り、
フリード長は大気ゆらぎ強度の高さ方向の積分値によって定義される量である。また、大気の構造が同じでも観測波長が長いほどフリード長は大きくなる。一般的に大気ゆらぎの強度を議論する上では可視光 (λ = 500nm) かつ天頂方向 (γ = 0) におけるフリード長を用いることが多い。天文観
測サイトにおけるフリード長は日本国内 (岡山) では 7cm 程度 [18]、マウナケア山頂では 13.5cm程度である [19]。

2. シーイング (seeing) θseeing
シーイングは、大気ゆらぎの影響によって広がった点像分布関数のサイズ (FWHM) によって大気ゆらぎの強度を評価するパラメータである。フリード長を用いて以下のようにかける。

θseeing = 0.98

λr0(2.20)

天文観測サイトにおけるシーイングは日本国内 (岡山) では 1.5 秒角程度 [18]、マウナケア山頂では0.75 秒角程度である [19]。
3. アイソプラナティック角 (isoplanatic angle) θ0
天球面上で、大気ゆらぎによる波面の乱れが同一であると見なせる角度領域をアイソプラナティッ
ク角と呼ぶ。以下の式で表される。

[

2.91k 2θ0 =(cosγ)8/3

]− 53

∫
h

5/3

2
CN
(h)dh

(2.21)

このパラメータによって、1 つのガイド星でどの程度の領域の補償ができるかを見積もることができる。

4. コヒーレンス時間 (coherence time) τ0
コヒーレンス時間は、大気ゆらぎによる波面の乱れが同一と見なせる時間間隔を表現するパラメー
タであり、フリード長を用いて以下の式で表される。

τ0 = 0.31

r0
V0(2.22)

ここで、V0 は大気ゆらぎの平均移動速度であり、

[∫
V0 =

2
V (h)5/3 CN
(h)dh
∫ 2
CN (h)dh] 35(2.23)である。コヒーレンス時間を計算するには大気ゆらぎ強度の高さ分布だけではなく大気ゆらぎの移
動速度の高さ分布も必要になるが、このパラメータによってどれほどの速度で補償の制御ループを回す必要があるかを見積もることができる。

第2章

大気ゆらぎの高さ分布の推定

202.2 大気ゆらぎの高さ分布の推定手法
大気ゆらぎの高さ分布の推定にはいくつか確立された手法が存在する。口径の大きな望遠鏡が登場し大気ゆらぎの影響が議論されるようになった時代には、各サイトの大気を定量的に評価するためにすでに必要だったからである。前節で見たように、アイソプラナティック角やコヒーレンス時間を計算するためには大気ゆらぎの高さ分布の情報が必要である。ここではいくつかの伝統的な手法を、その特徴や問題点とともに紹介してゆく。

2.2.1 バルーンを用いた直接測定
バルーンを飛ばすことによって直接的に大気ゆらぎの強度を測定する実験は、Barletti et al.1977[20]、Coulman et al.1995[21]、Azouit & Vernin2005[22]、Mchugh et al.2008[23] などによって行われている。
この手法のメリットは、直接的に測定できるので推定誤差が小さく、高さ方向に非常に細かい分解をした分布が得られることである。一方で、1 回の高さ分布測定に長い時間がかかること、天文観測と同時に行うことの技術的難しさから観測方向におけるリアルタイムな高さ分布推定には適さない手法である。

図 2.2 Azouit & Vernin2005[22] の気球実験によって得られたチリの Cerro Pachon における大気2
ゆらぎの高さ分布。横軸が大気ゆらぎ強度を表す CNの値で縦軸が標高を示している。

図 2.2 は、Azouit & Vernin2005[22] において得られたチリの Cerro Pachon における大気ゆらぎの高さ分布である。非常に細かい高さ方向の分解能で測定できている。また、この分解能で見てもゆらぎ強度
分布の概形は滑らかではなく、ある高度で急激に大気ゆらぎ強度が強くなる部分があることから大気ゆら第2章

大気ゆらぎの高さ分布の推定21

ぎが非常に薄い層状に存在していることが分かる。

2.2.2 2 つの波面センサーの相関を用いる方法
大気ゆらぎに関するリアルタイムな情報は補償光学系の波面センサーに各時刻の波面ゆらぎとして刻まれている。波面センサーデータから大気ゆらぎの高さ分布を復元する手法が SLOpe Detection And

Ranging(SLODAR; Wilson 2002[25]) である。この手法の基本的な原理の概念図が図 2.3 である。

図 2.3 SLODAR の概念図。2 つの波面センサーで小さな角度 θ 離れた 2 つの星をそれぞれ観測する。2 つの波面センサーデータの相互相関を取るとき、大気ゆらぎの高度に応じたずらし量に対して大気ゆらぎの強度に応じた相関のパワーが測定される。

図 2.3 に示すように、SLODAR では二重星や 2 つのレーザーガイド星といった小さな角度 θ だけ離れた 2 つの星を 2 つの波面センサーでそれぞれ観測する。この時、同一の大気ゆらぎによる影響を受けた
別々の天体からの光は、大気ゆらぎの高さに応じた距離 hθ だけ離れて開口に入る。したがって 2 つの波面センサーの測定量（星像の変位）の相互相関を取ることで大気ゆらぎの高さごとの情報を引き出すことができる。相関のパワーが大気ゆらぎ強度に比例することから、大気ゆらぎ強度の高さ分布を得ることができる。最も重要なのは、これらの手法は補償光学で用いられる波面センサーで行われるということである。観測と同じ方向の高さ分布推定が得られる、特別な装置を必要としないなどの点で非常に有用であ第2章

大気ゆらぎの高さ分布の推定

22

り、広く用いられている手法の 1 つである。
一方でデメリットも存在する。この手法によって得られる高さ方向の分解能は Wilson 2002[25] によ
ると、

δH =

D

(2.24)

nsub θ

である。ここで、D は望遠鏡の口径、nsub はマイクロレンズアレイによる望遠鏡の口径の分割数、θ は二重星の離角である。また、感度のある最も高いゆらぎ高度は、

Hmax = nsub δH =D
θ(2.25)と表される。例えばすばる望遠鏡の AO188 に用いられている波面センサーの場合は、

D = 8.2[m]nsub = 12(2.26)(2.27)(2.28)であるので、LTAO(θ ∼ 20”) の場合で δH = 7.0[km], Hmax = 85[km] であり、WFAO(θ = 10′ ) の場合で δH = 0.23[km], Hmax = 2.8[km] となる。したがって 1 章で述べた δH < 1[km] かつ Hmax > 20[km]の条件を満たすことはできない。このように、波面センサーの相関を用いる手法においては参照光源となるレーザーガイド星の離角が補償光学の要求から決まってしまっているため、大気ゆらぎの高さ分布に必要とされる高さ方向の分解能や最大の高度を得ることができなくなる。

2.2.3 MASS-DIMM
MASS-DIMM は、MASS(Multi Aperture Scintillation Sensor) と DIMM(Diﬀerential Image Motion
Monitor) という 2 つの手法を組み合わせた手法である。簡単にそれぞれの手法について説明すると、MASS は単一の星像の明るさゆらぎの情報から大気ゆらぎ強度の高さ分布を推定できる手法であり、DIMM は単一の星像の位置ゆらぎの情報から大気ゆらぎ強度の積分値 (フリード長やシーイング) を推定できる手法である。両者ともに単一の星の観測で推定ができることが大きなメリットである。MASS の手法では各高度における大気ゆらぎ強度がわかるが、∼0.5km よりも低い大気ゆらぎについては原理上感度が悪くなる。一方 DIMM の手法では大気ゆらぎ強度を高さ方向に足し合わせた値が求められる。したがって DIMM の結果から MASS の結果を差し引くことで地表近くの強度成分を求めることができ、両者の組み合わせによって地表層を含めた高さ分布が得られる。この 2 つの手法を 1 つの装置で実践できるシステムが Kornilov et al.2007 [26] によって開発されており、現在世界各地で観測サイト調査やシー
イングモニターとして使われている。図 2.4 に示すのが、実際の MASS-DIMM 装置の写真である。写真[A] からわかる通り装置自体は 20cm 程度の非常にコンパクトなものである。写真 [C] のように、小型の望遠鏡に設置して用いられる。

第2章

大気ゆらぎの高さ分布の推定

23

図 2.4 [A] 観測サイト調査やシーイングモニターとして広く用いられている MASS-DIMM 装置。

[B] 比較的安価に制作することができるため数多く作られ様々なサイトで用いられている。[C] 口径
20cm 程度の小さな望遠鏡に取り付けて運用される。非常にコンパクトな装置であることがわかる。

ここからはそれぞれの手法の原理の詳細な説明を行う。

MASS(Muti Aperture Scintillation Sensor) は、Kornilov et al. 2003[13]、Tokovinin et al. 2003a[27]によって提案された方法である。MASS では、大気ゆらぎに起因する星の明るさゆらぎの情報から各高度の大気ゆらぎの強度を推定する。まず、大気ゆらぎに起因する明るさゆらぎのパワースペクトラムは式

2.18 から、ΦI (fx , fy ) =

∑
layer

4 × Φχ =

∑
j

1.53f

−11/3

{

sin(πλhj f 2 )
λ

}2
2
CN
(hj )∆hj(2.29)であった。ここで j は大気層の数に対応するインデックスである。今、大気ゆらぎが 1 つの層のみから2なっていて、ゆらぎ強度が CN(hj )∆hj = 6.9 × 10−13 [m1/3 ]（観測波長 500nm でのフリード長が 10cmに対応する。）であるとし、大気ゆらぎの高度を変えたときの明るさゆらぎのパワースペクトラムを描くと、図 2.5 のようになる。ここから分かるのは、高い大気ゆらぎほど明るさゆらぎのパワーへの寄与が大きいということである。また、パワースペクトラムのピークとなる点は、ゆらぎ層の高度が高くなるほど低周波数側へ移動してゆくことがわかる。言い換えると、様々な空間スケールで星の明るさ変動を測定す
ることによってどの高さにどれくらいの強度を持った大気ゆらぎが存在するかを推定することが可能で
ある。

第2章

大気ゆらぎの高さ分布の推定

100

power

10 1

24

height of turbulent layer : 1 km
height of turbulent layer : 10 km
height of turbulent layer : 20 km

10 2
10 3
10 4
10 5 100

図 2.5

101
spatial frequency [m^-1]

102

大気ゆらぎに起因する星の明るさゆらぎのパワースペクトラム。ある高度に 1 層のみ大気ゆらぎがあった場合の結果を大気ゆらぎの高度ごとに色を変えて示している。

図 2.6 MASS 手法の概念図 [28]。円環状に区切られた 4 つの開口でそれぞれに測光を行う。これは明るさゆらぎのパワースペクトラムに対して、特定の空間周波数成分を取り出していることに相当
する。

第2章

大気ゆらぎの高さ分布の推定25

MASS はこの原理を利用して大気ゆらぎの高さ分布を推定する。MASS 手法の概念図を図 2.6 に示す。MASS では単一の星に対して、円環状に区切られた 4 つの開口でそれぞれに測光を行う。4 つの開口

A,B,C,D の直径は主鏡上での長さに換算してそれぞれ、2.0cm,3.7cm,7.0cm,13cm である。例えば、図2.6 の開口 A は 4 つの中で最も小さな開口であり、高い空間周波数成分を取り出す。一方で開口 D は空間的に離れた円環形の開口であり、低い空間周波数成分を取り出す役割をしている。実際の運用としては
図 2.4[A] に示されているように、円環状に切り分けられたミラーが異なる傾きで配置されている。それぞれのミラーで反射した光が別々の検出器に入り測光されるようになっている。
このとき数百-千 Hz の高速撮像を行うことによって大気ゆらぎによる星の明るさ変動データを取得する。この明るさ変動データから以下の式で定義されるシンチレーションインデックス (Scintillation

Index; SI) を計算する。
[

IA
SA = V ar< IA >](2.30)[IAIB,
SAB = SA + SB − 2Cov
< IA > < IB >](2.31)

ここで、I は星の明るさ (単位時間あたりの光子数) であり添字の A や B は円環状に区切られた開口を表す。<>, V ar[], Cov[] はそれぞれ時間変動データから得られる平均値、分散値、共分散値を表している。SA のような、1 つの開口で測定される明るさ変動の分散に対応する量を Normal Scintillation

Index(NSI)、SAB のような、2 つの開口で測定される明るさ変動の共分散に対応する量を Diﬀerential
Scintillation Index(DSI) と呼んで区別する。MASS では 4 つの開口 A,B,C,D を用いるため 4 種類のNSI および 6 種類の DSI が一度の観測データから計算される。明るさ変動の測定は 4 つしかないが、DSIによって与えられる情報を含めることで合計 10 個の独立な測定量を取り出すことができ、その分推定できる大気ゆらぎの高さも増える。実際には 0.5km,1km,2km,4km,8km,16km の 6 つの高さに対して推定を行うことが可能である。
一方で数学的事実として、ゆらぎの分散値はゆらぎのパワースペクトラムの積分値によって与えられ、2 種類のゆらぎの共分散値はゆらぎ同士のクロススペクトラムの積分値によって与えられる。星の明るさ
ゆらぎのパワースペクトラム ΦI は式 2.18 で求められていたので、これらを踏まえると SI は、

∫ ∫
SA =
∫ ∫
SAB =

ΦI (fx , fy ) × AA (fx , fy )dfx dfy(2.32)

ΦI (fx , fy ) × AAB (fx , fy )dfx dfy(2.33)と書くことができる。ここで、A(fx , fy ) は開口の形状 (ここでは円形や円環形) による周波数フィルターであり、以下のように表される。

AA (fx , fy ) = |F[PA (x, y)]|2(2.34)AAB (fx , fy ) = |F[PA (x, y) − PB (x, y)]|2(2.35)

ここで P (x, y) は各開口形状の開口関数 (開口部を 1、遮蔽部を 0 とした関数) であり F[] はフーリエ変換を表す。式 2.32-式 2.35 からわかるように、SI の値は明るさゆらぎのパワースペクトラムと開口形状の第2章

大気ゆらぎの高さ分布の推定

26

情報から計算できることがわかる。逆に、SI の値を観測から測定することができれば明るさゆらぎのパ
ワースペクトラムに含まれる各高度の大気ゆらぎ強度の情報を推定することができる。
定式化をしてゆく。式 2.32, 式 2.33 は同じ形式となっているためここからは添字 i のみを用いて開口を区別してゆく。i ∈ (A, B, C, D, AB, AC, AD, BC, BD, CD) といった具合である。式 2.18, 式 2.32から、

}2sin(πλhj f 2 )2CN(hj )∆hj × Ai (fx , fy )dfx dfy
λ
j
[∫ ∫
]
{
}2
∑
sin(πλhj f 2 )
−11/3
2
=
1.53f
× Ai (fx , fy )dfx dfy CN(hj )∆hj
λj∑=Wi,j Jj

Si =

∫ ∫ ∑

1.53f −11/3{(2.36)(2.37)(2.38)

j

ただし、

∫ ∫
Wi,j =

1.53f

−11/3

{

sin(πλhj f 2 )λ}2
× Ai (fx , fy )dfx dfy

2
Jj = CN(hj )∆hj(2.39)(2.40)となる。SI と各高度の大気ゆらぎ強度の関係式は連立 1 次方程式の形に書くことができることがわかった。行列の形で書くと、

SA SB 
 SC  
WA,0.5km
 SD  W
 SAB  
B,0.5km..
 SAC  = .
  SAD WCD,0.5km
 SBC 
 SBD SCD

WA,1km···


J0.5kmWA,16km  J1km 
J2km   J4km  J8km 
J16km


(2.41)となる。今 Jj は高さ hj の大気ゆらぎ層の強度に層の厚みをかけたものであり、このパラメータも大気ゆ
らぎ強度の指標としてよく用いられる。Wi,j は Jj から Scintillation Index Si へと変換する係数行列となっており、式 2.39 から分かるように開口形状の情報のみから求められる。MASS 手法は、様々な開口形状による星の明るさ観測から求めた Si と開口形状の情報から計算される Wi,j から、式 2.41 に表される連立 1 次方程式をといて各高度の大気ゆらぎ強度 Jj を求める手続きである。
次に DIMM の手法の説明を行う。DIMM(Diﬀerential Image Motion Monitor) は、Sarazin & Roddier

1990[29] によって体系化された方法である。DIMM は、大気ゆらぎに起因する星の位置ゆらぎの情報から大気ゆらぎ強度の高さ方向の積分値（フリード長やシーイング）を推定する。DIMM 手法の概念図を
図 2.7 に示す。

第2章

大気ゆらぎの高さ分布の推定

27

図 2.7 DIMM 手法の概念図 [30]。単一の星を一定距離離れた 2 つの開口で観測することによって 2つの星像を得る。これらの位置の時間変動を測定することによって、開口間距離に対応する空間周波数の変動成分を取り出すことができる。測定した変動のパワーを理論的なパワースペクトラムとフィッ
ティングすることによって大気ゆらぎ強度の測定が可能である。

DIMM では 2 つの開口で、ある単一の星を観測する。DIMM の測定量は検出器上に得られる 2 つの星像の相対的な位置変動である。MASS 同様に数百-千 Hz の高速撮像を行うことによって大気ゆらぎによる星の位置変動データを取得する。図 2.7 の左側の開口 A による星像の位置を (xA , yA )、右側の開口 B
による星像の位置を (xB , yB ) とするとき、

σx2 ≡< (xA − xB )2 > [rad2 ](2.42)

σy2 ≡< (yA − yB )2 > [rad2 ](2.43)が測定される。ここで <> は時間変動データから得られる平均値を表している。また、開口 A,B がx 方向に距離 d だけ離れており、開口面上の位置 (x, y) における x, y 方向の波面の傾きをそれぞれ
α(x, y), β(x, y) とするとき、これらの値は
σx2 =< (α(0, 0) − α(d, 0))2 >= Dα (d, 0)(2.44)

σy2 =< (β(0, 0) − β(d, 0) >= Dβ (d, 0)(2.45)とも書くことができる。DIMM の測定量は波面の傾きの構造関数の値であることがわかる。
一方で構造関数はパワースペクトラムからも導出することができる。まず、大気ゆらぎに起因する位置
ゆらぎのパワースペクトラムは式 2.16、2.17 から、

Φα (fx , fy ) =∑λ2 fx2 Φψ =Φβ (fx , fy ) =

layer

2
9.7 × 10−3 (2π)2 fx2 f −11/3 cos2 (πλhj f 2 )CN(hj )∆hj(2.46)

2
9.7 × 10−3 (2π)2 fy2 f −11/3 cos2 (πλhj f 2 )CN(hj )∆hj(2.47)

j

layer

∑

∑

λ2 fy2 Φψ =

∑
j

第2章

大気ゆらぎの高さ分布の推定28であった。ここで j は大気層の数に対応するインデックスである。このパワースペクトラムの式をフーリエ変換することで以下の位置ゆらぎの相関関数 Bα (x, y), Bβ (x, y) が得られる（Wiener-Khinchin の定理）。

λ2 ∂ 2Bψ (x, y)
4π 2 ∂x2
λ2 ∂ 2Bβ (x, y) = − 2 2 Bψ (x, y)4π ∂yBα (x, y) = −(2.48)(2.49)

さらに式 2.5 を用いることで、

λ2 ∂ 2Dψ (x, y)
8π 2 ∂x2
λ2 ∂ 2Dψ (x, y)Bβ (x, y) =8π 2 ∂y 2

Bα (x, y) =(2.50)(2.51)とかける。ここで、Dψ (x, y) は地上に伝搬後の位相ゆらぎの構造関数である。位相ゆらぎの情報は低周
波数側では伝搬の前後で変化しない（理由は図 2.8 を参照のこと）。

101

after 1 km propagation
after 10 km propagation
after 20 km propagation
before propagation

100

power

10 1
10 2
10 3
10 4
10 5 100
図 2.8

101
spatial frequency [m^-1]102

光の位相ゆらぎのパワースペクトラム。色の違いは伝搬距離の違い (=大気ゆらぎの高度の違い) を表している。5[m−1 ] よりも低空間周波数帯においては位相ゆらぎのパワーは伝搬距離によらな
いことがわかる。

したがって伝搬前の位相ゆらぎの構造関数 Dϕ (x, y) を用いて、λ2 ∂ 2Dϕ (x, y)
8π 2 ∂x2
λ2 ∂ 2Bβ (x, y) ∼Dϕ (x, y)
8π 2 ∂y 2

Bα (x, y) ∼(2.52)(2.53)

第2章

大気ゆらぎの高さ分布の推定29とかける。大気ゆらぎを通過することによって生じる位相ゆらぎの構造関数は Tatarski 1971[31] によって求められており、フリード長 r0 を用いて(√
Dϕ (x, y) = 6.88

x2 + y 2
r0)5/3(2.54)であるので、Bα , Bβ は、

[
]
1 2 2
2
2 −1/6
2 −7/6(x + y )− x (x + y )
3
[
]
1 2 2
2
2 −1/6
2 −5/3
2 −7/6(x + y )− y (x + y )Bβ (x, y) ∼ 0.145λ r03−5/3Bα (x, y) ∼ 0.145λ2 r0(2.55)(2.56)とかける。再び式 2.5 を用いれば DIMM の測定量は、

Dα (d, 0) = 2[Bα (0, 0) − Bα (d, 0)](2.57)Dβ (d, 0) = 2[Bβ (0, 0) − Bβ (d, 0)](2.58)と表される。
−5/3 −1/3Bα (d, 0) = 0.0968λ2 r0d(2.59)−5/3 −1/3Bβ (d, 0) = 0.145λ2 r0d(2.60)である。なお Bα (0, 0), Bβ (0, 0) の値は式 2.55,2.56 を見ると無限大に発散してしまうことがわかる。しかし実際の DIMM の開口は無限小の「点」にはなっておらず有限の開口サイズで平均した値となるので発散しない。Bα (0, 0), Bβ (0, 0) の実効的な値は Fried 1965[32], Fried 1975[33], Tatarski 1971[31] によって調べられており、−5/3

Bα (0, 0) = Bβ (0, 0) = 0.179λ2 r0

D−1/3(2.61)

である。最終的に、DIMM の測定量は[−5/3

σx2 = 0.358λ2 r0

D−1/3 1 − 0.541[−5/3 −1/3
σy2 = 0.358λ2 r0
D

1 − 0.810

(
(

d
D
d
D

)1/3 ](2.62))1/3 ](2.63)となる。DIMM 手法は、2 つの星像の相対的な位置変動の測定と式 2.62,2.63 からフリード長 r0 を求める手続きである。

第2章

大気ゆらぎの高さ分布の推定

302.2.4 各手法のまとめ
ここまでで紹介した大気ゆらぎの高さ分布を測定する既存の手法の特徴ついて、1 章で述べたトモグラ
フィック波面推定に要求される項目ごとに表 2.1 にまとめる。このように、既存の大気ゆらぎの高さ分布
高度範囲

高度分解能

即時性

観測と測定の方向

0-20km をカバー

∼1km

あり

一致

0-30km

∼10m

なし

不一致

波面センサーの相関を用いる方法

0-85km(LTAO)

∼7km

なし

一致

波面センサーの相関を用いる方法

0-2.8km(WFAO)∼0.23kmなし

一致

0-16km

6 種類の高度

あり

不一致

要求
バルーンを用いた直接測定

MASS-DIMM

表 2.1 大気ゆらぎの高さ分布推定に用いられる既存の手法。

推定手法ではトモグラフィック波面推定の事前プロファイルに要求される性能を満たすものはないことが
わかる。

31

第3章

波面センサーを用いた MASS-DIMM の提案
3.1 Shack-Hartmann 型波面センサーを用いた MASS-DIMM 法の提案
ここまで、1 章では複数のレーザーガイド星を用いた補償光学には大気ゆらぎの高さ分布についての事
前情報が必要であることを述べた。2 章では既存の測定手法を紹介し、それらの方法では事前プロファイ
ルに求められる条件を満たさないことを述べた。本章で述べるのは本研究で新たに検討した手法についてである。それは、Shack-Hartmann 型波面センサーを用いた MASS-DIMM 手法である。MASS-DIMM法のメリットは 1 つの星の観測からプロファイルを求めることで高い高度のゆらぎが推定できなくなる問
題が生じないことである。一方で専用の装置を必要とし、それに伴いトモグラフィ補償光学で推定すべき
大気の方向と異なった方向の事前情報を得てしまうというデメリットもあった。そこで波面センサーを用いて同手法を実践することができればこの問題を解決できると考えた。実際、MASS では様々な開口パ
ターンにおける明るさ変動の観測が必要となるが、これは波面センサーのサブ開口の選び方によって再現できる。さらに、従来の MASS では 4 つの開口パターンしか持たないためそれらに対応する 6 つの高度のプロファイルしか得られなかった。波面センサーのサブ開口の選び方によって従来よりも多い開口パ
ターンを取ることができるため高度方向の分解能も上がることが期待できる。また、DIMM では 2 つの開口による像の位置変動の観測が必要となるが、これもどの 2 つのサブ開口を選ぶかによって再現できる。波面センサーと MASS-DIMM 手法とを組み合わせることによって事前プロファイルの要求を満たす
ことができると考え、手法の検証を行った。ここからは波面センサー用いて MASS-DIMM を行うためにはどのような光学系のパラメータが必要になるかについて検討してゆく。

3.2 MASS を行うためにMASS 手法において様々な開口形状で明るさゆらぎを観測することは、大気ゆらぎの高度によって変わる特定の空間周波数成分を取り出していることに相当する。式 2.13 で示した明るさゆらぎのパワース
ペクトラムを再掲すると
2
Φχ (fx , fy ) = 9.7 × 10−3 k 2 f −11/3 sin2 (πλhf 2 )CN
∆h

(3.1)

第3章

波面センサーを用いた MASS-DIMM の提案

32である。したがって、パワースペクトラムのピークを取る空間周波数は2πλhfpeak=π2√

fpeak =

(3.2)
1
2λh

(3.3)

で与えられる。この周波数の逆数と大気ゆらぎの高さ h との関係をいくつかの観測波長 λ に対してプ
ロットしたのが図 3.1 である。縦軸が大気ゆらぎの高さ、横軸がパワースペクトルのピークを取る空間ス
ケールであり、その高さの大気ゆらぎが作る明るさゆらぎの代表的な相関長と言える。大気ゆらぎの典型的な存在領域である高度 0-20km を考えると、明るさゆらぎの相関長は 0-15cm に現れる。このことは実際の MASS 装置がそうであったように、波面センサーによる瞳の分割を細かくして数 cm 程度の小さな
空間スケールでサンプリングする必要があることを意味している。例えばサブ開口のサイズが 5[cm] の場
合は MASS が感度をもつ最も低い高度はおよそ 2.5[km]、サブ開口のサイズが 2.5[cm] の場合はおよそ

height of turbulent layer [km]1[km] である。

図 3.1

20.0
17.5
15.0
12.5
10.0
7.5
5.0
2.5
0.0

observation wavelength : 500 nm
observation wavelength : 600 nm
observation wavelength : 700 nm

0.0

2.5

5.0 7.5 10.0 12.5
1/(peak frequency) [cm]

15.0

17.5

大気ゆらぎの高度と、大気ゆらぎによって作られる明るさゆらぎのパワーが最大となる空間周波数との関係。大気ゆらぎの高度が高いほど明るさゆらぎの相関長は大きくなるが、最も高い大気ゆらぎ（高度 20km）を考えてもその相関長は 20cm 以下である。

一方で、小さな空間スケールでのサンプリングを行うと 1 つのサブ開口あたりに入ってくるガイド星からの光子の数が減ってしまい検出が難しくなる。そこで次にサブ開口の空間スケールに対して、どれくら
い明るいガイド星が必要となるのかについて考えてゆく。すばる望遠鏡の AO188 に用いられている波面センサーでは 8.2m の主鏡を 188 に分割しており（実際には等分割ではない）、その場合で mR < 13 の明
るい星がガイド星として用いられる。この明るさは、波面センサーにおいて光を分割した場合においても

第3章

波面センサーを用いた MASS-DIMM の提案

33

像の位置変動を測定するのに十分な SN 比が得られるために必要な要求となっている。したがってサブ開口の口径を D[cm] とした時、必要となる R バンド等級は、(
13 − 2.5log

D28002 /188)(3.4)と見積もることができる。これをプロットしたのが、図 3.2 である。サブ開口のサイズが 5.0[cm] の場
合にはおよそ 8 等のレーザーガイド星、2.5[cm] の場合にはおよそ 6 等のレーザーガイド星が必要となる。サブ開口のサイズが 2.5[cm] よりも小さくなってくると必要となるガイド星の明るさは急激に明るくなる。実際の補償光学においてレーザーによって作ることのできる星の明るさは、例えばすばる望遠鏡の現在のシステムの場合では 14 等程度 (レーザーガイド星システム運用開始当初は 11 等)[1] であるが、

ULTIMATE-SUBARU プロジェクトへ向けたレーザーのアップグレードによって 8-10 等まで上がる予定である。したがって、レーザーガイド星 + 波面センサーによる MASS を行う場合、その瞳分割数はおよそ 5[cm] よりも大きい必要がある。

required guide star's brightness

R-band magnitude

10
8
6
4
2
0
0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5 20.0
subaperture size [cm]
図 3.2

以上を踏まえるとすばる望遠鏡のアップグレードされたレーザーガイド星 + 波面センサーを用いてMASS を行うために必要な波面センサーの分割数は 5.0[cm] よりも粗いサンプリングであり、この時
2.5[km] よりも高い大気ゆらぎに感度をもつ MASS を実現することができる。

3.3 DIMM を行うために
DIMM ではその原理上、星像の位置ゆらぎの低空間周波数成分を観測する必要があった。より具体的には 20km の大気ゆらぎを考えた場合 5[m−1 ] よりも低い空間周波数において、位置ゆらぎのパワースペ

第3章

波面センサーを用いた MASS-DIMM の提案

34

クトラムが大気の位相ゆらぎのパワースペクトラムと一致することを用いる。したがって、およそ 20cm
よりも離れた 2 点の測定点が必要となる。DIMM を行う場合、必要となる望遠鏡の口径は 20cm よりも大きい必要がある。

35

第4章

東北大学 50cm 望遠鏡と波面センサーによる星の観測
4.1 光学系
３章で検討した波面センサー光学系で MASS-DIMM 手法が実践可能であることを検証するために、波面センサー光学系を設計し、東北大学 50cm 望遠鏡を用いて星の観測を行なった。本章では設計した波面
センサー光学系や観測の詳細について述べる。

4.1.1 波面センサー光学系
今回設計・製作した波面センサー光学系の構成を図 4.1 に示す。

図 4.1 設計・製作した波面センサーの構成。赤色の領域は幾何光学的な光線の束である。

このような光学系の配置は波面センサーとしては最も基本的なものとなっている。望遠鏡から出てきた
光をコリメータレンズで平行光にし、望遠鏡の瞳に対して共役な位置に配置されたマイクロレンズアレイ
によって多数の像を作る。最後にリレー光学系と呼ばれる 2 枚のレンズによって F 値を変えながら像をカメラの検出器に再結像するというものである。コリメータの焦点距離によってコリメート光のビームサ
イズが決まる。そのビーム径をマイクロレンズアレイ (MLA) で何分割するかによって瞳の分割サイズが決まる。さらにリレーレンズ系を構成する 2 枚のレンズの焦点距離の比によって検出器面上のどれほどの第4章

東北大学 50cm 望遠鏡と波面センサーによる星の観測

36

領域に像が広がり、1 つ 1 つの像がどれほどの大きさになるかが決まる。したがって今回の波面センサーを特徴付けるパラメータは• コリメータの焦点距離 fcol
• マイクロレンズアレイのピッチ pMLA
• マイクロレンズの焦点距離 fMLA
• リレーレンズ系の焦点距離の比 Rrelayである。これらのパラメータの値を決めるための条件について説明する。まず 3 章での結論から波面センサーによる MASS を実現するための瞳分割はサブ開口の主鏡上での直径 Dsubap が Dsubap > 5.0cm であった。しかし、この条件はすばる望遠鏡でのアップグレードされたレーザーガイド星を使用することが前提にある。今回の実験は東北大学にて行い、観測天体の明るさは自由に選ぶことができることを考慮し、より細かい瞳分割での MASS 実践を検証することにした。３章での検討結果から、観測天体の明るさの制限がない場合サブ開口のサイズが 2.5cm で 1km より高い大気ゆらぎの測定ができるため、今回の実験では、

Dsubap = 2.5[cm](4.1)とした。Dsubap は、望遠鏡の (有効) 焦点距離 ftel 、コリメータの焦点距離 fcol 、マイクロレンズアレイのピッチ pMLA を用いて、

Dsubap =

ftel
pMLA
fcol

(4.2)

とかける。必要なコリメータ焦点距離とマイクロレンズアレイピッチの比は、

fcol
ftel
=
pMLA
Dsubap(4.3)

で計算される。
また、リレーレンズによってスポット像同士の間隔 dMLA は、

dMLA = pMLA × Rrelay(4.4)に拡大される。DIMM を行うために必要な測定距離は ddimm = 20[cm] であったが、今回のサブ開口がDsubap = 2.5[cm] であることから ddimm /Dsubap = 8 個以上のサブ開口が検出器に収まらないといけな
い。検出器の物理的なサイズを Ldet × Ldet とすると DIMM を行うための条件は、

ddimm< LdetDsubap
Ldet Dsubap
pMLA Rrelay <ddimm

pMLA Rrelay ×

(4.5)
(4.6)

となる。
一方でスポット像の広がり dspot は観測波長 λ、マイクロレンズアレイの焦点距離 fMLA を用いて、dspot = 1.22λ

fMLA
Rrelay
pMLA

(4.7)

第4章

東北大学 50cm 望遠鏡と波面センサーによる星の観測37と表される。この値がスポット像同士の間隔 dMLA に比べて十分小さくないとスポット像の裾野同士が重なってしまい、それぞれのスポット像の明るさの測定が難しくなる。したがって、

dspot << dMLA

(4.8)
(4.9)

すなわち、

p2MLA
>> 1.22λfMLA(4.10)

も重要な条件である。
今回用いた望遠鏡、検出器、フィルターの仕様から、

ftel = 6025[mm](4.11)

Ldet = 3.12[mm](4.12)

λ = 600[nm](4.13)であるので、式 4.3,4.6,4.10 の条件を満たすような光学系パラメータとして、

fcol = 35[mm](4.14)

pMLA = 150[µm](4.15)fMLA = 4.67[mm](4.16)

Rrelay = 1.18

(4.17)とした。採用した光学系を表 4.3 に示す。
光学系

焦点距離 [mm]

ピッチ [mm]

型番

コリメータ

35

なし

Thorlabs 社製 AC254-035-A

光学フィルター

なし

なし

不明

マイクロレンズアレイ

4.67

0.15

Thorlabs 社製 MLA150-5C-M

リレーレンズ 1

75

なし

Thorlabs 社製 AC254-075-A1

リレーレンズ 2

88.3

なし

mks newport 社製 KBX061

表 4.1

備考

Bessell R

波面センサー光学系に用いたレンズの仕様。

これらの光学系の互いの距離については、光学シミュレーションソフトの ZEMAX を用いて検出器上での像のサイズが最も小さくなるように最適化して決めた。その情報をもとに実験室での組み上げ、調整を行なった。図 4.2 に最適化した光学系による ZEMAX 上でのシミュレーションスポットイメージを
示す。

第4章

東北大学 50cm 望遠鏡と波面センサーによる星の観測38

図 4.2 ZEMAX によってシミュレートされた検出器上でのスポット像のイメージ。正方形の領域が今回の検出器面を表している。

4.1.2 EMCCD
波面センサーの検出器には研究室で保有している EMCCD(Electron Multiplying Charge CoupledDevice) を用いた。以下 Robbins et al.2003[34] にしたがって説明する。通常の CCD カメラの感度は、読み出しノイズ（光電効果によって発生した電荷を電圧値として読み出す際に印加されるノイズ）によって制限される。読み出しノイズは読み出しのレートが高くなるほど大きくなる。従来の CCD では読み
出しレートを遅くすることによってその高い感度を実現している。一方で高速撮像には適していない。

EMCCD は従来の CCD よりも長い電荷転送部をもち、その中で 40-45[V] ほどの高い電圧をかけること
によって衝突電離を起こし電子数を増倍する。1 つの増倍レジスターにおいて衝突電離が起こる確率は 1
つの電子あたりおよそ 1.5% 程度であるが、増倍レジスターを何段階にも重ねることによって大きな像倍率が得られる。これによって少ない光子数でも実効的に低い読み出しノイズで測定できる。EMCCD による実効読み出しノイズは以下のように表される。

√
σef f =

F 2 (S + Sdark ) +2
σreadout
M2(4.18)

ここで S は天体や空からのシグナル、Sdark は検出器の暗電流によるシグナル、σreadout は読み出しノイズ、M は電子の増倍率、F は電子増倍時のノイズファクターである。
今回用いた EMCCD の各値について述べる。電子増倍率は過去の測定実験の結果から電子増倍用ク
ロックの電圧と検出器の温度、電子増倍率の関係が分かっておりこれを用いた。今回の観測では電子増倍用クロックの値は 43.00242[V]、検出器の温度は −30[℃] とし、このときの電子増倍率はおよそ 300 である。読み出しノイズについても過去の測定による σreadout = 11[ADU ] を用いた。測定されていなかった

第4章

東北大学 50cm 望遠鏡と波面センサーによる星の観測

39

ノイズファクターについては Robbins et al. 2003[34] をもとに理論的な値を用いた。これによると、増倍率 M および増倍レジスターの段階数が大きいとき F は√2 に収束する。また量子効率（光子数から電子数への変換ファクター）は 1 を仮定し、電圧値からカウント値 (ADU) への変換ファクターは過去の測定から 18.5[electron/ADU] を用いた。

4.2 東北大学 50cm 望遠鏡を用いた観測
観測は 2019 年 10 月 16 日に行なった。観測天体にはデネブ（mR = 1.14[35]）を用いた。500[frame/sec]
で 60 秒間の高速撮像を 1 回のデータ取得とし、1 時間の間に 8 回のデータ取得を行なった。観測時の光学系の様子を図 4.3[A] に示す。また、図 4.3[B],[C] は EMCCD カメラに取り付けられた波面センサー光学系の様子である。

図 4.3 [A]50cm 望遠鏡に取り付けられた波面センサー光学系と EMCCD カメラの様子。[B][C] 筒状のレンズホルダーに波面センサー光学系に用いられるレンズが並べられており、銀色の EMCCD カ
メラにねじ込み式に取り付けられるようにした。

40

第5章

新手法の実践および結果
ここからは、実際に取得した Shack-Hartmann 型波面センサーのデータを用いて新しく提案する手法を実践した過程および結果について述べてゆく。

5.1 星の変位・明るさの時間変動
取得したデータから星の位置と明るさの時間変動の情報を取り出した解析について説明する。図 5.1 左に示すのが 1 分間星を撮像して取得した 30000 フレームのシャックハルトマン画像をバイアスフレームを減算したのち平均化したものである。主鏡が光学系の瞳になっているため全体は円形であり中心部分は副鏡による遮蔽が見えている。全体的に左側にシフトしているのは望遠鏡と波面センサーからなる全体の光学系の光軸中心と検出器の中心位置とがずれてしまっていることが原因と考えられる。

1min(30000frames) mean

0

6000

6000

40

5000

40

5000

60

4000

60

4000

80

3000

80

3000

2000

100

1000

120
0

図 5.1

7000

20

20

40

60
80
y[pixel]100

120x[pixel]x[pixel]20

Region used for spot detection

0

7000

2000

100

1000

120
0

20

40

60
80
y[pixel]

100

120

0

左：1 分間の星の高速撮像によって取得した 30000 フレームの画像を平均化したもの。右：左の図に対して、解析に用いることができないスポット像の領域をマスクしたもの。この画像をもとにスポット像の位置のリファレンスを作成し、各フレームのスポット像の検出に用いた。

この平均画像から解析に用いることのできるスポット像を選択するために、ドーナツ型の外側と内側の縁、検出器によって部分的にスポット像がけられてしまっている箇所、検出器の感度が十分でない領域

第5章

新手法の実践および結果

41(y > 85 の部分) をマスクしたのが図 5.1 右である。この画像から明るいピクセルを探し出し、それぞれの位置で明るさ重心を計算することでスポット像の位置のリファレンスを作成した。さらに、スポット像のガウシアンフィットから求めたスポット像のエアリー直径が 4.45 ピクセルだったことから直径 5 ピク
セルの範囲をスポット像と定義し、全てのスポット像のカウント値を差し引いた画像の中央値をスカイ
バックグラウンドとして見積もった。
次に、スポット像の位置のリファレンスをもとに、30000 フレーム全ての画像に対してスポット像の検出を行なった。ここでもスポット像は明るさ重心位置を中心とする直径 5 ピクセルの円形領域と定義し、各スポット像の重心位置と総カウントの測定を行なった。
得られたスポット像のカウントの変動の原因は大気ゆらぎによるものだけでなく、薄雲による大気の透過率の変化なども考えられる。そこで、大気ゆらぎに起因するものかどうかを調べるために、スポット像のカウントのヒストグラムを調べたのが、図 5.2,5.3 である。大気ゆらぎによって星が明るさ変動をする
場合、明るさの分布が対数正規分布に従うことが知られている [36]。対数正規分布とは、ある確率変数の対数をとった値が正規分布に従うときにその確率変数が従う分布のことであり以下の式で表される。

)(
A(ln x − µ)2
f (x) = √
exp −
2σ 2
2πσx

(5.1)

ここで、A は定数倍、µ, σ は分布の形状を決めるパラメータであるが正規分布の場合とは異なり、分布の平均値、標準偏差値を表しているわけではない。図 5.2 のヒストグラムはある 1 つのスポット像のカウント値の分布を、図 5.3 のヒストグラムは全てのスポット像のカウント値の総和の分布を示している。それぞれの図における赤い線が対数正規分布でフィットした線であり、フィッティングパラメータは図の右上に記してある。いずれの場合も適当なパラメータによってフィットできていることがわかる。

Spot count histgram of one subapertureA = 5.7e + 7= 1.0e + 1= 4.7e 1

1750
1500
Number

1250
1000
750
500
250
0

0

25000 50000 75000 100000 125000 150000 175000
Spot count

図 5.2 ある 1 つのスポット像に着目し、そのスポット像が 1 分間に取るカウントの分布。

第5章

新手法の実践および結果

42

Spot count histgram of all subapertures
1000

A = 1.1e + 9= 1.5e + 1= 1.2e 1

Number

800

600

400

200

0

2500000 3000000 3500000 4000000 4500000 5000000 5500000 6000000
Spot count

図 5.3 全てのスポット像のカウントの総和が 1 分間に取る値の分布。

最後に、各スポット像の位置と明るさの統計情報として、時間変動の平均値、分散値、スポット像同士の共分散値を計算した。その際考慮したフォトンノイズと検出器の読み出しノイズについて記す。まず、
いずれもランダムノイズであるため平均値や共分散値を計算するときにはこれらは寄与しない。明るさ変動の分散を計算する際には、以下に示す成分を減算することで補正した。
22σnoise[ADU2 ] = CFM(Itot + Isky ) + σreadout(5.2)

まず、第 1 項がフォトンノイズ、第 2 項が読み出しノイズに相当する。各文字の意味は、• C: 検出器で決まる、電子数からカウント値への変換ファクター• F : EMCCD のノイズファクター• M : EMCCD の増倍率
• Itot : そのスポットの明るさ変動の平均値 [ADU]• Isky : その観測時刻のスカイバックグラウンドの値 [ADU]• σreadout : 検出器の読み出しノイズ [ADU]である。位置変動の分散を計算する際には Tokovinin 2002[30] を参考に、以下に示す成分を減算することで補正した。
2
σnoise[pixel2 ] = 0.5∑(0.70λ/D)2σ2+ readout(xi,j − < x >)2NSI2tot
window

こちらも第 1 項がフォトンノイズ、第 2 項が読み出しノイズに相当していて、各文字の意味は、

(5.3)

第5章

新手法の実践および結果

43• λ: 観測波長• D: サブ開口の口径• N: 光子数• S: ストレル比• σreadout : 検出器の読み出しノイズ [ADU]• Itot : そのスポットの明るさ変動の平均値 [ADU]• xi,j : スポット内の各ピクセルの位置 [pixel]• < x >: スポット中心の位置 [pixel]である。ストレル比については Tokovinin 2002[30] によると、以下のように見積もることができる。

Imax 4S=
Itot π

(

λCCD
D∆x

)2
(5.4)

この式において、• λCCD : CCD の感度が最も良い波長• ∆x: CCD のピクセルサイズ [rad]である。

5.2 MASS
前節で取り出したマイクロレンズアレイのスポット像の明るさ変動の情報を用いて、MASS 手法を実践した過程について述べてゆく。

5.2.1 開口パターンの取り方
MASS 手法において様々な開口形状を選択することは特定の空間周波数成分を取り出すことに相当する。これまで MASS-DIMM 装置として用いられてきた手法では円環形に 4 つの開口を取ることによって 4 種類の空間周波数成分に相当する明るさゆらぎの情報を取り出し、それらの相関を計算していた。今回シャックハルトマンセンサーを使うことで様々な開口パターンを取ることができる。ここでは選択した
開口パターンについて述べる。図 5.4 に示したのが今回定めた開口パターンの一例である。2 つのスポットによって定義される 2 点開口を考えることとした。特定の空間周波数成分を取り出すという目的があるため、ある特定の距離によって決まるような開口パターンが適していると考えられるためである。この定義に従うと、今回のシャックハルトマンデータにおいて取りうる開口形状のパターン数すなわち 50cm の主鏡上で取れる長さスケールのパターンは 105 通りとなる。

第5章

新手法の実践および結果

44

図 5.4 MASS の開口パターンは 2 つのスポットの取り方のパターンで決めた。今回のデータにおいて取りうる開口形状のパターン数すなわち距離 d の取り方は 105 通りとなった。

また、2 つの開口パターンから決まる Diﬀerential 成分について述べる。これまで用いられてきた方法
では開口が 4 種類しかないため開口同士の組み合わせは 6 種類にとどまるが、今回の方法では開口が 105
種類もあるため組み合わせの総数は非常に多くなる。計算時間などの問題から、組み合わせを取るのは 2
つの開口の中心が一致しているもののみに絞ることとした。図 5.5 に一例を示す。この場合、取りうる開口同士の組み合わせは 641 種類となった。したがって、本研究における MASS では 105+641=746 種類の開口パターンの情報を用いて大気ゆらぎの高さ分布を推定した。

図 5.5 Diﬀerential 成分を計算する開口パターンの組は、両者の中心位置が一致した開口同士とした。

第5章

新手法の実践および結果

45

5.2.2 Scintillation Index の計算方法
Scintillation Index(SI) を式 2.30, 式 2.31 を用いて計算してゆく。今回開口は 2 つのスポットの距離だ
けで決まるものと定義しているので、SI をスポット像の統計量 (平均、分散、共分散) で書き換えることができる。以下のようになる。

[]
IASA = V ar< IA >
[]V ar Ii
[ <Ii > ]=V ar Ii +Ij<Ii +Ij >{ V ar(I )(5.5)(A がスポット i からなる場合)(A がスポット i, j からなる場合)

(5.6)

i

=

<Ii >2
V ar(Ii )+V ar(Ij )+2Cov(Ii Ij )(<Ii >+<Ij >)2

(5.7)

]
IB
IA
,
(5.8)
SAB = SA + SB − 2Cov
< IA > < IB >
)(Ii +IjIk
SA + SB − 2Cov
(A がスポット i, j から B がスポット k からなる場合),( <Ii +Ij > <Ik > )=Ii +IjIk +IlSA + SB − 2Cov
(A がスポット i, j から B がスポット k, l からなる場合)<Ii +Ij > , <Ik +Il >[={Cov(I ,Ik )+Cov(Ij ,Ik )SA + SB − 2 (<Ii i>+<Ij >)<Ik >SA + SB − 2

Cov(Ii ,Ik )+Cov(Ij ,Ik )+Cov(Ii ,Il )+Cov(Ij ,Il )(<Ii >+<Ij >)(<Ik >+<Il >)(5.9)(5.10)

ここで、I は天体からのフォトンカウント値、i, j の添字はスポットの番号を表している。このように、各開口パターンの SI はスポット毎の明るさ変動の統計量から計算できる。
また、今回取得したのは 2 ミリ秒の積分時間で高速撮像したデータであるため、2 ミリ秒よりも短い時
間スケールの変動成分は見えなくなるという影響がある。したがって、理想的な積分時間 0 ミリ秒でのSI はより大きな値をとるはずである。Tokovinin et al. 2003 [27] では、積分時間 τ のデータから求めた
Scintillation Index(SI) と積分時間 2τ のデータから求めた SI の線型外挿から積分時間 0ms の SI を推定する手法が行われており、本研究もこれに従った。以下のような式で表される。

S0ms = 2Sτ ms − S2τ ms(5.11)

なお積分時間 2τ のデータは、積分時間 τ で取得したシャックハルトマン像の画像データセットから時間
的に隣同士にある 2 枚の画像を平均化することで実質的に積分時間 2τ の画像データセットを作成した。

5.2.3 明るさゆらぎの相関関数としての Scintillation Index
ここで、各時刻に取得した 1 分間 (=30000 フレーム) のデータから算出したシンチレーションインデッ
クス (SI) の振る舞いについての結果を示す。図 5.6 に示すのが NSI を、開口形パターン（開口を構成す
る 2 つのスポットの距離）の関数として表したものである。色の違いは観測時刻の違いになっている。

第5章

新手法の実践および結果

46

22:39
22:46
22:51
22:56
23:11
23:27
23:33
23:39
23:50

normal scintillation index

0.6
0.5
0.4
0.3
0.2
0.1

0

5 10 15 20 25 30 35
spot separation on primary mirror [cm]40

図 5.6 開口を構成するスポット間距離と星像の明るさ変動との関係。時刻によって大気乱流強度を反映して曲線が上下している一方、スポット間距離 15cm 程度よりも大きな距離で曲線が平坦になるという特徴は変わらない。

離れた 2 つのスポットが相関を持って明るさ変動を起こしているとすると、その総和である開口の明るさも大きな時間変動を示す。一方、2 スポットの明るさ変動が無相関であれば開口の明るさ変動も小さくなるはずである。つまり、NSI は星像の明るさゆらぎの相関関数に対応している。したがって図 5.6 は大気ゆらぎが明るさ変動に与える空間スケールを表していると見ることができる。観測時刻ごとに大気乱流強度が異なるため NSI の値も異なる。一方カーブの形状の違いは相似形的な変化であり、主鏡上でおよそ15cm のところから長い側ではカーブがフラットになるという特徴は変わらない。このことは大気ゆらぎによる星像の明るさ変動が相関を持つ空間スケールがおよそ 15cm 以下であることを示しており、図 3.1で示した、高さ 20km までの大気ゆらぎが寄与した場合の空間スケールと矛盾のない結果になっている。

5.2.4 Weighting Function
ここでは SI と大気ゆらぎの高さ分布とを結びつける係数行列である、Weighting Function(WF) の計算について述べる。WF の定義式を再掲すると、∫ ∫Wi,j =

1.53f

−11/3

{

sin(πλhj f 2 )λ}2
× Ai (fx , fy )dfx dfy

AA (fx , fy ) = |F[PA (x, y)]|2
AAB (fx , fy ) = |F[PA (x, y) − PB (x, y)]|(5.12)(5.13)
2(5.14)

第5章

新手法の実践および結果47であった。ただし、P (x, y) は各開口形状の開口関数 (開口部を 1、遮蔽部を 0 とした関数) であり F[] はフーリエ変換を表す。開口パターン i は 5.2.1 で選んだ 746 通りである。高さ方向 j は今回は 1-50km まで 2km ごとに 25 通り想定することにした。典型的な大気ゆらぎの存在する高さスケール 20km を大幅に超え、目標とする大気ゆらぎ分解能 ∆h ∼1km を下回る想定をしている理由を説明する。前節でも触れたように観測天体の高度角の変化によって見かけの大気ゆらぎの高さが変わる。例えば観測天体の高度角が θ の場合は、高度 20km にある大気ゆらぎの見かけの高さは 20/ sin θkm となる。今回観測したデネブの場合は時刻によっては θ ∼ 30[deg] 程度になるため、想定する大気ゆらぎの高さは 20km の 2 倍よりも大きい 50km とし高さ分解能は要求の 2 倍の 2km に設定した。
したがって WF 行列は 746 × 25 の行列となる。この行列の各要素を数値計算により求めた。今回用
いた波面センサーのサブアパーチャサイズは直径 25[mm] の円形であるため、1pixel を 2.5[mm] とする

210 × 210 のサンプリング で各開口関数 P (x, y) を表現した。これをフーリエ変換することによって得られる空間周波数領域のサンプリング は 1pixel あたり 0.39[m−1 ] であり式 5.12 の積分を計算するのに十分な細かさであることを確認した。
計算して得られた WF の一例を示す。図 5.7 は Normal Weighting Function をプロットしたもので、横軸に大気ゆらぎの高さ、縦軸に WF の値を示している。開口パターンの違い (すなわち開口をなす 2 つのスポット間距離) を色の違いによって示している。まず初めにわかるのは、どの線も大気ゆらぎの高度に対して単調増加になっている点である。これは、高いゆらぎ層によって作られる明るさゆらぎの方がより大きなパワーを持っていることに対応している。また、2 つのスポット間距離が短いほど WF の値が大きくなっていることは、同じ大気ゆらぎのもとでは 2 つのスポット間距離が短い開口の方が明るさ変動が大きくなることを意味している。このことは、図 2.5 において高空間周波数側に行くほど様々な高度の大気ゆらぎからの寄与があり、トータルのゆらぎのパワーが大きくなることに対応している。図 5.6 で見たように、実際に得られた SI にもこの性質が現れている。また 2 つのスポット間距離が長いほど、ゆらぎ層の高度に対して下に凸型のグラフ形状となっていることは、2 つのスポット間距離が長い (すなわちゆらぎの低空間周波数成分を観測する) 開口ほど高い高度の大気ゆらぎからの寄与が大きいことと対応している。

新手法の実践および結果

log(Weight value [m 1/3])

第5章

48

1012

1011

0

2spots' separation : 0.0 [cm]2spots' separation : 2.5 [cm]2spots' separation : 5.0 [cm]2spots' separation : 7.5 [cm]2spots' separation : 10.0 [cm]
10
20
30
40
50
Height of turbulent layer [km]
図 5.7 Normal Weighting Function

図 5.8 は Diﬀerential Weighting Function をプロットしたものである。Normal Weighting Func-tion(NWF) と 異 な り 、あ る 高 度 に weight の ピ ー ク を 持 つ よ う な 形 状 に な っ て い る 。Kornilov etal.2003[13] によると、このピークの位置は 2 つの開口パターンのうち空間スケールが小さい方の直√√径が ∼ λh となる高さにあたる。 λh は、高さ h の大気ゆらぎが作る明るさゆらぎの代表的な相関長である。NWF との概形の大きな違いから NWF とは相補的な関係になっており、連立方程式の独立性を上げるはたらきをしていると言える。

新手法の実践および結果

log(Weight value [m 1/3])

第5章

49

1011

0

10
20
30
40
Height of turbulent layer [km]50

図 5.8 Diﬀerential Weighting Function

5.2.5 行列の逆解き手法
観測から求めた Scintillation Index から大気ゆらぎの高さ分布を推定する解析的計算は以下の様な連
立 1 次方程式となる。

⃗ = W J⃗S

(5.15)⃗ は各開口パターンに対する Scintillation Index を並べたベクトル、W は Weighting Function
ただし S
行列、J⃗ は各高さでの大気ゆらぎ強度を並べたベクトルである。今、746 通りの開口パターンで得られた情報をもとに、25 通りの高さ方向の分布を求めようとしているため過決定方程式となっている。通常このような行列を解く手法としては最小二乗解を求める手法が用いられる。以下の式で与えられるカイ二乗
値を最小化する手法である。

⃗ =χ2 (J)∑ (Si − W Ji )2i

σi2(5.16)

ここで、Si , σi は測定から求められる SI の値とそのノイズを表している。全体を J⃗ の関数と見たとき、これを最小化する J⃗ が最小二乗解である。最小二乗解の解析解は Press et al.1992[37] によって既に求められており、⃗J⃗leastsquare = (W T W )−1 W S(5.17)

によって表される。しかし今回の問題では WF の形状にも表れている通り、連立する式は相関が強い（すなわち独立性が低い）。このような問題で式 5.17 にしたがって解くと負の乱流強度を含んだ解を得ること

第5章

新手法の実践および結果

50がほとんどであり、物理的に意味のある解を得ることができない。そこで今回は反復法によって式 5.16を最小化する解を探索した。これは Tokovinin et al.2003[27] においても行われている方法である。具体的には、BFGS 法 (Broyden-Fletcher-Goldfarb-Shanno algorithm) を採用した。BFGS 法は非線形最適化問題に対するよく知られた解法である。与えられた評価関数の 2 回微分 (ヘッセ行列) までの情報を計算して評価関数の停留点を探してゆくニュートン法の仲間である。ただしヘッセ行列の計算には計算量の問題から近似的解法を用いており、この点から準ニュートン法と呼ばれる方法に分類されている。BFGS
法が準ニュートン法の中で最も普及している方法であるためこの手法を選んだ。また、反復法を用いることによって解の探索範囲を制限し、正値の解を得ることができるのもメリットである。他の観測サイトにおいて測定された大気ゆらぎの高さ分布 [38][39] から典型的な大気ゆらぎ強度の値の範囲を考慮し、全てのゆらぎ層の高さ j に対して

−17 < log(Jj [m1/3 ]) < −11(5.18)という制約を加えた。一方で反復法を用いることで生じるデメリットとしてはローカルミニマムの問題がある。与えた評価関数が J⃗ に対して複雑な形状をしている場合、局所的な停留点がいくつもあることによって反復の初期値ごとに異なる解を求めてしまう問題である。そこで、ランダムな初期値から最小化の過程を行うことを何度も繰り返し (今回は 1000 回)、得られた多数の解に対してその評価値が良いものか
ら 10% に入るものだけを抽出し、その平均とばらつきを評価することにした。

5.2.6 大気ゆらぎの高さ分布得られた大気ゆらぎの高さ分布の推定結果を図 5.9 に示す。9 つのグラフは 9 つの観測時刻に対応して
2
dh[m1/3 ] を取っている。時刻とおり、各グラフでは縦軸に地表からの高度を、横軸に大気ゆらぎ強度 CN

ともに高さ分布の形状が変わってゆくのに対して、どの観測時刻の結果に対しても最も地表に近い層で
22CNdh ∼ 10−12 [m1/3 ] と最もゆらぎが強いこと、高度 ∼ 18km 以上において CN
dh < 10−17 [m1/3 ] とゆらぎが無視できるほど弱いことなどの特徴は変わらない。これらの結果から、２章で述べた大気ゆらぎの高さ分布の典型的構造は、東北大学のサイトにおいても成立していることがわかった。また同時に、今回の推定手法によってデータセットによらない高さ分布の推定ができていることが支持される結果となった。一方で観測時刻によっては、推定した高さ分布につく強度のエラーが大きい部分がある。これは今回の測定データを説明する大気ゆらぎの高さ分布の可能性がいくつもあることを意味しており、手法が持つ本質的な高さ方向の分解能が、想定した 2km よりも粗いことを示している。

新手法の実践および結果

Observation time : 22:39

30
25

20

20

20

15

Altitude[km]

25

15
10

10

5

5

5

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

Observation time : 22:56

30

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

Observation time : 23:11

30
25

20

20

20
Altitude[km]25Altitude[km]

25

15

15
10

10

5

5

5

30

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

Observation time : 23:33

30

Observation time : 23:39

30

25

25

20

20

20
Altitude[km]

25

15

15
10

10

5

5

5

102

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

102

Observation time : 23:50

15

10

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

Observation time : 23:27

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

102

15

10

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

Observation time : 22:51

15

10

30

Altitude[km]Observation time : 22:46

25

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

Altitude[km]30Altitude[km]Altitude[km]3051Altitude[km]第5章

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

図 5.9 今回実践した新しい手法によって得られた大気ゆらぎの高さ分布。縦軸が地表面からの高度、横軸が大気屈折率の構造定数にゆらぎ層の厚みをかけた値になっている。地表付近に最も強い乱流層があり、高度 20km 以上ではエラーの範囲で十分大気ゆらぎが弱いことがわかる。

第5章

新手法の実践および結果

52

また、今回の推定結果のフィッティングの適合度についての結果を図 5.10 に示す。いずれの図も観測
時刻 22 時 39 分の推定結果に対するものとなっており、左の図は推定の評価関数として用いたカイ二乗
値の値の分布を示したものである。モンテカルロ的に行なった 1000 回の試行のヒストグラムになっている。今回は上位 10% を用いているため推定の reduced カイ二乗値は ∼ 5 であることがわかる。一方で右の図は観測によって求められた Scintillation Index(SI) と、得られた高さ分布から求められる SI との相対残差を開口パターンのヒストグラムとして示したものである。相対残差は以下の式に従って計算した。

Relative Residual =SIobserved − SImodelled
SImodelled(5.19)

相対残差の値は全ての開口パターンに対して 7% 以内に収まっていることがわかる。これら 2 つの結果か
らも、観測データを十分に説明できる推定結果になっていることがわかる。また、他の時刻の結果についてもカイ二乗値は 4 < χ2 < 6 の範囲に、相対残差は 10% 以内であった。

50
40
NumberNumber

102

101

30
20
10

100

101

102
Reduced chi square value

0 6

4

2
0
2
4
SI relative residual [%]6

図 5.10 推定した高さ分布のフィッティングの適合度を調べた結果。左：1000 回の最適化試行に対す
るカイ二乗値の分布。最終的な高さ分布推定は全体の上位 10% すなわちカイ二乗値 ∼ 5 の結果を用いている。右：推定された高さ分布と大気モデルから計算される SI と観測された SI の相対誤差を開口パターンのヒストグラムとして表したもの。全ての開口パターンに対して観測された SI を 10% 以
下の残差の範囲で説明できる推定結果となっている。

次に、各高度における大気ゆらぎ強度の時間変化を調べた結果を図 5.11 に示す。図 5.9 に示した高さ
分布を 0-5km,5-10km,10-15km,15-20km,20-25km,25km 以上の 6 つに分け、それぞれの高度帯におけるゆらぎ強度の時間変化を示している。5km ごとに分けた各高度帯における大気ゆらぎ強度のオーダーは

5-10km 帯では時間とともに徐々に強くなってゆく一方、それ以外の高度帯ではほとんど変化していない。
20km 以上の高度帯における強度はそれ以下の強度に比べて 1 桁-2 桁小さく、十分無視できることがわ
かる。

第5章

新手法の実践および結果

53

Cn^2dh/1.0e-14 [m^1/3]

102

0-5km
5-10km
10-15km15-20km20-25km
25km-

101
100
10 1
10 2
10 3
22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

図 5.11 各高さの大気ゆらぎ強度の時間変動。図 5.9 に対して 5km ごとに積分した大気ゆらぎ強度
を、時間の関数として表している。

また、得られたゆらぎ強度の高さ分布を積分することによって、フリード長やシーイング、アイソプラ
ナティック角を算出することができる。図 5.12 にそれらの値の時間変化を示す。観測時刻 23 時 11 分を境に大気の状態が悪くなっていることがわかる。この傾向は、図 5.6 において観測時刻 23 時 11 分以降のSI の値が急激に大きくなっていることと対応していると考えられる。また、図 5.11 と組み合わせて考えると、5-10km の大気ゆらぎ強度の増加に伴って、パワーの強い 0-5km の大気ゆらぎ強度が強くなったことが原因であると推察できる。

6.5

Fried parameter

seeing value

6.0
5.0

[arcsec]

[cm]

5.5
4.5
4.0
3.5
22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

3.0
2.8
2.6
2.4
2.2
2.0
1.8
1.6
22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

図 5.12 大気ゆらぎ強度の積分値に対応する物理量の時間変化。左がフリード長、右がシーイングでありいずれも 500nm における天頂方向の値になっている。

第5章

新手法の実践および結果

54

5.3 DIMM
ここからは、5.1 節において取り出したマイクロレンズアレイのスポット像位置の変動情報を用いてDIMM 手法を行なった結果を示す。

5.3.1 開口の取り方
DIMM 手法では 2 つの開口で得られる星像からその相対的な位置の時間変動を測定する。これまで
MASS-DIMM 装置として用いられてきた手法ではある決まった距離だけ離れた 2 つの開口を用いていた。波面センサーを用いた今回の解析では、様々な開口間距離に対応する 2 つのスポット像を自由に選択することができる。取りうる開口間距離の分だけ測定点が増えるため、従来の手法に比べて測定のランダ
ム誤差による影響が少なくできると考えられる。

5.3.2 Diﬀerential Image Motion の計算
次に、5.1 で取得した各スポット像の位置変動の統計量（平均・分散・共分散）から 2 つの星像の相対的な位置変動 (Diﬀerential Image Motion) を計算する流れを説明する。式 2.42,2.43 から 2 つの方向における Diﬀerential Image Motion は、

σx2 ≡< (xi − xj )2 >=< x2i > + < x2j > −2 < xi xj >(5.20)= V ar[xi ] + V ar[xj ] − 2Cov[xi , xj ]σy2 ≡< (yi − yj )2 >=< yi2 > + < yj2 > −2 < yi yj >(5.21)= V ar[yi ] + V ar[yj ] − 2Cov[yi , yj ]と表される。但し i, j はスポットのインデックスであり、今 x 方向を 2 つの開口が離れる方向（longitu-

dinal 方向）としている。5.1 節で取り出したスポット像ごとの位置変動の統計情報をこれらの式に代入することによって計算できる。

5.3.3 構造関数の理論曲線とのフィット測定した Diﬀerential Image Motion から、式 2.62,2.63 を用いてフリード長 r0 を求める過程について述べる。波面センサーデータを用いることで、測定量である σx2 , σy2 は様々な開口間距離に対して取得される。したがって式 2.62,2.63 と測定量とを最小二乗フィットすることによって r0 の値を求めた。
図 5.13 に、フィッティングの様子を示す。横軸に 2 つの開口間の距離を、縦軸に σx , σy の値を表している。青系の色が開口の離れている方向 (longitudinal direction) すなわち σx 、赤系の色がそれと直交す
る方向 (transverse direction) すなわち σy を示している。実際には 2.5[cm] 以上の開口間距離に対して網羅的に Diﬀerential Image Motion が計算されるが、フィッティングは開口間距離が 25cm 以上離れている場合のデータ点に対してのみ行なっている。これは DIMM の原理上、空間スケールの大きい (すな
わち空間周波数スケールの小さい) 位置ゆらぎのパワーには全ての高度の大気ゆらぎの影響が同等に含まれており、大気ゆらぎ全体のパワーを反映しているためである。

第5章

新手法の実践および結果

55

square root of structure function [rad]0.000005
0.000004
0.000003
0.000002
0.000001
0

10

20
30
baseline length [cm]longitudinal,2ms
transverse,2ms
longitudinal,4ms
transverse,4ms
40

図 5.13 Diﬀerential Image Motion の測定値と式 2.62,2.63 とのフィッティングの様子。エラーの範囲でよくフィットできていることがわかる。

5.3.4 積分時間の補正
DIMM の積分時間補正について述べる。MASS についても述べたとおり、今回のデータは 2 ミリ秒の積分時間での高速撮像であり、その積分時間よりも短い時間スケールの変動成分が見えていないことが考えられる。理想的には無限小の積分時間であるべきであり、その補正方法についてはいくつかの手法が挙げられている [29][30]。今回は、シーイング値が積分時間 τ に対して exp(−aτ ) の関係で変化していくものと仮定して、以下の式によって積分時間の補正を行なった。

θseeing,0ms =(θseeing,2ms )2
θseeing,4ms(5.22)

この式はフリード長に対して書き換えた場合も同様の形式になる。

r0,seeing,0ms =(r0,seeing,2ms )2
r0,seeing,4ms(5.23)

5.3.5 フリード長とシーイングの推定結果
DIMM によって求めた、フリード長やシーイングといった大気ゆらぎの高さ分布の積分値に対応する物理量の推定結果について述べる。図 5.14 は、各時刻での測定データに対して DIMM の手法によって求めた観測波長 500nm におけるフリード長を表している。色の違いは 2 つの開口が離れている方向(longitudinal direciton) とそれに直交する方向 (transverse direction) の星像ゆらぎから算出したフリー第5章

新手法の実践および結果

56

ド長であることを示している。両者の値自体には大きいところで 1.5cm ほどの差異が見られるものの、Fried parameter

6.75
6.50
6.25
6.00
5.75
5.50
5.25
5.00

seeing value
DIMM-longitudinal
DIMM-transverse

2.0

DIMM-longitudinal
DIMM-transverse

1.9[arcsec][cm]

振る舞いの傾向はよく合っている。

1.8
1.7
1.6

22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

1.5

Observation time

22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

図 5.14 DIMM から求めた大気ゆらぎ強度の積分値に対応する物理量の時間変化。左がフリード長、右がシーイングでありいずれも 500nm における値になっている。

5.4 MASS と DIMM の結果の対応関係
MASS と DIMM はそれぞれ、星像の明るさゆらぎと位置ゆらぎの情報を用いて大気ゆらぎの推定を行う手法である。星像の明るさは光の振幅、位置は光の位相にそれぞれ相当するため、測定する物理量が独立であることから MASS と DIMM は手法として独立だと言える。そこで、MASS によって求めたフリード長とシーイングを DIMM によって求めたそれらと比較したのが図 5.15 である。時間変動の全体的な振る舞いは両者で同様になっている。大気の状態が観測時間の前半では徐々に良くなっていき 23 時

11 分を過ぎたあたりで急激に悪くなり、その後また次第に良くなっていったという流れは有意に結論づけられる。一方でフリード長やシーイングの絶対値には乖離が見られる。またそれぞれの手法の原理を考えると、MASS は高度 ∼0km の地表層には感度を持たない一方で DIMM は全ての高度に対して感度を持つはずである。したがって一般的に最もゆらぎの強い地表層を含む DIMM の方が悪い大気状態を示すことが予想される。今回の結果は逆の関係になっている。国内の他サイトのシーイング値を考えると、MASS による結果は比較的妥当な値となっているのに対して、DIMM は大気状態を良い方向に過大評価しているように見える。

第5章

新手法の実践および結果

57

Fried parameter
6.5
6.0

2.6

DIMM-longitudinal
DIMM-transverse
MASS

2.4

5.5

[arcsec][cm]

seeing value
DIMM-longitudinal
DIMM-transverse
MASS

5.0

2.2
2.0

4.5

1.8

4.0

1.6
22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

22:39 22:46 22:51 22:56 23:11 23:27 23:33 23:39 23:50

Observation time

図 5.15 フリード長およびシーイングの時間変化を MASS による結果と DIMM による結果とで比較したもの。

58

第6章

議論
6.1 MASS の結果に関して
6.1.1 大気ゆらぎの高さ分布の妥当性
今回の推定手法で得られた大気ゆらぎの高さ分布の妥当性を議論する。まずプロファイルの概形については 5 章でも述べたように、地表層のゆらぎ強度が最も強いこと、高度 10-15km にかけて強度の局所的なピークがあること、高度 20km よりも高い位置にはゆらぎはほとんどないことなどの特徴が見られた。
これらの特徴は２章において述べたような典型的な大気構造と矛盾しないものとなっている。
次に他のサイトとも結果を比較してみる。本来異なる観測サイト間ではその緯度や地理的な構造が異なるため、大気ゆらぎの高さ分布は直接比較できるものではない。しかし今回測定した仙台市上空の大気についてはこれまでに取得された情報が他に存在しないため、他の観測サイトの結果とおおよその特徴のみを比べることとした。実際、大気の大局的な構造が大きく異なるということは考えにくいため、大気構造を特徴付けるパラメータ同士は比較できないもののプロファイルの概形や各高度にどれくらいの割合の大気ゆらぎ強度が存在しているかを比較することができると考えられる。図 6.1 に示すのは Ono

et al.2016[38] において測定されたマウナケア（標高 4000m）における大気ゆらぎの高さ分布である。
SLODAR と呼ばれる手法 (本論文では波面センサーの相関を用いる手法として前述) によって求められたものであり、黒の実線で描かれているものが Ono et al.2016 における結果となっている。また、6.2 に示すのは Farley et al.2018[40] において測定されたパラナル (標高 2600m) における大気ゆらぎ分布であり、Stereo-SCIDAR と呼ばれる手法 (本論文では波面センサーの相関を用いる手法として前述) によって求められたものである。いずれの結果も観測サイトの標高から測った大気までの距離を大気の高度として算出しているため、マウナケアであればおよそ 4km、パラナルであればおよそ 2.6km のオフセットがあることに注意が必要である。一方で地表付近において地面との摩擦により大きな乱流構造が生まれるという特徴は標高によらず変わらない。したがっておよそ高度 1km 以下の値は地面との摩擦による強い乱流成分を、それより高い高度における値は実際の大気の乱流強度を算出していると見ることができる。これらと今回の結果との比較からも、得られた高さ分布の概形には共通する部分が多く見られる。特に Farley

et al. 2018 とは海抜 3-4km 付近で急激に乱流強度が弱くなるという、細かいプロファイルの形状もよく一致している。

第6章

議論

59

図 6.1 Ono et al.2016[38] によって調べられたマウナケアにおける大気ゆらぎのプロファイル。本研
究で求めた図とは、横軸はリニアスケールとなっていることに注意が必要である。

図 6.2 Farley et al.2018[40] によって調べられたパラナルにおける大気ゆらぎのプロファイル。本研
究で求めた図とは横軸・縦軸が入れ替わっていることに注意が必要である。

6.1.2 大気ゆらぎの強度の推定エラー
本研究で推定した大気ゆらぎの高さ分布は全体の傾向として典型的な大気構造を反映したプロファイルになっており、かつ他の時刻における観測データからも同様の傾向が見えていた。一方で観測時刻によっては、ゆらぎ強度につくエラーが大きくなっている箇所がある。このエラーは今回の手法の場合、反復法の初期値の違いによって異なる解に行き着いていることが原因である。ここでは反復法を用いた行列の逆

第6章

議論

60

解き手法の妥当性について考察を行う。
まず、反復法の初期値をランダムに選ぶモンテカルロ計算の試行回数について検討する。この回数が小さい場合、多くの試行においてグローバルミニマムにたどり着かないということが起こりうる。これによって様々なローカルミニマムの解がカイ二乗値の上位 10% に入り込み、エラーを大きくしている可能性がある。一方で試行回数が大きいほど計算に膨大な時間がかかり、事前情報更新のタイムスケールに間に合わないという状況になってしまう。そこで今回の 1000 回という試行回数が妥当なものであるかを判断するために、様々な試行回数における推定を行なった。図 6.3 に示すのがその結果である。図 6.3 のそ
れぞれの線は、試行の回数 N を N = 500, N = 2000 と変えて行った結果を今回の N = 1000 のものと2
重ねてプロットしたものである。比較的乱流強度の高い CN
dh > 10−14 [m1/3 ] におけるプロファイルの概形を比較した時、N = 500 と N = 1000 の結果の間にはずれが見られるが N = 1000 と N = 2000 の結果はよく一致している。このことから、今回の N = 1000 の試行回数によって十分な解の空間の探索ができていると言える。

N=500

30

N=1000

30

25

25

20

20

20

15

Altitude[km]25Altitude[km]Altitude[km]

30

15

15

10

10

10

5

5

5

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

N=20000 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

図 6.3 ランダムに初期値を選ぶ回数を変えた場合の大気ゆらぎの高さ分布の比較。1000 回の試行で十分な解の探索ができていることがわかる。

また図 6.4 は、試行回数 N を変えて行った結果のカイ二乗値と相対誤差の値をヒストグラムとして重ねたものである。

第6章

議論

61

N=2000N=1000
N=500

102
101

N=2000N=1000N=500

80
NumberNumber

103

60
40
20

100

101
Reduced chi square value

0 10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5
SI relative residual [%]図 6.4 想定する大気ゆらぎ層の数を変えた場合の SI のカイ 2 乗値と相対残差の比較。

次に、高さ方向の分解能について検討する。今回の 2.5[cm] の開口を用いた MASS では、感度を持つ
最低の大気ゆらぎ高度は 1km 以下になるものの高度の幅については 1km 以下を達成できていない可能性がある。このような場合、最適化が過決定な状態なため様々な解の候補が生じており、エラーを大きくしている可能性が考えられる。そこで想定する大気層の数を変えて同様の最適化を行なった。図 6.5 は想定する大気層の数を 17 層、50 層と変えて行なった結果を今回のものと重ねてプロットしたものである。大気層の数を変えて行なった場合も層の間隔は等間隔とした。天体の高度角の影響を考えると、17 層、25
層、50 層における高さ方向の分割はそれぞれ 2.2km、1.4km、0.72km に対応する。図 6.5 からわかるのは、想定する大気ゆらぎ層の数によって推定エラーが大きく変わることである。50 層の結果ではほとんどの高度のゆらぎ強度に対して上限値しか与えられないのに対して、17 層の結果ではエラーはかなり小さくなっている。

17 layers

30

25 layers

30

25

25

20

20

20

15

Altitude[km]25Altitude[km]Altitude[km]

30

15

15

10

10

10

5

5

5

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

50 layers

0 10 3 10 2 10 1 100 101
Cn2dh/10 14[m1/3]

102

図 6.5 大気層の数を変えた場合の大気ゆらぎの高さ分布の比較。推定エラーは高度方向の分解能によ
る縮退によって生じていることが分かる。

第6章

議論

62

図 6.6 は、大気層の数を変えて行った結果のカイ二乗値と相対誤差の値をヒストグラムとして重ねたものである。全ての場合で十分にカイ二乗値の小さい解が得られていること、測定の SI を 10% の誤差で説明できる解になっていることがわかる。図 6.5 と組み合わせて考えることによって、大気層の数が多すぎると様々な解によって測定が説明できてしまうことがわかる。今回の手法においては高さ方向の分解能と推定エラーとがトレードオフの関係にあることがわかった。ここで重要なのは、トモグラフィ計算の事前情報として許容されるエラーの大きさとの比較であるが、これについてはまだ評価されておらず今後数値
シミュレーションなどによる評価が課題となる。

103

17 layers
25 layers
50 layers

NumberNumber102
101

17 layers
25 layers
50 layers

80
60
40
20

100

102
101
Reduced chi square value
図 6.6

0

10.0 7.5 5.0 2.5 0.0 2.5 5.0 7.5 10.0
SI relative residual [%]大気層の数を変えた場合の SI のカイ 2 乗値と相対残差の比較。

最後に、大気ゆらぎ強度の推定結果につくエラーの原因としてはモデルとのずれも考えられる。MASS-

DIMM では星像の位置や明るさのパワースペクトラムを計算する際に 3 つの仮定をしている。それは大気乱流がコルモゴロフモデルに従うこと、大気中の光の伝搬がフレネル伝搬で記述できること、弱乱流近似（乱流が比較的弱く各店の位相の平均値との差が 1rad よりも十分小さくなること）であった。これらのうちのいずれかが今回の東北大学のサイトでは満たされていない可能性がある。実際にモデルとの適合度を示すカイ 2 乗値の値は、全ての観測時刻の結果に対して 4 < χ2 < 6 であり 1 以下になるような測定
データセットはなかった。推定において十分な解空間の探索ができていることを考慮すると、このカイ 2乗値はモデルとの本質的なずれを表している可能性が考えられる。

6.2 DIMM の結果に関して
今回の DIMM 手法では、星像の重心位置の情報を用いることでフリード長やシーイングを測定したがそれらの値は MASS 手法による結果とはコンシステントにはならなかった。日本国内におけるシーイン
グ測定の結果を踏まえると、今回の DIMM 手法の結果が想定されるよりも非常に良い大気状態を示していると考えられる。実際に、シャックハルトマン波面センサーで得られたスポット像をガウシアンフィッ
ティングすることによって見積もられるシーイングサイズは 4 秒角程度となっており、この値は MASSの結果とも整合する。今回の像の重心検出手法では、大気ゆらぎによる像の位置変動を正確に測定できていない可能性が考えられるため、重心検出手法の再検討を含め改善の余地が残る結果となった。

第6章

議論

63

6.3 トモグラフィの事前情報取得手法としての機能性
今回新たに検証した手法は、Fusco & Costile2010[11] によって得られたシミュレーションの結果をもとに、トモグラフィ手法を用いた補償光学に必要となる大気ゆらぎの事前情報取得を目指したものであった。重要となる指標は数十分程度のタイムスケールで情報更新ができる即時性、観測や補償と同じ方向の大気ゆらぎの推定ができること、大気ゆらぎの存在する高度 20km までのゆらぎ強度をおよそ 1km の分解能で推定できることであった。これらを踏まえて今回の手法が実際に運用できるものであるかどうか議論してゆく。まず、今回の手法の計算アルゴリズムでは反復法を用いた高さ分布推定位の段階に多くの時間を要する。試行回数 N < 1000 程度であれば計算全体として数十分の範囲に収まる計算量である。

6.1.2 での結果からも N = 1000 の試行回数で十分な解の評価ができていたため即時性に関する部分は問題にならないと考える。また既存の MASS-DIMM とは異なり、波面センサーを用いた手法実践を実現したことによって観測や補償と完全に同一な方向の大気ゆらぎの推定ができることがわかった。さらに今回の手法は単一の星の観測で実践できるため、感度をもつ高度に上限がないことも大きなメリットである。一方で MASS が感度を持つ高度の下限は今回の実験では 1km 程度だが、実際の運用において観測天体（レーザーガイド星）が暗いことを考えると高度の下限は 2.5km 程度となる。DIMM によって高度方向の強度積分値が得られたとしても 0-2.5km のゆらぎをひとまとめにした推定しかできない。すなわち

SLODAR や SCIDAR とは反対に、低高度での分解能が悪くなるということが分かった。また、高さ方向の分解能は推定のエラーとトレードオフの関係にあり、定量的に評価するにはどれほどの推定エラーがトモグラフィ計算に影響するのかについての検討が必要であることが分かった。

6.4 今後の展望
今後はまず、今回の DIMM の結果が MASS の結果と整合しなかった点について手法の再検討を行う。
スポット像のサイズから予想される実際のシーイングサイズよりも小さく見積もられる結果となっていることから、位置変動を捉えるための重心検出方法の改善が必要と考えられる。
また前節でも述べたように、トモグラフィ計算の事前情報としてはどれほどの精度が求められるのかについてのシミュレーション等が重要である。その結果によって今回の推定の高さ方向の分解能が決定し、

δh ∼ 1km を満たせているのかどうかの議論が可能となる。
実際の運用を考えたときに低高度に感度がないという結論であったが、これについては SLODAR やSCIDAR と組み合わせることによって解決できるかもしれない。SLODAR,SCIDAR も波面センサーを用いる手法であるため今回検証した方法と共存が可能と考えられる。このような組み合わせについてより定量的な検討も進めていきたい。

64

第7章

結論
すばる望遠鏡をはじめとする 8m 級の望遠鏡において、可視光領域・広視野での補償光学を実現するには複数のレーザーガイド星を用いたトモグラフィ補償光学が必要である。また近い将来、トモグラフィの技術は 30m 級の望遠鏡の補償光学における基盤技術となることが予想される。一方で、トモグラフィ計
算には大気乱流の構造に関する事前情報が大きな影響を与える。本研究では、大気乱流強度の高さ分布を推定する手法についての検討を行なった。
手法を評価する上で重要と考えたのは、数十分程度の時間分解能で推定できること、観測と同じ方向における推定ができること、高さ方向の分解能および範囲である。そこで望遠鏡のサイト調査にこれまで広く使われてきた MASS-DIMM の手法を補償光学に用いられる波面センサーで行うことを考えた。波面センサーの情報を用いることで観測と同一方向のリアルタイムな推定が可能となる。また波面センサーのサ
ブ開口の選び方によって様々な開口パターンを作ることができることは、従来の MASS-DIMM に無い高
さ方向の分解能が達成できると考えた。
実際に波面センサー光学系を設計し 50cm 望遠鏡に取り付けて星の明るさ・位置変動データを取得・解析を行なった。結果として、典型的な大気構造を反映した高さ分布を得ることができ、シーイングの値の時間変化は独立な推定法である MASS と DIMM において同様の傾向が示された。一方で、DIMM の手法から得られる結果が大気状態を良い方向に過大評価する結果となり、地表層の推定ができなかった。また今回の手法の高さ方向の分解能を定量的に評価するためには、トモグラフィ計算としての事前情報推定の精度を今後検討してゆく必要があることもわかった。今後はこれらの課題に取り組むことによって本手法の最終的な評価を進めてゆきたい。さらに SLODAR や SCIDAR と今回の MASS 手法を組み合わせる手法についても検討することで、今後のトモグラフィ補償光学の時代において最適な大気ゆらぎの高さ
分布推定法を提案したい。

65

謝辞
修士課程の 2 年間で多くの方々にお世話になりました。
まず、指導教員の秋山正幸教授に心より感謝申し上げます。秋山さんには学部生の期間を含めて 3 年間、研究や勉学に関する指導や助言をいただきました。会議や出張等で忙しい時でも時間を割いて、丁寧に相談にのってくださいました。また多くの学会や研究会に参加する機会を与えていただいたことは、非常に貴重な経験となりました。来年度からも引き続き博士課程学生としてお世話になります。研究者としてまだまだ未熟者ですが、どうぞこれからもよろしくお願いいたします。
また、国立天文台の大屋真准教授、美濃和陽典准教授、大野良人助教に深く感謝申し上げます。毎週のミーティングでいただいた多くのアドバイスが研究を進める上で非常に役立ちました。大屋さんには補償光学研究会などでも多くの助言をいただきました。美濃和さん、大野さんには、ハワイでの光学作業の際に大変お世話になりました。本当にありがとうございました。
秋山研の皆様をはじめとする、天文学教室の皆様に感謝申し上げます。ゼミ等でいただいた質問によって私自身の理解が深まった場面がたくさんありました。本当にありがとうございました。
最後に、これまでの学生生活を支え、これからの博士課程への進学を快く受け入れてくれた家族に心から感謝いたします。

66

参考文献
[1] Performance of subaru ao188. https://subarutelescope.org/Observing/Instruments/AO/
performance.html.
[2] Yosuke Minowa, Yutaka Hayano, Shin Oya, Makoto Watanabe, Masayuki Hattori, Olivier
Guyon, Sebastian Egner, Yoshihiko Saito, Meguro Ito, Hideki Takami, et al. Performanceof subaru adaptive optics system ao188. In Adaptive Optics Systems II, volume 7736, page77363N. International Society for Optics and Photonics, 2010.
[3] T Sean Ross. Limitations and applicability of the maréchal approximation. Applied optics,
48(10):1812–1818, 2009.
[4] Warren Skidmore, Sebastian Els, Tony Travouillon, Reed Riddle, Matthias Schöck, Edison
Bustos, Juan Seguel, and David Walker. Thirty meter telescope site testing v: seeing andisoplanatic angle. Publications of the Astronomical Society of the Pacific, 121(884):1151, 2009.
[5] R Genzel, NM Förster Schreiber, D Rosario, P Lang, D Lutz, E Wisnioski, E Wuyts, S Wuyts,
K Bandara, R Bender, et al. Evidence for wide-spread active galactic nucleus-driven outflows
in the most massive z 1-2 star-forming galaxies. The Astrophysical Journal, 796(1):7, 2014.
[6] M Akiyama, Y Minowa, N Kobayashi, K Ohta, and I Iwata. Adaptive optics imaging of lyman
break galaxies as progenitors of spheroids in the local universe. Proceedings of the International
Astronomical Union, 3(S245):447–450, 2007.
[7] NM Förster Schreiber, Alvio Renzini, Chiara Mancini, Reinhard Genzel, Nicolas Bouché, Giovanni Cresci, Erin KS Hicks, Simon J Lilly, Yingjie Peng, Andreas Burkert, et al. The sins/zcsinf survey of z 2 galaxy kinematics: Sinfoni adaptive optics–assisted data and kiloparsec-scale
emission-line properties. The Astrophysical Journal Supplement Series, 238(2):21, 2018.
[8] Thayne Currie, Ruslan Belikov, Olivier Guyon, N Jeremy Kasdin, Christian Marois, Mark
Marley, Kerri Cahoy, Michael McElwain, Eduardo Bendek, Marc Kuchner, et al. Using groundbased telescopes to mature key technologies and advance science for future nasa exoplanet direct
imaging missions. arXiv preprint arXiv:1803.05453, 2018.
[9] Seppo Mattila, Petri Väisänen, Duncan Farrah, Andreas Efstathiou, WPS Meikle, Tomas
Dahlén, Claes Fransson, Paulina Lira, Peter Lundqvist, Göran Östlin, et al. Adaptive optics discovery of supernova 2004ip in the nuclear regions of the luminous infrared galaxy iras
18293–3413. The Astrophysical Journal Letters, 659(1):L9, 2007.

参考文献

67

[10] Rainer Schödel, T Ott, R Genzel, R Hofmann, M Lehnert, A Eckart, N Mouawad, T Alexander,
MJ Reid, R Lenzen, et al. A star in a 15.2-year orbit around the supermassive black hole atthe centre of the milky way. Nature, 419(6908):694–696, 2002.
[11] Thierry Fusco and Anne Costille. Impact of cn2 profile structure on wide-field ao performance.
In Adaptive Optics Systems II, volume 7736, page 77360J. International Society for Optics andPhotonics, 2010.
[12] A Costille and T Fusco. Impact of cn2 profile on tomographic reconstruction performance:
application to e-elt wide field ao systems. In Adaptive Optics Systems III, volume 8447, page
844757. International Society for Optics and Photonics, 2012.
[13] Victor Kornilov, Andrei A Tokovinin, Olga Vozyakova, Andrei Zaitsev, Nicolai Shatsky, Serguei F Potanin, and Marc S Sarazin. Mass: a monitor of the vertical turbulence distribution.
In Adaptive Optical System Technologies II, volume 4839, pages 837–845. International Society
for Optics and Photonics, 2003.
[14] François Roddier. V the eﬀects of atmospheric turbulence in optical astronomy. In Progress in
optics, volume 19, pages 281–376. Elsevier, 1981.
[15] Andrei Nikolaevitch Kolmogorov. Local turbulence structure in incompressible fluids at very
high reynolds numbers. In Dokl. Akad. Nauk SSSR, volume 30, 1941.
[16] Alexander Mikhailovich Obukhov. Structure of temperature field in turbulent flow. Technicalreport, AIR FORCE SYSTEMS COMMAND WRIGHT-PATTERSON AFB OH FOREIGNTECHNOLOGY DIVISION, 1970.
[17] John W Hardy. Adaptive optics for astronomical telescopes, volume 16. Oxford University Presson Demand, 1998.
[18] Bun ei Sato, Eiji Kambe, Yoichi Takeda, Hideyuki Izumiura, and Hiroyasu Ando. Developmentof iodine cells for the subaru hds and the okayama hides: Ii. new software for precise radial
velocity measurements. Publications of the Astronomical Society of Japan, 54(6):873–882, 2002.
[19] M Schöck, S Els, R Riddle, W Skidmore, T Travouillon, R Blum, E Bustos, G Chanan, SG Djorgovski, P Gillett, et al. Thirty meter telescope site testing i: overview. Publications of the
Astronomical Society of the Pacific, 121(878):384, 2009.
[20] R Barletti, G Ceppatelli, L Paterno, A Righini, and N Speroni. Astronomical site testingwith balloon borne radiosondes-results about atmospheric turbulence, solar seeing and stellar
scintillation. Astronomy and Astrophysics, 54:649–659, 1977.
[21] CE Coulman, J Vernin, and A Fuchs. Optical seeing―mechanism of formation of thin turbulent
laminae in the atmosphere. Applied optics, 34(24):5461–5474, 1995.
[22] M Azouit and J Vernin. Optical turbulence profiling with balloons relevant to astronomy andatmospheric physics. Publications of the Astronomical Society of the Pacific, 117(831):536, 2005.
[23] John P McHugh, George Y Jumper, and Mark Chun. Balloon thermosonde measurements overmauna kea and comparison with seeing measurements. Publications of the Astronomical Societyof the Pacific, 120(874):1318, 2008.

参考文献

68

[24] Remy Avila, Jean Vernin, and Elena Masciadri. Whole atmospheric-turbulence profiling withgeneralized scidar. Applied Optics, 36(30):7898–7905, 1997.
[25] Richard W Wilson. Slodar: measuring optical turbulence altitude with a shack–hartmann
wavefront sensor. Monthly Notices of the Royal Astronomical Society, 337(1):103–108, 2002.
[26] V Kornilov, A Tokovinin, N Shatsky, O Voziakova, S Potanin, and B Safonov. Combinedmass–dimm instruments for atmospheric turbulence studies. Monthly Notices of the Royal Astronomical Society, 382(3):1268–1278, 2007.
[27] A Tokovinin, V Kornilov, N Shatsky, and O Voziakova. Restoration of turbulence profile fromscintillation indices. Monthly Notices of the Royal Astronomical Society, 343(3):891–899, 2003.
[28] Andrei Tokovinin. Turbulence profiles from the scintillation of stars, planets and moon. In
Workshop on Astronomical Site Evaluation (Eds. I. Cruz-González, J. Echevarrı́a & D. Hiriart),
volume 31, pages 61–70, 2007.
[29] Marc Sarazin and F Roddier. The eso diﬀerential image motion monitor. Astronomy andAstrophysics, 227:294–300, 1990.
[30] Andrei Tokovinin. From diﬀerential image motion to seeing. Publications of the Astronomical
Society of the Pacific, 114(800):1156, 2002.
[31] VI Tatarski. The eﬀects of the turbulent atmosphere on wave propagation. US Department ofCommerce, 1971.
[32] David L Fried.

Statistics of a geometric representation of wavefront distortion.

JoSA,

55(11):1427–1435, 1965.
[33] David L Fried. Diﬀerential angle of arrival: theory, evaluation, and measurement feasibility.
Radio Science, 10(1):71–76, 1975.
[34] Mark Stanford Robbins and Benjamin James Hadwen. The noise performance of electron multiplying charge-coupled devices. IEEE transactions on electron devices, 50(5):1227–1232, 2003.
[35] JR Ducati. Vizier online data catalog: Catalogue of stellar photometry in johnson’s 11-color
system. VizieR Online Data Catalog, 2237, 2002.
[36] Xiaoming Zhu and Joseph M Kahn. Free-space optical communication through atmosphericturbulence channels. IEEE Transactions on communications, 50(8):1293–1300, 2002.
[37] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery. Numerical
recipes in c++. The art of scientific computing, 2:1002, 1992.
[38] Yoshito H Ono, Carlos M Correia, Dave R Andersen, Olivier Lardière, Shin Oya, Masayuki
Akiyama, Kate Jackson, and Colin Bradley. Statistics of turbulence parameters at maunakea
using the multiple wavefront sensor data of raven. Monthly Notices of the Royal Astronomical
Society, 465(4):4931–4941, 2016.
[39] OJD Farley, J Osborn, T Morris, T Fusco, B Neichel, C Correia, and RW Wilson. Identifying
optical turbulence profiles for realistic tomographic error in adaptive optics. Monthly Noticesof the Royal Astronomical Society, 2019.
[40] OJD Farley, J Osborn, T Morris, M Sarazin, T Butterley, MJ Townson, P Jia, and RW Wilson.

参考文献

69

Representative optical turbulence profiles for eso paranal by hierarchical clustering. Monthly
Notices of the Royal Astronomical Society, 481(3):4030–4037, 2018.

70

付録 A

大気ゆらぎの統計的取り扱い
A.1 ゆらぎと統計量
A.1.1 ゆらぎ (Fructuation)ある物理量 A に対して、そのゆらぎ ξA は以下で定義される。

ξA (⃗x) = A(⃗x)− < A(⃗x) >(A.1)

ここで、< A > は考える空間の範囲内での平均である。

A.1.2 自己相関関数 (Auto Correlation Function)
ある位置 ⃗
x における物理量とそこから ⃗r だけ離れた位置 ⃗x + ⃗r における物理量との間にどれだけの相関があるかを表す量として、自己相関関数が考えられる。物理量 A の自己相関関数を次で定義する。

BA (⃗x) =< A(⃗r)A(⃗r + ⃗x) >(A.2)

したがって、ゆらぎの自己相関関数は以下でかける。

BξA (⃗x) =< ξA (⃗r)ξA (⃗r + ⃗x) >(A.3)⃗x = ⃗0 とした時、これは A の分散に一致する。
2
BξA (⃗0) =< ξA (⃗r)2 >= σA(A.4)

また、無限に離れた二点の間に相関はないものとして

BξA (∞) = 0(A.5)とできる。ゆらぎの自己相関関数を分散で規格化して表現することもある。

ρξA (⃗x) =BξA (⃗x)
BξA (⃗x)=
2⃗
σA
BξA (0)

この表現をすると、ゆらぎの自己相関関数の値は −1 ∼ +1 の範囲を取ることになる。

(A.6)

付録 A 大気ゆらぎの統計的取り扱い

71

A.1.3 相互相関関数 (Cross Corelation Function)
自己相関関数と同様に、物理量 A, B の相互相関関数は以下のように定義される。

BA,B (⃗x) =< A(⃗r)B(⃗r + ⃗x) >(A.7)

したがって、ゆらぎの相互相関関数は

BξA,ξB (⃗x) =< ξA (⃗r)ξB (⃗r + ⃗x) >

(A.8)

であり、⃗
x = ⃗0 の時、これは A, B の共分散に一致する。

BξA,ξB (⃗0) =< ξA (⃗r)ξB (⃗r) >= Cov(A, B)(A.9)

ゆらぎの相互相関関数を共分散で規格化して表現することもある。

ρξA,ξB (⃗x) =BξA,ξB (⃗x)
BξA,ξB (⃗x)=⃗
Cov(A, B)
BξA,ξB (0)(A.10)

この表現をすると、ゆらぎの相互相関関数の値は −1 ∼ +1 の範囲を取ることになる。

A.1.4 自己構造関数 (Auto Structure Function)
物理量の相関を表す量として構造関数を考えることもできる。物理量 A の自己構造関数を次で定義する。

DA (⃗x) =< (A(⃗r) − A(⃗r + ⃗x))2 >(A.11)DA (⃗x) =< A2 (⃗r) + A2 (⃗r + ⃗x) − 2A(⃗r)A(⃗r + ⃗x) >(A.12)

すなわち、= 2(BA (⃗0) − BA (⃗x))(A.13)

ゆらぎの自己構造関数は以下でかける。

DξA (⃗x) =< (ξA (⃗r) − ξA (⃗r + ⃗x))2 >(A.14)=< (A(⃗r)− < A > −A(⃗r + ⃗x)+ < A >)2 >(A.15)= DA (⃗x)(A.16)

つまり、物理量 A のゆらぎの自己構造関数は A 自身の自己構造関数に一致する。つまり、平均値が時間変動するような物理量に対して構造関数は有用である。また、BξA (∞) = 0 であることを用いると、DξA (∞) = 2BξA (⃗0) なのでBξA (⃗x) =が成り立つ。

1
1
DξA (∞) − DξA (⃗x)
2
2

(A.17)

付録 A 大気ゆらぎの統計的取り扱い

72

A.1.5 相互構造関数 (Cross Structure Function)
自己構造関数と同様に、物理量 A, B の相互構造関数は以下のように定義される。

DA,B (⃗x) =< (A(⃗r) − B(⃗r + ⃗x))2 >(A.18)

DA,B (⃗x) =< A2 (⃗r) + B 2 (⃗r + ⃗x) − 2A(⃗r)B(⃗r + ⃗x) >(A.19)

すなわち、= BA (⃗0) + BB (⃗0) − 2BA,B (⃗x)

(A.20)

ゆらぎの相互構造関数は
22
DξA,ξB (⃗x) = σA
+ σB
− 2BξA,ξB (⃗x)

(A.21)

2
2
DξA,ξB (⃗0) = σA
+ σB
− 2Cov(A, B)(A.22)であり ⃗x = ⃗0 の時である。

A.1.6 パワースペクトラム (Power Spectrum)
ある物理量のパワースペクトラムはその物理量のフーリエ変換の大きさの二乗で定義される。

WA (f⃗) = |F(A(⃗x))|2(A.23)

パワースペクトラムと自己相関関数との間には以下のような関係がある。(Wiener-Khinchine theorem)

BA (⃗x) = F(WA (f⃗))(A.24)

したがって、パワースペクトラムを用いて物理量 A の分散を表現することもできて、∫
2σA= BξA (⃗0) = F(WξA (f⃗))|⃗x=⃗0 =WξA (f⃗)df⃗(A.25)となる。

A.2 大気ゆらぎによる影響
無限遠にある天体からやってくる波長 λ の光の波面が、高さ h ∼ δh にある薄い大気乱流層によって乱されて我々の元に届く状況を考える。

付録 A 大気ゆらぎの統計的取り扱い

73

A.2.1 光の状態の記述
光の状態は複素振幅 Ψ であらわす。高さ z 、その高さの平面内での位置 ⃗x における波の複素振幅を

Ψ = Ψz (⃗x)

(A.26)とかくことにする。ここで無限遠において複素振幅が 1（振幅が 1、位相が 0）となるような規格化をすると、

Ψ∞ (⃗x) = 1(A.27)

Ψh+dh (⃗x) = 1(A.28)とかける。また、乱流層を通ることで位相が変化するため、

Ψh (⃗x) = exp[iϕ(⃗x)]∫ h+dhϕ(⃗x) = kn(⃗x, z)dz

(A.29)(A.30)hとかける。n は屈折率のゆらぎである。乱流層を抜けた後の光については、フレネル回折の効果により、

(
)]
x2
1
exp iπΨ0 (⃗x) = [Ψh (⃗x)] ∗
iλh
λh
[

(A.31)とかける。記号 ∗ は畳み込み積分を表す。以上をまとめると、

Ψ∞ (⃗x) = 1(A.32)

Ψh+dh (⃗x) = 1(A.33)Ψh (⃗x) = exp[iϕ(⃗x)]
[
(
)]
1
x2
Ψ0 (⃗x) = Ψh (⃗x) ∗exp iπiλhλh(A.34)(A.35)である。

A.2.2 様々な物理量のゆらぎと統計量
天体からの光に関する物理量とそのゆらぎの統計量を列挙する。

1. 乱流層の屈折率 @z = hξN (⃗x) ≡ N (⃗x)− < N >(A.36)

2
DξN (⃗x) = CN(h)x2/3(A.37)

2WξN (f⃗) = 9.7 × 10−3 CN(h)f −11/3(A.38)(A.39)

ξN は定義、DN は Obukhov’s low から、WN はコルモゴロフモデルから決まっている。

付録 A 大気ゆらぎの統計的取り扱い

74

2. 乱流層中の光路 @z = h(ただし、基準は z = h + dh)∫ h+dhl(⃗x) ≡

N dh

(A.40)

ξN dh(A.41)h∫ h+dhξl (⃗x) =h2Dξl (⃗x) = 2.91CN
dhx5/3(A.42)

2Wξl (f⃗) = dh × WξN = 9.7 × 10−3 CN(h)dhf −11/3(A.43)＜導出＞
まず、Bξl (⃗x) ≡< ξl (⃗r)ξl (⃗r + ⃗x) >∫ h+dh ∫ h+dh=dz
dz ′ < ξN (⃗r, z)ξN (⃗r + ⃗x, z ′ ) >h∫ h+dh=dz

dζBξN (⃗x, ζ)(A.46)

h−z
∫ ∞

dz
h

(A.45)h∫ h+dh−zh∫ h+dh∼(A.44)

∫ ∞

= dh
−∞

−∞

dζBξN (⃗x, ζ)dζBξN (⃗x, ζ)(A.47)(A.48)であるので、Dξl (⃗x) ≡ 2(Bξl (⃗0) − Bξl (⃗x))∫= 2dh dζ[BξN (⃗0, ζ) − BξN (⃗x, ζ)]∫= dh dζ[DξN (⃗x, ζ) − DξN (⃗0, ζ)]∫2= dh dζCN(ζ)[(x2 + ζ 2 )1/3 − ζ 2/3 ]∫2∼ dhCN(h) dζ[(x2 + ζ 2 )1/3 − ζ 2/3 ]2= 2.91CN(h)dhx5/3(A.49)(A.50)(A.51)(A.52)(A.53)(A.54)であり、Wξl (f⃗) = F −1 [Bξl ]∫ ∞= dhdζF −1 [BξN (⃗x, ζ)]−∞∫ ∞dζBξN (f⃗, ζ)= dh(A.55)(A.56)(A.57)−∞= dhWξN (f⃗, 0)である。

(A.58)

付録 A 大気ゆらぎの統計的取り扱い

75

3. 光の位相 @z = h(ただし、基準は z = h + dh)∫ h+dhϕ(⃗x) ≡ kl = k

N dh(A.59)h∫ h+dhξϕ (⃗x) = kξl = k

ξN dh
h

2

(A.60)(

2

2
Dξϕ (⃗x) = k Dξl = 2.91k CN
dhx5/3 = 6.88

x
r0

) 53

2
Wξϕ (f⃗) = k 2 Wξl = 9.7 × 10−3 k 2 CN(h)dhf −11/3(A.61)(A.62)

4. 光の複素振幅 @z = hΨ(⃗x) = exp[iξϕ ])(1BΨ (⃗x) = exp − Dξϕ
2
(
(
))
1
DΨ (⃗x) = 2 1 − exp − Dξϕ2(A.63)(A.64)(A.65)＜導出＞
まず、BΨ (⃗x) ≡< Ψ(⃗r)Ψ∗ (⃗r + ⃗x) >(A.66)=< exp[i(ξϕ (⃗r) − ξϕ (⃗r + ⃗x))] >(A.67)≡< exp[iG(⃗x)] >(A.68)∞ n∑i< Gn >
n!
n=0
]
∞ [
∑
(−1)m
i(−1)m
2m
2m+1
=
<G >+
<G
>
(2m)!
(2m + 1)!
m=0

=

∞
∑
(−1)m

< G2m >
(2m)!
m=0
(
)m
∞
∑
1
1
2
=
− <G >
m!
2
m=0
(
)m
∞
∑
1
1
2
=
− < (ξϕ (⃗r) − ξϕ (⃗r + ⃗x)) >
m!
2
m=0
]
[
1
= exp − Dξϕ (⃗x)
2=

(A.69)(A.70)(A.71)(A.72)(A.73)(A.74)であるので、DΨ (⃗x) = 2[BΨ (⃗0) − BΨ (⃗x)](())
1= 2 1 − exp − Dξϕ (⃗x)
2である。

(A.75)(A.76)

付録 A 大気ゆらぎの統計的取り扱い

76

5. 光の複素振幅 @z = 0(ただし、位相の基準は z = h、振幅は 1 に規格化)()
1x2
ψ(⃗x) = 1 + ξϕ ∗exp iπ
λh
λh
(
)
2
1
x
ξψ (⃗x) = ξϕ ∗exp iπλhλh⃗Wξψ (f ) = Wξϕ(A.77)(A.78)(A.79)

Bψ (⃗x) = BΨ(A.80)

Dψ (⃗x) = DΨ

(A.81)＜導出＞

z = 0 での計算をするにあたって、一つの大きな近似をする。それは弱乱流近似と呼ばれ ξϕ (⃗x) <<1 というものである。この近似は天文観測において、天頂方向についてはよく成り立つ近似である。（天頂角が 60 度を越えると成り立たない [Young]）また、位相の基準を z = h に取り直す。

(
)]
x2
1
exp iπψ(⃗x) = Ψ(⃗x) ∗
iλh
λh
[
(
)]
1
x2
= exp[iξϕ ] ∗exp iπ
iλh
λh
[
(
)]
1
x2
∼ (1 + iξϕ ) ∗exp iπ
iλh
λh
[
(
)]
[
(
)]
2
1
1
x
x2
=1∗
exp iπ+ iξϕ ∗
exp iπ
iλh
λh
iλh
λh
(
)]
[
2
x
1
exp iπ= 1 + ξϕ ∗λh
λh[

(A.82)(A.83)(A.84)(A.85)(A.86)

ただし最後の変形で、任意の関数 g(x) に対して

1 ∗ g(x) = F −1 [F[1 ∗ g(x)]]=F−1[F[1] × F[g(x)]]= F −1 [δ(f ) × g̃(f )]=F−1[δ(f )]=1(A.87)(A.88)(A.89)(A.90)(A.91)であることを用いた。すると、

Bψ (⃗x) =< ψ(⃗r)ψ ∗ (⃗r + ⃗x) >

[

(

)] [
(
)]
1
x
1
x2
∗
=< Ψ(⃗r)Ψ (⃗r + ⃗x) > ∗exp iπ∗exp −iπ
iλh
λh
−iλh
λh
∗
=< Ψ(⃗r)Ψ (⃗r + ⃗x) > ∗δ(⃗x)(A.92)

2= BΨ

(A.93)(A.94)(A.95)であるので、Dψ = DΨ
である。パワースペクトラムについては後で証明する。

(A.96)

付録 A 大気ゆらぎの統計的取り扱い

77

6. 光の実数振幅と位相 @z = 0
[ 2]
1
x
cos π
λh
λh
[ 2]
1
x
θ(⃗x) ≡ Arg[ψ] = ξϕ ∗
sin π
λh
λh
[ 2]
1
x
ξχ (⃗x) = ξϕ ∗
cos π
λh
λh
[ 2]
1
x
ξθ (⃗x) = ξϕ ∗
sin π
λh
λh
2
⃗
Wξχ (f ) = Wξϕ × sin (πλhf 2 )χ(⃗x) ≡ |ψ| = 1 + ξϕ ∗(A.97)(A.98)(A.99)(A.100)(A.101)Wξθ (f⃗) = Wξϕ × cos2 (πλhf 2 )(A.102)

ψ = 1 + ξψ = χexp[iθ](A.103)

1 + Re[ξψ ] = χcos[θ](A.104)Im[ξψ ] = χsin[θ](A.105)

＜導出＞

と置くと、

となるので、

χ = ((1 + Re[ξψ ])2 + (Im[ξψ ])2 )1/2

(A.106)= (1 + 2Re[ξψ ] + |ξψ | )(A.107)∼ (1 + 2Re[ξψ ])(A.108)∼ 1 + Re[ξψ ](A.109)

2 1/2

1/2

[ 2]
1
x
= 1 + ξϕ ∗
cos π
λh
λh
(
)
−1
θ = tanIm[ξψ ](1 + Re[ξψ ])−1
∼ tan

−1

∼ tan

−1

(A.111)(Im[ξψ ](1 − Re[ξψ ]))(A.112)(Im[ξψ ])(A.113)∼ Im[ξψ ]= ξϕ ∗(A.110)

[

]
2

1
x
sin π
λh
λh

(A.114)(A.115)である。またパワースペクトラムについては、畳み込みのフーリエ変換に対して、F T (f ∗ g) = F T (f ) × F T (g)(A.116)

付録 A 大気ゆらぎの統計的取り扱い

78が成り立つことから、複素振幅、振幅、位相のゆらぎのパワースペクトラムはそれぞれ、

(

[
]) 2
1
x2
exp iπλhλh([ 2 ]) 2
1xWξχ (f⃗) = Wξϕ (f⃗) × F T
cos π
λh
λh
[ 2 ]) 2(x
1Wξθ (f⃗) = Wξϕ (f⃗) × F T
sin π
λh
λh

Wξψ (f⃗) = Wξϕ (f⃗) × F T(A.117)(A.118)(A.119)と表せる。ここで、

[
])
1
x2
exp ±iπ
FT
λh
λh
[
]
∫
∫
1
x2 + y 2= dx dy exp ±iπ
exp [−2πi(fx x + fy y)]
λh
λh
(∫
[
])2
x2
1
dxexp ±iπ
− 2πifx x
=
λh
λh
(∫
[
])2
1
π
2
2
=
dxexp ±i (x ∓ λhfx ) ∓ iπλhfx
λh
λh
(∫
)2
[
[
]
π 2]
1
2
2
exp ∓iπλh(fx + fy )dxexp ±i x
=
λh
λh
(√ )2
[
]
1
λh
=
exp ∓iπλhf 2λh∓i(= ±iexp[∓iπλhf 2 ](A.121)(A.122)(A.123)(A.124)(A.125)(A.126)= sin(πλhf ) ± icos(πλhf )2(A.120)2(A.127)であることから、

(

[
])
x2
1
exp iπ= iexp[−iπλhf 2 ]
FT
λh
λh
(
[ 2 ])
1x
FT
cos π= sin(πλhf 2 )λhλh([ 2 ])
1x
FT
sin π= cos(πλhf 2 )λh
λh(A.128)(A.129)(A.130)であり、ゆらぎのパワースペクトラムはWξψ (f⃗) = Wξϕ (f⃗)(A.131)2Wξχ (f⃗) = Wξϕ (f⃗) × sin (πλhf )(A.132)Wξθ (f⃗) = Wξϕ (f⃗) × cos (πλhf )(A.133)Wξψ (f⃗) = Wξχ (f⃗) + Wξθ (f⃗)(A.134)222と書くことができる。
ここで、重要な関係式として

付録 A 大気ゆらぎの統計的取り扱い

79が得られた。
また、ここでは証明は省くが多層の場合の結果のみを載せておく。多層の場合は各層からの寄与の線型結合になる。したがってWξψ (f⃗) = 9.7 × 10−3 k 2 f −11/3Wξχ (f⃗) = 9.7 × 10−3 k 2 f −11/3Wξθ (f⃗) = 9.7 × 10−3 k 2 f −11/3

∫
∫
∫

2
CN
(h)dh

(A.135)2CN(h)dhsin2 (πλhf 2 )(A.136)2CN(h)dhcos2 (πλhf 2 )(A.137)という関係式が成り立つ。また、高周波数側では f −11/3 の影響でパワーが小さくなるため、いま
低周波数側だけ考えて πλhf 2 << 1 であるような近似（near field approximation という）をすると、Wξχ (f⃗) ∼ 9.7 × 10−3 k 2 f −11/3∫2CN(h)dh(πλhf 2 )2
∫
−3
4 1/3
2
= 9.7 × 10 × 4π f
h2 CN
(h)dh
∫
2
Wξθ (f⃗) ∼ 9.7 × 10−3 k 2 f −11/3 CN(h)dh= Wξψ (f⃗)(A.138)(A.139)(A.140)(A.141)

7. 光の明るさ（明るさの平均値に対する相対値）@z = 0I(⃗x) ∼ 1 + 2ξχ

[ 2]
x
1
cos π
ξI (⃗x) = 2ξχ = 2kξl (⃗x) ∗λhλhξ˜I (f⃗) = 2ξ˜χ = 2k ξ˜l (f⃗) × sin(πλhf 2 )WξI (f⃗) = 4Wξχ = 4k 2 Wξl (f⃗)sin2 (πλhf 2 )()2sin(πλhf 2 )−11/3
2= 1.53fCN(h)dh
λ(A.142)(A.143)(A.144)(A.145)(A.146)＜導出＞
天体からの光の強度（intensity）i(⃗x) は実数振幅の 2 乗に比例するので比例定数を α として、

i(⃗x) = αχ2 (⃗x)

(A.147)i(⃗x) = α(1 + ξχ )2(A.148)∼ α(1 + 2ξχ )(A.149)とかける。χ = 1 + ξχ であるから、

アンサンブル平均（時間平均）をとれば、< ξχ >= 0 なので、< i(⃗x) >∼ α

(A.150)

付録 A 大気ゆらぎの統計的取り扱い

80

したがって、i(⃗x)= 1 + 2ξχ
< i(⃗x) >(A.151)

左辺を I(⃗x) とおいているので上述のような式になる。

8. 波面の傾き @z = 0
x 方向の傾きを α、y 方向の傾きを β とする。
(
[ 2 ])λ ∂λ ∂1xα(⃗x) = −θ(⃗x) = −
ξϕ ∗
sin π
2π ∂x
2π ∂x
λh
λh
(
[ 2 ])λ ∂λ ∂1xξθ (⃗x) = −
ξϕ ∗
sin π
ξα (⃗x) = −
2π ∂x
2π ∂x
λh
λh
2
2
Wξα (f⃗) = λ f Wξθ (f⃗)
x
2

λ ∂2
λ2 ∂ 2λ2 ∂ 2
B
(⃗
x
)
=
D
(⃗
x
)
∼
Dξϕ (⃗x)ξθξθ
4π 2 ∂x2
8π 2 ∂x[28π 2 ∂x2
]
λ2
5
1 2 2
−5/3
2
2 −1/6
2 −7/6=× 6.88 × r0
× × (x + y )− x (x + y )
8π 233Dξα (⃗x) = 2(Bξα (⃗0) − Bξα (⃗x))
Bξα (⃗x) = −(A.152)(A.153)(A.154)(A.155)(A.156)(A.157)