CPS/WTK惑星大気放射計算ミニワークショップ
神戸大学CPS 2018.3.9 13:05-14:05「あかつきIR2データと放射計算による金星雲頂構造の推定」

佐藤隆雄
ISAS/JAXA
1© ISAS/JAXA

放射伝達モデル2

放射伝達モデル
 理論 [平行平面大気を仮定]
放射伝達方程式
1. 光の消散 2𝜋

2. 多重散乱1

太陽反射光： 1-3
熱放射： 1,2,4𝑑𝐼 𝜏, 𝜇, 𝜙𝜛
𝜇= 𝐼 𝜏, 𝜇, 𝜙 −න න 𝑃 𝜏, 𝜇, 𝜙, 𝜇′ , 𝜙 ′ 𝐼 𝜏, 𝜇′ , 𝜙 ′ 𝑑𝜇′ 𝑑𝜙 ′𝑑𝜏
4𝜋 0 −1𝜛𝜏−𝑃 𝜏, 𝜇, 𝜙, −𝜇0 , −𝜙0 𝜋𝐹0 𝑒𝑥𝑝 −− 1−𝜛 𝐵 𝑇4𝜋𝜇0

3. 太陽光の一次散乱

4. 熱放射

𝐼: 放射輝度, 𝜏: 光学的厚さ, 𝜛: 一次散乱アルベド, 𝑃: 散乱位相関数,𝜋𝐹0 : 太陽直達光, 𝜇, 𝜙 : 光の入射(射出)方向, 𝐵 𝑇 : プランク関数

𝐼 𝜏, 𝜇, 𝜙 のままでは数値計算が困難⇒放射輝度𝐼は𝜏だけの関数に．
𝜙： フーリエ級数展開
𝐼 𝜏, 𝜇, 𝜙 = ෍𝑚

2 − 𝛿0𝑚 𝐼𝑚 𝜏, 𝜇 𝑐𝑜𝑠𝑚 𝜙 − 𝜙′𝜇： ガウス・ルジャンドル求積法で離散化

𝐼𝑚 𝜏, 𝜇 = ෍ 𝐼𝑚,𝑖 𝜏 𝜇𝑖 𝑤𝑖𝑖3

放射伝達モデル
 ”doubling and adding method” [e.g., Goody and Young, 1983]-散乱大気の放射伝達過程を扱う代表的な手法の一つ
多重散乱の取り扱い方： 層を次々に合成 (C層=A層+B層)(𝑅: 反射行列， 𝑇: 透過行列)
行列計算
1回散乱＋多重散乱を経た反射光

R 𝐶 T𝐶R A TAR∗A TA∗R B TBR∗𝐵 T𝐵∗1本の反射光として記述する

A層
C層
B層
R∗𝐵 T𝐵∗

放射伝達モデル
 元の層まで、上下の層を足し合わせていくのがadding法，足し合わせる層が同一の光学的特性をもつとき、特にdoubling法と呼ぶ。

 では初期層はどうとるか？
1回散乱程度の薄い層 𝜏0 ≈ 10−6 まで層を分割して出発。
∵解析的に解が計算できる (この場合は，𝑅 = 𝑅∗ ，𝑇 = 𝑇 ∗ )𝜛𝑅 𝜏0 , 𝜇, 𝜙, 𝜇′, 𝜙′ =𝑃 𝜏, 𝜇, 𝜙, −𝜇′ , 𝜙 ′4 𝜇 + 𝜇′1 − 𝑒𝑥𝑝 −𝜏01 1+𝜇 𝜇′𝜛𝜏0𝜏0𝑃 𝜏, −𝜇, 𝜙, −𝜇′ , 𝜙 ′ 𝑒𝑥𝑝 −− 𝑒𝑥𝑝 − ′4 𝜇 − 𝜇′𝜇𝜇𝑇 𝜏0 , 𝜇, 𝜙, 𝜇′, 𝜙′ =𝜛𝜏0𝜏′ , 𝜙 ′ 𝑒𝑥𝑝 − 0𝑃
𝜏,−𝜇′,𝜙,−𝜇4𝜇′2𝜇′

