つくばにおける FTIR 観測による微量気体高度分布の解析：CO 発生源の推定

つくばにおける FTIR 観測による微量気体高度分布

の解析：CO 発生源の推定

著者

鈴木 大将

学位授与機関

Tohoku University

修士論文

つくばにおける FTIR 観測による微量気体高度分布の解析

：CO 発生源の推定

東北大学大学院 環境科学研究科

先端環境創成学専攻 文化環境学コース

鈴木 大将

令和 3 年

指導教員

中島 英彰 客員教授

審査員

（〇：主査）

〇 中島 英彰 客員教授

1 土屋 範芳 教授

2 町田 敏暢 客員教授

3 村田 功 准教授

iii

要旨

一酸化炭素（CO）は主に化石燃料の燃焼やバイオマス燃焼, 炭化水素の酸化によって生成 される. CO は対流圏オゾン発生の前駆物質であるとともに，OH ラジカルとの反応によっ て様々な大気微量成分の存在量に影響を与えるため大気化学において極めて重要な物質で ある．大気中では数週間～２ヶ月程度の寿命を持ち，その存在量や発生源には地域的な偏 りがある．近年 CO の全球的な濃度変化に大きな影響を与えているアジア地域ではより正 確な変動の理解が必要であり，日本での観測はアジア地域の CO の理解に貢献できる．つ くば（36.05°N, 140.12°E）では、2001 年 5 月から高分解能フーリエ変換赤外分光計 （FTIR）による大気微量成分の観測を行っている．FTIR はその高波数分解能 （0.0035cm-1）を生かし，観測した吸収スペクトルの線形から存在量の高度分布のリトリ ーバルが可能である．本研究では，つくば上空の CO の発生源および発生地域を明らかに することを目的として，CO およびシアン化水素（HCN）の高度分布を用いた CO 発生源 の推定を行った．大気中の HCN は数年程度の寿命を持ち，主にバイオマス燃焼によって 生成される．本解析において HCN をバイオマス燃焼のプロキシとして用いた． 2010 年 4 月－2019 年 5 月の期間について，つくばの FTIR 観測スペクトルによる CO および HCN の高度分布の導出を行った．観測スペクトルの解析にはロジャーズの最適推 定法 [Rodgers, 2000] を基に開発されたインバージョン解析プログラム”SFIT4”を使用し た．また同期間について，CO，HCN 高度分布の相関解析，火災データ，粒跡線解析を用 いて，つくば上空の CO の発生源および発生地域の推定を行った．発生地域の推定には人 工衛星 Terra/Aqua 搭載の MODIS センサによる火災データと環境研 METEX による後方 粒跡線解析の結果を使用した． リトリーバルで得られた CO の高度分布ついて高度 0-5km, 5-18km の二層に分解し，同 じ日時，高度に対応する HCN との相関解析を行った．解析の結果，春季の CO と HCN のパーシャルカラムにおいて，高度 5-18km におけるそれぞれの値の間に有意な相関関係 が確認された. また高度 0-5km においては有意な相関関係が認められなかった．相関解析 から得られた CO/HCN 比から，高度 5-18km 春季について，大気中の CO パーシャルカ ラムにおけるバイオマス燃焼起源 CO の存在量の割合を求めた．高度 5-18km の 3-5 月に おいてバイオマス燃焼起源の CO の割合は平均で 70％程度であると推測する. この値は， つくば上空の大気に牧草や残渣の燃焼による CO が含まれている場合，過剰に見積もられ ている可能性がある．火災データと後方粒跡線解析の複合解析の結果, 有意な相関関係が 見られた春季の高度 5-18km において，中国北東部，ユーラシア大陸の緯度 50-60°N 帯 の地域の林野火災を起源とする CO である可能性が示された．

iii

目次

第 1 章 序論 1.1 一酸化炭素(CO) ... 1 1.2 CO の生成と消失 ... 1 1.3 CO の濃度分布および排出量の経年変化 ... 3 1.4 シアン化水素(HCN) ... 6 1.5 本研究の目的 ... 6 第 2 章 手法 2.1 フーリエ変換型赤外分光計 ... 7 2.1.1 特徴 ... 7 2.1.2 原理 ... 7 2.1.3 装置 ... 11 2.1.4 データセット ... 12 2.2 解析方法 ... 13 2.2.1 スペクトルの形状 ... 13 2.2.2 高度分布導出の概要 ... 15 2.2.3 高度分布の導出手法 ... 16 2.2.4 Averaging Kernel ... 19

2.2.5 Degree of Freedom for Signal ... 20

2.2.6 入力パラメータ ... 21 2.3 CO，HCN のスペクトルフィッティング ... 22 2.3.1 リトリーバルパラメータ ... 22 2.3.2 スペクトルフィッティングの結果 ... 22 2.3.3 つくば CO，HCN 高度分布の結果 ... 25 2.4 METEX による粒跡線解析 ... 27 2.4.1 粒跡線計算ツール METEX ... 27 2.4.2 粒跡線の計算 ... 27 2.5 MODIS による火災データ ... 28 2.5.1 MODIS センサ ... 28 2.5.2 火災強度 FRP ... 28 2.5.3 データセット ... 28

iii 第 3 章 つくば上空の CO 発生源の推定 3.1 解析する二層の高度範囲の決定 ... 30 3.1.1 解析手法の概要 ... 30 3.1.2 下層の高度範囲の決定 ... 30 3.1.3 上層の高度範囲の決定 ... 32 3.2 CO－HCN 相関解析 ... 32 3.2.1 上層 5－18km ... 32 3.2.2 下層 0－5km ... 33 3.2.3 Biomass-CO パーシャルカラム量の推察 ... 34 第 4 章 つくば上空の CO 発生地域の推定 4.1 つくばを起点とする粒跡線の特徴 ... 38 4.2 火災データの特徴 ... 41 4.3 粒跡線と火災データの複合解析 ... 44 第 5 章 結論 ... 48 謝辞 ... 49 付録 ... 50 参考文献 ... 64

1

第 1 章 序論

1.1 一酸化炭素(CO)

一酸化炭素(CO)は，対流圏オゾン(O3)，メタン(CH4)のような温室効果ガスの前駆物質 であることから，対流圏大気化学や気候変動において重要な役割を持つ物質である．CO は 主に化石燃料の不完全燃焼やバイオマス燃焼，大気中の炭化水素の酸化によって生成され， OH ラジカルとの反応によって消失する．CO の大気中での寿命は数週間から二ヶ月程度で あり，多様な発生源を持つことから全球的に濃度は不均一な物質である．

1.2 CO の生成と消失

CO の発生源は地域や季節によって多様であり，化石燃料の不完全燃焼やバイオマス燃焼， メタン等の炭化水素の酸化によって生成される．バイオマス燃焼とは，化石燃料以外の燃料 の燃焼を指し，動植物等を起源とする生物由来の物質を燃やすことを表す．具体的には，森 林火災，農業活動，残渣の燃焼等があり，自然起源と人為起源の両方を含んでいる．表 1.1.1 にグローバルな CO の収支を示す．対流圏発生源は，化石燃料の燃焼や熱帯での農業活動 (バイオマス燃焼)による人為発生源が大部分を占める．日本における人為起源 CO につい て，燃料施設の対策が行き届いていることから自動車による排出が主である．しかし，燃焼 効率の向上やガソリンエンジンと電気モーターとの併用により，地表大気中の CO 濃度は 減少傾向にある． 表 1.1.1 グローバルな CO の収支 [D. J. ジェイコブ，訳：近藤 (2002)] 推定値の範囲（Tg CO yr-1_） 発生源 1800 - 2700 化石燃料の燃焼/産業活動 300 - 550 バイオマス燃焼 300 - 700 植生 60 - 160 海洋 20 - 200 メタンの酸化 400 - 1000 他の炭化水素の酸化 200 - 600 消失源 2100 - 3000 対流圏 OH による酸化 1400 - 2600 成層圏 ～100 土壌への取り込み 250 - 640

2 図 1.1.1 日本における CO の濃度変化 [環境省，ホームページ] 一般局：特定の発生源の影響を受けない場所での測定 自排局：自動車排気ガスの影響を強く受ける場所での測定 CO の消失源は，主に対流圏での OH ラジカルとの酸化反応である．反応式は以下であ る．M は反応の触媒となる物質を示す．

𝐶𝑂 + 𝑂𝐻 → 𝐶𝑂

₂

+ 𝐻

(1.2.1)

𝐻 + 𝑂

₂

+ 𝑀 → 𝐻𝑂

₂

+ 𝑀

(1.2.2)

OH ラジカルは O3の光乖離によって生じる気体であり，以下の反応式によって生成される． hνは紫外線，O(1_{D)は励起状態の酸素原子を表す．}

𝑂

₃

+ ℎ𝜐 → 𝑂

₂

+ 𝑂(

1

𝐷)

(1.2.3)

𝑂(

1

𝐷) + 𝑀 → 𝑂 + 𝑀

(1.2.4)

𝑂(

1

𝐷) + 𝐻

₂

𝑂 → 2𝑂𝐻

(1.2.5)

