JP2019007836A - 温度測定方法及び装置 - Google Patents

Abstract

【課題】熱赤外分光放射計を用いた建築分野における気温分布の逆推定をする温度測定方法を提供する。【解決手段】建築物又は都市空間より放出される赤外線を観測し、赤外線観測値を分析する。気温をｘ、赤外線観測値をｙ、空気の温度分布に関する事前確率分布をＰ（ｘ）、赤外線観測値の観測誤差の確率密度関数をＰ（ｙ｜ｘ）、温度分布に関する事後分布をＰ（ｘ｜ｙ）とし、Ｐ（ｘ）及びＰ（ｙ｜ｘ）に基づいて、Ｐ（ｘ｜ｙ）を最大とする事後気温分布を推定する。【選択図】図７

Description

本発明は、建築物又は都市空間における気温分布を、熱赤外分光放射計を用いて測定する温度測定方法及び装置に関する。

建築空間や都市空間における気温分布を測定する場合、気温を測定したい場所にそれぞれ温度計を設置して測定するという方法がある。そのような方法によって、気温分布を測定するには多くの温度計が必要であり、実測の制約によっては、温度計を設置することが困難な場合もある。また、温度計は熱容量を持つため時定数が問題となる場合もあり、さらに測定原理から放射の影響による誤差が混入しやすいといった問題点が挙げられる。

特開昭６３−１００３３９号公報（特許文献１）には、熱容量の小さい材料で形成された検知板を、温度分布を測定すべき空間に設置し、その検知板の温度を測定する赤外線放射温度計によって、建物内部の温度分布を測定する装置が開示されている。

また、特開平１０−３１８８４４号公報（特許文献２）には、人工衛星から海面を撮影して得た熱赤外画像データと撮影時の大気状態情報とに基づいて、海面温度を推定する温度推定装置が開示されている。更に、リモートセンシング分野においては、人工衛星や地上から、分光放射計によって気温の鉛直分布を観測する手法が実用化されている。例えば、Kaplan,L.D.;Opt.Soc.Am.:49(10),1959,p.1004.（非特許文献１）において、Kaplan(1959)は気体の吸収係数が異なる複数の波長を観測することで気温の鉛直分布を得られる可能性があることが開示されている。また、Wark,D.Q.;Hilleary,D.T.;Science.165(3899),1969,p.1256.（非特許文献２）には、人工衛星により気温の鉛直分布の観測ができることが開示されている。また、特開２００７−３３０８号公報（特許文献３）には、航空機等に搭載された熱赤外センサの観測値により、地表面温度を精度よく推定する地表面温度推定方法及びそのためのプログラムが開示されている。さらに、大気の影響を評価する放射伝達計算に基づいて観測対象領域の地表面温度を推定する際、大気パラメータを含む放射伝達方程式を用いて地表面輝度温度を推定する方法が開示されている。このように、リモートセンシング分野においては、気温の鉛直分布（地表の上下方向に沿った気温分布）を測定することが開示されている。

ここで、従来の放射伝達過程の概要と観測方程式の定式化について説明する。

放射伝達過程については、放射は空気中を伝播する間に大気中の物質により吸収や散乱をうけて減衰すること、他の方向からの散乱や空気からの放射により増幅することが一般に知られている。その放射伝達過程を記述した放射伝達方程式の一般解は、境界面が黒体で、局所熱力学平衡状態とすると、観測方程式として、下記数１のように表せる。

数１について説明すると、数１の左辺は分光放射計に入射するエネルギー、数１における右辺第１項は境界面からの寄与、第２項は空気からの寄与を表している。なお、
であり、熱赤外光は波長が可視光線より長く、散乱の影響が殆どないので、散乱を無視できる。このため、光の放射伝達方程式による
は、数１のように近似することができる。

図１は、境界面からの放射が観測機（例えば熱赤外分光放射計）に入射する様子を示す模式図である。図１では、ｎ層に分割された空気層からの観測機に対する寄与を表しており、観測機に近い空気層ほど、観測機への寄与は大きいことを矢印の太さによって表している。ここで、数１における第２項の積分をn 層に離散化した空気層として計算する。そして、効率的なプログラミング及び理解の簡単化のため、観測機で測定した観測値をｙ、各空気層における気温の真値（真の気温分布）をｘとし、さらに誤差εを明示的に表現し、数１の観測方程式が下記数２のように簡略に表せると約束する。なお、数１のT(s’)が数２のｘ，数１のIv（s）が数２のｙに対応する。

各空気層における気温（真の気温分布）ｘが、観測値ｙに対して寄与する程度については、数３のように表すことができ、ヤコビアン（Jacobian matrix）と呼ばれる。

観測値ｙの数値計算をするため、数２について真の気温分布ｘの近傍ｘ₀において、数３のヤコビアンを用いて、数２をテイラー展開することができる。数２のテイラー展開は下記数４として表せる。