光学的厚さ，一次散乱アルベド，散乱位相関数：入力パラメータとして必要

5放射伝達モデル
 放射計算コード• 散乱大気中での太陽放射 (Sato et al., 2013)• 散乱大気中での熱放射• 灰色大気中での熱放射 (Sato et al., 2014)➢ 入力モジュールを変更すると，分子・エアロゾルの光学的特性を計算しLine-byLineを実行．
 数値計算上考慮していること• Renormalization of scattering phase function (Hansen, 1971)1න 𝑃 𝜏, 𝜇, 𝜙, 𝜇′ , 𝜙 ′ 𝑑Θ = 14𝜋 4𝜋• Second order scattering in homogeneous layer (Hovenier et al., 2004)
初期層において，二次散乱まで解析的に求める➢ 初期層の厚さを厚くできる，特殊ケース 𝜛 = 1 でも精度を維持できる．

• Separation of single scattering from azimuthal Fourier series of reflection andtransmission matrices (de Haan et al., 1987)
波長より粒径が大きいエアロゾルを考慮する場合に効果あり

𝑅 𝜏0 , 𝜇, 𝜙, 𝜇′, 𝜙′= 𝑅1 𝜏0 , 𝜇, 𝜙, 𝜇′, 𝜙′ + ෍ 2 − 𝛿0𝑚 𝑅 𝑚 𝜏0 , 𝜇, 𝜇′ − 𝑅1𝑚 𝜏0 , 𝜇, 𝜇′ 𝑐𝑜𝑠𝑚 𝜙 − 𝜙′𝑚6放射伝達モデル
 放射計算コードの精度検証：金星版Rstarとの比較

図．(左) Akatsuki/IR2 2.02 μm，(右) Akatsuki/LIRの場合の精度検証．

IR2, LIRのどのチャンネルでも0.5 % (*) 以下で整合* Rstarの雲の光学テーブルが粗かったということでまだ改善の余地あり． 7

「あかつき」IR2昼面データを用いた
金星雲頂構造の推定

 最終ゴール：• ピクセルごとに反射太陽光の明るさを雲頂高度に焼き直したい．

8

はじめに雲層雲層図．(左) PV probe entry siteにおける金星雲の高度分布 [Knollenberg and Hunten, 1980],(中) PV probe entry siteにおける入力太陽光フラックスの高度分布 [Tomasko et al., 1980],(右) PV/OCPPの紫外画像 [Rossow et al., 1980].

 厚い硫酸雲 (50-70 km) [τ ~ 20-40@可視光]• 入射太陽光の~76%が雲層で反射される [Crisp and Titov, 1997]．
• 大気・地表で吸収する太陽光エネルギーの~50%が (主に紫外未知吸収物質によって) 雲層 (>64 km)で吸収される [Tomasko et al., 1980]．
金星における雲は，大気化学，大気運動，熱収支を考える上で重要な役割を果たしている．

9

金星探査機「あかつき」

金星周回軌道から雲の上層から下層まで透視する3次元的なリモートセンシングによって金星の気候のしくみに迫る• スーパーローテーションはな
ぜ起こるのか• 大気は上下，南北方向にどう循環しているのか• 金星全体を隙間なく覆う雲はどう作られるのか• 雷は起こるのか• 活火山はあるか10

「あかつき」の科学機器

中間赤外カメラ (LIR)

雷・大気光カメラ (LAC)

紫外イメージャ (UVI)
1 μmカメラ (IR1)2 μmカメラ (IR2)超高安定発信器 (USO)「あかつき」の諸元• 1.04 m×1.45 m×1.40 m• 502 kg (内訳：バス機器56%，燃料38%，科学機器7%)
11

「あかつき」 IR2が目指す科学 Data set
2.02 μm• Internal released version

2.26, 2.32 μm
1.735 μm

図．夜面2波長 (1.735 μmと2.26 μm) を用いた金星下層雲の疑似画像.