上の反応式のように，O3の光乖離には紫外線が必須であり，OH の生成には水蒸気が必要 である．そのため OH ラジカルは，日射が強く，海洋からの水蒸気の供給量が増加する夏 季に多く生成される．結果として，夏季は対流圏 CO が減少し，冬季に増加する季節変化を 示す．また，この OH ラジカルとの反応を通して，同じく OH を主な消失源とする CH4の 大気中の寿命に影響を与えている[Gaubert et al., 2017]． OH ラジカルによる CO の消失(式(1.2.1)，式(1.2.2))の際に生成された HO2は，以下の反 応により OH を再生成し，NO2を生成する．

𝐻𝑂

₂

+ 𝑁𝑂 → 𝑂𝐻 + 𝑁𝑂

₂

(1.2.6)

さらに，NO2は光乖離によって NO と O3を再生成する．反応式は以下の式で表される．

𝑁𝑂

₂

+ ℎ𝜐 → 𝑁𝑂 + 𝑂

(1.2.7)

𝑂 + 𝑂

₂

+ 𝑀 → 𝑂

₃

+ 𝑀

(1.2.8)

式(1.2.1)，(1.2.2)，(1.2.6)，(1.2.7)，(1.2.8)による反応の結果を正味として以下の式で表す

3 ことができる．

正味：𝐶𝑂 + 2𝑂

₂

→ 𝐶𝑂

₂

+ 𝑂

₃

(1.2.9)

以上のように，CO は温室効果ガスであるメタンや対流圏オゾンの存在量に影響を与えてお り，大気化学や気候変動において重要な物質である．

1.3 CO の濃度分布および排出量の経年変化

図 1.3.1 に WDCGG (The World Data Centre for Greenhouse Gases)のデータをもとに緯 度帯別に平均した大気中の CO 月平均濃度の経年変化[気象庁]を示す．1900 年頃をピーク に，現在まで全球的な濃度の減少傾向が確認されている．さらに 1.1.2 で述べたように，冬 季から春季にかけて濃度が高く，夏季には濃度が低くなる季節変動がみられる．また，北半 球中高緯度で濃度が高く，南半球に向かうにつれて濃度が低くなっており，CO の主な放出 源が北半球中高緯度にあり，大気中への水蒸気の供給が多く OH ラジカルの量が多い熱帯 海洋上に運ばれるにつれて消滅することを示唆している． 図 1.3.1 緯度帯別に平均した大気中の一酸化炭素月平均濃度の経年変化 [気象庁 website] 図 1.3.2 に Worden et al. (2013)の衛星観測による地域別 CO の経年変化を示す．2000-2012 年における MOPITT の観測による気柱全量の 12 カ月移動平均の経年変化が地域別 に色分けされている．どの地域についても明瞭な減少傾向が見られる中，経済成長にある中

4

国(赤線)やインド(紫線)においても CO 濃度の減少傾向が見られている．産業活動による二 酸化炭素(CO2)の排出量が年々増加する両国において，燃料の燃焼効率の改善等の対策が

CO 濃度の減少に反映されているのが確認できる．

図 1.3.2 CO 気柱量経年変化の地域別トレンド [Worden et al., 2013]

Kurokawa and Ohara (2020) に よ る ア ジ ア 地 域 排 出 イ ン ベ ン ト リ REAS(Regional Emission inventory in ASia)では，アジア地域における燃料消費，自動車走行量，人口等の 統計データ，排出係数，規制動向などもとに，CO を含む大気汚染物質，温室効果ガスの発 生源別・燃料別・地域別の人為起源排出量の算出が可能である．図 1.3.3 (a)に REAS によ るアジアの地域別 CO 排出量の経年変化，(b)に地域分けの定義を示す．人為起源 CO につ いて，中国(水色)において 2005 年以降から排出量の減少傾向がみられ，Worden et al. (2013) による中国の CO 濃度の減少傾向と整合的な結果が得られていることが分かる．中国やイ ンドを含むアジア地域での CO の動向は，今後の全球的な CO の濃度変化の重要な要素の 一つである．

5

図 1.3.3 (a) アジア地域人為起源 CO 放出量の経年変化 (b) (a)の地域分け定義 [Kurokawa and Ohara, 2020]

SEA(South East Asia)：(b)の国名青色の地域 OEA(Other East Asia)：(b)の国名緑色の地域 OSA(Other South Asia)：(b)の国名茶色の地域

(a)

6

1.4 シアン化水素(HCN)

シアン化水素(HCN)とは，主にバイオマス燃焼を発生源とし，バイオマス燃焼のトレー サーとして用いられる物質である．Li et al. (2009) では，HCN の発生源とその放出量につ いて，バイオマス燃焼から 0.4 – 3.2 Tg N yr-1_{，家庭用バイオマス燃料の燃焼から 0.2 Tg N} yr-1_{と推定している．また，Bertschi et al. (2003) は HCN の発生源とその放出量について，} 世界の化石燃料の燃焼源からの放出量を 0.04 Tg N yr-1_{と推定しており，バイオマス燃料の} 燃焼による放出量と比較すると無視できるほど小さい．よって本研究では，つくば上空で観 測された HCN はすべてバイオマス燃焼起源であるとみなし，HCN をバイオマス燃焼起源 CO のプロキシとして用いる．

1.5 本研究の目的

アジア地域において，実際の CO 濃度の変化に対して，人為起源 CO 排出量の変化が与 える影響は確かに大きいが，アジア域の CO の動向をより正確に理解するためには，自然 起源の CO の影響を明らかにすることが望ましい．また，アジア域の風下に位置する日本 での CO の観測、解析はアジア域の CO の理解に大きく貢献できる．よって本研究では， CO の自然起源および人為起源の発生源であるバイオマス燃焼に着目し，つくばの FTIR で 観測された CO および HCN の高度分布を用いて，つくば上空のバイオマス燃焼起源 CO の 存在量を推定することを目的とする．さらに，粒跡線解析と火災マップを用いることで，つ くば上空の CO の具体的な発生方法と発生地域を推定することを目指す．

7

第 2 章 手法

2.1 フーリエ変換型赤外分光計

2.1.1 特徴

フーリエ変換型赤外分光（通称 FTIR：Fourier Transform Infrared Spectrometer）とは， 分子が赤外の波長領域に吸収を持つこと利用して試料を分析する分析装置である．マイケ ルソン干渉計を用いた分光計で，光の光路差により起こる干渉光のインターフェログラム を測定しフーリエ変換することで強度スペクトルを得る．一度の測定で広い波長領域のス ペクトルを測定できること，また，それを高い波数分解能で測定できることが FTIR 観測の 利点である．

2.1.2 原理

図 2.1.1 にマイケルソン干渉計の原理図を示す． 図 2.1.1 マイケルソン干渉計 光源からマイケルソン干渉計に強度𝐼₀の単色光が入射する場合を考える．干渉計の可動鏡を 動かしたときの 2 光束間の光路差を図 2.2.1 中の𝑑を用いて2𝑑 = 𝐷とし，この波長における 半透鏡の反射率，透過率を 𝑅₀，𝑇₀，単色光の波数を𝜐₀とするとき，観測される干渉波形は，

𝐼

_𝑜𝑏𝑠

= 2𝑅

₀

𝑇

₀

𝐼

₀

(1 + (cos 2𝜋𝜐

0

𝐷))

(2.1.1)

と表せる．ここで，半透鏡による吸収がなく，入射前と反射（透過後）の光エネルギーが保 存されるものと仮定する．その場合，最大効率は𝑅₀= 𝑇0= 0.5の時であるため，最大効率と 仮定して式(2.1.1)を変形すると

8

𝐼

_𝑜𝑏𝑠

=

𝐼

0

2

(1 + (cos 2𝜋𝜐

0

𝐷))

(2.1.2)

と書き換えることができる．以上は単色光の場合である．ここで，光源が連続光を出して居 る時の干渉波形を，光源からの光放射のスペクトル強度 B(ν)を用いて一般化すると

𝐼

_𝑜𝑏𝑠

(𝐷) = ∫

𝐵(𝜐)

2

(1 + (cos 2𝜋𝜐𝐷))𝑑𝜐

∞ 0

(2.1.3)

と表せる．また 𝐷 → ∞ としたとき，余弦関数は𝐷の微小な変化に対して±1 の間を激しく 振動し，平均すると 0 になる．すなわち

𝐼

_𝑜𝑏𝑠

(∞) = ∫

1

2

∞ 0

𝐼

_𝑜𝑏𝑠

(0)

(2.1.4)

この式(2.1.4)は，干渉出力波形（光源から出力された光のエネルギー）の直流成分（半透鏡 を透過し光源方向へと向かう光のエネルギー）を表す．このことは 𝐷 = 0 のとき，エネル ギーの半分が光源側へ戻ることからも理解できる．式(2.1.3)から式(2.1.4)を引くことで，干 渉光のインターフェログラムを得ることができる．インターフェログラム𝑃(𝐷)は，

𝑃(𝐷) = 𝐼

_𝑜𝑏𝑠

(𝐷) −

1

2

𝐼

𝑜𝑏𝑠

(0)

= ∫

𝐵(𝜐)

2

∞ 0

