(57)【特許請求の範囲】
【請求項１】
時間的かつ空間的な欠落を含まない画像を生成する画像生成システムであって、
プログラムを実行するプロセッサと、前記プロセッサが実行するプログラムを格納するメモリとを備え、
前記プロセッサは、
地上の同じ領域を撮影した分解能が異なる第１の画像と第２の画像とを用いて、前記第１の画像と前記第２の画像とを変換するためのスペクトル補正パラメータを算出し、
撮影時刻が異なる二つの第２の画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出し、
前記時間変化補正パラメータを用いて、画像の欠落を含む既存の第２の画像から任意の時刻の第２の画像を生成し、
前記スペクトル補正パラメータを用いて、前記任意の時刻の第２の画像から、当該任意の時刻の第１の画像であって、画像の欠落を含む画像を生成し、
前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記任意の時刻と異なる時刻に撮影された第２の画像から、前記任意の時刻の第１の画像のうち欠落している領域の画像を生成し、
前記生成された任意の時刻における画像の欠落を含む第１の画像と、前記生成された任意の時刻の第１の画像のうち欠落している領域の画像とを合成することによって、前記欠落がない任意の時刻の第１の画像を生成する画像生成システム。
【請求項２】
時間的かつ空間的な欠落を含まない画像を生成する画像生成システムであって、
プログラムを実行するプロセッサと、前記プロセッサが実行するプログラムを格納するメモリとを備え、
前記プロセッサは、
地上の同じ領域を撮影した分解能が異なる第１の画像と第２の画像とを、地上制御点を用いて位置を合わせ、
スペクトル混合モデルを用いて、前記第１の画像から、前記第２の画像と同じ分解能の画像を生成し、
前記位置合わせに使用された第２の画像のスペクトルと、前記スペクトル混合モデルを用いて生成された第２の画像のスペクトルとの差によって、前記第１の画像と前記第２の画像とを変換するためのスペクトル補正パラメータを算出し、
撮影時刻が異なる二つの第２の画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出し、
前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記第１の画像と撮影時刻が異なる第２の画像から、前記第１の画像のうち欠落している画像を生成する画像生成システム。
【請求項３】
請求項１又は２に記載の画像生成システムであって、
前記第１の画像は、高分解能画像であって、
前記第２の画像は、前記第１の画像より低い分解能の低分解能画像であって、
前記プロセッサは、
地上の同じ領域を撮影した前記高分解能画像と前記低分解能画像とを用いて、前記高分解能画像と前記低分解能画像とを変換するためのスペクトル補正パラメータを算出し、
撮影時刻が異なる二つの低分解能画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出し、
前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記高分解能画像と撮影時刻が異なる低分解能画像から、前記高分解能画像のうち欠落している画像を生成する画像生成システム。
【請求項４】
請求項１又は２に記載の画像生成システムであって、
前記プロセッサは、
撮影時刻が異なる第２の画像の画素を比較して時間変化パラメータを算出し、
前記第２の画像において、雲がかかっている領域の近くの画素を用いて、欠落している領域の画素を推定し、
前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記推定された画素を取り込んだ第２の画像から、前記第１の画像のうち欠落している画像を生成する画像生成システム。
【請求項５】
請求項４に記載の画像生成システムであって、
前記プロセッサは、
推定すべき画素の近傍に所定の領域を定め、
前記所定の領域に含まれる画素と、前記推定すべき画素との間の距離の逆数によって重み係数を定め、
前記所定の領域に含まれる画素のスペクトルと前記定められた重み係数とを乗じた値を加算することによって、前記雲がかかっている領域の画素を推定する画像生成システム。
【請求項６】
計算機を用いて、時間的かつ空間的な欠落を含まない画像を生成する画像生成方法であって、
前記計算機は、プログラムを実行するプロセッサと、前記プロセッサが実行するプログラムを格納するメモリとを有し、
前記方法は、
前記プロセッサが、地上の同じ領域を撮影した分解能が異なる第１の画像と第２の画像とを用いて、前記第１の画像と前記第２の画像とを変換するためのスペクトル補正パラメータを算出するステップと、
前記プロセッサが、撮影時刻が異なる二つの第２の画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出するステップと、
前記プロセッサが、前記時間変化補正パラメータを用いて、画像の欠落を含む既存の第２の画像から任意の時刻の第２の画像を生成するステップと、
前記プロセッサが、前記スペクトル補正パラメータを用いて、前記任意の時刻の第２の画像から、当該任意の時刻の第１の画像であって、画像の欠落を含む画像を生成するステップと、
前記プロセッサが、前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記任意の時刻と異なる時刻に撮影された第２の画像から、前記任意の時刻の第１の画像のうち欠落している領域の画像を生成するステップと、
前記プロセッサは、前記生成された任意の時刻における画像の欠落を含む第１の画像と、前記生成された任意の時刻の第１の画像のうち欠落している領域の画像とを合成することによって、前記欠落がない任意の時刻の第１の画像を生成するステップと、を含む画像生成方法。
【請求項７】
計算機を用いて、時間的かつ空間的な欠落を含まない画像を生成する画像生成方法であって、
前記計算機は、プログラムを実行するプロセッサと、前記プロセッサが実行するプログラムを格納するメモリとを有し、
前記方法は、
前記プロセッサが、地上の同じ領域を撮影した分解能が異なる第１の画像と第２の画像とを用いて、前記第１の画像と前記第２の画像とを変換するためのスペクトル補正パラメータを算出するステップと、
前記プロセッサが、撮影時刻が異なる二つの第２の画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出するステップと、
前記プロセッサが、前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記第１の画像と撮影時刻が異なる第２の画像から、前記第１の画像のうち欠落している画像を生成するステップと、を含み、
前記スペクトル補正パラメータを算出するステップでは、前記プロセッサは、
地上制御点を用いて、前記第１の画像と前記第２の画像との位置を合わせ、
スペクトル混合モデルを用いて、前記第１の画像から、前記第２の画像と同じ分解能の画像を生成し、
前記位置合わせに使用された第２の画像のスペクトルと、前記スペクトル混合モデルを用いて生成された第２の画像のスペクトルとの差によって、前記スペクトル補正パラメータを算出する画像生成方法。
【請求項８】
請求項６又は７に記載の画像生成方法であって、
前記第１の画像は、高分解能画像であって、
前記第２の画像は、前記第１の画像より低い分解能の低分解能画像であって、
前記スペクトル補正パラメータを算出するステップでは、前記プロセッサは、地上の同じ領域を撮影した前記高分解能画像と前記低分解能画像とを用いて、前記高分解能画像と前記低分解能画像とを変換するためのスペクトル補正パラメータを算出し、
前記時間変化補正パラメータを算出するステップでは、前記プロセッサは、撮影時刻が異なる二つの低分解能画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出し、
前記画像を生成するステップでは、前記プロセッサは、前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記高分解能画像と撮影時刻が異なる低分解能画像から、前記高分解能画像のうち欠落している画像を生成する画像生成方法。
【請求項９】
請求項６又は７に記載の画像生成方法であって、
前記画像を生成するステップでは、前記プロセッサは、
撮影時刻が異なる第２の画像の画素を比較して時間変化パラメータを算出し、
前記第２の画像において、雲がかかっている領域の近くの画素を用いて、欠落している領域の画素を推定し、
前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記推定された画素を取り込んだ第２の画像から、前記第１の画像のうち欠落している画像を生成する画像生成方法。
【請求項１０】
請求項９に記載の画像生成方法であって、
前記プロセッサは、推定すべき画素の近傍に所定の領域を定めるステップと、
前記プロセッサは、前記所定の領域に含まれる画素と、前記推定すべき画素との間の距離の逆数によって重み係数を定めるステップと、
前記プロセッサは、前記所定の領域に含まれる画素のスペクトルと前記定められた重み係数とを乗じた値を加算することによって、前記雲がかかっている領域の画素を推定するステップとを含む画像生成方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像を生成する画像生成システムに関する。
【背景技術】
【０００２】
人工衛星や航空機などで空中から撮影した画像において、雲がかかっている地域など遮蔽物で隠されている領域では地上の事物が撮影されない。特に、特定の期間に観測された情報を必要とするアプリケーションがある。例えば、農業分野では、作物が生長する期間の情報が生育管理や収量予測に重要である。
【０００３】
中〜高分解能画像は、画質が十分で、作物の詳細な情報を提供することから、農業分野で広く利用されている。しかし、中〜高分解能画像は、低い頻度（多くの場合、２週間以上）でしか同じ地域が観測されず、さらに、天候（例えば、雲）の影響を受ける。この観測時間を観測時間窓（acquisition window period）と称する。このため、有用な衛星画像を入手することは困難である。
【０００４】
また、衛星の再来周期は衛星が運用される軌道によって定まる。中〜高分解能画像を撮影する衛星は、中〜低高度軌道で運用される。しかし、技術的困難性から、衛星の再来頻度の増加は難しい。例えば、衛星の数を増やせば、衛星の再来頻度を増やすことができるが、コストが高く、実現が困難である。
【０００５】
一方、高々度軌道で運用される観測衛星の再来頻度は多い（例えば、１日）。また、高々度衛星は、詳細な空間情報を含まない低分解能画像を高い時間頻度で提供する。低分解能画像は、低コストで、自由に利用できる場合が多い。
【０００６】
雲の影響がない（又は、少ない）衛星画像は、非常に有用であるが、利用できる機会が少ない。このため、画像から雲を除去する技術が必要とされている。従来技術として以下の二つがある。一つ目は、雲で覆われている領域を、同じ領域の雲がない他の画像で置き換える方法である。しかし、雲で覆われている領域と雲がない画像とでは時間的な情報が一致しない。二つ目は、雲で覆われている領域を、同じ画像の周囲の領域の情報を用いて補完する、すなわち、他の領域の雲がない画像で置き換える方法である。しかし、雲で覆われている領域と雲がない画像とでは空間的な情報が一致しない。
【０００７】
この技術分野の背景技術として非特許文献１及び２がある。非特許文献１には、厚い雲によって画像が取得できないことは、地上観測衛星に共通する問題であることから、近隣の類似する画素による補間によって雲を除去する方法が開示される。また、非特許文献２には、スペクトルアンミキシングを適用することによって、雲で覆われた領域のサブピクセルから雲による影響を分離する方法が開示される。
【先行技術文献】
【非特許文献】
【０００８】
【非特許文献１】Yanfeng Gu 外５名、Multiple-kernel learning-based unmixing algorithm for estimation of cloud fractions with MODIS and CloudSat data、Geoscience and Remote Sensing Symposium (IGARSS), 2012 IEEE International、pp1785-1788、２０１２年７月
【非特許文献２】Xiaolin Zhu 外３名、A Modified Neighborhood Similar Pixel Interpolator Approach for Removing Thick Clouds in Landsat Images、Geoscience and Remote Sensing Letters, IEEE (Volume:9, Issue:3)、pp521-525、２０１２年５月
【発明の概要】
【発明が解決しようとする課題】
【０００９】
雲がかかっていることによって利用できない情報を回復するための鍵は、情報が豊富で信頼できる手がかりを見出すことである。前述した一つ目の方法では、異なる時刻の画像を利用するが、対象地域における時間的な変化が考慮されていない。また、二つ目の方法では、対象地域の情報は周囲の地域と類似しているが、対象地域の実際の情報は全く含まれていない。
【００１０】
このため、本発明は、時間的かつ空間的な手がかりを考慮して、欠落している情報を導出する方法を提案する。
【課題を解決するための手段】
【００１１】
本願において開示される発明の代表的な一例を示せば以下の通りである。すなわち、時間的かつ空間的な欠落を含まない画像を生成する画像生成システムであって、プログラムを実行するプロセッサと、前記プロセッサが実行するプログラムを格納するメモリとを備え、前記プロセッサは、地上の同じ領域を撮影した分解能が異なる第１の画像と第２の画像とを用いて、前記第１の画像と前記第２の画像とを変換するためのスペクトル補正パラメータを算出し、撮影時刻が異なる二つの第２の画像を用いて、時間的な変化を補正するための時間変化補正パラメータを算出し、前記時間変化補正パラメータを用いて、画像の欠落を含む既存の第２の画像から任意の時刻の第２の画像を生成し、前記スペクトル補正パラメータを用いて、前記任意の時刻の第２の画像から、当該任意の時刻の第１の画像であって、画像の欠落を含む画像を生成し、前記スペクトル補正パラメータ及び前記時間変化補正パラメータを用いて、前記任意の時刻と異なる時刻に撮影された第２の画像から、前記任意の時刻の第１の画像のうち欠落している領域の画像を生成し、前記生成された任意の時刻における画像の欠落を含む第１の画像と、前記生成された任意の時刻の第１の画像のうち欠落している領域の画像とを合成することによって、前記欠落がない任意の時刻の第１の画像を生成する。
【発明の効果】
【００１２】
本発明の代表的な実施の形態によれば、雲がかかった地域の情報を正確に導出することができる。前述した以外の課題、構成及び効果は、以下の実施例の説明により明らかにされる。
【図面の簡単な説明】
【００１３】
【図１】本発明の実施例の画像生成方法の概要を説明する図である。
【図２Ａ】本実施例の画像間で変換するための式を示す図である。
【図２Ｂ】本実施例の画像間で変換するための式を示す図である。
【図３】既存の雲がかかっている領域Ｂ１０（ｔ）を推定する方法のフローチャートである。
【図４】スペクトル補正パラメータｄ１を算出する処理（Ｓ１）の詳細のフローチャートである。
【図５】時間変化補正パラメータｄ２を算出する処理（Ｓ２）の詳細のフローチャートである。
【図６】ステップＳ２３の処理を説明する図である。
【図７】任意の時刻の衛星画像Ａ１０（ｔｊ）を推定する方法のフローチャートである。
【図８】画像Ａ１０（ｔｊ）を作成する処理（Ｓ５）の詳細のフローチャートである。
【図９】本実施例の画像生成システムの論理的な構成を示すブロック図である。
【図１０】本実施例の画像生成システムの物理的な構成を示すブロック図である。
【発明を実施するための形態】
【００１４】
まず、本発明の実施例の概要を説明する。
【００１５】
本実施例は、遮蔽物で覆われた領域を含まない画像（例えば、雲がかかっていない衛星画像）を生成する計算機システムに関し、時系列の他の衛星画像から、衛星画像中で雲がかかっている領域で欠落している地上の事物からの反射光の情報を回復する。特に、重要な観測タイミングにおいて、雲で覆われた影響を除去することができる。このため、衛星画像の有用性を向上する。
【００１６】
本実施例の主要な点は、雲がかかっていない領域を用いて時間的変化を導出し、空間的類似性に従って雲がかかっている領域に適用する。具体的には、画像の時間的な変化を求めるために時系列の衛星画像を用いて、求められた時間的な変化を処理対象の領域に適用し、雲がかかっていない低分解能画像（雲がかかっている高分解能画像と同じ時刻で）を合成する。そして、低分解能画像を高分解能画像に投影し、雲がかかっている領域で必要な情報を生成する。
【００１７】
図１は、本発明の実施例の画像生成方法の概要を説明する図である。
【００１８】
本実施例では、高分解能の衛星画像１及び低分解能の衛星画像２を用いる。前述したように、高分解能画像は、撮影コストが高く、低頻度にしか撮影されない画像であり、例えば、低高度衛星や航空機が１月に１回程度撮影する画像である。一方、低分解能画像は、撮影コストが低く、高頻度に撮影されている画像であり、高高度衛星が１日に数回撮影する画像である。
【００１９】
本実施例について人工衛星が撮影した衛星画像を用いて説明するが、本発明は航空機が撮影した画像にも適用できる。また、衛星の一時的な故障によって撮影できなかった領域の画像を回復するためにも使用することができる。
【００２０】
あるタイミングｔにおいて撮影された衛星画像１及び衛星画像２では、領域Ｂに雲がかかっており地上の事物（例えば、建物、草原の植生など）を撮影することができなかった。また、タイミングｔｉにおいて撮影された衛星画像２では、領域Ｂには雲がかかっておらず地上の事物を撮影することができた。さらに、タイミングｔｉは衛星画像１の撮影周期ではなく、高分解能画像は撮影されていない。
【００２１】
本発明は、このような状況において、画像列１の衛星画像Ａ１０（ｔ）から画像列２の衛星画像Ａ２０（ｔ）へのスペクトルの変化ｄ１と、画像列２の衛星画像Ａ２０（ｔｉ）からＡ２０（ｔ）への時間的変化ｄ２とを用いて、シミュレーションによって時刻ｔにおける雲を除去した画像Ｂ１０（ｔ）を生成することができる。
【００２２】
さらに、本発明の方法を用いると、シミュレーションによって得られた任意の時刻ｔｊの画像列２の衛星画像Ａ２０（ｔｊ）を用いて、画像列１の衛星画像Ａ１０（ｔ）からＡ２０（ｔ）へのスペクトルの変化ｄ１と、画像列２の衛星画像Ａ２０（ｔ）からＡ２０（ｔｊ）への時間的変化ｄ３とを用いて、時刻ｔｊにおける高分解能画像Ａ１０（ｔｊ）をシミュレーションによって生成することができる。
【００２３】
すなわち、本発明の方法を用いることによって、将来も含む任意の時刻の衛星画像を精度よく求めることができる。
【００２４】
図２Ａは、本実施例の画像間で変換するための式を示す図である。
【００２５】
数式（１）は、画像列１の衛星画像Ａ１と画像列２の衛星画像Ａ２とがスペクトル補正パラメータｄ１で変換できることを示す。同じ時刻ｔで撮影された衛星画像列１と衛星画像２とは、同じ地上の事物を異なるセンサによって異なる分解能で撮影したものである。このため、各画像を撮影するときに使用されたセンサの周波数特性や感度の違いによって、撮影される画像が異なる。このため、このセンサによる画像の違いをスペクトル補正パラメータｄ１によって変換する。
【００２６】
数式（２）は、画像列２の衛星画像Ａ２（ｔ）とのＡ２（ｔｉ）とが時間変化補正パラメータｄ２で変換できることを示す。異なる時刻ｔとｔｉとで撮影された衛星画像２は、時間的変化ｄ２によって変換することができる。画像列１の二つの衛星画像を用いて時間変化補正パラメータｄ２を算出してもよい。なお、時間的変化ｄ２は画像が撮影された時間差ｔｉ−ｔの関数である。
【００２７】
数式（３）は、数式（１）及び（２）を用いて、異なる時刻に撮影された衛星画像列１と衛星画像２とを変換するための数式である。すなわち、画像列２の衛星画像Ａ２（ｔ）とのＡ２（ｔｉ）とが差分ｄ２で変換でき、画像列１の衛星画像Ａ１と画像列２の衛星画像Ａ２とがスペクトル差分ｄ１で変換できる。このため、画像列２の衛星画像Ａ２（ｔｉ）と画像列１の衛星画像Ａ１（ｔ）とがスペクトル差分ｄ１及び時間変化補正パラメータｄ２で変換できる。
【００２８】
図２Ｂは、本実施例で分解能が異なる画像間での画素の変換方法を説明する図である。
【００２９】
衛星画像１の分解能Ｒ１が衛星画像２の分解能Ｒ２より小さい場合、及び、衛星画像１の分解能Ｒ１と衛星画像２の分解能Ｒ２とが同じ場合、衛星画像２の画素から衛星画像１の画素を決めることができる。
【００３０】
一方、衛星画像１の分解能Ｒ１が衛星画像２の分解能Ｒ２より大きい場合、画像２の画素は画像１の画素が混合したものとなる。このため、衛星画像２の画素から、衛星画像１の画素を推定する。例えば、高分解能の衛星画像１の複数の画素の特徴量と低分解能の衛星画像２の一つの画素の特徴量とを対応付けたデータを作成し、低分解能の衛星画像２の画素から高分解能の衛星画像１の画素を推定してもよい。
【００３１】
図３は、既存の雲がかかっている領域Ｂ１０（ｔ）を推定する方法のフローチャートである。
【００３２】
まず、時刻ｔの画像列１の衛星画像Ａ１と画像列２の衛星画像Ａ２との雲がかかっていない領域を用いてスペクトル補正パラメータｄ１を算出する（Ｓ１）。ステップＳ１の処理の詳細は図４を用いて後述する。
【００３３】
次に、撮影時刻が異なる二つの衛星画像Ａ２を用いて時間変化補正パラメータｄ２を算出する（Ｓ２）。ステップＳ２において、関数ｆ２は地上の事物の種類毎の時間変化モデルである。ステップＳ２の処理の詳細は図５を用いて後述する。なお、画像列１の二つの衛星画像Ａ１を用いて時間変化補正パラメータｄ２を算出してもよい。
【００３４】
最後に、補正パラメータｄ１及びｄ２を用いて、異なる時刻及び異なる分解能の雲なし画像Ｂ２０（ｔｉ）から、雲がかかっている領域の画素を算出し、雲なし画像Ｂ１０（ｔ）を作成する（Ｓ３）。
【００３５】
図４は、スペクトル補正パラメータｄ１を算出する処理（Ｓ１）の詳細のフローチャートである。
【００３６】
まず、ステップＳ１では、時刻ｔの画像Ａ１（ｔ）と画像Ａ２（ｔ）の空間的な位置を合わせる（Ｓ１１）。これは、地上に設けられたＧＣＰ（クロスポイント、エッジポイントなどの地上制御点）を二つの画像上で特定し、特定された地上制御点の位置を合わせることによって、二つの画像の位置を合わせる。
【００３７】
次に、画像Ａ１（ｔ）及びＡ２（ｔ）の大気補正を行う（Ｓ１２）。大気補正によって、地上観測機器から得られた大気の情報を用いて放射伝達を計算し、大気による散乱や吸収の影響を除去することができる。
【００３８】
その後、スペクトル混合モデルを用いて、画像Ａ１（ｔ）の画素を、画像Ａ２（ｔ）と同じ分解能の画像Ａ２（ｔ）’を作成する（Ｓ１３）。具体的には、ＬＤＡを用いて、高分解能画像の分布モデルと低分解能画像の分布モデルとを生成し、低分解能画像の１画素の特徴量と高分解能画像の複数画素の特徴量とを対応させた特徴関連モデルを構築する。そして、構築された特徴関連モデルを用いて、低分解能画像の画素から高分解能画像の画素との関係を推定する。
【００３９】
そして、画像Ａ２（ｔ）の画素のスペクトルと生成された画像Ａ２（ｔ）’の画素のスペクトルとの差によって、スペクトル補正パラメータｄ１を算出する（Ｓ１４）。スペクトル補正パラメータｄ１は、低分解能画像Ａ２のスペクトル毎に算出される。
【００４０】
図５は、時間変化補正パラメータｄ２を算出する処理（Ｓ２）の詳細のフローチャートである。
【００４１】
まず、ステップＳ２では、時刻ｔｉの雲なしの衛星画像列２の画像Ａ２０（ｔｉ）を取得する（Ｓ２１）。
【００４２】
次に、雲がかかっていない領域で、時刻ｔｉと時刻ｔの画素から時間変化パラメータを算出する（Ｓ２２）。具体的には、画像Ａ２０（ｔ）の画素のスペクトルと画像Ａ２０（ｔｉ）の画素のスペクトルとの差を算出する。これによって、画素のスペクトル毎にスペクトルの時間変化モデルを算出することができる。これは、同じスペクトルの画素は同じ地上の事物を撮影していると推定され、同じ地上の事物は時間の経過によってスペクトルが同じ変化をすると推定されるからである。数式において、ｆ２は、画素のスペクトルから時間変化補正パラメータｄ２を求める関数である。なお、ｄ２は、画像Ａ２０（ｔｉ）の画素のスペクトルによって定まるパラメータである。
【００４３】
雲がかかっている領域の近くの画素を用いて、雲がかかっている領域の画素Ｐｉを推定する（Ｓ２３）。
【００４４】
図６は、ステップＳ２３の処理を説明する図である。
【００４５】
推定すべき画素Ｐｉを中心とした所定の大きさの円を作成する。作成された円には、ｎ個の画素が含まれる。また、画素Ｐから画像Ａ２０（ｔ）のｊ番目の画素Ａ２０（ｊ）までの距離の逆数を重み係数ｋｊとする。ｊは、作成された円に含まれる画素を制御するためのパラメータである。そして、作成された円に含まれる画素Ａ２０（ｊ）のスペクトルに距離の逆数の重み係数を乗じて、それらの和を算出する。この和を画素数ｎで除することによって、推定すべき画素Ｐｉのスペクトルを算出することができる。
【００４６】
このように、近くの画素が表す地上の事物は類似している可能性が高いことを利用して、すなわち、画素間の距離の逆数ｋｊによって重み付けをして、雲がかかっている領域の画素を空間的類似性を用いて推定することができる。
【００４７】
なお、雲がかかった領域の地上の事物が、地上の観測などによって分かっている場合、事物の種類毎に定められた時間変化モデルｆ２を用いて、画素毎に変化を算出してもよい。例えば、画素のスペクトルから事物が小麦であることが分かっている場合、小麦の生育状態とスペクトルの変化の関係を示す生育モデルを用いて、雲がかかっている領域の画素を推定する。
【００４８】
次に、雲がかかっている領域における時間変化を補正するためのパラメータｄβを、時間変化補正パラメータｄ２及び推定された画素Ｐｉを用いて算出する（Ｓ２４）。具体的には、全ての画素について、時間変化補正パラメータｄ２と画素Ｐｉのスペクトルを乗じて、それらの和を算出する。すなわち、関数ｆｄは、全てのｉについて、Ｐｉ×ｄ２の和を算出する関数である。
【００４９】
これによって、ｄ２、すなわち異なる時刻に撮影された二つの衛星画像Ａ２０（ｔ）とＡ２０（ｔｉ）とをパラメータとした時間変化補正パラメータｄβ（ｄ２）を算出することができる。
【００５０】
以上、本発明の画像生成方法について、遮蔽物の例を雲として、雲がかかっている領域の推定を例にして説明したが、本発明は、雲以外の遮蔽物の影響によって撮影されていない領域を含む画像にも適用することができる。
【００５１】
次に、任意の時刻の衛星画像を推定する方法について説明する。
【００５２】
図７は、任意の時刻の衛星画像Ａ１０（ｔｊ）を推定する方法のフローチャートである。
【００５３】
まず、画像列１の時刻ｔの衛星画像Ａ１と画像列２の時刻ｔの衛星画像Ａ２との雲がかかっていない領域を用いてスペクトル補正パラメータｄ１を算出する（Ｓ１）。ステップＳ１の処理の詳細は図４を用いて説明した。
【００５４】
次に、所望の時刻ｔｊの画像列２の画像Ａ２０（ｔｊ）を推定する（Ｓ４）。この推定には、既存の方法を用いることができる。例えば、図３のステップＳ２で算出された時間変化補正パラメータｄ２を、時刻差ｔｊ−ｔに適用して、画像Ａ２０（ｔ）から画像Ａ２０（ｔｊ）を算出することができる。
【００５５】
その後、撮影時刻が異なる二つの衛星画像Ａ２を用いて時間変化補正パラメータｄ３を算出する（Ｓ５）。なお、ステップＳ５の処理は、前述したステップＳ２の処理（図３、図５）と同じであり、関数ｆ２は地上の事物の種類毎の時間変化モデルである。
【００５６】
なお、時間変化補正パラメータｄ３は、画像列２の既存の二つの画像Ａ２から求め、さらに時刻差によって変化したものでもよい。例えば、図３のステップＳ２で算出された時間変化補正パラメータｄ２を、時刻差ｔｊ−ｔに適用して、時間変化補正パラメータｄ３を算出してもよい。
【００５７】
最後に、補正パラメータｄ３を用いて、時刻ｔの画像Ａ１０（ｔ）から時刻ｔｊの画像の画素を算出し、画像Ａ１０（ｔｊ）を作成する。また、補正パラメータｄ１及びｄ３を用いて、時刻ｔの画像Ｂ２０（ｔｊ）から、雲がかかっている領域の画素を算出し、雲なし画像Ｂ１０（ｔｊ）を作成する（Ｓ６）。
【００５８】
図８は、画像Ａ１０（ｔｊ）を作成する処理（Ｓ６）の詳細のフローチャートである。
【００５９】
まず、画像列１において、時間変化補正パラメータｄ２を用いて、画像Ａ１０（ｔ）から、時刻ｔｊの画像Ａ１０（ｔｊ）を作成する（Ｓ６１）。しかし、この段階で作成された画像Ａ１０（ｔｊ）には、元の画像Ａ１０（ｔ）と同様に、雲がかかっている領域を含む。
【００６０】
次に、補正パラメータｄ１及びｄ３を用いて、雲なし画像Ｂ２０（ｔｉ）から、雲がかかっている領域の画素を算出し、異なる分解能の雲なし画像Ｂ１０（ｔｊ）を作成する（Ｓ６２）。
【００６１】
その後、作成された雲なし画像Ｂ１０（ｔｊ）を画像Ａ１０（ｔｊ）に取り込むことによって、雲がかかっていない画像Ａ１０（ｔｊ）を作成する（Ｓ６３）。
【００６２】
図９は、本実施例の画像生成システムの論理的な構成を示すブロック図である。
【００６３】
本実施例の画像生成システムは、入力部１１０、表示部１２０、データ記録部１３０、演算部１４０、データ出力部１５０及び記憶部１６０を有する計算機である。入力部１１０、表示部１２０、データ記録部１３０、データ出力部１５０及び記憶部１６０は、演算部１４０を介して（又は、相互にバスによって）接続される。
【００６４】
入力部１１０は、画像入力部１１１及び地図情報入力部１１２を有する。画像入力部１１１は、高分解能画像及び低分解能画像が入力される装置であり、地図情報入力部１１２は、地図情報が入力される装置である。画像入力部１１１及び地図情報入力部１１２は、例えば、光ディスクドライブ、ＵＳＢインターフェース等のデータの入力を受け付けるデバイスや、キーボード、タッチパネル、マウスなどのヒューマンインターフェースによって構成される。なお、画像入力部１１１と地図情報入力部１１２とは、同じ入力装置によって構成しても、異なる入力装置によって構成してもよい。
【００６５】
表示部１２０は、画像表示部１２１及び地図情報表示部１２２を有する。画像表示部１２１は、処理される画像を表示する表示装置である。地図情報表示部１２２は、入力された地図情報を表示する表示装置であり、地図中のクロスポイント、エッジポイントなどの地上基準点（ＧＣＰ）を指定するために用いられる。なお、画像表示部１２１と地図情報表示部１２２とは、同じ表示装置によって構成しても、異なる表示装置によって構成してもよい。
【００６６】
データ記録部１３０は、この画像生成システムによって処理される画像データを格納する不揮発性記憶装置であり、例えば、ハードディスク装置、不揮発性メモリ等によって構成される。演算部１４０は、プロセッサを含み、プログラムを実行することによって本画像生成システムで行われる処理を実行する。
【００６７】
データ出力部１５０は、この画像生成システムによって処理された結果を出力する装置で、例えば、プリンタ、プロッタ等によって構成される。記憶部１６０は、演算部１４０によって実行されるプログラムを格納する記憶装置であり、例えば、ハードディスク装置、不揮発性メモリ等によって構成される。
【００６８】
図１０は、本実施例の画像生成システムの物理的な構成を示すブロック図である。
【００６９】
本実施例の画像生成システムは、演算部１０、記憶装置２０、通信インターフェース３０及び媒体ドライバ４０を有する。
【００７０】
演算部１０は、プログラムを実行するプロセッサ（ＣＰＵ）１０１、不揮発性の記憶素子であるＲＯＭ１０２及び揮発性の記憶素子であるＲＡＭ１０３を有する。ＲＯＭ１０２は、不変のプログラム（例えば、ＢＩＯＳ）などを格納する。ＲＡＭ１０３は、記憶装置２０に格納されたプログラム及びプログラムの実行時に使用されるデータを一時的に格納する。
【００７１】
演算部１００が実行するプログラムは、リムーバブルメディア（ＣＤ−ＲＯＭ、フラッシュメモリなど）又はネットワークを介して計算機に提供され、非一時的記憶媒体である記憶装置に格納される。このため、画像生成システムを構成する計算機は、リムーバブルメディアからデータを読み込むインターフェースを有するとよい。
【００７２】
記憶装置２０は、例えば、磁気記憶装置、フラッシュメモリ等の大容量かつ不揮発性の記憶装置であり、プロセッサ１０１によって実行されるプログラム及びプログラムの実行時に使用されるデータを格納する。すなわち、プロセッサ１０１によって実行されるプログラムは、記憶装置２０から読み出されて、ＲＡＭ１０３にロードされて、プロセッサ１０１によって実行される。
【００７３】
通信インターフェース３０は、所定のプロトコルに従って、他の装置との通信を制御するネットワークインターフェース装置である。媒体ドライバ４０は、画像生成システムに導入されるプログラムやデータが格納された記録媒体５０を読むためのインターフェース（例えば、光ディスクドライブ、ＵＳＢポート）である。
【００７４】
本実施例の画像生成システムは、物理的に一つの計算機上で、又は、論理的又は物理的に構成された複数の計算機上で構成される計算機システムであり、同一の計算機上で別個のスレッドで動作してもよく、複数の物理的計算機資源上に構築された仮想計算機上で動作してもよい。例えば、空中写真を提供するセンシング事業者が、この画像生成システムをクラウド環境で提供してもよい。
【００７５】
本実施例の画像生成システムは、スタンドアロン形式のコンピュータに実装しても、クライアント・サーバ型の計算機システムに実装してもよい。クライアント・サーバ型の計算機システムに実装する場合、サーバが演算処理を実行し、クライアントが入力データを受け付け、演算結果を出力する。
【００７６】
以上に説明したように、本発明の実施例によると、正確かつ高信頼性で雲などの遮蔽物で覆われた画像を回復することができる。特に、空間的変化及び時間的変化に多くの情報を取り込んで、雲で覆われた領域の画像を高精度で回復することができる。また、高分解能画像中の雲の影響を低減し、短時間間隔で高分解能画像を得ることができるので、高分解能画像の有用性を向上することができ、高分解能画像の取得コストを低減することができる。
【００７７】
また、雲の影響を低減するだけでなく、将来も含む任意の時刻の衛星画像を精度よく求めることができる。
【００７８】
このようにして、高効率かつ低コストで高分解能画像を生成する画像生成システムを構築することができる。
【００７９】
なお、本発明は前述した実施例に限定されるものではなく、添付した特許請求の範囲の趣旨内における様々な変形例及び同等の構成が含まれる。例えば、前述した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに本発明は限定されない。また、ある実施例の構成の一部を他の実施例の構成に置き換えてもよい。また、ある実施例の構成に他の実施例の構成を加えてもよい。また、各実施例の構成の一部について、他の構成の追加・削除・置換をしてもよい。
【００８０】
また、前述した各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等により、ハードウェアで実現してもよく、プロセッサがそれぞれの機能を実現するプログラムを解釈し実行することにより、ソフトウェアで実現してもよい。
【００８１】
各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリ、ハードディスク、ＳＳＤ（Solid State Drive）等の記憶装置、又は、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に格納することができる。
【００８２】
また、制御線や情報線は説明上必要と考えられるものを示しており、実装上必要な全ての制御線や情報線を示しているとは限らない。実際には、ほとんど全ての構成が相互に接続されていると考えてよい。
【符号の説明】
【００８３】
１０  演算部
２０  記憶装置
３０  通信インターフェース
４０  媒体ドライバ
５０  記録媒体
１００  演算部
１０１  プロセッサ
１０２  ＲＯＭ
１０３  ＲＡＭ
１１０  入力部
１１１  画像入力部
１１２  地図情報入力部
１２０  表示部
１２１  画像表示部
１２２  地図情報表示部
１３０  データ記録部
１４０  演算部
１５０  データ出力部
１６０  記憶部
【図１】
0006351496.tif 000002
【図２Ａ】
0006351496.tif 000003
【図２Ｂ】
0006351496.tif 000004
【図３】
0006351496.tif 000005
【図４】
0006351496.tif 000006
【図５】
0006351496.tif 000007
【図６】
0006351496.tif 000008
【図７】
0006351496.tif 000009
【図８】
0006351496.tif 000010
【図９】
0006351496.tif 000011
【図１０】
0006351496.tif 000012