図．観測波長と感度高度の関係 [Taylor, 1998].
図．IR2の外観 [Satoh et al., 2016]．

• 2.02 μm： CO2吸収帯における反射太陽光➢ 雲頂高度の可視化，金星昼面のみ観測可• 1.735, 2.26, 2.32 μm： 大気の窓，雲層下から射出される熱放射➢ 雲の光学的厚さ，CO分布の可視化，金星夜面のみ観測可12

「あかつき」 IR2の観測記録
 Data set• Internal released version
周回軌道
後初観測

定常観測
開始

昼面観測
シーズン
(4月–7月)夜面観測
シーズン(8月–10月)
IR-AE
異常停止(12月9日)図．2015–2016年のIR2観測記録．

13

先行研究に学ぶ「あかつき」データ解析へのヒント
朝側境界

夕側境界

649.4 nm
1553.7 nm
図．(左) VEX/VIRTIS熱放射データに基づく雲頂高度の緯度分布 [Haus et al., 2014],(右) VEX/SOIR太陽掩蔽データに基づく雲頂高度の緯度分布 [Luginin et al., 2016].

 欧州金星探査機Venus Express• 赤道域における平均雲頂高度はほぼ一定 [Ignatiev et al., 2009; Haus et al., 2014]• 平均雲頂高度の顕著なローカルタイム依存性はない [Luginin et al., 2016]➢ 平均的な雲モデル (ベースモデル) は，観測された赤道域平均の放射輝度の太陽位相角依存性 (位相曲線) を一様の雲頂高度で再現する必要がある．
14

観測データ 使用したデータ• 2.02 μmで取得した93画像(軌道番号12–16, 2016/4/4–5/25)• 公開データ (VCO-V-IR2-3-CDR-V1.0)• デコンボリューションと低周波フラット成分を追加処理• 各画像について，下記条件を満たすピクセルの放射輝度を平均✓ 𝜇 ≥ 0.1
✓ 𝜇0 ≥ 0.1✓ −30° ≤ 𝑙𝑎𝑡𝑖𝑡𝑢𝑑𝑒 ≤ 30°

図．衛星直下点における太陽位相角に対する(上) 視直径(中) 衛星直下点におけるローカルタイム(下) 赤道域平均の放射輝度

Orbit 12

Orbit 15

Orbit 13

Orbit 16

Orbit 1415

放射伝達モデル：入力パラメータ

 入力パラメータ
図．(左上) VIRAモデルより赤道域の温度プロファイル [Seiff et al., 1985](右上) 2.02 μm解析のために考慮した大気分子の高度プロファイル [Marcq et al.,2005; Haus et al., 2010](左下) 雲モデル [Haus et al., 2015].
16放射伝達モデル：入力パラメータ
Molecules

表．大気分子の吸収断面積を計算するためのパラメータ(主にHaus et al. [2013]に基づく).
DatabaseLine shape
Chi factor

Cutoff

H2O

HITEMP2010

Voigt

1.0

125 cm-1

CO2

HITEMP2010 (dayside)
HITEMP-Venus (nightside)

Voigt

Tonkov et al.
(1986)

125 cm-1

COHITRAN2014

Voigt

1.0125 cm-1

SO2HITRAN2014

Voigt

1.0125 cm-1

HF

HITRAN2014

Voigt

1.0125 cm-1

HCl

HITRAN2014

Voigt

1.0125 cm-1

OCS

HITRAN2014

Voigt

1.0125 cm-1

N2HITRAN2014

Voigt

1.0

125 cm-1

表．IR2の観測波長におけるCO2連続吸収係数(2.02 μm 以外はHaus et al. [2013]に基づく)
Spectral channel

1.735

2.02

2.26 and 2.32

Coefficient coefficients[cm-1 amagat-2]

4.1x10-9

0.043x10-917

HITEMPデータベースの比較
HITEMP比較：CO2• 比較対象： 02_1000k, private, HITEMP2010• 10-35 (cm-1/molecule cm-2) よりも弱い線吸収データは使用しない．
表．各データベースのCO2線吸収