(cos 2𝜋𝜐𝐷)𝑑𝜐

(2.1.5)

と求めることができる．𝑃(𝐷)は𝐷に関する偶関数になっており，実際には負の波数は無いが 存在し得るものと仮定すると，𝐵(𝜐) = 𝐵(−𝜐)となる．よって式(2.1.5)は偶関数の公式より，

𝑃(𝐷) = ∫

𝐵(𝜐)

2

∞ −∞

exp(𝑖2𝜋𝜐𝐷) 𝑑𝜐

(2.1.6)

と書くことができる． 式(2.1.6)について，𝑃(𝐷)は𝐵(𝜐)をフーリエ変換したものであるから， 逆フーリエ変換を行い，スペクトル𝐵(𝜐)について求めると

𝐵(𝜐) = ∫

𝑃(𝜐)

2

∞ −∞

exp(−𝑖2𝜋𝜐𝐷) 𝑑𝜐

= ∫ 𝑃(𝐷)

∞ 0

cos(2𝜋𝜐𝐷) 𝑑𝐷

(2.1.7)

となる．実際の観測では，𝐷に関して連続的な点を測ることはできないが，光路差Δ𝐷ごとの 𝑃(𝐷) = 𝐼𝑜𝑏𝑠(𝐷) − 1 2𝐼𝑜𝑏𝑠(0)を測り，式(2.1.7)を計算することでスペクトルを得ることができ る．すなわち

9

𝐵(𝜐) = ∑

{𝐼

𝑜𝑏𝑠

(𝑛∆𝐷) −

1

2

𝐼

𝑜𝑏𝑑

(0)} cos(2𝜋𝜐𝑛Δ𝐷)Δ𝐷

∞ 𝑛=0

(2.1.8)

と表すことができる．以上が FTIR の原理である．しかし，式(2.1.8)では 𝐷 は ∞ ~ − ∞ まで定義されているが，実際の観測では光路差𝐷を無限大にとることはできない．実際の観 測で可動鏡の動かせる範囲は

𝐷 ≤ |𝐷

_𝑚𝑎𝑥

|

(2.1.9)

となる．したがって観測されるスペクトルは 式(2.1.8) で示された理想的な真のスペクト ルとは異なる．実測スペクトル𝐵𝑝(𝜐) は

𝐵

_𝑝

(𝜐) = ∫

𝑃(𝐷)

2

𝐷_𝑚𝑎𝑥 −𝐷_𝑚𝑎𝑥

exp(−𝑖2𝜋𝜐𝐷) 𝑑𝐷

= ∫

𝑃(𝐷)

2

∞ −∞

𝐹(𝐷) exp(−𝑖2𝜋𝜐𝐷) 𝑑𝐷

(2.1.10)