水平方向から観測機に入射されるエネルギーを計算した例を図２に示し、ヤコビアン（Jacobian matrix）を計算した例を図３に示す。図３に示されるように、ヤコビアンは境界面と観測機との距離の増加するにつれて単調に減少する。つまり、観測機から離れた位置の空気からの寄与が小さくなる。このようなヤコビアンの特性は、人工衛星からの鉛直方向の観測と異なり、水平方向の観測に特有なものである。

気温分布を求めることは、数２の観測値ｙから真の気温分布ｘを算出することである。その前提は、気体（空気）の濃度が既知であること、境界面の温度及び反射率が既知であることである。

そのような前提があれば、空気をｎ層に分割すること、及び観測チャンネルをｍ個とすることによって、未知数ｎ個を有するｍ個の連立方程式を得ることができる。ｍ個のチャンネルで観測できる観測機（熱赤外分光放射計）を用いて測定することによって、連立方程式の数ｍが、ｎ以上（ｎ＜ｍ）となるようにする。その連立方程式を解くことによって、各空気層における気温（気温分布ｘ）を算出することができる。

特開昭６３−１００３３９号公報 特開平１０−３１８８４４号公報 特開２００７−３３０８号公報

Kaplan,L.D.;Opt.Soc.Am.:49(10),1959,p.1004. Wark,D.Q.;Hilleary,D.T.;Science.165(3899),1969,p.1256.

しかしながら、前記数２及び数４が示すように、実際の観測値ｙは誤差εが付加される。そして、解は、誤差εに対して鋭敏に応答する。このような鋭敏性が問題となり、実際の観測においては、十分に意味のある信憑性のある連立方程式の解が得られないことが知られている。

このような背景で、リモートセンシング分野において、地表に対して水平方向の気温分布、及び境界面と観測機との距離が気温分布の測定に関する開示は見られない。そして、地表に対して鉛直方向の気温分布の人工衛星からの測定では、気体成分が高度の上昇につれ減少するという性質を利用することができる。その一方で、気温の水平分布を測定する場合には、そのような性質を利用することはできない。

本発明は上述のような事情に基づいてなされたものであり、本発明の目的は、建築分野において地表面に対して主に水平方向の気温分布を、熱分光放射計を用いて逆推定をする温度測定方法及び装置を提供することにある。

本発明に係る温度測定方法は、建築物又は都市空間より放出される赤外線を観測し、赤外線観測値を分析することにより、前記建築物又は前記都市空間における空気の温度分布を推定することにより達成される。

また、本発明に係る温度測定方法の上記目的は、気温をｘ、前記赤外線観測値をｙ、前記温度分布に関する事前確率分布をＰ（ｘ）、前記赤外線観測値の観測誤差の確率密度関数をＰ（ｙ｜ｘ）、前記温度分布に関する事後分布をＰ（ｘ｜ｙ）とし、前記Ｐ（ｘ）及び前記Ｐ（ｙ｜ｘ）に基づいて、前記Ｐ（ｘ｜ｙ）を最大とする事後気温分布を推定することにより、或いは前記Ｐ（ｘ）、前記Ｐ（ｙ｜ｘ）及び前記Ｐ（ｘ｜ｙ）の関係は、式（１）により表され、
前記（１）式の右辺を最大とする前記事後気温分布を推定することにより、或いはガウス・ニュートン法を用いて、前記事後気温分布を最大にするＭＡＰ解を算出することより、より効果的に達成される。

また、本発明に係る温度測定装置の上記目的は、建築物又は都市空間より放出される赤外線を観測する観測手段と、赤外線観測値に基づいて、前記建築物又は前記都市空間における空気の温度分布を推定する推定手段とを有することにより達成される。

また、本発明に係る温度測定装置の上記目的は、前記推定手段は、気温をｘ、前記赤外線観測値をｙ、前記温度分布に関する事前確率分布をＰ（ｘ）、前記赤外線観測値の観測誤差の確率密度関数をＰ（ｙ｜ｘ）、前記温度分布に関する事後分布をＰ（ｘ｜ｙ）とし、前記Ｐ（ｘ）及び前記Ｐ（ｙ｜ｘ）に基づいて、前記Ｐ（ｘ｜ｙ）を最大とする事後気温分布を推定することにより、或いは前記Ｐ（ｘ）、前記Ｐ（ｙ｜ｘ）及び前記Ｐ（ｘ｜ｙ）の関係は、式（１）により表され、
前記推定手段は、前記（１）式の右辺を最大とする前記事後気温分布を推定することにより、或いはガウス・ニュートン法を用いて、前記事後気温分布を最大にするＭＡＰ解を算出することにより、より効果的に達成される。

本発明によれば、観測機（例えば、熱赤外分光放射計）を用いて、建築物、建築空間又は都市空間における気温分布の逆推定をすることができる。また、建築空間や都市空間における気温の平均値を得ることにも有効であり、代表性のある気温を得ることもできる。また、この方法は多くの温度計を必要とせず、さらに温度計の熱容量や放射による誤差の影響などを考慮しなくてよい。