02_1000k

privateHITEMP2010

IR1波数域(9,000-12,000 cm-1)

167

50,769187,434

IR2波数域(4,000-6000 cm-1)

26,756

390,185

402,507

LIR波数域(600-1600 cm-1)

34,952486,249

593,978➢ 出版の順番通り，どの波数域でもHITEMP2010が最も多くの線吸収データを有している．しかし，・・・HITEMPデータベースの比較 IR2
02_1000k, private, HITEMP2010

図．IR2波長帯におけるCO2線吸収• Private: HITEMP2010にない比較的強い線吸収データが1.735,2.320 μmフィルター付近に存在する．

HITEMPデータベースの比較 IR2: 1.735 μm

相対誤差：-33.7 %図．(左) データベースの違いによる計算スペクトルの違い(右) 地上観測スペクトル (Bézard et al., 1990)

大気窓の短波長側が全く合わない．
➢ 観測スペクトルを見る限り，HITEMP2010は全く合わない印象．

Haus雲モデルを使用した場合の計算スペクトル表．4種類のモード粒子のモード半径と分散 [Pollack et al., 1993](log-normalサイズ分布を仮定).
Mode 1 Mode 2 Mode 2’

Mode 3𝑟ҧ 𝜇𝑚0.30

1.00

1.403.65

硫酸雲の複素屈折率(H2SO4:H2O=75%:25%)

𝜎𝑙𝑜𝑔

1.56

1.29

1.23

1.28

[Palmer and Williams, 1975] Haus雲モデルを使用した場合の計算
スペクトル例• LBL計算• Adding-doubling method

本研究での雲頂高度の定義：
2.02 μmで大気上端から雲の光学的厚さが1となる高度

CO2

図．IR2昼面データのための計算
スペクトル例(波長間隔：0.005 cm-1，平均の太陽金星距離を使用)．
21

解析方法
1.

Mode 1とMode 2の数密度分布𝑁 𝑧 を簡素化
• 𝑁𝑚𝑜𝑑𝑒1 𝑧 = 𝛼𝑁𝑐 𝑒𝑥𝑝 − 𝑧 − 66 Τ𝐻• 𝑁𝑚𝑜𝑑𝑒2 𝑧 = 𝑁𝑐 𝑒𝑥𝑝 − 𝑧 − 66 Τ𝐻100• 𝑁𝑐 = 1ൗ 𝛼𝜎1 + 𝜎2 ‫ 𝑝𝑥𝑒 𝑧׬‬− 𝑧 − 66 Τ𝐻 𝑑𝑧𝑐2.
3.

4.

8次元パラメータ空間 𝛼, 𝑟1 , 𝑟2 , 𝐻, 𝑧𝑐 , 𝜇, 𝜇0 , Δ𝜙 の各組合せに対して，フィルター平均輝度を計算し参照テーブルを作成．
5次元パラメータ空間 𝛼, 𝑟1 , 𝑟2 , 𝐻, 𝑧𝑐 の各組合せに対して，各画像の疑似放射輝度マップを，手順2で得たテーブルと各ピクセルの観測ジオメトリ 𝜇, 𝜇0 , Δ𝜙 を基に内挿により求める．
4次元パラメータ空間 𝛼, 𝑟1 , 𝑟2 , 𝐻 の各組合せに対して，観測で得られた赤道域の位相曲線を最も再現する𝑧𝑐 を最小二乗法により求め，𝜒 2 の最小値を計算．
𝑖𝜒2 = ෍𝑖

𝐼𝑜𝑏𝑠 − 𝐼𝑠𝑖𝑚𝜎𝑜𝑏𝑠 𝑖𝑖 2𝜎𝑜𝑏𝑠 = 1(見積もり中のため)

𝛼: Mode1とMode2の混合比，𝑁𝑐 : 雲頂高度𝑧𝑐 となるための数密度，𝐻: 雲のスケールハイ
ト，𝜎1 , 𝜎2 : Mode 1とMode 2のフィルター平均消散断面積，𝑟1 , 𝑟2 : Mode 1とMode 2のモード半径22解析方法
 任意の観測ジオメトリ 𝜇, 𝜇0 , Δ𝜙 における放射輝度を参照テーブルから求める際の精度検証• 参照テーブルの計算範囲