𝐹(𝐷) = {

1 |𝐷| ≤ 𝐷

𝑚𝑎𝑥

0 |𝐷| ≥ 𝐷

_𝑚𝑎𝑥

(2.1.11)

となる．インターフェログラム𝑃(𝐷)に方形関数 𝐹(𝐷)を掛けることで，𝑃(𝐷)について𝐷 ≤ |𝐷𝑚𝑎𝑥| の範囲を切り出している．上式のような表現によって逆フーリエ変換が可能になる ので，

𝑃(𝐷)

2

𝐹(𝐷) = ∫ 𝐵

𝑝 ∞ −∞

(𝜐)𝑒𝑥𝑝(𝑖2𝜋𝜐𝐷)𝑑𝜐

(2.1.12)

と書くことができる．また，フーリエ変換が方形関数𝐹(𝐷)になるような関数𝑓(𝜐)で表すと，

𝐹(𝐷) = ∫ 𝑓(𝜐)𝑒𝑥𝑝(𝑖2𝜋𝜐𝐷)𝑑𝜐

∞ −∞

(2.1.13)

となる．𝑃(𝐷)は，真のスペクトル𝐵𝑝(𝜐)を用いて式(2.1.6)のように表すことができる．式 (2.1.12)に式(2.1.6)と式(2.1.13)を代入すると

∫ 𝐵

_𝑝 ∞ −∞

(𝜐) exp(𝑖2𝜋𝜐𝐷) 𝑑𝜐 =

1

2

𝑃(𝐷)𝐹(𝐷)

=

1

2

∫ 𝐵(𝜐) exp(𝑖2𝜋𝜐𝐷) 𝑑𝜐 ∫ 𝑓(𝜐) exp(𝑖2𝜋𝜐𝐷) 𝑑𝜐

∞ −∞ ∞ −∞

=

1

2

∫ 𝐵(𝜐)

∞ −∞

∗ 𝑓(𝜐) exp(𝑖2𝜋𝜐𝐷) 𝑑𝜐

(2.1.14)

10 と表せる．＊はコンボリューション（重量積分）の定理を用いたことを表す．両辺は=で結 ばれているので積分記号を外すと

𝐵

_𝑝

(𝜐) =

1

2

𝐵(𝜐) ∗ 𝑓(𝜐)

(2.1.15)

となる．すなわち，実測スペクトル𝐵𝑝(𝜐)は真のスペクトル𝐵(𝜐)と関数𝑓(𝜐)とのコンボリュ ーションで表現される．𝑓(𝜐)をフーリエ変換したものが𝐹(𝐷)である．𝐹(𝐷)は 𝐷 = 0 の軸に 対して対称で2𝐷𝑚𝑎𝑥の全幅を持つ高さが 1 の矩形波であるから，𝑓(𝜐)は𝑠𝑖𝑛𝑐関数であり

𝑓(𝜐) = 2𝐷

_𝑚𝑎𝑥

sin 2𝜋𝜐𝐷

𝑚𝑎𝑥

2𝜋𝜐𝐷

_𝑚𝑎𝑥

= 2𝐷

𝑚𝑎𝑥

𝑠𝑖𝑛𝑐(2𝜐𝐷

𝑚𝑎𝑥

)

(2.1.16)

と書かれる．ただし

{

sin(𝑥) =

sin 𝜋𝑥

𝜋𝑥

𝑠𝑖𝑛𝑐(0) = 1, 𝑠𝑖𝑛𝑐(𝑛) = 0, 𝑛:

整数

∫ 𝑠𝑖𝑛𝑐(𝑥)𝑑𝑥

∞ −∞

(2.1.17)

と定義される．式(2.1.16)は，ある波数𝜐に主ピークを持つ曲線で，フーリエ変換型赤外分光 計の装置と呼ばれる．装置関数は分光器の分解能を表しており，式(2.1.16)の中央ピーク位 置と 𝑠𝑖𝑛𝑐(2𝜐𝐷𝑚𝑎𝑥) が最初に 0 になる位置までの距離を分光器の波数分解能として定義す る．この定義は𝑅𝑎𝑦𝑙𝑒𝑖𝑔ℎ基準といわれ，

2𝜐𝐷

_𝑚𝑎𝑥

= 1

(2.1.18)

を満たす 𝜐 が波数分解能 Δ𝜐 であるから

Δ𝜐 =

1

2𝐷

_𝑚𝑎𝑥

(2.1.19)

と定義できる．また，強度半値全幅を波数分解能とする𝑇𝑎𝑦𝑙𝑜𝑟基準で表す場合は，

Δ𝜐 =

0.607

𝐷

_𝑚𝑎𝑥

(2.1.20)

11

2.1.3 装置

本研究では，つくばの国立環境研究所内(緯度 36.1°N，経度 140.1°E)の大気微量成 分スペクトル観測室に設置される FTIR (Bruker IFS125HR)によって観測されたデータを 用いて解析を行った．Bruker IFS125HR では 2010 年 4 月～現在まで観測が行われてい る．表 2.1.1 に FTIR の諸元を示す． 図 2.1.2 に装置のブロック図を示す．観測装置は太陽追尾器，光学系，FTIR 本体，コン ピュータからなるり，光源に太陽光を用いている．屋上の太陽追尾器から平面鏡を通して 太陽光を室内にある FTIR 本体に導入する．FTIR 本体に入射する太陽光の波数領域を制 限するために，大気の窓領域に対応した光学フィルタを使用している．表 2.1.2 に光学フ ィルタの波数領域を示す．本研究では，フィルタ#2 とフィルタ#4 を用いている．検出器 には，中間赤外から近赤外域（4000～10000cm-1_{）に高い感度を有するインジウムアンチ} モン（InSb）検出器を用いており，熱雑音を抑えるために液体窒素で冷却している． 表．2.1.1 FTIR の諸元 図 2.1.2 装置のブロック図 [国立環境研究所特別研究報告書, 2003]

型名

最高波数分解能

最大光路差

観測波数範囲

IFS125HR

0.0035cm

-1

_257cm

_500-4350cm

-1

12 表 2.1.2 光学フィルタの波数領域 波数範囲（cm-1_） フィルタ#1 フィルタ#2 フィルタ#3 フィルタ#4 フィルタ#5 フィルタ#6 4000-4350 2800-3700 2400-3200 2000-2600 1700-2200 500-1350

2.1.4 データセット

本研究では，2010 年 4 月 6 日～2019 年 5 月 30 日の期間内につくばで観測された CO， HCN の高度分布を解析に用いた．表 2.1.3 に解析に用いたつくばにおける CO（フィルタ #4）の毎月の観測データ数，表 2.1.4 に HCN（フィルタ#2）の毎月の観測データ数を示 す．一日に同じフィルタを用いて複数回の観測が行われている日に関しては，日平均値を 解析に用いるため，一回の観測としてカウントしている． 表 2.1.3 観測データ数（CO） 月 1 2 3 4 5 6 7 8 9 10 11 12 計 2010 5 4 6 3 4 2 0 7 7 38 2011 8 8 4 3 3 1 1 5 2 1 6 3 45 2012 9 8 7 5 3 3 2 4 2 4 5 5 57 2013 2 1 5 3 6 0 4 5 4 2 8 1 41 2014 11 4 5 7 3 2 1 3 3 3 4 5 51 2015 6 0 6 4 3 1 3 1 1 5 2 5 37 2016 5 3 6 4 2 3 3 2 0 1 5 4 38 2017 6 9 6 6 4 7 5 2 5 6 5 10 71 2018 4 4 2 4 2 1 2 2 3 4 4 5 37 2019 8 3 5 5 5 計 59 40 46 46 35 24 24 28 22 26 46 45

13 表 2.1.4 観測データ数（HCN） 月 1 2 3 4 5 6 7 8 9 10 11 12 計 2010 4 4 7 5 0 3 0 7 7 37 2011 8 8 4 3 3 1 1 2 2 1 6 4 43 2012 10 8 7 5 3 3 2 4 2 4 5 5 58 2013 2 1 6 3 6 0 4 4 5 2 8 1 42 2014 11 4 5 7 3 2 0 1 3 3 4 5 48 2015 6 1 7 4 3 1 1 1 1 5 2 5 37 2016 5 5 8 7 3 4 3 1 0 2 6 7 51 2017 7 9 6 6 4 7 3 0 5 6 5 10 68 2018 4 4 2 4 2 2 0 3 3 4 5 5 38 2019 8 2 5 5 5 25 計 61 42 50 48 36 27 19 16 24 27 48 49

2.2 解析方法

2.2.1 スペクトルの形状

光のある波長を 𝜐 とし，大気上端での太陽光強度を 𝐼₀(𝜐) とすると，地上で観測される光 強度 𝐼(𝜐) は，

𝐼(𝜐) = 𝐼

₀

(𝜐) exp(−𝜏(𝜐))

(2.2.1)

と表せる．𝜏(𝜐) は大気の光学的厚みである．ここで簡単のために吸収する大気成分が 1 種 類だけであるし，平行平面大気を仮定すると，

𝜏(𝜐) = ∫ 𝜎(𝜐, 𝑧)𝜌(𝑧)

1

cos 𝜃

∞ 0

𝑑𝑧

(2.2.2)

と書ける．ここで 𝜎(𝑧, 𝜐)は吸収断面積，𝜎(𝑧)は高度𝑧における吸収成分の数密度，𝜃は太陽 天頂角を表す．式(2.2.1)に式(2.2.2)を代入すると，

𝐼(𝜐) = 𝐼

₀

(𝜐)𝑒𝑥𝑝 [− ∫ 𝜌(𝑧)𝜎(𝜐, 𝑧)

1

cos 𝜃

𝑑𝑧

∞ 0

]

(2.2.3)

となり，上記の式によって，地上で観測される吸収スペクトルが表現される．最終的には， 高度𝑧における吸収成分の数密度を求める．

14 吸収断面積𝜎(𝜐, 𝑧)は吸収線の強度𝑆 とその吸収線の波数方向への広がり𝐻によって決ま っており，

𝜎(𝜐, 𝑧) = 𝑆𝐻

(2.2.4)

となる．まず，吸収線の強度𝑆 について，

𝑆 = 𝑆

₀

[

ℎ𝑐

𝑘

𝐸" (

1

𝑇

₀

−

1

𝑇

)]

(2.2.5)

で表される．𝑆0, ℎ, 𝐸”, 𝑚はそれぞれ，𝑇0における線強度，プランク定数，lower stage energy，

吸収線強度の温度依存係数（吸収成分の分子構造の違いによる吸収の差を調節するための 係数，直線系分子の場合 𝑚 = 3 2⁄ ，対称性分子の場合 𝑚 = 1）を表す． 次に吸収線の波数方向への拡がり𝐻 について，大気成分が太陽光を吸収する場合，その 吸収はある固有の波数のみの光だけではなく，その光の中心波数の周囲にまで及ぶ．このよ うに，吸収線が中心波数の周囲に拡がりを持つ原因として以下の三つが上げられる． (1) 分子の励起状態の寿命による拡がり (2) 分子間の衝突や相互作用による拡がり：ローレンツ線形 (3) 分子の熱運動による拡がり：ドップラー線形 しかし大気中において，(1)の分子の励起状態の寿命は，分子間の衝突の時間よりはるかに 長い．(1)の拡がりは(2)の拡がりに対して限りなく小さいため通常は考慮しなくてよい．そ のため(2)のローレンツ線形および(3)のドップラー線形について考える． ローレンツ幅（ローレンツ線形の拡がり幅）𝛼𝐿は，

𝛼

_𝐿

(𝑃, 𝑇) = 𝛼

_𝐿

(𝑃

₀

, 𝑇

₀

)

𝑃

𝑃

₀

[

𝑇

₀

𝑇(𝑧)

]

𝑛

(2.2.6)

ドップラー幅（ドップラー線形の拡がり幅）𝛼𝐷は，

𝛼

_𝐷

(𝑇) =

𝜐

0

𝑐

[

2𝑘𝑇(𝑧)

𝑚

]

1 2

(2.2.7)

でそれぞれ与えられる．𝑃0, 𝑇0 はそれぞれ標準状態の気温(= 296𝐾), 気圧(= 1013.25ℎ𝑃𝑎) で，𝛼𝐿(𝑃0, 𝑇0)は標準状態におけるローレンツ幅，𝑧は高度である．𝜐0は吸収線の中心波数， 𝑐は光速，𝑘はボルツマン定数，𝑛はローレンツ幅の温度依存係数，𝑚は吸収成分の分子質量 である．式(2.2.6)，式(2.2.7)より，ローレンツ幅は気圧と気温に依存し，ドップラー幅は気 温に依存することが確認できる．大気成分による吸収で得られる吸収線の拡がりは，以上の

15 ローレンツ線形とドップラー線形の両方を重ね合わせた形状となり，これを近似した吸収 線の拡がりをフォイト線形という．地球大気の場合，成層圏以下の高度では，高度差による 気温の変化よりも気圧の変化が卓越するため気圧に依存するローレンツ幅が線形に大きく 影響し，逆に中間圏以上の高度では，気圧の変化よりも気温の変化が卓越するためドップラ ー幅による線形が支配的となる．このフォイト線形𝐻は，式(2.2.6)，式(2.2.7)の二つの吸収 線形の畳み込みで与えられる． 𝐻 (𝛼𝐿 𝛼𝐷 ,𝜐 − 𝜐0 𝛼𝐷 ) = 1 𝜋3⁄2 𝛼𝐿 𝛼𝐷 ∫ 1 (𝜐′_{− 𝜐} 0)2+ 𝛼𝐿2 𝑒𝑥𝑝 (−(𝜐 − 𝜐 ′₎2 𝛼𝐷2 ) ∞ −∞ 𝑑𝜐′ _(2.2.8)

𝑆0, 𝐸", 𝜐0, 𝑛は，TCCON(The Total Carbon Column Observing Network)による吸収線パラメ

ータ ATM16 [Toon et al. 2015]を使用している．

式(2.2.5)，式(2.2.8)より，式(2.2.4)の𝜎(𝜐, 𝑧)は気温と気圧に依存していることが分かる．よ って気温および気圧の高度分布が分かれば計算により𝜎(𝜐, 𝑧)を求めることができる．以上よ り，𝐼(𝜐)を観測できればインバージョンにより数密度𝜌(𝑧)を求めることができる．

2.2.2 高度分布導出の概要

スペクトルフィッティングによる高度分布導出の概要を図 2.2.1 に示す．2.2.1 で述べた ように，スペクトルの吸収線の線幅は気圧に依存しており，地上で観測されるスペクトル は，各気圧（高度）下に対応する吸収スペクトルの重ね合わせである．よって観測スペク トルの吸収線幅を細かく区切り，最適な波数幅から得られる情報を各高度に代表させるこ とで，より真値に近い高度分布を解として求めることができる．観測スペクトルから高度 分布を求めるためには，対象成分のアプリオリ（初期高度分布）を仮定する必要がある． そこで，仮定した成分のアプリオリが実際に分布する場合のスペクトル，すなわち計算ス ペクトルを，気象データ（気温，気圧），吸収線パラメータから放射伝達方程式より算出 し，計算スペクトルと観測スペクトルの差が最小となるようにインバージョン法を用いて フィッティングを行うことで，実際の高度分布を導入する．本研究では，地表（0.03km） から大気上端(106.00km)まで 48 層に分解し，層ごとの吸収量から成分の数密度を求めて いる．

16

図 2.2.1 高度分布導出の原理 [国立環境研究所特別研究報告書, 2003]

2.2.3 高度分布の導出手法

高度分布の導出には，SFIT4 というスペクトルフィッティングプログラムを用いて解析 を行った．SFIT4 は，NASA の Langley Research Center (LaRC)と National Institute of Water and Atmosphere (NIWA)による SFIT2 [e. g. Rinsland et al., 1998]を基に開発されて いる．リトリーバルの方法として，Optimal Estimation Method（通称 Rodgers 法

[Rodgers, 2000]）を用いたインバージョンにより高度分布を求めている． 以下，Rodgers 法の概要を述べる． 観測スペクトル𝑦と求めたい高度分布𝑥の関係をモデル化すると，次式のように書くことが できる．

𝑦 = 𝐹(𝑥, 𝑏)

(2.2.9)

𝐹はフォワードモデルと呼ばれる．𝑏は𝑥以外のすべてのフォワードモデルパラメータであ る．さらに，観測データには必ずノイズが含まれるので，それを誤差𝜀とすると式(2.2.9) は，

𝑦 = 𝐹(𝑥, 𝑏) + 𝜀

(2.2.10)

となる．𝑥は連続関数であるため，有限の波数分解能によって観測されるスペクトルから 真の𝑥を推定することはできない．そのため大気を有限の層数に分け，層ごとに近似され

17 た最適解𝑥̂を求めた．基本的な解法は最小二乗法である．フォワードモデルが𝑥について線 形可能な場合，

𝑦 = 𝐾𝑥 + 𝜀

(2.2.11)

𝑦 = [

𝑦

₁

𝑦

2

⋮

𝑦

𝑛

] ,

𝜀 = [

𝜀

₁

𝜀

₂

⋮

𝜀

_𝑛

] ,

𝑥 = [

𝑥

₁

𝑥

₂

⋮

𝑥

_𝑚

] ,

𝐾 = [

𝐾

11

𝐾

12

𝐾

₂₁

𝐾

₂₂

…

…

𝐾

_1𝑚

𝐾

_2𝑚

⋮ ⋮

⋱

⋮

𝐾

_𝑛1

𝐾

_𝑛2

⋯ 𝐾

_𝑛𝑚

]

(2.2.12)

とする．式(2.2.12)について，𝐾 =𝜕𝐹 𝜕𝑥であり，これは荷重積分関数(weighting function)で ある． 𝑛はある波数領域を観測装置の波数分解能で分解できるだけの数を表し，𝑚は有限 の数に分解した大気の層数となる． 誤差𝜀の二乗和が最小になる解を求めたいため．式(2.2.11)を𝜀について書き直すと

𝜀 = 𝑦 − 𝐾𝑥

(2.2.13)

となる．誤差の二乗和を𝑓とすると

𝑓 = 𝜀

₁2

+ 𝜀

₂2

+ ⋯ + 𝜀

_𝑛2

= 𝜀

𝑇

𝜀

= (𝑦 − 𝐾𝑥)

𝑇

(𝑦 − 𝐾𝑥)

}

(2.2.14)

と書ける．𝑇は転置行列を表している．𝑓が最小となるのは，𝑥で微分した式が０になると きであるので，式(2.2.13)について，𝜕𝑓 𝜕𝑥= 0とすると，

𝜕𝑓

𝜕𝑥

= −2𝑦

𝑇

_{𝐾 + 2𝑥}

𝑇

_(𝐾

𝑇

_{𝐾) = 0}

_(2.2.15)

𝑥

𝑇

(𝐾

𝑇

_{𝐾) = 𝑦}

𝑇

_𝐾

_(2.2.16)

と表すことができる．両辺の転置をとり，(𝐴𝐵)𝑇_{= 𝐴}𝑇_𝐵𝑇_{であり，(𝐴}𝑇_{𝐴)が対称行列である} ことから式(2.2.17)は

(𝐾

𝑇

_{𝐾)𝑥 = 𝐾}

𝑇

_𝑦

_(2.2.17)

となり，これを𝑥について解くと

𝑥̂ = (𝐾

𝑇

𝐾)

−1

𝐾

𝑇

𝑦

(2.2.18)

このようにして最小二乗解を求めることができる． 実際には観測誤差が観測ごとに含まれるので，最小二乗法に対して誤差の標準偏差の逆数 で重み付けをする．式(2.2.10)で重みを導入した式(2.2.14)を表すと，

𝑓 = [𝑦 − 𝐹(𝑥, 𝑏)]

𝑇

𝑆

_𝑒−1

[𝑦 − 𝐹(𝑥, 𝑏)]

(2.2.19)

となる．𝑆𝑒はノイズ共分散行列である．この𝑓が最小となるような解を求める． しかしこのままだと解が不安定になりやすい．例えば，極端に大きなノイズが乗ってしま った観測データによって解が発散してしまう場合などである．観測された𝑦には誤差𝜀が含

18 まれており，そのために求めたい𝑥が一意に決まらない場合がある．この問題を解決する ために，アプリオリを用いた拘束条件を最小二乗法によって導入し，解を収束させる．式 (2.2.19)を用いて表すと

𝑓 = [𝑦 − 𝐹(𝑥, 𝑏)]

𝑇

𝑆

𝑒−1

[𝑦 − 𝐹(𝑥, 𝑏)] + (𝑥 − 𝑥

𝑎

)

𝑇

𝑆

𝑎−1

(𝑥 − 𝑥

𝑎

)

(2.2.20)

𝑥𝑎はアプリオリであり，𝑆𝑎はその重み付けをするアプリオリの共分散行列である．式 (2.2.20)が最小になる解の推定法を Rodgers 法という．Rodgers 法を用いる場合重要なの が𝑆𝑒, 𝑆𝑎である．これらは各高度における重み係数を対角成分に収めた行列である．高度 毎に情報が独立している場合はそのまま対角行列となるが、通常では高度情報は完全に分 離できないため、非対角要素には相関係数（共分散係数）が入ることになり、共分散行列 となる。 SFIT4 を使用する際，𝑆𝑒が大きく（観測ノイズが小さく），𝑆𝑎の不確実性が大き い（アプリオリによる拘束が小さい）と式(2.2.20)の式は第一項が卓越し，𝑥 は観測値がよ り反映されたものとなる．逆に𝑆𝑒が小さい（観測ノイズが大きい）ため，𝑆𝑎の確実性を大 きくすると，𝑥 はアプリオリがより反映されたものとなる．第一項と第二項のバランスを とるための𝑆𝑒と𝑆𝑎の調節が必要である． 式(2.2.14)と同様に式(2.2.20)の線形近似が可能であるなら

𝑓 = [𝑦 − 𝐾

_𝑥

]

𝑇

_𝑆

𝑒−1

[𝑦 − 𝐾

𝑥

] + (𝑥 − 𝑥

𝑎

)

𝑇

𝑆

𝑎−1

(𝑥 − 𝑥

𝑎

)

(2.2.21)

となり．𝑥で微分した式から直接解が得られる．

𝑥̂ = (𝐾

𝑇

𝑆

_𝑒−1

𝐾 + 𝑆

_𝑎−1

)

−1

(𝐾

𝑇

𝑆

_𝑒−1

𝑦 + 𝑆

_𝑎−1

𝑥

_𝑎

)

= 𝑥

_𝑎

+ (𝐾

𝑇

𝑆

_𝑒−1

𝐾 + 𝑆

_𝑎−1

)

−1

𝐾

𝑇

𝑆

_𝑒−1

(𝑦 − 𝐾𝑥

_𝑎

)

}