熱赤外分光放射計に入射するエネルギーを示す模式図である。 水平方向から観測機（熱赤外分光放射計）に入射されるエネルギーを計算した例である。 ヤコビアン（Jacobian matrix）を計算した例を示す特性図である。 本発明の構成例を示すブロック図である。 本発明の分光放射計の構成を示すブロック図である。 本発明の演算部の構成を示すブロック図である。 感度解析の具体的な動作を示すフローチャートである。 標準ケースの模式図である。 標準ケースのＡＫ（アベレ−ジングカーネル）及びＡＫareaを示す図である。 標準ケースのＭＡＰ解を示す図である。 事前確率分布が１１±５℃である場合のＭＡＰ解を示す図である。 気温が一様でない場合（４〜５層目の気温が低下する場合）のＭＡＰ解を示す図である。 観測機と境界面との距離が１０ｍの場合のＭＡＰ解を示す図である。 観測機と境界面との距離が１０００ｍの場合のＭＡＰ解を示す図である。 空気層の分割数を３にした（３層に離散化の）場合のＭＡＰ解を示す図である。 空気層の分割数を１０にした（１０層に離散化の）場合のＭＡＰ解を示す図である。 観測機の分解能（ＦＷＨＭ）を標準ケースの１／２倍である３０［cm^-1］）の場合のＭＡＰ解を示す図である。 観測機の分解能（ＦＷＨＭ）を標準ケースの３倍である５［cm^-1］）の場合のＭＡＰ解を示す図である。 水蒸気に−５％ｒｈの推定誤差がある場合のＭＡＰ解を示す図である。 水蒸気に＋５％ｒｈの推定誤差がある場合のＭＡＰ解を示す図である。 屋外実験の空間の平面図を示す図である。 屋内実験の空間の平面図を示す図である。 実験に用いた分光放射計の写真である。 実験に用いた平面黒体炉の写真である。 屋内実験におけるＭＡＰ解と実測温度とを示す図である。 屋外実験におけるＭＡＰ解と実測温度とを示す図である。 屋内実験におけるＡＫ（アベレ−ジングカーネル）及びＡＫareaを示す図である。 屋外実験におけるＡＫ（アベレ−ジングカーネル）、及びＡＫareaを示す図である。

本発明においては、上述のような誤差による鋭敏性を解消するために、最大事後確率推定法（以下、ＭＡＰ（Maximum A Posteriori)法とする）を用いる。なお、ＭＡＰ法は、ベイズの定理を用いて事後確率分布を推定する方法である。

以下、建築空間において、分光放射計を用いて地表面に対し水平方向の気温分布を逆推定する方法の原理について説明する。

ここで先ず、ベイズの定理を簡単に説明する。

気温ｘ及び観測値（データ：観測結果）ｙの確率は、それぞれＰ（ｘ）及びＰ（ｙ）と表わされる。特に気温ｘの確率Ｐ（ｘ）は、気温ｘに関する事前確率分布を表している。また、観測値ｙという条件下における気温ｘ（事象ｘ）の条件付確率はＰ（ｘ｜ｙ）と表され、事後確率分布と呼ばれる。そして、気温ｘ（事象ｘ）という条件下における観測値ｙの条件付確率はＰ（ｙ｜ｘ）と表され、観測誤差の確率密度関数若しくは尤度と呼ばれる。なお、確率Ｐ（ｘ）及びＰ（ｙ｜ｘ）については、本発明では、正規分布を仮定している。また、本発明の実施形態のＭＡＰ法で用いられるベイズの公式は、数５のように表される。

なお、Ｐ（ｙ）は、観測データｙが生起する確率であり、ＭＡＰ推定においては、計算する必要はない。即ち、観測値ｙが定数であることを考慮し、数５の右辺定数項Ｐ（ｙ）を無視すると、数６のようなベイズの定理（「事後確率分布∝尤度×事前確率分布」）になる。

ＭＡＰ法においては、数６のベイズの定理によって与えられた事後確率分布Ｐ（ｘ｜ｙ）が最大となるような気温分布ｘをＭＡＰ解x_MAPとして求める。

また、ＭＡＰ解x_MAPを求める場合には、ＭＡＰ法は、観測前に知りうる、気温ｘに関する事前確率分布Ｐ（ｘ）（確率的な情報）を用いる。

ＭＡＰ法においては、Ｐ（ｙ｜ｘ）Ｐ（ｘ）が最大となるようなｘを求め、見つかったｘがＭＡＰ解x_MAPである。

一方、事後確率分布Ｐ（ｘ｜ｙ）を最大にするｘは、対数をとったｌｏｇ（Ｐ（ｘ｜ｙ））を最大にするｘと一致する。