0.067 ≤ 𝜇 ≤ 10.067 ≤ 𝜇0 ≤ 10° ≤ Δ𝜙 ≤ 180°
➢ 合計10,164 (22×22×21) 点• スプライン補間✓ 精度検証• ランダムに発生させた1,000点の疑似観測ジオメトリに対して，A) 計算間隔：0.005 cm-1
ジオメトリ：直接B) 計算間隔：0.005 cm-1
ジオメトリ：内挿C) 計算間隔：0.1 cm-1
ジオメトリ：内挿B−AA
0.00±0.02 (%)

𝐶−AA
0.3±0.5 (%)

23
図．参照テーブルの内挿に関する精度検証．

結果• 各パラメータの感度を調べるため，9種類のシナリオについて計算した．
Mode2: スケール 雲頂高度:
𝑟2 (μm)𝑧𝑐 (km)ハイト:
𝐻 (km)

Mode 1とMode 2の混合比: 𝛼*Mode 1:𝑟1 (μm)*,
𝜎1 *

ベスト
フィット
χ2

Case1

0-

0.014935

1.20

6.070.29

Case2#

1.45(0.30, 1.56)

0.014937

1.25

6.570.09

Case35(0.30, 1.56)

0.0152701.30

6.569.82

Case450(0.30, 1.56)

0.018324

1.45

6.568.30

Case5500(0.30, 1.56)

0.035290

1.50

5.067.41

Case6

1.45(0.085, 1.45)

0.014935

1.20

6.070.29

Case75(0.085, 1.45)

0.014933

1.20

6.070.29

Case850(0.085, 1.45)

0.014934

1.20

6.070.28

Case9500(0.085, 1.45)

0.014823

1.20

6.570.08* 固定値.
(0.30, 1.56): Pollack et al. [1993]，(0.085, 1.45): Luginin et al. [2016]# Case2がHaus雲モデルに最も近いモデル.

24

結果: Case2 (標準モデル)
図．(上) 𝑟2 と𝐻に対する𝜒 2 と雲頂高度𝑧𝑐 のマップ(下) 観測とベストフィットモデルの位相曲線の比較．

Orbit 12

Orbit 15

Orbit 13Orbit 16

Orbit 14

 ベストフィットパラメータ• Mode 2のモード半径𝑟2 = 1.25 μm• 雲のスケールハイト
𝐻 = 6.5 km• 雲頂高度𝑧𝑐 = 70.09 km➢ ベストフィットモデルは観測の位相曲線をよく再現している．
25

結果: Case2 (標準モデル)
 感度チェック
𝑟2 : 1.25 μm, 𝐻: 3.0 km
𝑟2 : 1.25 μm, 𝐻: 5.0 km
𝑟2 : 1.25 μm, 𝐻: 6.5 km (ベストフィット)𝑟2 : 1.25 μm, 𝐻: 7.0 km
𝑟2 : 1.25 μm, 𝐻: 9.0 km

𝑟2 : 1.00 μm, 𝐻: 6.5 km
𝑟2 : 1.25 μm, 𝐻: 6.5 km (ベストフィット)𝑟2 : 1.50 μm, 𝐻: 6.5 km

図．(上) 雲のスケールハイトに対する位相曲線変化(下) Mode 2のモード半径に対する位相曲線変化．

• 雲のスケールハイト 𝐻➢ 高高度情報を含む中～高位相角 (>~80˚) データに感度あり➢ 上層雲はかなり広がった層構造をもつ．
• Mode 2のモード粒子 𝑟2➢ 位相角 (10-40˚) の位相曲線の形に感度あり

26結果: Mode 1 (Pollack) の場合
 Mode 1 (Pollack) を使った場合のMode 1とMode 2の混合比依存性

α=1.45α=1.45
α=5
α=50
α=500