(2.2.22)

𝜀 = 𝑦 − 𝐾𝑥なので， 𝐺 = (𝐾𝑇_𝑆 𝑒−1𝐾 + 𝑆𝑎−1) −1 𝐾𝑇_𝑆 𝑒−1とすると，式(2.2.22)は以下のように 書き直すことができる．

𝑥̂ = 𝑥

_𝑎

+ 𝐺(𝑦 − 𝐾𝑥)

= 𝑥

_𝑎

+ 𝐺[𝜀 + 𝐾(𝑥 − 𝑥

_𝑎

)]

}

(2.2.23)

つまり、𝐾で高度分布の次元からスペクトルの次元へ変更して、観測スペクトルデータと 重ね合わせ，𝐺で観測データと初期値の重み付けをすると同時に再び高度分布の次元へ戻 し、解を導出している．

19

2.2.4 Averaging Kernel

導出可能な高度やその分解能の指標としてAveraging Kernel (A)がある．図 2.2.2 にそ の例を示す．図中の青曲線の本数は計算した高度分布の数，つまり分解した大気の総数分 ひかれている．Averaging Kernel を確認することによって，各層の最適解高度分布が，ど の高度の情報から再現されたものなのかを可視化することができる．つまり，Averaging kernel を調べることで高度情報が取り出せる高度範囲と観測の実質的な高度分解能を知る ことができる．A は

𝐴 = 𝐺𝐾

(2.2.24)

で定義される．𝐴を用いて式(2.2.23)を表すと，

𝑥̂ = 𝑥

_𝑎

+ 𝐴(𝑥 − 𝑥

_𝑎

) + 𝐺𝜀