また、尤度Ｐ（ｙ｜ｘ）及びＰ（ｘ）に関する確率分布は、全て正規分布を仮定しているため、例えば、尤度Ｐ（ｙ｜ｘ）は数７のように表せる。なお、特に、観測誤差に関する共分散行列をＳε、ＭＡＰ解の分散をＳ_MAP、事前確率分布の分散をＳ_aと表す。

したがって、尤度Ｐ（ｙ｜ｘ）に関して、数８のようになる。また、事前確率分布Ｐ（ｘ）に関して、数９のようになる。

そして、数８及び数９におけるconstは定数である。そのため、Ｐ（ｙ｜ｘ）Ｐ（ｘ）は、数１０のように表される。

本発明における目的は、Ｐ（ｙ｜ｘ）Ｐ（ｘ）を最大にするｘを見つけることである。それは数１０の右辺である、
を最小にするようなｘを見つけることで達成される。そのためには、一般的にガウス・ニュートン法を用いることができる。

一方で、ＭＡＰ解x_MAPには事前確率分布からの影響が含まれるため、真値との差のみで評価することは危険である。そこで、事前確率分布P(x)がＭＡＰ解x_MAPに与える影響を評価するため、後述するように、Root Mean Squared Error（以下、ＲＭＳＥとする）、アベレージングカーネル（以下ＡＫ、又は「Ａ」と表記する）、ＡＫarea、及びDegree Of Freedom for Signal（以下、ＤＯＦＳとする）という指標を用いる。（ｘ_aは事前確率分布P(x)の期待値（平均）を表す。）
まず、アベレージングカーネルについて説明する。ＭＡＰ解x_MAPを、真値（真の気温分布）ｘ、アベレージングカーネルＡ、事前確率分布の期待値ｘ_a、単位行列Iを用いて表現すると、数１１のように表せる。

Ａ≡ＧＫ
により定義されるＡをアベレージングカーネル（Averaging kernel）と呼ぶ。

数１１によれば、ＭＡＰ解x_MAPとは、真の気温分布ｘと事前確率分布ｘ_aとを、それぞれアベレージングカーネルＡと、単位行列ＩからアベレージングカーネルＡを引いたＩ−Ａによって重み付け平均したものであること示している。つまり、Ａが単位行列に近いほど、ＭＡＰ解x_MAPは正しく気温分布を反映しているということである。ＡをＡ≡ＧＫにより計算することで、ＭＡＰ解x_MAPを評価できる。なお、共分散行列Ｓ（分散共分散行列とも言う。）とは，分散（散らばり具合を表す指標）の概念を多次元確率変数に拡張して行列としたものである。また、共分散行列は変数（事前確率分布、観測値）の散らばり具合や相関という情報を集約したものである。

次に、ＲＭＳＥ及びＤＯＦＳについて説明する。

ＲＭＳＥは、ＭＡＰ解x_MAPと真の気温分布ｘとに基づいて計算された２乗平均平方根誤差である。また、アベレージングカーネルＡのトレース tr(Ａ)は、ＤＯＦＳと呼ばれる数で、観測によって明らかになった未知数の数を表す。また、アベレージングカーネルＡを行方向に足したものはＡＫareaと呼ばれ、各位置における観測データが寄与する程度を示すものである。これらは、ＭＡＰ解x_MAPの有効性の指標である。

以下に、本発明の実施の形態を、図面を参照して詳細に説明する。

図４は、本発明の温度測定装置１０の構成例を示すブロック図である。温度測定装置１０は、熱源１１（例えば、黒体炉）からの放射される赤外線の分光を観測する分光放射計１２、及び観測値に基づいて温度分布を演算する演算部１３で構成される。

そして、図５は、本発明の分光放射計１２の構成を示すブロック図である。分光放射計１２は、紫外から遠赤外までの広い波長域の光線の分光測定が可能である。図５は、分光放射計１２を入射部１２０、集光部１２１、ハーフミラー１２２、ＣＣＤカメラ１２３、チョッパ１２５、内部黒体１２６、連続可変フィルタ（以下、ＣＶＦとする）１２７及び検出部１２８を用いて構成した例である。

まず、入射部１２０を介して、黒体炉１１から放射される光線が、集光部１２１に入射される。集光部１２１は、光線を集光し、ハーフミラー１２２に光線を出射する。ハーフミラー１２２を通過した光線は、ＣＣＤカメラ１２３に入射される。ＣＣＤカメラ１２３は、入射した光線の映像を映像信号に変換する。なお、映像信号はコントローラ１２４に出力され、コントローラ１２４は、映像信号を後述する表示部１３ｂに入力できるデジタル信号に変換するようにしてもよい。

また、ハーフミラー１２２によって反射された光線は、チョッパ１２５に入射される。そして、チョッパ１２５が閉じて該光線を遮断した際には、内部黒体１２６の放射が、後述するＣＶＦ１２７を介して、検出部１２８に出射される。一方、チョッパ２５が開いて該光線を通した際には、集光部１２１からの該光線がＣＶＦ１２７を介して、検出部１２８に入射される。