α=5

ほぼ同じ

α=50α=500

図．(上) 𝑟2 と𝐻に対する𝜒 2 マップのα依存性(下) 観測とベストフィットモデルの位相曲線の比較．

Mode 1 (Pollack) を増やすと，フィットは悪くなる➢ αが小さい (α≤5) 場合は，IR2データのみから混合比を決めるのは不可能．
27

結果: Mode 1 (Luginin) の場合
 Mode 1 (Luginin) を使った場合のMode 1とMode 2の混合比依存性

α=1.45α=1.45
α=5α=50
α=500α=5

ほぼ同じ

α=50α=500

図．(上) 𝑟2 と𝐻に対する𝜒 2 マップのα依存性(下) 観測とベストフィットモデルの位相曲線の比較．

Mode 1 (Luginin) をかなり増やしても フィットは変わらない．
➢ IR2データのみから混合比を決めるのは不可能．
28

議論，先行研究との比較
 Mode 1のモード半径とMode 1とMode 2の混合比への感度• 消散断面積 (2 μm) (m2/particle): Mode 2 (i.e., 𝑟2 =1.0 μm, 𝜎2 =1.29) = 1.0×10-11
Mode 1 (Pollack) = 2.7×10-13, Mode 1 (Luginin) = 2.3×10-16➢ Mode 1 (Luginin): Mode 2に対して十分無視できるレベル➢ αが小さい (α≤5) 場合は，IR2昼面データでMode 1の光学特性は制約できない.
 Mode 2のモード半径，雲のスケールハイト，雲頂高度• Rossi et al. [2015]: VEX/SPICAV, 近赤外昼面偏光分光→ 𝑟2 ~0.85 μm• Luginin et al. [2016]: VEX/SOIR, 可視近赤外太陽掩蔽分光→ 𝑟2 =0.66 μm• Lee et al. [2012] VEX/VIRTIS+VeRa, 近赤外夜面分光→ H=2.2-5.4 km• Luginin et al. [2016]: VEX/SOIR,可視近赤外太陽掩蔽分光→ H=3-4 km✓ 典型的な (既知として使われてきた) Mode 2のモード半径 (𝑟2 =1 μm) を使うと，スケールハイトは5 kmになる.
• Cottini et al. [2012]: VEX/VIRTIS-H, 近赤外昼面分光→ 𝑧𝑐 =69.5 km at 2.5 μm• Haus et al. [2014]: VEX/VIRTIS, 近赤外夜面分光→ 𝑧𝑐 =71 km at 1 μm
• Fedrova et al. [2016]: VEX/SPICAV, 近赤外昼面分光→ 𝑧𝑐 =70.2 km at 1.48 μm✓ 導出した雲頂高度 (𝑧𝑐 ~70 km) は整合的，Mode 1に対するモデル依存度も低い．
29

まとめと今後の課題
 まとめ• 各画像に対して雲頂高度マップを求めるため，「あかつき」IR2の2.02 μmの昼面画像 (2016年4月–5月) を放射伝達計算を基に解析している．
• 第一段階として，赤道域の平均放射輝度の位相曲線を再現するベースとなる雲モデルを導出した．
• 標準モデル (Case2) の場合，𝑟2 =1.25 μm, 𝐻=6.5 km, 𝑧𝑐 =70.09 kmの組み合わせでベストフィットが得られた．雲頂高度は先行研究と整合的だが，Mode 2のモード半径や雲のスケールハイトは先行研究よりも大きい結果となった．Mode 1の光学特性や混合比は，IR2昼面データ単独では制約できないことが分かった．
 今後の課題• 観測誤差 (バイアス，ランダム) の見積もり➢ 解の許容範囲の決定に重要• 上記ベースモデルを基に，雲の光学的厚さ (＝雲頂高度) のみ を調節することで，各画像の雲頂高度マップを作成する• IR2 (雲頂高度) と UVI (未知吸収物質) の模様の比較• その他のパラメータ (温度プロファイル) に対する雲モデルの感度検証30