(2.2.25)

となる．この式から分かるように，最適解𝑥 は真値𝑥と初期値𝑥𝑎の重ね合わせであり，𝐴は それにかける係数である．最適解を得たい高度に対し，𝐴のピークがその高度に近く鋭い ほど，得られた最適解はその高度の真値を強く反映している．逆にピークがなだらかであ ったり、最適解を得たい高度からずれていたりすると、得られた最適解はほかの高度の影 響を強く受けていることになり、その高度の情報が取り出せていない。𝐴はどの高度に感 度があるかを表しており、𝐴の半値幅が概ね高度分解能に相当する． 式(2.2.25)を

𝑥̂ = 𝐴𝑥 + (𝐼 − 𝐴)𝑥

_𝑎

+ 𝐺𝜀

(2.2.26)

とする．𝐼は単位行列である．𝐴が単位行列に近いと第 2 項が 0 に近くなるので，最適解が 真の値により近いことがわかる．𝐴が 0 行列に近いと最適解は初期値により拘束され、初 期値の寄与が大きくなる．

20

図 2.2.2 2010 年 4 月～2019 年 5 月の全解析から得られた CO Averaging Kernels の例

2.2.5 Degree of Freedom for Signal

地表から大気上端までを高さとし，単位面積を底面とする気柱に含まれる分子量をトー タルカラム(気柱全量)と呼ぶのに対し，大気の任意の高度間に含まれる分子量をパーシャ ルカラムと呼ぶ．FTIR 観測による高度分布の導出では，大気を複数の層に分解し，層ご との濃度を計算するため，気柱全量に加えて任意の高度間のパーシャルカラムを求めるこ とができる．しかし，各層の最適解高度分布を得るために使われるスペクトル情報の採用 率は高度によって偏りがあるため，求めた高度分布は，高度によって観測スペクトルから 得られる情報量に差かある．その指標となるのが Degree of Freedom for Signal (DOFS)で ある．ある任意の高度間において，下端から上端までの DOFS の積分値が１を越えていれ ば，その高度は独立した一つの層として扱うことができるだけの情報が得られていると判 断できる．図 2.2.3 に例として CO の DOFS を地表から積分したプロファイルを示す．高 度 8.01km 以下，また高度 8.01-15.96km では，DOFs がそれぞれ 1.710，0.965 であり， 独立した層としての十分な情報を観測スペクトルから得られていることが確認できる．こ のことは，図 2.2.2 のアベレージングカネルが地表から高度約 20km の間にピークを持つ ことからも理解できる．もし地表から大気上端までの DOFs の積分が 1 程度にしか満たな い場合，得られた高度分布は，高度を分解しパーシャルカラムを得るのに必要な情報量を 含んでいないことを意味する．

21

図 2.2.3 DOFS プロファイル（2010 年 4 月～2019 年 5 月，つくば，CO）

2.2.6 入力パラメータ

解析には，入力パラメータとして，気温気圧プロファイル，吸収線パラメータ，対象成 分の初期高度プロファイルが必要である．気温気圧プロファイルには、アメリカ国立環境 予測センター（National Centers for Environmental Prediction）による NCEP 再解析デー タを用いた．計算スペクトルの作成に必要な吸収線パラメータには，TCCON(The Total Carbon Column Observing Network)による ATM16 [Toon et al. 2015]を使用している． ATM16 により，吸収線の中心波数𝑣0，標準状態におけるローレンツ幅𝑎𝐿0，ローレンツ幅

の温度依存係数𝑛，標準状態における線強度𝑆0，Lower state energy E”といったパラメータ

を全て得ることが出来る．各成分の初期高度プロファイルには，気象モデル WACCM （Whole Atmosphere Community Climate Model）によるつくばにおける 40 年平均値を用 いた．WACCM の開発は，NCAR Community Earth System Model （CESM）を共通の ツールとして使用して，上層大気モデリング，中層大気モデリング，対流圏モデリングを 統合することで構成されている[Marsh et al ., 2013]．

22

2.3 CO，HCN のスペクトルフィッティング

2.3.1 リトリーバルパラメータ

表 2.3.1 に CO および HCN 高度分布のリトリーバルパラメータを示す．より正確な高 度分布を導出するため，RMS（Root Mean Square residual：残差の二乗平均平方根）と DOFS の値を指標とし，より真値に近い高度分布か導出可能な観測データを使用してい る．CO の観測データについて，全観測データ数 755 のうち，RMS 値が 1.0 を超えるも の，DOFS 値が 2.5 を下回るもの計 18 の観測データを解析から除外した．また HCN の観 測データについて，全観測データ数 801 のうち，RMS 値が 0.4 を超えるもの，DOFS 値 が 1.5 を下回るもの計 123 の観測データを解析から除外した． 表 2.3.1 CO，HCN のリトリーバルパラメータ Species CO HCN Algorithm SFIT4

A priori profiles WACCM v6

Micro-Windows(cm-1₎

MW1: 2657.7-2058.0 MW1: 3268.06-3268.4 MW2: 2069.56-2069.76 MW2: 3287.1-3287.3

MW3: 2157.5-2159.15 MW3: 3299.4-3299.595 Retrieved species CO2, OCS, H2O, N2O H2O, C2H2, CO2, O3

Retrieval constraint RMS ≤ 1.0 RMS ≤ 0.4 DOFS ≥ 2.5 DOFS ≥ 1.5

2.3.2 スペクトルフィッティングの結果

図 2.3.1 に CO の解析を行う各 3 つの波数領域（MW：Micro-Window)における，CO2, OCS, H2O, N2O の吸収線と，スペクトルフィッティングとその残差の結果を示す．また， 図 2.3.2 に HCN の解析を行う各 3 つの波数領域における，H2O, C2H2, CO2, O3の吸収線 と，スペクトルフィッティングの結果を示す．各図中，下部のグラフがスペクトルフィッテ ィングの結果と各リトリーバル成分の吸収線を示しており，青線が観測スペクトル，オレン ジ線がフィッティングを行った計算スペクトルの結果を表している．各図中，上部のグラフ が観測スペクトルと計算スペクトルのフィッティング残差の結果を表し，値が小さいほど フィッティングが良い結果で，単位は％である．フィッティング残差の二乗平均が RMS と なる．

23

図 3.2.1 CO のフィッティングスペクトル結果

24

図 3.2.2 HCN のフィッティングスペクトル結果

25

2.3.3 つくば CO，HCN 高度分布の結果

図 2.3.3 に導出したつくばにおける CO および HCN 高度分布の平均(2010.4.6～ 2019.5.30)を示す．黒線がリトリーバルによって得られた高度分布の平均を示し，赤線が初 期値として用いた WACCM の高度分布の平均を示している．つくば CO の高度分布につい て，地表付近のみに濃度のピークを持つ初期高度分布とは異なり，地表付近と高度 8km 付 近の二か所に濃度のピークを持つ高度分布が得られた．HCN の高度分布について，初期高 度分布は高度 10km 付近に濃度のピークが見られていたが，導出された高度分布では，地表 付近にピークが確認される結果となった． 図 2.3.4 に CO および HCN 高度分布の平均(2010.4.6～2019.5.30)のアベレージングカネ ルを示す．CO のアベレージングカネルについて，高度分布で確認された地表付近，高度 8km 付近を含む高度 0-18km にピークが見られ，高度 18km 付近までは解析を行うための 十分な情報が観測スペクトルから得られていることが確認された．HCN のアベレージング カネルについても高度 0-20km にピークが見られ，CO の解析に対応する高度において十分 真値を反映した高度分布が得られていることが確認された． 図 2.3.3 CO，HCN 高度分布の平均 (2010.4.6～2019.5.30) (a)CO, (b)HCN

26 図 2.3.4 2010.4.6～2019.5.30 の全解析で得られた高度分布のアベレージングカネル (a) CO, (b)HCN

(a)

(a)

(a)

(b)

27

2.4 METEX による粒跡線解析

2.4.1 粒跡線計算ツール METEX

大気成分の移流経路を解析する手段として，国立環境研究所地球環境研究センター (NIES-CGER : National Institute for Environmental Studies – Global Environmental Research)で開発された粒跡線計算ツール METEX (Meteorogical Date Explorer)[Zeng et al., 2010]を使用した．METEX では，大気の塊を質点としてみなし，気象データからその 質点の移動軌跡（粒跡線）を求めることができる．本研究ではアメリカ環境予測センター （NCEP）の再解析データを気象データとして用いて計算された粒跡線データを用いた．

2.4.2 粒跡線の計算

対象の期間内における FTIR によるすべての CO 観測日時に対して，つくばの各高度を起 点とする後方粒跡線の計算を用いて行った．表 2.4.1 に本研究における粒跡線解析の諸元 を示す．実際には表中の起点高度に加えて高度 10，12，14，16km を起点とする粒跡線の 計算もおこなったが，大気がジェット気流に乗り粒跡線が地表付近に到達しなかったこと から，地表に発生源を持つ CO の起源地域の推定には適さないと判断し，解析から除い た． 表 2.4.1 粒跡線の諸元

期間

2010. 04. 06 ～ 2019. 05. 30

長さ