ＣＶＦ１２７は、円周上に分光特性が異なるフィルタを配置されているものである。ＣＶＦ１２７は、所定の軸の回りで回転するように構成されている。例えば、ＣＶＦ１２７は、図示していないがモータによって回転するようにしても良く、モータの回転角度に応じて、光線は分光特性の異なるフィルタを通過することができる。そうすることによって、ＣＶＦ１２７を１回転させるごと、観測値ｙ（例えば、分光データ）を取得することができる。

コントローラ１２４には、チョッパ１２５、ＣＶＦ１２７及び検出部１２８との間で制御線及びデータ線が配設され、制御又はデータの入出力することができる。これにより、例えばチョッパ１２５及びＣＶＦ１２７それぞれの回転角速度は、コントローラ１２４を用いて制御することができる。また、検出部１２８は光線の輝度値に応じたデジタル信号に変換することができるし、コントローラ１２４によって、検出部１２８が変換をするタイミング、さらにデジタル信号の形式を選択することができる。

なお、分光放射計１２における、光学系（入射部１２０、集光部１２１、ハーフミラー１２２）、ＣＶＦ１２７、検出部１２８の特性を変更することにより、測定波長の範囲、または測定する視野角を変更することができる。また、ＣＶＦ１２７を任意の波長で固定して使用すると、その波長の放射量の経時変化が測定することができる。以上の説明したように分光放射計１２は、一般的な分光光度計（例えば、ＦＴＩＲ）と構成が異なり、測定試料を装置内部に入れるタイプではない。このため、大きな測定物、加熱・冷却した測定物、または遠方の対象物を計測することができる特徴がある。

また、分光放射計１２で測定した観測値ｙ（例えば、分光データ）は、専用のソフトウェアで様々な解析することができる。例として、ＣＳＶ形式で出力することができるため、一般的の表計算ソフト及び数値計算を用いて解析することも可能である。

次に、図６は、本発明のＭＡＰ法を用いた気温分布の逆推定を行う演算部１３の構成を示すブロック図である。

演算部１３は、ＣＰＵ１３ａ及び表示部１３ｂで構成され、ＣＰＵ１３ａは、放射伝達計算部１、ヤコビアン記憶部２、制御部３、入力部４、最適解処理部５、平均二乗誤差算出部６、事前確率分布記憶部７で構成されている。各構成との関係、機能及び動作を順に説明する。

先ず、放射伝達計算のためのパラメータ（距離、境界面の温度、離散化（分割数）のための層の幅（分割幅）、気体（空気）の濃度）が、放射伝達計算部１に入力される。

次に、放射伝達計算部１は、パラメータに基づいて、ＭＯＤｅｒａｔｅ ｒｅｓｏｌｕｔｉｏｎ ＴＲＡＮｓｍｉｓｓｉｏｎ（以下、ＭＯＤＴＲＡＮとする）を用いて、放射伝達計算を行い、ヤコビアンＫを算出する。

そして、放射伝達計算部１は、ヤコビアンＫをヤコビアン記憶部２に出力し、ヤコビアン記憶部２はヤコビアンＫを記憶する。

また、制御部３は、入力部４、最適解処理部５、平均二乗誤差算出部６のそれぞれに対して、制御信号を入出力して制御する。

具体的には、入力部４を介した制御部３からの制御信号に応じたヤコビアン記憶部２は、ヤコビアンＫを入力部４に出力する。一方、分光放射計１２からの観測値ｙ（観測データ）は、入力部４に入力され、観測値ｙは、演算部１３（例えば、コンピュータ）内部での処理に適したデータ形式に変換される。また、制御部３からの制御信号に応じて、事前確率分布記憶部７から事前確率分布Ｐ（ｘ）が最適解処理部５に入力される。更に、制御部３からの制御信号に応じて、入力部４から、観測値ｙ及びヤコビアンＫが最適解処理部５に入力される。

そして、制御部３からの制御信号に応じて、最適解処理部５は、事前確率分布Ｐ（ｘ）、観測値ｙ及びヤコビアンＫに基づいて、ガウス・ニュートン法によって、最適解処理を行う。制御部３は制御信号を最適解処理部５に所定の回数出力し、ＭＡＰ解x_MAPを計算させる。そして、所定の回数の繰り返し計算によって、事後確率分布Ｐ（ｘ｜ｙ）が最大となるようなＭＡＰ解x_MAP及び分散を計算することができる。

そして、最適解処理部５は、所定の複数回の繰り返し計算が終了したことを示す制御信号を制御部３に出力すると、制御部３は、最適解処理部５に対して平均二乗誤差算出部６へＭＡＰ解x_MAP及び分散S_MAPを出力するように指示する制御信号を出力する。そして、平均二乗誤差算出部６は、ＭＡＰ解x_MAPに基づいて、ＲＭＳＥ及びＤＯＦＳの計算を実行する。

