地表面温度推定方法及びそのためのプログラム
【課題】 航空機等に搭載された熱赤外センサの観測値より、地表面温度を精度よく推定する。
【解決手段】 航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、熱赤外センサの観測条件を入力し、放射伝達プログラムの放射伝達計算結果から、観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される大気パラメータを得、得られた大気パラメータを用いて放射伝達方程式を演算して地表面輝度温度を得、そして、得られた地表面輝度温度と放射率から観測対象領域の地表面温度を推定することを特徴とする。
【解決手段】 航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、熱赤外センサの観測条件を入力し、放射伝達プログラムの放射伝達計算結果から、観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される大気パラメータを得、得られた大気パラメータを用いて放射伝達方程式を演算して地表面輝度温度を得、そして、得られた地表面輝度温度と放射率から観測対象領域の地表面温度を推定することを特徴とする。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、航空機等に搭載されたセンサによる観測値から地表面温度を推定する地表面温度推定方法に関する。
【背景技術】
【０００２】
従来、電磁波を用いて対象を非接触で計測するリモートセンンシングといわれる技術がよく知られている。特に、地球科学では、人工衛星や航空機などに搭載されたセンサで地球表面を観測する技術を指し、紫外域、可視域、赤外域、マイクロ波域と様々な電磁波が使用されている。使用する波長によって得られる情報は異なるが、熱赤外域（８−１２μｍの波長域）では温度の情報が得られる。
【０００３】
熱赤外域におけるリモートセンシングの観測値は、地表面の温度と放射率、大気の水蒸気・二酸化炭素・オゾンなどの吸収、エアロゾルの散乱、さらに大気自身からの放射など、様々な要素を含んでいる。
【０００４】
このように、人工衛星や航空機などに搭載の熱赤外センサで観測される温度には、観測する高度までの大気の影響が含まれるため、実際の地表面温度よりも低い値となるのが一般的であり、必ずしも地表面温度を正確に表していない。そこで、大気の影響を考慮して地表面温度を推定する必要がある。
【０００５】
大気の影響を考慮して地表面温度を推定する方法として、例えば、大気を複数のモデルに分類し、そのモデルを用いて大気の影響を評価して地表面温度を推定する方法が提案されている（例えば、特許文献１参照。）。この特許文献1による方法では、光の放射輝度データを導出するための放射伝達方式を熱赤外光の属性情報により簡略化し、人工衛星より海面温度を推定している。
【特許文献１】特開平１０−３１８８４４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、特許文献１に記載された方法では、放射伝達方程式の簡略化の過程で、大気の放射量を無視できる量としており、高温多湿などで大気の放射が大きい場合にはその影響を無視できず、誤差が大きくなる。
【０００７】
また、上記方法では大気モデルとして中緯度夏モデルを使用しているが、例えば日本国内の陸上で使用する場合、季節や場所により、大気の状態は大きく異なる。一例として、大気パラメータ（透過率、光路輝度、天空照度）を推定する放射伝達プログラムＭＯＤＴＲＡＮに標準で用意されているモデル（以下、「デフォルトモデル」とする。）による温度の鉛直分布を、図６Ａ，Ｂに示す。図６Ａは、各モデルにおける温度の鉛直分布を表し、図６Ｂは、図６Ａの要部を拡大したものである。この６モデルだけでは、季節や場所に応じた大気の状態を決め細かく反映することはできない。なお、ＭＯＤＴＲＡＮは、米国空軍ＡＦＲＬ（Air Force Research Laboratory）により開発されたプログラムである。
【０００８】
また、表面温度の推定には地表面における放射率を与えることが必要であるが、この値は観測対象領域の被覆により異なる。水面の放射率は、８−１２μｍの波長域ではほぼ１でよいが、例えば地表面などで放射率が１とできないとき、あるいはその変動が大きいときには誤差が大きく適用することができない。
【０００９】
本発明は斯かる点に鑑みてなされたものであり、航空機等に搭載された熱赤外センサの観測値より、地表面温度を精度よく推定することを目的とする。
【課題を解決するための手段】
【００１０】
上記課題を解決するため、本発明は、航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、熱赤外センサの観測条件を入力し、放射伝達プログラムの放射伝達計算結果から、観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される大気パラメータを得、得られた大気パラメータを用いて放射伝達方程式を演算して地表面輝度温度を得、そして、得られた地表面輝度温度と放射率から観測対象領域の地表面温度を推定することを特徴とする。
なお、上記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、また高層気象データは、少なくとも気圧、温度、湿度、高度から成る。
また、熱赤外センサの観測条件は、少なくとも観測対象領域の標高、航空機等の飛行高度、観測波長、観測角から成る。
【００１１】
上記構成によれば、上記大気の影響を評価する（大気パラメータを計算する）ときに、放射伝達モデルのデフォルトモデルを使用せず、必要な大気圧、温度、湿度等の鉛直プロファイルを観測時刻・観測対象領域における高層気象データや現地地上気象データより与えている。このため、観測時期や場所、大気の状態などによる誤差に影響されない。また、熱赤外センサの観測角を入力値とすることにより、観測対象領域の中心部と端部における大気の影響の差を反映することができる。これらのことから精度のよい大気パラメータが得られる。
【００１２】
また本発明は、上述した発明において、土地利用データを基に観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、観測対象領域の地表面温度を推定する構成とする。
【００１３】
上記構成とした場合、地表面の組成（被覆）を反映して放射率を決定することができ、精度のよい地表面輝度温度を得ることができる。
【発明の効果】
【００１４】
本発明によれば、観測対象上空から観測した熱赤外センサの観測値より、観測時期や場所、大気の状態などによらず、地表面温度を精度よく推定することができる。また、地表面の組成（被覆）を反映した精度のよい地表面温度を得ることができる。
【発明を実施するための最良の形態】
【００１５】
以下、航空機に熱赤外センサを搭載したリモートセンシングの例について説明する。図１は、リモートセンシングによる実際の地表観測における放射メカニズムを示したものである。図１に示す航空機４に搭載された熱赤外センサでは、地表からの放射１の他に、大気の上向き放射（光路輝度）２、および大気下方放射（天空照度）の地表面反射成分３が加算されて観測される。
【００１６】
大気中を通過する電磁波は、大気中に含まれる様々な物質によって吸収、散乱される。その度合いは電磁波の波長に強く依存する。この大気の吸収、散乱の中には、成層圏オゾンによる紫外線の吸収、二酸化炭素や水蒸気による熱赤外線の吸収、空気分子やエアロゾルによる散乱などがある。このような大気の影響のために、地表面の観測を目的とするリモートセンシングでは「大気の窓」と呼ばれる比較的大気による吸収、散乱の影響の少ない波長域が用いられる。
熱赤外域と呼ばれる８−１２μｍの波長域は「大気の窓」の１つであり、様々なリモートセンシングに使用されているが、大気の影響が完全にないわけではないので、定量的に計測を行うためには、大気の影響を補正する必要がある。
【００１７】
航空機に搭載の熱赤外センサによる観測値は、大気上端（飛行高度）における温度であるため、大気の吸収や散乱の影響を取り除く必要がある（大気補正）。大気補正手法として、航空機での観測と同期して地上観測を行い、両者を比較する方法、放射伝達計算に基づき、大気の影響を評価する方法がある。本発明は、後者の放射伝達計算に基づいて熱赤外センサの観測値の大気補正を行なうものである。
【００１８】
図２は、熱赤外センサによる観測値から地表面温度を推定するまでの手順を示す図である。観測値から地表面温度を推定するためには、以下の３つの過程を経る。
１．輝度校正
２．大気補正
３．温度−放射率分離
【００１９】
輝度校正とは、熱赤外センサが観測したエネルギーを表す電気信号(ＤＮ値)を物理量である輝度温度（もしくは放射輝度）に変換することである。一般的には、標準（黒体）を用いた実験などから、熱赤外センサ機種毎に固有の変換式が得られている。輝度校正より得られた輝度温度は、大気上端（航空機高度)における輝度温度である。この値には、大気の吸収や散乱、放射などの影響が含まれる。
【００２０】
大気補正とは、上述したような大気の影響を取り除き、地表面における輝度温度（もしくは放射輝度）を推定することである。得られた地表面輝度温度は、温度と放射率の関数で表される。
【００２１】
温度−放射率分離とは、地表面輝度温度から地表面温度Ｔｓと地表面放射率εを求めることをいう。
【００２２】
熱赤外センサによって波長λで観測される大気上端の放射輝度と地表面からの放射輝度（熱放射）は、以下の放射伝達方程式（１）、および地表面での熱放射式（２）でそれぞれ表される。
【数１】
【数２】
ここで、Ｉλ（Ｚ）、Ｉλ（Ｚ０）はそれぞれ大気上端高度Ｚでの放射輝度、地表面Ｚ０からの放射輝度を表している。τλ（Ｚ，Ｚ０）は高度Ｚ，Ｚ０間の波長λにおける透過率、Ｉ↑λ、Ｆ↓λ、ελはそれぞれ波長λにおける光路輝度、天空照度および地表の放射率、Ｔｓは地表面温度、Ｂλは波長λにおけるプランク関数である。なお、観測角の表記は省略している。
【００２３】
このような大気上端と地上間の大気の影響を考慮するために、大気補正手法が開発されてきた。代表的な大気補正手法には、単バンドアルゴリズム、昼夜アルゴリズム、差分吸収アルゴリズムなどがある。ここでは、単バンドアルゴリズムを例に説明する。
単バンドアルゴリズムは、観測時の大気情報を基に放射伝達計算を行う方法で、大気補正後に温度−放射率分離を行なうことにより地表面温度を得る。具体的には、式（１）に大気パラメータのうち透過率、光路輝度を与えて地表面の熱放射を求めた後、式（２）に残りの大気パラメータである天空照度及び地表面放射率を与えて地表面温度を求める。この場合、大気圧、温度、湿度といった鉛直プロファイル（情報）の与え方が、地表面温度の推定精度に大きく影響する。
【００２４】
以上より、航空機に搭載の熱赤外センサより得られたデータから精度よく地表面温度を推定するためには、大気補正、および温度−放射率分離手法の開発が重要である。
本発明においては、放射伝達プログラム（放射伝達モデル、あるいは放射伝達コードともいう。）によって、大気の影響を評価し、放射伝達方程式を用いて地表面の輝度温度を推定する。上記大気の影響を評価する（大気パラメータを計算する）ときに、放射伝達プログラムのデフォルトモデルを使用せず、必要な大気圧、温度、湿度等の鉛直プロファイルを観測時刻・観測対象領域における衛星データや現地気象データより与える。
【００２５】
図３は、大気補正、および温度−放射率分離を含む地表面温度推定処理のフローチャートである。なお、放射伝達計算（Ｓ５）を行う放射伝達プログラムに対するデータの入出力、および放射伝達計算以外の処理は、例えばコンピュータ装置が所定のプログラムを読み込んで、それを実行することにより行われる。プログラムは、通常は、コンピュータ装置と一体に組み込まれて使用されるが、コンピュータ装置が読み取り可能な記録媒体、例えばＣＤ−ＲＯＭ（コンパクトディスク型ＲＯＭ）等に記録され、地表面温度推定処理実行時にコンピュータ装置の記憶装置にインストールされるものであってもよい。
【００２６】
（大気補正）
本例では、単バンドアルゴリズムを使用し、放射伝達方程式を解く際に必要な大気パラメータ（透過率、光路輝度、天空照度）は、放射伝達プログラムＭＯＤＴＲＡＮにより導出する。本例においては、大気パラメータの計算（放射伝達計算）にＭＯＤＴＲＡＮを使用するが、これに限られるものではなく、観測時刻・観測対象領域における鉛直プロファイルから大気パラメータを導出できれば、他の放射伝達プログラムでも構わない。
【００２７】
上記放射伝達計算に使用する大気の鉛直プロファイル中の温度、湿度は、既往の６つのデフォルトモデル（熱帯、中緯度夏、中緯度冬、高緯度夏、高緯度冬、US Standard）では観測時刻や観測対象領域により大きく異なる大気状態を反映することができず、算出される地表面輝度温度の誤差も大きくなる。
【００２８】
そこで、観測対象領域の鉛直プロファイルとして、地上気象データは各地の測候所で観測される地上の気象データ（Ｓ２）より、高層気象データは全球解析データ（Ｓ３）より与える。エアロゾルやオゾンなど気体の量については、それほど大きな影響がないので、後述するデフォルト値を与えることとする。
【００２９】
なお、地上データについては、各測候所のデータに重み付けを行い内挿するようにしてもよい。例えば、観測対象領域と、隣接する測候所との距離に応じて重み付けを行い、重み付けされた気象データを放射伝達プログラムに入力する。
【００３０】
全球解析データとは、全球大気の３次元格子点毎に種々の気象パラメータを持つデータセットのことで、ラジオゾンデ、人工衛星、船舶、航空機、地上測候所などによる３次元気象観測データを取り組んで作成される。例えばＮＣＥＰ(National Centers for Environmental Prediction)データは、経度・緯度分解能は２．５°、鉛直方向は地表から２．７ｈＰａ間隔で２８層に区切られた格子点毎に、高度、温度、相対湿度、風速などの気象パラメータを有する。観測は、００Ｚ，０６Ｚ，１２Ｚ，１８Ｚの１日４回行われる（ＺはＵＴＣ時刻）。
【００３１】
さらに、上記地上気象データ（Ｓ２）および高層気象データ（Ｓ３）に加えて、観測対象領域の地表面の標高、航空機の飛行高度、赤外線センサの波長および観測角（あるいは入射角）のデータ（Ｓ４）を与える。航空機の場合、観測角の影響が大きい。例えば、高度７００ｋｍ上空から観測する人工衛星に搭載したセンサの場合、最大観測角は高々１０°であるので、影響は小さい。一方航空機搭載センサの場合、観測角が２５°を越えるものもある。
【００３２】
上述した各種データ（Ｓ２），（Ｓ３），（Ｓ４）と、熱赤外センサの観測輝度温度データ（Ｓ１）を、放射伝達プログラムＭＯＤＴＲＡＮに入力して放射伝達計算を行い（Ｓ５）、第１の大気パラメータ（透過率、光路輝度、天空照度）を得る（Ｓ６）。
【００３３】
ここで、熱赤外センサの感度特性（例えば計測可能な波長域の両端で感度が低い等）を補正するために、上記Ｓ６で得られた大気パラメータに対し、熱赤外センサの波長（バンド）に応じた重み付けを行う。図４に大気透過率の計算結果の一例を、図５に光路輝度の計算結果の一例を示す。図４，図５に示すように、大気透過率および光路輝度ともに波長によって計算値が異なり、特に８−９μｍ域では大気による吸収の影響が大きい。これらの計算結果に対して波長毎の重み付けを行い、重み付け後の大気パラメータ（透過率、光路輝度）を算出する（Ｓ７）。天空照度についても同様に行なう。
【００３４】
そして、重み付けされた大気パラメータのうち透過率、光路輝度を放射伝達方程式（１）に代入して演算し（Ｓ８）、地表面輝度温度を計算する（Ｓ９）。
【００３５】
なお、上記地上気象データおよび高層気象データ（全球解析データ）を用いた各大気の鉛直プロファイルは、一例として下記のように与える。
１）地上の気圧、気温、湿度は、１時間おきに観測される地上気象データのうち、航空機搭載の熱赤外センサが観測した時刻に最も近い時刻、観測対象領域に最も近い測候所のデータを与える。
２）高層の高度、気圧、気温、湿度は、６時間おきに観測されるＮＣＥＰデータのうち、航空機搭載の熱赤外センサが観測時刻に最も近い時刻、観測対象領域に最も近い格子点のデータを与える。
３）エアロゾルやオゾンなど気体の鉛直プロファイルは、国内で観測の場合、中緯度夏モデルか中緯度冬モデルのデフォルトモデルを与える(５−１０月：夏モデル、１１−４月：冬モデル)。
【００３６】
（温度−放射率分離）
上記大気補正を行った後、地表面温度を推定するためには、温度−放射率分離の処理が必要となる。地表面での熱放射式（２）より、地表面の輝度温度から地表面温度と放射率を求めることは、未知数が２つに対して方程式数が１つと少ないため、厳密には不可能である。そこで、放射率の値を仮定することにより、地表面温度を推定する（Ｓ１０）。
【００３７】
ここで、上記放射率は、地表面の組成（被覆）により異なる。精度よく地表面温度を推定するためには、地表面の組成に応じた放射率を与える必要がある。そこで、土地利用データなどを基に、土地被覆を種類別に数カテゴリーに分類し、それぞれのカテゴリー毎に所定の放射率を与え、地表面温度を推定する。
【００３８】
さらに、例えば国土地理院の１０ｍメッシュ土地利用データと熱赤外センサで得られる画像データを重ね、１画素毎に土地被覆のカテゴリーを与え、放射率を設定するようにしてもよい。このようにした場合、よりきめ細かい地表面温度の推定が可能になる。
【００３９】
以上の一連の方法により 地表面温度を得ることができるが、必ずしも地表面輝度温度から地表面温度を求める必要はなく、例えばユーザの目的に応じて、地表面輝度温度か地表面温度のどちらかを表示するようにしてもよい。
【００４０】
以上述べたように、本手法では放射伝達計算によって大気の影響を評価するので、観測時に地上での地表面温度を測定する必要がない。
【００４１】
また、上記の実施形態例によれば、観測対象領域の側に気象測候所があり、地上の気圧、気温、相対湿度等のデータが得られる国内外全ての領域に適用可能である。
【００４２】
また、上記実施形態例では、放射伝達プログラムに与える鉛直プロファイルにデフォルトモデルを用いず、地上気象データおよび高層気象データの観測データを使用するとともに、航空機の観測角などの観測条件を放射伝達プログラムに入力している。このため、地表面温度を精度よく推定することができる。
【００４３】
また、同一の地表面温度推定方法を様々な時期や地域に適用可能なため、観測時刻や観測対象領域が異なったデータ間の比較が可能である。
【００４４】
さらに、土地利用データと熱赤外センサで得られる画像データを重ね、例えば１画素毎に土地被覆のカテゴリーを与えて放射率を設定するようにした場合、よりきめ細かく地表面温度の推定が行なえる。
【００４５】
なお、本発明は、上述した各実施の形態例に限定されるものではなく、例えば人工衛星に搭載の熱赤外センサを用いて観測するなど、その他本発明の要旨を逸脱しない範囲において、種々の変形、変更が可能であることは勿論である。
【図面の簡単な説明】
【００４６】
【図１】リモートセンシングによる実際の地表観測における放射メカニズムを示した図である。
【図２】本発明の一実施形態例に係る観測値から地表面温度を推定するまでの手 順を示す図である。
【図３】本発明の一実施形態例に係る地表面温度推定方法を示すフローチャートである。
【図４】本発明の一実施形態例に係る大気透過率の計算結果の一例を示す図である。
【図５】本発明の一実施形態例に係る光路輝度の計算結果の一例を示す図である。
【図６】Ａは本発明の一実施形態例に係るＭＯＤＴＲＡＮのデフォルトモデルによる温度の鉛直分布、ＢはＡの要部の拡大図である。
【符号の説明】
【００４７】
１…地表からの放射成分、２…大気の上向き放射成分（光路輝度）、３…大気下方放射（天空照度）の地表面反射成分、４…航空機
【技術分野】
【０００１】
本発明は、航空機等に搭載されたセンサによる観測値から地表面温度を推定する地表面温度推定方法に関する。
【背景技術】
【０００２】
従来、電磁波を用いて対象を非接触で計測するリモートセンンシングといわれる技術がよく知られている。特に、地球科学では、人工衛星や航空機などに搭載されたセンサで地球表面を観測する技術を指し、紫外域、可視域、赤外域、マイクロ波域と様々な電磁波が使用されている。使用する波長によって得られる情報は異なるが、熱赤外域（８−１２μｍの波長域）では温度の情報が得られる。
【０００３】
熱赤外域におけるリモートセンシングの観測値は、地表面の温度と放射率、大気の水蒸気・二酸化炭素・オゾンなどの吸収、エアロゾルの散乱、さらに大気自身からの放射など、様々な要素を含んでいる。
【０００４】
このように、人工衛星や航空機などに搭載の熱赤外センサで観測される温度には、観測する高度までの大気の影響が含まれるため、実際の地表面温度よりも低い値となるのが一般的であり、必ずしも地表面温度を正確に表していない。そこで、大気の影響を考慮して地表面温度を推定する必要がある。
【０００５】
大気の影響を考慮して地表面温度を推定する方法として、例えば、大気を複数のモデルに分類し、そのモデルを用いて大気の影響を評価して地表面温度を推定する方法が提案されている（例えば、特許文献１参照。）。この特許文献1による方法では、光の放射輝度データを導出するための放射伝達方式を熱赤外光の属性情報により簡略化し、人工衛星より海面温度を推定している。
【特許文献１】特開平１０−３１８８４４号公報
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、特許文献１に記載された方法では、放射伝達方程式の簡略化の過程で、大気の放射量を無視できる量としており、高温多湿などで大気の放射が大きい場合にはその影響を無視できず、誤差が大きくなる。
【０００７】
また、上記方法では大気モデルとして中緯度夏モデルを使用しているが、例えば日本国内の陸上で使用する場合、季節や場所により、大気の状態は大きく異なる。一例として、大気パラメータ（透過率、光路輝度、天空照度）を推定する放射伝達プログラムＭＯＤＴＲＡＮに標準で用意されているモデル（以下、「デフォルトモデル」とする。）による温度の鉛直分布を、図６Ａ，Ｂに示す。図６Ａは、各モデルにおける温度の鉛直分布を表し、図６Ｂは、図６Ａの要部を拡大したものである。この６モデルだけでは、季節や場所に応じた大気の状態を決め細かく反映することはできない。なお、ＭＯＤＴＲＡＮは、米国空軍ＡＦＲＬ（Air Force Research Laboratory）により開発されたプログラムである。
【０００８】
また、表面温度の推定には地表面における放射率を与えることが必要であるが、この値は観測対象領域の被覆により異なる。水面の放射率は、８−１２μｍの波長域ではほぼ１でよいが、例えば地表面などで放射率が１とできないとき、あるいはその変動が大きいときには誤差が大きく適用することができない。
【０００９】
本発明は斯かる点に鑑みてなされたものであり、航空機等に搭載された熱赤外センサの観測値より、地表面温度を精度よく推定することを目的とする。
【課題を解決するための手段】
【００１０】
上記課題を解決するため、本発明は、航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、熱赤外センサの観測条件を入力し、放射伝達プログラムの放射伝達計算結果から、観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される大気パラメータを得、得られた大気パラメータを用いて放射伝達方程式を演算して地表面輝度温度を得、そして、得られた地表面輝度温度と放射率から観測対象領域の地表面温度を推定することを特徴とする。
なお、上記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、また高層気象データは、少なくとも気圧、温度、湿度、高度から成る。
また、熱赤外センサの観測条件は、少なくとも観測対象領域の標高、航空機等の飛行高度、観測波長、観測角から成る。
【００１１】
上記構成によれば、上記大気の影響を評価する（大気パラメータを計算する）ときに、放射伝達モデルのデフォルトモデルを使用せず、必要な大気圧、温度、湿度等の鉛直プロファイルを観測時刻・観測対象領域における高層気象データや現地地上気象データより与えている。このため、観測時期や場所、大気の状態などによる誤差に影響されない。また、熱赤外センサの観測角を入力値とすることにより、観測対象領域の中心部と端部における大気の影響の差を反映することができる。これらのことから精度のよい大気パラメータが得られる。
【００１２】
また本発明は、上述した発明において、土地利用データを基に観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、観測対象領域の地表面温度を推定する構成とする。
【００１３】
上記構成とした場合、地表面の組成（被覆）を反映して放射率を決定することができ、精度のよい地表面輝度温度を得ることができる。
【発明の効果】
【００１４】
本発明によれば、観測対象上空から観測した熱赤外センサの観測値より、観測時期や場所、大気の状態などによらず、地表面温度を精度よく推定することができる。また、地表面の組成（被覆）を反映した精度のよい地表面温度を得ることができる。
【発明を実施するための最良の形態】
【００１５】
以下、航空機に熱赤外センサを搭載したリモートセンシングの例について説明する。図１は、リモートセンシングによる実際の地表観測における放射メカニズムを示したものである。図１に示す航空機４に搭載された熱赤外センサでは、地表からの放射１の他に、大気の上向き放射（光路輝度）２、および大気下方放射（天空照度）の地表面反射成分３が加算されて観測される。
【００１６】
大気中を通過する電磁波は、大気中に含まれる様々な物質によって吸収、散乱される。その度合いは電磁波の波長に強く依存する。この大気の吸収、散乱の中には、成層圏オゾンによる紫外線の吸収、二酸化炭素や水蒸気による熱赤外線の吸収、空気分子やエアロゾルによる散乱などがある。このような大気の影響のために、地表面の観測を目的とするリモートセンシングでは「大気の窓」と呼ばれる比較的大気による吸収、散乱の影響の少ない波長域が用いられる。
熱赤外域と呼ばれる８−１２μｍの波長域は「大気の窓」の１つであり、様々なリモートセンシングに使用されているが、大気の影響が完全にないわけではないので、定量的に計測を行うためには、大気の影響を補正する必要がある。
【００１７】
航空機に搭載の熱赤外センサによる観測値は、大気上端（飛行高度）における温度であるため、大気の吸収や散乱の影響を取り除く必要がある（大気補正）。大気補正手法として、航空機での観測と同期して地上観測を行い、両者を比較する方法、放射伝達計算に基づき、大気の影響を評価する方法がある。本発明は、後者の放射伝達計算に基づいて熱赤外センサの観測値の大気補正を行なうものである。
【００１８】
図２は、熱赤外センサによる観測値から地表面温度を推定するまでの手順を示す図である。観測値から地表面温度を推定するためには、以下の３つの過程を経る。
１．輝度校正
２．大気補正
３．温度−放射率分離
【００１９】
輝度校正とは、熱赤外センサが観測したエネルギーを表す電気信号(ＤＮ値)を物理量である輝度温度（もしくは放射輝度）に変換することである。一般的には、標準（黒体）を用いた実験などから、熱赤外センサ機種毎に固有の変換式が得られている。輝度校正より得られた輝度温度は、大気上端（航空機高度)における輝度温度である。この値には、大気の吸収や散乱、放射などの影響が含まれる。
【００２０】
大気補正とは、上述したような大気の影響を取り除き、地表面における輝度温度（もしくは放射輝度）を推定することである。得られた地表面輝度温度は、温度と放射率の関数で表される。
【００２１】
温度−放射率分離とは、地表面輝度温度から地表面温度Ｔｓと地表面放射率εを求めることをいう。
【００２２】
熱赤外センサによって波長λで観測される大気上端の放射輝度と地表面からの放射輝度（熱放射）は、以下の放射伝達方程式（１）、および地表面での熱放射式（２）でそれぞれ表される。
【数１】
【数２】
ここで、Ｉλ（Ｚ）、Ｉλ（Ｚ０）はそれぞれ大気上端高度Ｚでの放射輝度、地表面Ｚ０からの放射輝度を表している。τλ（Ｚ，Ｚ０）は高度Ｚ，Ｚ０間の波長λにおける透過率、Ｉ↑λ、Ｆ↓λ、ελはそれぞれ波長λにおける光路輝度、天空照度および地表の放射率、Ｔｓは地表面温度、Ｂλは波長λにおけるプランク関数である。なお、観測角の表記は省略している。
【００２３】
このような大気上端と地上間の大気の影響を考慮するために、大気補正手法が開発されてきた。代表的な大気補正手法には、単バンドアルゴリズム、昼夜アルゴリズム、差分吸収アルゴリズムなどがある。ここでは、単バンドアルゴリズムを例に説明する。
単バンドアルゴリズムは、観測時の大気情報を基に放射伝達計算を行う方法で、大気補正後に温度−放射率分離を行なうことにより地表面温度を得る。具体的には、式（１）に大気パラメータのうち透過率、光路輝度を与えて地表面の熱放射を求めた後、式（２）に残りの大気パラメータである天空照度及び地表面放射率を与えて地表面温度を求める。この場合、大気圧、温度、湿度といった鉛直プロファイル（情報）の与え方が、地表面温度の推定精度に大きく影響する。
【００２４】
以上より、航空機に搭載の熱赤外センサより得られたデータから精度よく地表面温度を推定するためには、大気補正、および温度−放射率分離手法の開発が重要である。
本発明においては、放射伝達プログラム（放射伝達モデル、あるいは放射伝達コードともいう。）によって、大気の影響を評価し、放射伝達方程式を用いて地表面の輝度温度を推定する。上記大気の影響を評価する（大気パラメータを計算する）ときに、放射伝達プログラムのデフォルトモデルを使用せず、必要な大気圧、温度、湿度等の鉛直プロファイルを観測時刻・観測対象領域における衛星データや現地気象データより与える。
【００２５】
図３は、大気補正、および温度−放射率分離を含む地表面温度推定処理のフローチャートである。なお、放射伝達計算（Ｓ５）を行う放射伝達プログラムに対するデータの入出力、および放射伝達計算以外の処理は、例えばコンピュータ装置が所定のプログラムを読み込んで、それを実行することにより行われる。プログラムは、通常は、コンピュータ装置と一体に組み込まれて使用されるが、コンピュータ装置が読み取り可能な記録媒体、例えばＣＤ−ＲＯＭ（コンパクトディスク型ＲＯＭ）等に記録され、地表面温度推定処理実行時にコンピュータ装置の記憶装置にインストールされるものであってもよい。
【００２６】
（大気補正）
本例では、単バンドアルゴリズムを使用し、放射伝達方程式を解く際に必要な大気パラメータ（透過率、光路輝度、天空照度）は、放射伝達プログラムＭＯＤＴＲＡＮにより導出する。本例においては、大気パラメータの計算（放射伝達計算）にＭＯＤＴＲＡＮを使用するが、これに限られるものではなく、観測時刻・観測対象領域における鉛直プロファイルから大気パラメータを導出できれば、他の放射伝達プログラムでも構わない。
【００２７】
上記放射伝達計算に使用する大気の鉛直プロファイル中の温度、湿度は、既往の６つのデフォルトモデル（熱帯、中緯度夏、中緯度冬、高緯度夏、高緯度冬、US Standard）では観測時刻や観測対象領域により大きく異なる大気状態を反映することができず、算出される地表面輝度温度の誤差も大きくなる。
【００２８】
そこで、観測対象領域の鉛直プロファイルとして、地上気象データは各地の測候所で観測される地上の気象データ（Ｓ２）より、高層気象データは全球解析データ（Ｓ３）より与える。エアロゾルやオゾンなど気体の量については、それほど大きな影響がないので、後述するデフォルト値を与えることとする。
【００２９】
なお、地上データについては、各測候所のデータに重み付けを行い内挿するようにしてもよい。例えば、観測対象領域と、隣接する測候所との距離に応じて重み付けを行い、重み付けされた気象データを放射伝達プログラムに入力する。
【００３０】
全球解析データとは、全球大気の３次元格子点毎に種々の気象パラメータを持つデータセットのことで、ラジオゾンデ、人工衛星、船舶、航空機、地上測候所などによる３次元気象観測データを取り組んで作成される。例えばＮＣＥＰ(National Centers for Environmental Prediction)データは、経度・緯度分解能は２．５°、鉛直方向は地表から２．７ｈＰａ間隔で２８層に区切られた格子点毎に、高度、温度、相対湿度、風速などの気象パラメータを有する。観測は、００Ｚ，０６Ｚ，１２Ｚ，１８Ｚの１日４回行われる（ＺはＵＴＣ時刻）。
【００３１】
さらに、上記地上気象データ（Ｓ２）および高層気象データ（Ｓ３）に加えて、観測対象領域の地表面の標高、航空機の飛行高度、赤外線センサの波長および観測角（あるいは入射角）のデータ（Ｓ４）を与える。航空機の場合、観測角の影響が大きい。例えば、高度７００ｋｍ上空から観測する人工衛星に搭載したセンサの場合、最大観測角は高々１０°であるので、影響は小さい。一方航空機搭載センサの場合、観測角が２５°を越えるものもある。
【００３２】
上述した各種データ（Ｓ２），（Ｓ３），（Ｓ４）と、熱赤外センサの観測輝度温度データ（Ｓ１）を、放射伝達プログラムＭＯＤＴＲＡＮに入力して放射伝達計算を行い（Ｓ５）、第１の大気パラメータ（透過率、光路輝度、天空照度）を得る（Ｓ６）。
【００３３】
ここで、熱赤外センサの感度特性（例えば計測可能な波長域の両端で感度が低い等）を補正するために、上記Ｓ６で得られた大気パラメータに対し、熱赤外センサの波長（バンド）に応じた重み付けを行う。図４に大気透過率の計算結果の一例を、図５に光路輝度の計算結果の一例を示す。図４，図５に示すように、大気透過率および光路輝度ともに波長によって計算値が異なり、特に８−９μｍ域では大気による吸収の影響が大きい。これらの計算結果に対して波長毎の重み付けを行い、重み付け後の大気パラメータ（透過率、光路輝度）を算出する（Ｓ７）。天空照度についても同様に行なう。
【００３４】
そして、重み付けされた大気パラメータのうち透過率、光路輝度を放射伝達方程式（１）に代入して演算し（Ｓ８）、地表面輝度温度を計算する（Ｓ９）。
【００３５】
なお、上記地上気象データおよび高層気象データ（全球解析データ）を用いた各大気の鉛直プロファイルは、一例として下記のように与える。
１）地上の気圧、気温、湿度は、１時間おきに観測される地上気象データのうち、航空機搭載の熱赤外センサが観測した時刻に最も近い時刻、観測対象領域に最も近い測候所のデータを与える。
２）高層の高度、気圧、気温、湿度は、６時間おきに観測されるＮＣＥＰデータのうち、航空機搭載の熱赤外センサが観測時刻に最も近い時刻、観測対象領域に最も近い格子点のデータを与える。
３）エアロゾルやオゾンなど気体の鉛直プロファイルは、国内で観測の場合、中緯度夏モデルか中緯度冬モデルのデフォルトモデルを与える(５−１０月：夏モデル、１１−４月：冬モデル)。
【００３６】
（温度−放射率分離）
上記大気補正を行った後、地表面温度を推定するためには、温度−放射率分離の処理が必要となる。地表面での熱放射式（２）より、地表面の輝度温度から地表面温度と放射率を求めることは、未知数が２つに対して方程式数が１つと少ないため、厳密には不可能である。そこで、放射率の値を仮定することにより、地表面温度を推定する（Ｓ１０）。
【００３７】
ここで、上記放射率は、地表面の組成（被覆）により異なる。精度よく地表面温度を推定するためには、地表面の組成に応じた放射率を与える必要がある。そこで、土地利用データなどを基に、土地被覆を種類別に数カテゴリーに分類し、それぞれのカテゴリー毎に所定の放射率を与え、地表面温度を推定する。
【００３８】
さらに、例えば国土地理院の１０ｍメッシュ土地利用データと熱赤外センサで得られる画像データを重ね、１画素毎に土地被覆のカテゴリーを与え、放射率を設定するようにしてもよい。このようにした場合、よりきめ細かい地表面温度の推定が可能になる。
【００３９】
以上の一連の方法により 地表面温度を得ることができるが、必ずしも地表面輝度温度から地表面温度を求める必要はなく、例えばユーザの目的に応じて、地表面輝度温度か地表面温度のどちらかを表示するようにしてもよい。
【００４０】
以上述べたように、本手法では放射伝達計算によって大気の影響を評価するので、観測時に地上での地表面温度を測定する必要がない。
【００４１】
また、上記の実施形態例によれば、観測対象領域の側に気象測候所があり、地上の気圧、気温、相対湿度等のデータが得られる国内外全ての領域に適用可能である。
【００４２】
また、上記実施形態例では、放射伝達プログラムに与える鉛直プロファイルにデフォルトモデルを用いず、地上気象データおよび高層気象データの観測データを使用するとともに、航空機の観測角などの観測条件を放射伝達プログラムに入力している。このため、地表面温度を精度よく推定することができる。
【００４３】
また、同一の地表面温度推定方法を様々な時期や地域に適用可能なため、観測時刻や観測対象領域が異なったデータ間の比較が可能である。
【００４４】
さらに、土地利用データと熱赤外センサで得られる画像データを重ね、例えば１画素毎に土地被覆のカテゴリーを与えて放射率を設定するようにした場合、よりきめ細かく地表面温度の推定が行なえる。
【００４５】
なお、本発明は、上述した各実施の形態例に限定されるものではなく、例えば人工衛星に搭載の熱赤外センサを用いて観測するなど、その他本発明の要旨を逸脱しない範囲において、種々の変形、変更が可能であることは勿論である。
【図面の簡単な説明】
【００４６】
【図１】リモートセンシングによる実際の地表観測における放射メカニズムを示した図である。
【図２】本発明の一実施形態例に係る観測値から地表面温度を推定するまでの手 順を示す図である。
【図３】本発明の一実施形態例に係る地表面温度推定方法を示すフローチャートである。
【図４】本発明の一実施形態例に係る大気透過率の計算結果の一例を示す図である。
【図５】本発明の一実施形態例に係る光路輝度の計算結果の一例を示す図である。
【図６】Ａは本発明の一実施形態例に係るＭＯＤＴＲＡＮのデフォルトモデルによる温度の鉛直分布、ＢはＡの要部の拡大図である。
【符号の説明】
【００４７】
１…地表からの放射成分、２…大気の上向き放射成分（光路輝度）、３…大気下方放射（天空照度）の地表面反射成分、４…航空機
【特許請求の範囲】
【請求項１】
航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、
前記放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、前記熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、前記熱赤外センサの観測条件を入力する処理と、
前記放射伝達プログラムの放射伝達計算結果から、前記観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される前記大気パラメータを得る処理と、
得られた大気パラメータを用いて前記放射伝達方程式を演算し、前記地表面輝度温度を得る処理と、
得られた地表面輝度温度と放射率から前記観測対象領域の地表面温度を推定する処理から成る
ことを特徴とする地表面温度推定方法。
【請求項２】
前記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、
前記高層気象データは、少なくとも気圧、温度、湿度、高度から成る
ことを特徴とする請求項１に記載の地表面温度推定方法。
【請求項３】
前記熱赤外センサの観測条件は、少なくとも前記観測対象領域の標高、前記航空機等の飛行高度、観測波長、観測角から成る
ことを特徴とする請求項１または２に記載の地表面温度推定方法。
【請求項４】
前記放射伝達プログラムの放射伝達計算により得られた大気パラメータに対して、前記熱赤外センサの観測波長に応じた重み付けを行なった上で、前記放射伝達方程式の演算を行なう
ことを特徴とする請求項１〜３のいずれかに記載の地表面温度推定方法。
【請求項５】
土地利用データを基に前記観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、前記観測対象領域の地表面温度を推定する
ことを特徴とする請求項１〜４のいずれかに記載の地表面温度推定方法。
【請求項６】
航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定するプログラムであって、
前記放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、前記熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、前記熱赤外センサの観測条件を入力する手順と、
前記放射伝達プログラムの放射伝達計算結果から、前記観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される前記大気パラメータを得る手順と、
得られた大気パラメータを用いて前記放射伝達方程式を演算し、前記地表面輝度温度を得る手順と、
得られた地表面輝度温度と放射率から前記観測対象領域の地表面温度を推定する手順と
をコンピュータに実行させるための地表面温度推定プログラム。
【請求項７】
前記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、
前記高層気象データは、少なくとも気圧、温度、湿度、高度から成る
ことを特徴とする請求項６に記載の地表面温度推定プログラム。
【請求項８】
前記熱赤外センサの観測条件は、少なくとも前記観測対象領域の標高、前記航空機等の飛行高度、観測波長、観測角から成る
ことを特徴とする請求項６または請求項７に記載の地表面温度推定プログラム。
【請求項９】
前記放射伝達プログラムの放射伝達計算により得られた大気パラメータに対して、前記熱赤外センサの観測波長に応じた重み付けを行なった上で、前記放射伝達方程式の演算を行なう
ことを特徴とする請求項６〜８のいずれかに記載の地表面温度推定プログラム。
【請求項１０】
土地利用データを基に前記観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、前記観測対象領域の地表面温度を推定する
ことを特徴とする請求項６〜９のいずれかに記載の地表面温度推定プログラム。
【請求項１】
航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する方法であって、
前記放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、前記熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、前記熱赤外センサの観測条件を入力する処理と、
前記放射伝達プログラムの放射伝達計算結果から、前記観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される前記大気パラメータを得る処理と、
得られた大気パラメータを用いて前記放射伝達方程式を演算し、前記地表面輝度温度を得る処理と、
得られた地表面輝度温度と放射率から前記観測対象領域の地表面温度を推定する処理から成る
ことを特徴とする地表面温度推定方法。
【請求項２】
前記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、
前記高層気象データは、少なくとも気圧、温度、湿度、高度から成る
ことを特徴とする請求項１に記載の地表面温度推定方法。
【請求項３】
前記熱赤外センサの観測条件は、少なくとも前記観測対象領域の標高、前記航空機等の飛行高度、観測波長、観測角から成る
ことを特徴とする請求項１または２に記載の地表面温度推定方法。
【請求項４】
前記放射伝達プログラムの放射伝達計算により得られた大気パラメータに対して、前記熱赤外センサの観測波長に応じた重み付けを行なった上で、前記放射伝達方程式の演算を行なう
ことを特徴とする請求項１〜３のいずれかに記載の地表面温度推定方法。
【請求項５】
土地利用データを基に前記観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、前記観測対象領域の地表面温度を推定する
ことを特徴とする請求項１〜４のいずれかに記載の地表面温度推定方法。
【請求項６】
航空機等に搭載された熱赤外センサを用いて得られた観測データから、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定するプログラムであって、
前記放射伝達計算を行なう放射伝達プログラムに対し、少なくとも、前記熱赤外センサの観測データ、観測時の観測対象領域の状態を反映する地上気象データと高層気象データ、前記熱赤外センサの観測条件を入力する手順と、
前記放射伝達プログラムの放射伝達計算結果から、前記観測対象領域の地表面輝度温度を導出するための放射伝達方程式に使用される前記大気パラメータを得る手順と、
得られた大気パラメータを用いて前記放射伝達方程式を演算し、前記地表面輝度温度を得る手順と、
得られた地表面輝度温度と放射率から前記観測対象領域の地表面温度を推定する手順と
をコンピュータに実行させるための地表面温度推定プログラム。
【請求項７】
前記地上気象データは、少なくとも気圧、温度、湿度、標高から成り、
前記高層気象データは、少なくとも気圧、温度、湿度、高度から成る
ことを特徴とする請求項６に記載の地表面温度推定プログラム。
【請求項８】
前記熱赤外センサの観測条件は、少なくとも前記観測対象領域の標高、前記航空機等の飛行高度、観測波長、観測角から成る
ことを特徴とする請求項６または請求項７に記載の地表面温度推定プログラム。
【請求項９】
前記放射伝達プログラムの放射伝達計算により得られた大気パラメータに対して、前記熱赤外センサの観測波長に応じた重み付けを行なった上で、前記放射伝達方程式の演算を行なう
ことを特徴とする請求項６〜８のいずれかに記載の地表面温度推定プログラム。
【請求項１０】
土地利用データを基に前記観測データの土地被覆を所定種類に分類し、同種の土地被覆には同一の放射率を与え、前記観測対象領域の地表面温度を推定する
ことを特徴とする請求項６〜９のいずれかに記載の地表面温度推定プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２００７−３３０８（Ｐ２００７−３３０８Ａ）
【公開日】平成１９年１月１１日（２００７．１．１１）
【国際特許分類】
【出願番号】特願２００５−１８２６４４（Ｐ２００５−１８２６４４）
【出願日】平成１７年６月２２日（２００５．６．２２）
【出願人】（０００１３５７７１）株式会社パスコ (102)
【Ｆターム（参考）】
【公開日】平成１９年１月１１日（２００７．１．１１）
【国際特許分類】
【出願日】平成１７年６月２２日（２００５．６．２２）
【出願人】（０００１３５７７１）株式会社パスコ (102)
【Ｆターム（参考）】
[ Back to top ]