120 hours（5 days）Backward

起点経緯

36.05°N，140.12°E

起点高度

0 , 0.5 , 1, 2 , 4 , 6 , 8 km

起点日時

CO のスペクトルが観測された日時

範囲

10 ~ 80°N，0 ~150°E

28

2.5 MODIS による火災データ

2.5.1 MODIS センサ

本研究では，NASA（National Aeronautics and Space Administration）の FIRMS（Fire Information for Resource Management System）が提供している Terra/Aqua 衛星（the MODerate Resolution Imaging Spectroradiometer）センサにより観測された火災データの緯度 経度，日時，火災強度を用いた．実際に使用したデータは NASA，FIRMS のウェブサイトから ダウンロードしたものである．MODIS センサは Terra，Aqua 両衛星に搭載され，同一地点を 昼と夜一日に 1～2 回観測する．火災データの水平分解能は 1km×1km であり，データの期間 は 2000 年から現在までである．

2.5.2 火災強度 FRP

火災強度 FRP（Fire Radiative Power）とは，衛星から観測される林野火災の強度を表す物理 量であり，火災からの放射量（[J/s]または[MW]）で表せる．FRP が高い火災では，放出され た煙が高高度に達し，境界層を越えて自由対流圏まで輸送されることがあり，FRP と放出され る煙の到達高度の関連性が報告されている[Martin et al., 2010]． FRP は以下の式で概算される．

𝐹𝑅𝑃 ≈

𝐴

𝑝𝑖𝑥

𝜎

𝑎𝜏

₄

(𝐿

4

− Γ

4

)

𝐿4は 4μm 帯の放射輝度，Γ4は 4μm のバックグラウンドの放射輝度，𝐴𝑝𝑖𝑥は MODIS のピクセ ルの面積，𝜎はシュテファン-ボルツマン定数（5.6704 × 10−8_𝑊𝑚−2_𝐾−4_），𝜏 4は 4μm 帯の大気 透過率，𝑎はセンサー固有の経験的定数（3.0 × 10–9 _{W m}–2 _sr–1 _μm–1 _K–4_{）を表している[Giglio et} al., 2016]．本研究では，FRP を指標とし，発生した成分が長距離の輸送が可能な火災データの みを解析に用いている．

2.5.3 データセット

本研究では，緯度 10 - 80°N，経度 0 - 150°E の範囲内で発生した火災のデータを用いて いる．その範囲を図 2.5.1 に示す．火災データを用いてつくば上空の CO の発生源を推定し たいため，火災によって発生した煙が長距離輸送の可能な大気の高度まで到達したデータ を用いたい．図 2.5.2 に煙の高さと FRP の関係を示す．Martin et al. (2010) では，MODIS による FRP と，火災による煙の高さの相関関係を報告しており，FRP と煙の高さの関係が

29 回帰直線として得られている．本研究ではこの結果を利用し，火災による煙の高さが大気の 境界層（地表から高度 1000km）を越えるような FRP の火災データを扱うものとする．計 算の結果，煙の高さが 1000m に相当する FRP の大きさは 320MW であった．よって本研 究では，FRP が 320MW 以上の火災データのみを使用した． 図 2.4.1 火災データの解析範囲

30

第 3 章 つくば上空の CO 発生源の推定

3.1 解析する二層の高度範囲の決定

3.1.1 解析手法の概要

本研究では，つくばにおける FTIR 観測で得られた CO および HCN の高度分布を用い て，つくば上空の CO の発生源とその存在量を推定することを目指す．ユーラシア大陸の 風下に位置する日本で観測される大気成分の気柱全量は，観測地周辺の国内に発生源を持 つものだけではなく，アジア州等の日本から距離のある国外に発生源を持つ成分が含まれ ていることが期待できる．長距離輸送によってつくばに到達する CO および HCN の特徴 は，地表付近の大気よりも，より高い高度の存在量に反映されやすいと考えられる．よって 本研究では，解析を行う高度範囲を上層と下層の二層に分解し，FTIR 観測による高度分解 性能を活かし，つくば上空の CO および HCN の高度分布から，本研究で定義する地表付近 および上空の２つの高度領域を抜き出し，二層のパーシャルカラム量の相関解析の結果を 比較することでつくば上空の CO 発生源の推定をおこなった．

3.1.2 下層の高度範囲の決定

解析する二層の高度範囲について説明する．まず下層となる高度範囲の決定について，本 研究で求めた CO 高度分布の分解能に基づいて決定した．図 3.1.2 に本研究で計算した CO および HCN 高度分布の DOFS それぞれを示す．第 2 章，図 2.3.4 中の CO，HCN の Averaging Kernels はともに地表から約 20km にピークを持つことから，CO，HCN の高度 分布は，地表―20km の高度範囲において真値を反映できていることが確認された．計算さ れた CO 高度分布を複数に分解して解析を行う際，分解されたそれぞれの高度範囲がリト リーバルの際に観測スペクトルから十分な情報を得られている必要があるため，二層それ ぞれの高度範囲において DOFS の値が１を上回ることが望ましい．図 3.1.1 中の高度 5km において，CO の DOFS 値は 1.220 となり 1 を上回る，また高度 5ｋｍはＣＯ高度分布の地 表付近のピークを包含できるため，本研究では，下層の高度範囲を高度 0－5km と定義す る．

31

図 3.1.1 2010 年 4 月～2019 年 5 月の全解析で得られた DOFS プロファイル (a)CO, (b)HCN

(a)

32

3.1.3 上層の高度範囲の決定

次に上層となる高度範囲の決定について，上層の高度範囲は本研究で求めた CO 高度分 布の線形の特徴から決定した．第 2 章，図 2.3.3 に本研究で求めた CO および HCN 高度分 布の平均を示している．黒線および赤線はそれぞれリトリーバルで得られた高度分布とア プリオリの平均を表している．本研究の対象期間において CO，HCN の FTIR 観測によっ てそれぞれ，755（うち 18 を解析から除外），801（うち 123 を解析から除外）観測データ が得られており，黒線は期間内すべての観測スペクトルから計算された高度分布の平均を 示している．リトリーバルで得られた CO 高度分布について，地表面と上空およそ 8km の 二か所に混合比のピークが確認された．上層の高度範囲には，この二つの混合比のピークの うち，上空のピークの特性を解析に反映させたい．このピークによるカーブは高度 4km と 18km に極小を持つ．よって上層の高度範囲を 4－18km とするのが理想であるが，下層と の領域の重複を避け，本研究では上層の高度範囲を高度 5－18km と定義する．下層の高度 範囲において，HCN のアベレージングカネルの DOFS 値は 0.579 であり 1 を上回ってい ないが，ＣＯ高度分布における二か所の混合比のピークの特性を各層の解析に反映させる ために，本研究では高度 0－5ｋｍを下層の高度範囲として解析を行った．

3.2 CO－HCN 相関解析

3.2.1 上層 5－18km

つくば上空の大気に含まれる CO とバイオマス燃焼の関係性を明らかにするために，高 度 5－18km および高度 0－5km において，CO および HCN パーシャルカラム量の相関解 析を行った．始めに高度 5－18km の結果について，図 3.2.1 に，上層における季節別 CO， HCN の相関を示す．3,4,5 月を春季，6,7,8 月を夏季，9,10,11 月を秋季，12,1,2 月を冬季と 季節分けを行い，CO および HCN について解析期間内すべての観測日のパーシャルカラム 量日平均値をプロットしている．回帰分析の結果，春季，夏季，秋季，冬季について，決定 係数（R2_{）はそれぞれ 0.59，0.40，0.54，0.38 となり，P 値はそれぞれ 1.24×10}-25_{, 2.86×} 10-7_{, 3.05×10}-17_{, 5.53×10}-18_{となった．春季と秋季において，夏季と冬季より比較的大きな} 決定係数が確認された．表 3.2.1 に CO，HCN 相関解析における各月の回帰分析の決定係 数を示す．各月も良い相関関係を示す中，3,4,5 月は一貫して高い決定係数を示し，春季で の高い決定係数と整合的な結果を示した．一方で 10,11 月は良い相関関係が確認されず，秋 季の高い決定係数とは不調和な結果となった．また 7,8 月は比較的良い相関関係を示した が，夏季での比較的低い決定係数とは不整合な結果となった．

33 図 3.2.1 高度 5-18km，季節別 CO-HCN パーシャルカラム相関 (a)3-5 月，(b)6-8 月，(c)9-11 月，(d)12-2 月 表 3.2.1 高度 5-18km，CO，HCN パーシャルカラム相関の月別決定係数

3.2.2 下層 0－5km

次に高度 0－5km の結果について，図 3.2.2 に下層における季節別 CO-HCN の相関を示 す．上層と同じく 3,4,5 月を春季，6,7,8 月を夏季，9,10,11 月を秋季，12,1,2 月を冬季と季 節分けを行い，CO および HCN について解析期間内すべての観測日のパーシャルカラム量 日平均値をプロットしている．回帰分析の結果，春季，夏季，秋季，冬季について，決定係 数はそれぞれ 0.13，0.40，0.10，0.20 となり，P 値はそれぞれ 3.30×10-5_{, 3.09×10}-7_{, 1.60} ×10-3_{, 1.80×10}-8_{となった．春季，秋季，冬季において，上層とは対象的に，CO，HCN パ} ーシャルカラムの間に有意な相関関係は確認されず，P 値においても上層より大きな値を示 した．また夏季の決定係数，P 値は上層と同程度の値を示した．表 3.2.2 に CO-HCN 相関 解析における各月の回帰分析の決定係数を示す．すべての月において，上層よりも小さい決 定係数を示した．特に，春季の月において上層の決定係数を大きく下回る値が確認された．