最後に、ＭＡＰ解x_MAP、分散、ＲＭＳＥ及びＤＯＦＳという統計データが、演算部１３の表示部１３ｂ（例えば、コンピュータのモニタ）に表示される。

次に、ＭＡＰ法を用いた数値シミュレーションによる逆推定の精度検証と感度解析の方法とを説明する。

感度解析の具体的な動作を図７に示すフローチャートを用いて説明する。先ず、パラメータ（距離、境界面の温度、離散化（分割数）のための層の幅（分割幅）、真の気温分布及び気体（空気）の濃度）がシステムに入力される（ステップＳ１０）。放射伝達過程の計算を行うプログラムであるＭＯＤＴＲＡＮを用いて、放射伝達計算を行い、真の観測値及びヤコビアンＫを算出し記憶する（ステップＳ２０）。なお、ＭＯＤＴＲＡＮは、波数に対する大気の透過率及び放射輝度を計算する。次に、真の観測値ｙ及びヤコビアンＫをファイルなどから読み出す（ステップＳ３０）。次に、真の観測値ｙに対し正規分布に従うランダム誤差を伴った観測値ｙ＋ε及びヤコビアンＫが入力される（ステップＳ４０）。ここで、気温に関する事前確率分布（事前情報）P(x)を読み込む（ステップS４１）。事前確率分布と誤差が加算された観測値に対して、ガウス・ニュートン法を用いて、事後確率分布Ｐ（ｘ｜ｙ）が最大となるようなＭＡＰ解x_MAPを計算する（ステップS５０）。気温に関するＭＡＰ解x_MAP及び分散S_MAPを出力する（ステップ６０）。ＭＡＰ解x_MAPを１０００回計算したか否か判定する（ステップＳ７０）。繰り返し計算した回数が１０００回未満であれば、ステップＳ４０に戻り、繰り返し計算した回数が１０００回であれば、ステップＳ８０に進む。その計算結果に基づいて、ＭＡＰ解x_MAPについてのＲＭＳＥ（ＲＭＳ_error）及びＡＫのトレースであるＤＯＦＳを計算する（ステップＳ８０）。

次に、数値シミュレーションによって、建築スケールでの水平気温分布逆推定の実現可能性を明らかにする。また、逆推定に影響を及ぼすパラメータの変化と、ＭＡＰ解x_MAPに与える影響との関係を調査するために感度解析をした結果を説明する。そこで、ＭＡＰ法を用いた数値シミュレーションの対象は、表１及び図８に示すように、観測機から境界面までの距離を１００ｍ、気温を１５℃一様に設定された街区スケール（以下、標準ケースとする。）を想定する。また、表２に示すように、標準ケースの事前確率分布は１７±５℃とする。また、標準ケースのＡＫ及びＡＫareaを図９に示す。図９において、1層目のＡＫareaは、“1”を超えているため、観測による寄与が大きいと考えられる。観測機からの距離が離れるほど、ＡＫareaが単調減少しており，ＭＡＰ解x_MAPに対する事前確率分布Ｐ（ｘ）の寄与が増加することが示される。

また、標準ケースにおけるＭＡＰ法の入力条件を、下記表２に示す。

そして、標準ケースにおけるＭＡＰ解x_MAPを図１０に示す。標準ケースにおける1層目のＭＡＰ解x_MAPは真値に近かった。しかし、観測機からの距離が大きくなるにつれ，ＭＡＰ解x_MAPは事前確率分布に近づく。１層目におけるＭＡＰ解x_MAPの標準偏差は最小である。一方、２層目のＭＡＰ解x_MAPの標準偏差は最大であり、さらに距離が大きくなるにつれて、２層目以降のＭＡＰ解x_MAPの標準偏差は減少する。このような傾向によって、気温の水平分布の逆推定は特徴付けられる。結果的に、ＭＡＰ解x_MAPのＲＭＳＥである１．２±０．６℃は、事前確率分布のＲＭＳＥである２．０℃より小さいから、標準ケースにおけるＲＭＳＥは減少している。結論的に、気温分布の逆推定に対して、ＭＡＰ法が有効であることが示された。なお、標準ケースにおける、ＲＭＳＥ、ＤＯＦＳは、それぞれ１．２、１．７であり、後述するパラメータ（条件）の影響の調査の際に基準となる数値である。

次に、気温分布の逆推定に対するパラメータ（条件）の影響を調査する。パラメータ（条件）としては、事前確率分布、一様でない気温（真の気温分布が一様でない。）、距離（境界面と観測機との間の距離）、層分割ｎ（空気層の分割数）、観測機の分解能（以下、ＦＷＨＭ（Full Width at Half Maximum）［cm^-1］）という）、水蒸気の推定誤差が挙げられる。

先ず、事前確率分布を１７±５℃から１１±５℃に変化させた場合にけるＭＡＰ解x_MAPは図１１に示すように、観測値近傍での精度は大きく変化しない。

