増田 耕一、伊藤 孝士**、楳田 貴郁**
東京大学 理学部 地球物理学教室、**学生
Kooiti MASUDA, Takashi ITO and Takafumi UMEDA
Geophysical Institute, Faculty of Science, University of Tokyo
地表面熱収支の各項の分布は等値線図による地図集の形で出版されている。そこではさまざまな投影図法が使われているため、相互に比較したり、統合して利用することが困難だった。本研究では、いくつかの地図集から、等値線の情報を一つの地図の上に合わせて表示するための試験的システムを作った。実用的利用までにはまだ考慮しなければならない問題点が多数ある。
地表面の熱収支(エネルギー収支)をBudyko (1971)に従って定式化すれば、次のようになる。
R = LE + P + A (1)(図1a参照)。ここで、各項は面積あたり時間あたりのエネルギー流量(エネルギーフラックス密度)であり、Rは正味放射(太陽放射・地球放射を合わせ正味下向き)、Pは顕熱輸送(熱伝導と対流による地表面から大気への熱輸送)、Aは地表面からより深い層(土壌・海水)への熱輸送、LE(L: 水の質量あたりの気化熱、E: 蒸発量)は水の蒸発による潜熱輸送である。その他の項(雪氷の融解熱など)は多くの場合値が小さいので無視されている。地表面の水収支は、図1bのように土壌の層を考えて
r = E + ΔW + f (2)と書くのがふつうである。rは降水量、Eは蒸発量、ΔWは土壌水分貯留量の時間変化、fは流出量である。蒸発量Eを通じて、熱収支と水収支は深く結びついている。
人間活動の影響を取り除いた自然植生を考える上で、第１に考えるべきものは気候条件であろうが、その場合の説明変数として、地表面熱収支・水収支の各項をとることが適当であることが、いくつかの研究によって示唆されている。Budyko(1971)は正味放射およびそれと降水量との比を植生帯と、Uchijima・Seino(文献省略)はこれらのパラメタと自然植生の光合成生産量を関連づけた。Lieth, Boxら(文献省略)は、蒸発量と植生帯および光合成生産量を関連づけた。
ここではまず、熱収支に的をしぼり、これまでに出版されている分布図の情報を統合して活用することを目標として作業を開始した。
われわれは地球上の任意の位置での熱収支項の値、あるいは任意の面積での集計量を知りたい。しかし地表面熱収支の各項を直接的な方法で観測している地点は少ない。(筑波大学水理実験センターはその一つである。) したがって、熱収支項の分布は、通常の気象観測地点で観測されている、気温、地表面温度(海上)、湿度、風速、雲量、降水量(陸上)などのデータを活用して推定しなければならない。熱収支気候図はいずれもこのような方法で作られている。
ここで使われる方法の概略を図2aに示す。{1}{2}{3}が入力、{7}が出力である。少数の特別観測点では、比較的直接的な方法で、求めたい物理量の値が観測されている{1}。これと並行して、特別観測点・通常観測点の双方で、通常の気象変数が観測される{2}{3}。次に、理論または経験則を用いて、通常観測量から求めたい物理量への変換式の形を構築する{4}。この変換式のパラメタは、{1}{2}の特別観測の結果を用いて決定される{5}。{3}から入力される通常観測点のデータを{4}{5}で得られた半経験的変換式に当てはめれば{6}、任意の通常観測点での値が得られる{7}。なお、パラメタの値が空間的に変わると考えられる場合には、特別観測点のそれぞれで得られたパラメタの値を空間的に内挿してから通常観測点で式を適用する。
熱収支気候図は、本来地球上の面的分布を表現するものであるが、伝統的に等値線(コンター)による表現がとられてきた。等値線とは値が一定の位置を連ねた線である。等値線による面的データの表現は、ポリゴン形式の一種と考えることもできる。2つの等値線にはさまれた領域に、いくつ以上いくつ未満という「値」が対応しているからである。なお、熱収支の場合、陸上と海上で不連続があるため、ポリゴンを形成するのは等値線のほか海岸線も考慮しなければならない。
本研究でめざしたことは、等値線図の形で与えられたデータから、なんらかの空間内挿法によって任意の位置の値を推定し、使いやすい形式(たとえば緯度・経度等間隔のメッシュ形式)のディジタルデータとすることである。この問題は、いくつかに分割できる。
全世界の熱収支気候図集として有名なものは、Budyko (1963)である。ただし、この地図帳では、両極地方や、中国の一部に等値線がひかれていないところがあり、また海上の情報は少ないようである。そこで、海上については、Oberhuber (1988)で、中国については高・陸(1981)で置きかえることにした。Oberhuberがデータ不足として結果を示していない南半球高緯度の海上については、ソ連の南極アトラス(ATLAS, 1966)を利用した。ただしこれに含まれているのは１月と７月の正味放射だけである。なお、最近世界の各大陸の熱収支項の分布を示したHenning (1989)の地図帳が出版されたが、本研究では参照したにとどめ、ディジタイズは行なわなかった。
本研究で参照した熱収支気候図には、次のようなさまざまな投影法が使われている。(図3[HTML版では省略]参照。ただし一部は縮小された図を他の文献から引用したが、文献名の詳細は省略した。)
| Budyko (1963) | 世界 | 修正Winkel図法 (推定、後述) |
| Oberhuber (1988) | 海洋表面 | 方眼図法 |
| 高・陸 (1981) | 中国 | 単純円錐図法(推定) |
| ATLAS (1966) | 南極域 | ステレオ図法 |
| Henning (1989) | 世界の陸上 | モルワイデ図法 |
方眼図法(正距円筒図法ともいう)は、経度・緯度をそのままx, yにとって表示したものであり、等面積でも等角でもないが、データ処理上は最も単純な投影法である。ここで統合した結果を出力する図の形式としては、単位緯度と単位経度の長さが同じになる(つまり緯線・経線の方眼が正方形になる)方眼図法を採用した。一方、Oberhuberの図では緯度・経度1°あたりの長さが同じでない(方眼が長方形である)。
「推定」というのは、地図帳自体に図法の記述がなく、われわれが試験的に地図を作成して重ねることにより推定したものである。特に、Budyko (1963, 1971)に使われている地図投影法については、いずれの文献にも説明がない。最も多くの投影法の具体例をのせているSnyder and Voxland(1989)の本を見ても、一致するものは見つからなかった。そこで、この本にあるうちで最も近い Winkel-Tripel図法の式を修正して、Budykoの図と重なるように経験的にパラメタを決めた。これを仮に修正Winkel図法と呼ぶ。福井(1985)の世界地図に採用された図法も同じと思われる。
作業手順を図4に示す。計算機プログラムの名前がはいっているが、説明は省略する。基本的には、杉谷(1990、第10章、円錐図法の場合)と同じである。
上記[1]の、地図上の等値線の位置の読みとりは、ディジタイザを用いて手で等値線をなぞることによって行なった。データが大量になる場合は、スキャナで地図を画像として読みこみ、計算機によって等値線の位置を抽出したほうがよいと思われるが、他に経緯線や海岸線もある図から等値線を抽出するソフトウェアを組むことは大事業となると思われたので、今年度は避けた。
上記[2]は地図投影変換の逆変換である。これは図2bのように考えて処理した。経験式の決定と適用という意味で、考え方は図2aと似ている。
経度 x 緯度 y を地図投影変換式 f, gに代入して得られた地図上の位置を(ut, vt)とする。
ut = f(x, y)
(3)
vt = g(x, y)
一方、ディジタイザから得られた地図上の座標を(u, v)とすると、これは(ut, vt)との間に
ut = au + bv + c
(4)
vt = du + ev + f
という関係があると仮定する{4}。(4)はアフィン変換と呼ばれ、衛星データの幾何補正などによく使われている。この式の係数は、ディジタイザ上の地図の置き方にも依存する。
(4)式の係数を求めるためには、理論値(ut, vt){1}とディジタイザによる測定値(u, v){2}の両方がわかっている点が必要である。この基準点としては地図上の経線・緯線の交点を利用する。基準点を6個使えば6個の係数を求めることができるが、測定誤差があることも考慮し、数十個の基準点をとって最小二乗法で係数 a, b, c, d, e, f を決定した{5}。最小二乗法のやり方としては、(ut, vt) それぞれの残差二乗和が最小になるようにする方法と、(ut, vt)を両方合わせた残差二乗和が最小になるようにする方法があり、結果にも若干の相違が生じる。ここでは前者のやり方を採用したが、積極的根拠はない。 (4)式の係数が決まればこれを適用して{6}、ディジタイズされた任意の位置で(u, v){3}からの規格化された座標値(ut, vt){7}が求められる。
残された問題は(ut, vt)から(x, y)への、投影変換式(3)の逆変換である。方眼図法、単純円錐図法、ステレオ図法の逆変換式は、解析的(数学的)に簡単に書きくだせた。しかし修正Winkel図法は、(本来のWinkel-Tripel図法でも同様であるが)、式が複雑であり、逆変換式を見つけることはできなかった。そこでわれわれは、これを非線形連立微分方程式と見て、Newton-Raphson法により数値的に解くことにした。f, gのx, yに関する導関数を求めることはむずかしくない。
作業にはパソコンPC-9801VX (CPU=80286 10MHz, 80287つき)を使用した。実際に読みとり、逆変換して作成した図は、1月と7月の正味放射R、顕熱P、潜熱LEの合わせて6枚である。正味放射は4つ、他の変数は南極アトラスを除く3つの出典を使った。一つの出典の一つの図から約1万点の位置をディジタイズするのにかかった時間は30分前後であった。また、逆変換の計算をして緯度・経度の情報にする時間は、Oberhuber, 高・陸、ATLASのそれぞれについては1--3分程度だったが、Budykoの図については2時間かかった。4つの出典から逆変換した位置をまとめて作図するのには約20分かかった。したがって、図5[HTML版では省略]に示す1月の正味放射分布図を作るのにかかった時間は、プログラミングを除いて、約4時間半ということになる。このうちBudykoの図のNewton法を使った逆変換に非常に時間がかかっている。
等値線図の情報を利用するための実用的システムをめざすためには、さらに上記[3][4]の処理方法も考えなければならないが、[2]の段階でも改良が必要である。Newton法による逆変換では時間がかかりすぎる。地図を小部分に分けて逆変換を線形(アフィン変換)で近似する方法や、東大地震研究所の地図システム(大縮尺の地形図などを想定したもの)で採用しているように、逆変換を最良多項式で近似する方法も検討したい。
その他、実際に作業してみて気づいた細かい注意点を列挙する。
松山洋(東大・地理・学生)、沼田直美(東大・地球物理・学生)両氏の協力を得た。