34 上層と下層の CO，HCN 相関の回帰分析の結果，春季において，上層では一貫して有意 な相関がみられ，下層では一貫して有意な相関が確認されなかった．この結果は，春季のつ くば上空において上層と下層で CO の発生源が明確に異なることを示す． 図 3.2.2 高度 0-5km，季節別 CO-HCN パーシャルカラム相関 (a)3-5 月，(b)6-8 月，(c)9-11 月，(d)12-2 月 表 3.2.2 高度 0-5km，CO，HCN パーシャルカラム相関の月別決定係数

3.2.3 Biomass-CO パーシャルカラム量の推察

第 1 章の 1.1.4 で述べた通り，HCN の発生源の大半はバイオマス燃焼である．よって本 研究では，つくば上空の HCN はすべてバイオマス燃焼起源であると仮定し，HCN をバイ オマス燃焼のプロキシとして解析に用いた．回帰分析において有意な直線的相関が認めら れる場合，回帰直線の傾きを二成分の存在量比としてとらえることができる．春季のつくば 上空の高度 5-18km において CO，HCN に有意な相関が確認されたことから，本研究では 回帰直線の HCN 増加量に対する CO 増加量を，HCN に対するバイオマス燃焼起源 CO の 存在量比( CO/HCN )Biomassとしてとらえ，以下の式から 3,4,5 月および春季におけるバイオ マス燃焼起源 CO，非バイオマス燃焼起源 CO を算出した．

35

[Biomass-CO] = [Biomass-HCN] × ( CO/HCN )Biomass

[Others-CO] = [All-CO] － [Biomass-CO]

[Biomass-CO]: 高度 5-18km におけるバイオマス燃焼起源 CO パーシャルカラム [Biomass-HCN]: 高度 5-18km における HCN パーシャルカラム

( CO/HCN )Biomass: 高度 5-18km における CO，HCN 相関の回帰直線の傾き

[Others-CO]: 高度 5-18km における非バイオマス燃焼起源 CO パーシャルカラム [All-CO]: 高度 5-18km における CO パーシャルカラム 高度 5-18km における 3,4,5 月および春季の CO，HCN 相関を図 3.2.3 に示す．3 月，4 月，5 月，春季における(CO/HCN)Biomassの値はそれぞれ 202.9，225.0，260.7，206.4 とな った．また表 3.2.3 に，これらの結果と Akagi et al. (2011)における植生別森林火災による CO，HCH 放出量比の比較を示す．高度 5-18km つくばにおける 3,4,5 月および春季の (CO/HCN)Biomass値は，熱帯林火災による CO, HCN 発生量比と近い値を示した．またアジ

ア州の多くの国が属する亜寒帯林および温帯林での森林火災による CO，HCN 発生量比と は整合的な結果が得られなかった．

36

表 3.2.3 植生別森林火災による CO/HCH 放出量比の比較

表 3.2.4 および図 3.2.4 に算出したバイオマス燃焼起源 CO および非バイオマス燃焼起源 CO のパーシャルカラム量の結果を示す．(a)はすべての月の( CO/HCN )Biomassとして，す

べ て の 月 に お い て 春 季 に お け る 回 帰 直 線 の 傾 き （ =206.4 ） の 値 を 使 用 し ， (b) は ( CO/HCN )Biomassとして，各月の回帰直線の傾き(Mar=202.9，Apr=225.0，May=260.7，

Spr=206.4)の値を使用している．春季つくば高度 5-18km におけるバイオマス燃焼起源 CO について，(a)の存在量および割合は各月で 4.60～6.29×1017 _[molec/m2_{]，65～75 % の変}

動幅を持ち，(b)の存在量および割合は各月で 4.52～7.94×1017 _[molec/m2_{]，64～95% の}

変動幅が確認された．各月の( CO/HCN )Biomass 値を使用した(b)の方が，春季通算の

( CO/HCN )Biomass 値を用いた(a)よりも，4,5 月において高いバイオマス燃焼起源 CO 存在

量を示した．

また，表 3.2.3 より，牧草や残渣の燃焼による CO/HCH 放出量比は，他の植生の燃焼よ りも大きくなっている．よってつくば上空の大気に牧草や残渣の燃焼による CO が含まれ ている場合，本研究で求めたバイオマス燃焼起源 CO の割合の値は，過剰に見積もられて いる可能性がある．

37

表 3.2.4 高度 5-18km，バイオマス燃焼起源 CO のパーシャルカラム量

図 3.2.4 高度 5-18km，バイオマス燃焼起源 CO のパーシャルカラム量 (a)すべての月の( CO/HCN )Biomass に，春季の回帰直線の傾きの値を使用

38

第 4 章 つくば上空の CO 発生地域の推定

4.1 つくばを起点とする粒跡線の特徴

本研究では，つくば(緯度 36.1°N，経度 140.1°E)において CO の FTIR 観測が行われ た日時に対して，つくばの各高度(0, 0.5, 1, 2, 4, 6, 8 km)を起点とし，120 時間（5 日間）遡 る粒跡線を計算した．ここでは，計算した粒跡線を地図上にプロットし，その画像から確認 された粒跡線の季節別，高度別特徴を述べる． つくばを起点とする粒跡線は，基本的のその経路をユーラシア大陸上に持つことが多く， 春季，秋季，冬季の粒跡線には共通して経路に一貫性が見られた．例として図 4.1.1 に解析 期間内における，つくばを起点とする高度 0km, 2km, 8km の 4 月すべての粒跡線を示す（他 の月の粒跡線は付録に示す）．高度 0km を起点とする粒跡線はつくばから北西方向に伸び， その終点は中国北東部やロシア南部付近の地域に集中している．ここから，起点の高度が上 がるにつれて，経路方向を北西から西方向に移し，終点までの距離が長くなる様子が確認さ れた． 夏季（6, 7, 8, 9 月）の粒跡線は，他の季節とは異なる特徴を示した．例として図 4.1.2 に つくばを起点とする高度 0km, 2km, 8km の 8 月すべての粒跡線を示す（他の月の粒跡線は 付録に示す）．夏季の粒跡線は，経路の方向に規則性がなく，様々な方向から，大きな蛇行 を繰り返すようにつくばに到達する様子が確認された．また夏季においても，起点の高度が 大きくなるにつれて，終点までの距離が長くなる様子が確認された．

39 図 4.1.1 高度別つくばを起点とする 4 月の粒跡線 (2010~2019) 上から高度 0km，2km，8km を起点とする粒跡線．色は年の違いを表す．

Altitude 0m

Altitude 2000m

Altitude 8000m

40 図 4.1.2 高度別つくばを起点とする 8 月の粒跡線 (2010~2019) 上から高度 0km，2km，18km を起点とする粒跡線．色は年の違いを表す．

Altitude 0m

Altitude 2000m

Altitude 8000m

41

4.2 火災データの特徴

表 4.2.1 に，2010 年 4 月 6 月～2019 年 5 月 30 日の期間において，解析範囲内で発生し た火災の年，月別の発生件数を示す．本研究で用いるデータはすべて FRP が 320MW 以上 の火災である．月別の火災発生数には大きな違いが見られた．春季の火災発生数が最も多く， 3, 4, 5 月の合計で全体件数の 48%を占めた．次に夏季の発生数が多く，6, 7, 8 月の合計で 全体件数の 38%であった．秋季（9, 10, 11 月）と冬季（12, 1, 2 月）の発生数はそれぞれ全 体の 9%, 5%となり，本研究の範囲地域（緯度 10 - 80°N，経度 0 - 150°E）では，火災の 大半が春季と夏季に発生しているという結果が得られた． 表 4.2.1 火災の年，月別の発生件数 (2010.4.6~2019.5.30, 10-80°N, 0-150°E ) 図 4.2.1 に各季節の火災発生地域の例として 4, 7, 10, 1 月の火災マップを示す（他の月は 付録で示す）発生数が最も多かった春季では．中国の北東部，ユーラシア大陸の緯度 50-60° N 帯，東南アジアに火災が集中していた．夏季においては，ロシア東部，中央アジアで多く の火災が確認された．秋季は中央アジア，アフリカの緯度 10°N 付近で火災が見られ，冬 季は同じくアフリカで火災が見られたがアジア域での火災はほとんど確認されなかった．

本研究で用いた FRP が 320MW 以上の火災データは，MODIS Fire Location Product (MCD14ML)によって，ほぼすべての火災が植生火災であると推定された．

42

April

July

43

図 4.2.1 月別火災マップ (2010.4.6~2019.5.30)

上から 4 月，7 月，10 月，1 月の火災．色は年の違いを表す．