次に、真の気温が一様でない気温分布（３層目の気温が相対的に高く、１〜２層目及び４〜５層目の気温が相対的に低い気温分布）に対して逆推定を行った結果、図１２に示すように、１〜３層目の気温の上昇は逆推定することができたものの、４〜５層目における気温の低下を推定することができなかった。このような結果から、本発明の気温分布の逆推定の方法は、観測機に近づくにつれ、事前確率分布の影響が小さくなり、推定された気温の精度が向上する傾向にあることが示された。

次に、境界面と観測機との間の距離を１０ｍ及び１０００ｍに設定した空気層に対して、気温分布の逆推定を行った結果をそれぞれ図１３及び図１４に示す。図１３及び図１４が示すように、いずれの場合でも、ＲＭＳＥ及びＤＯＦＳは大きく変化しなかった。したがって、本発明の気温分布の逆推定の方法は、建築スケール及び街区スケールのいずれにおいても同様の精度が得られることが示された。

次に、空気層の分割数を３にした空気層に対して、気温分布の逆推定を行った結果を図１５に示す。図１５が示すように、ＤＯＦＳは１．７であった。これに対し、ＲＭＳＥは１．３±０．８℃であり、標準ケースにおけるＲＭＳＥである１．２±０．６℃に比べ、わずかに悪化した。また、空気層の分割数を１０にした空気層に対して、気温分布の逆推定を行った結果を図１６に示す。図１６が示すように、ＤＯＦＳは１．７であった。これに対し、ＲＭＳＥは１．１℃へとわずかに改善した。このような結果については、ＤＯＦＳが向上すると観測値の寄与が大きくなる。このため、観測値のランダム誤差εの影響も大きくなったためと考えられる。

次に、観測機の分解能（ＦＷＨＭ）を１／２倍の３０［cm^-1］にし、気温分布の逆推定を行った結果を図１７に示す。図１７が示すように、ＤＯＦＳは、１.６に減じ、ＲＭＳＥは１．２±０．６℃のまま変化しなかった。また、観測機の分解能（ＦＷＨＭ）を３倍の５［cm^-1］にし、気温分布の逆推定を行った結果を図１８に示す。図１８が示すように、ＤＯＦＳは、２．１に上昇し、ＲＭＳＥは１．１±０．５℃になって改善できた。また、観測機から離れた空気層における気温の推定の精度は、平均的に向上できた。結果的に、観測機の分解能を上げることによって、精度は向上できた。

次に、水蒸気の相対湿度ｒｈ（relative humidity）に推定誤差がある空気層に対して、気温分布の逆推定を行った結果を図１９及び図２０に示す。前述のようなヤコビアンを作成する際に，水蒸気に−５％ｒｈの推定誤差がある場合、図１９が示すように、ＲＭＳＥが、標準ケースにおける１．２±０．６℃から２．０±０．９℃に変化（悪化）し、偏った推定結果となる傾向が示された。また、水蒸気の相対湿度ｒｈに＋５％ｒｈの推定誤差がある場合、図２０が示すように、ＲＭＳＥは１．５±０．８℃に変化（悪化）した。

まとめると、実施形態における気温分布の逆推定に関する感度解析においては、ＭＡＰ解x_MAPのＲＭＳＥは、事前確率分布のＲＭＳＥより改善することが示された。

次に、実施形態における建築スケールの水平気温分布の逆推定について説明する。

本発明の実施形態においては、温度が異なる（温度分布が存在する）２室を対象として、水平気温分布の逆推定を行った（屋内実験）。さらに、屋内外の温度差を有する建築空間を対象として、水平気温分布の逆推定を行った（屋外実験）。屋内実験の空間、屋外実験の空間のそれぞれの平面図を図２１及び図２２に示す。また、実験に用いた分光放射計及び平面黒体炉の仕様をそれぞれ表３及び表４に、また図２３及び図２４にそれぞれの写真を示す。屋内実験及び屋外実験の境界面の黒体炉の温度は、５０℃に設定した。

また、ＭＡＰ法に用いた条件を表５に示す。

温湿度計で実際に測定した値より２℃高い事前確率分布（Case 1）、及び温湿度計で実際に測定した値より２℃低い事前確率分布（Case 2）を用いて、屋内実験及び屋外実験を行った。なお、ヤコビアン（Jacobian matrix）計算のため、水蒸気濃度として温湿度計で測定した値が用いられた。また、無用なノイズを避けられるように、全チャンネルの内から相互情報量の多い１０チャンネルが、ＭＡＰ解x_MAPを計算するために用いられた。

屋内実験及び屋外実験におけるＭＡＰ解x_MAP、実測温度（温度計実測値）の分布及びＡＫareaを図２５〜２８に示す。図２５〜２８に示すように、屋内実験及び屋外実験において、気温の分布を再現することができた。なお、２層目以降では事前確率分布が顕著に影響した。また、図２７が示すように、１層目のＡＫareaは約０．８であるため、観測値がＭＡＰ解x_MAPに対して支配的に寄与した。一方で、図２７が示すように、４層目のＡＫareaは約０．２であるため、事前確率分布がＭＡＰ解x_MAPに対して支配的に寄与している。事前確率分布とＭＡＰ解x_MAPのＲＭＳＥを表６に示す。表６が示すように、屋内実験のＣａｓｅ２を除き、事前確率分布のＲＭＳＥより、ＭＡＰ解x_MAPのＲＭＳＥは、低減（改善）することができた。なお、事前確率分布のＲＭＳＥが２．０でない場合があるのは、位置を線形補間したためと考えられる。

以上の説明によって、本発明によれば、建築空間における観測対象となる気温について観測前に分かっている確率的な情報（事前確率分布）に基づいて、分光放射計を用いた水平気温分布を逆推定することが可能である。

なお、前述した本発明の実施形態は、何れも本発明を実施するにあたっての具体化の例（例えば、空気層の分割数、真の気温分布が一様であるか否か、境界面と観測機との距離、観測機の分解能、空気中の水蒸気の比重量等）を示したものに過ぎず、これらによって本発明の技術的範囲が限定的に解釈されてはならないものである。即ち、本発明は、その技術思想又はその主要な特徴から逸脱することなく、様々な形で実施することができる。

また、Ｐ（ｘ）、Ｐ（ｙ｜ｘ）については、本発明の実施形態においては、正規分布を仮定しているが、正規分布に限られたものではなく、例えば、ガンマ（γ）分布、ポアソン分布を用いても良い。ただしその場合はガウス・ニュートン法を用いることはできないので、マルコフ連鎖モンテカルロ法を用いる必要がある。

また、非線形最適化の手法として、ガウス・ニュートン法を挙げたが、例えばレーベンバーグ・マルカート法等を計算量と精度との兼ね合いに応じて用いてもよい。

１ 放射伝達計算部
２ ヤコビアン記憶部
３ 制御部
４ 入力部
５ 最適解処理部
６ 平均二乗誤差算出部
７ 事前確率分布記憶部
１０ 温度測定装置
１１ 熱源（黒体炉）
１２ 分光放射計
１３ 演算部
１３ａ ＣＰＵ
１３ｂ 表示部
１２０ 入射部
１２１ 集光部
１２２ ハーフミラー
１２３ ＣＣＤカメラ
１２４ コントローラ
１２５ チョッパ
１２６ 内部黒体
１２７ 連続可変フィルタ（ＣＶＦ）
１２８ 検出部

Claims (8)

1. 建築物又は都市空間より放出される赤外線を観測し、赤外線観測値を分析することにより、前記建築物又は前記都市空間における空気の温度分布を推定することを特徴とする温度測定方法。
2. 気温をｘ、前記赤外線観測値をｙ、前記温度分布に関する事前確率分布をＰ（ｘ）、前記赤外線観測値の観測誤差の確率密度関数をＰ（ｙ｜ｘ）、前記温度分布に関する事後分布をＰ（ｘ｜ｙ）とし、
   前記Ｐ（ｘ）及び前記Ｐ（ｙ｜ｘ）に基づいて、前記Ｐ（ｘ｜ｙ）を最大とする事後気温分布を推定する請求項１に記載の温度測定方法。
3. 前記Ｐ（ｘ）、前記Ｐ（ｙ｜ｘ）及び前記Ｐ（ｘ｜ｙ）の関係は、式（１）により表され、
   前記（１）式の右辺を最大とする前記事後気温分布を推定する請求項２に記載の温度測定方法。
4. ガウス・ニュートン法を用いて、前記事後気温分布を最大にするＭＡＰ解を算出する請求項３に記載の温度測定方法。
5. 建築物又は都市空間より放出される赤外線を観測する観測手段と、赤外線観測値に基づいて、前記建築物又は前記都市空間における空気の温度分布を推定する推定手段とを有することを特徴とする温度測定装置。
6. 前記推定手段は、
   気温をｘ、前記赤外線観測値をｙ、前記温度分布に関する事前確率分布をＰ（ｘ）、前記赤外線観測値の観測誤差の確率密度関数をＰ（ｙ｜ｘ）、前記温度分布に関する事後分布をＰ（ｘ｜ｙ）とし、
   前記Ｐ（ｘ）及び前記Ｐ（ｙ｜ｘ）に基づいて、前記Ｐ（ｘ｜ｙ）を最大とする事後気温分布を推定する請求項５に記載の温度測定装置。
7. 前記Ｐ（ｘ）、前記Ｐ（ｙ｜ｘ）及び前記Ｐ（ｘ｜ｙ）の関係は、式（１）により表され、
   前記推定手段は、前記（１）式の右辺を最大とする前記事後気温分布を推定する請求項６に記載の温度測定装置。
8. ガウス・ニュートン法を用いて、前記事後気温分布を最大にするＭＡＰ解を算出する請求項７に記載の温度測定装置。