JP2003323611A6

JP2003323611A6 - 衛星撮影画像のオルソ補正処理方法

Info

Publication number: JP2003323611A6
Application number: JP2003069813A
Authority: JP
Inventors: 文伸古村; 浩史上田; 隆司樋野; 高敏小平
Original assignee: Hitachi Software Engineering Co Ltd
Current assignee: Hitachi Software Engineering Co Ltd
Priority date: 2002-04-24
Filing date: 2003-03-14
Publication date: 2005-04-14
Anticipated expiration: 2023-03-14

Abstract

【課題】地上基準点がなくても適用でき、画像全体で誤差の発散がなく、計算量が過大とならず、地表面の起伏に起因する倒れ込み歪を補正可能とする。
【解決手段】基準回転楕円体地球モデル２３の海抜標高ゼロの基準面２６ａに対して、複数の異なる標高の基準面２６ｂ〜２６ｄを設定し、補正画像５を等間隔の格子でブロック１，２，３，……に分割し、これらブロック毎に、かつ各基準面２６ａ〜２６ｄ毎に、バイリニア関数による歪みのモデル式を求める。さらに、補正画像５の各画素位置Ａ’について、該当するブロックの標高データDEMの値を直接挟む上下の標高の２つの基準面のモデル式を用いて、夫々の標高値に対応する観測画像の画素位置を計算し、標高データDEMの値により、該２つの画素位置を線形補間して、所望の観測画像の画素位置を求める。
【選択図】図１

Description

【０００１】
【発明の属する技術分野】
本発明は、衛星から撮影した地上画像の歪みを補正する画像処理方法に係り、特に、地表の起伏に起因する歪みを補正するオルソ補正処理方法に関する。
【０００２】
【従来の技術】
始めに、衛星から撮影した画像（衛星画像）とその補正の原理について説明する。
【０００３】
図６は衛星による画像撮影の説明図であって、１は衛星、２は地表、３は地表観測線である。
【０００４】
同図において、衛星１は、下方に向けて搭載したセンサ（図示せず）で地表２からの電磁波を検知することにより、地表２を撮影する。かかるセンサが衛星用センサとして一般的なＣＣＤなどのラインセンサである場合、かかるセンサは、それを構成する１列のセル（センサ素子）の配列方向が衛星１の進行方向Ｌに直交するようにして、配置されており、かかるセンサによって同時に観測される地表２上の範囲が地表観測線３である。かかる地表観測線３は衛星１とともにその進行方向Ｌに移動し、この移動とともに各セルの検知データが所定のサンプリング周期で取り込まれる。かかる検知データの取込みタイミングは全セル同時である。この地表観測線３全体のデータが衛星画像の１ライン分に相当するものであり、従って、１回のデータ取込み（１サンプリング）で撮影画像の１ラインに相当するデータを取得することになる。このようにして、衛星１の軌道上の進行に伴い、地表観測線３による一定の幅（Swath）で地表２が連続観測される。
【０００５】
図７はこのようにして得られた衛星画像の補正の原理を示す図であり、４が観測画像、５が補正画像である。
【０００６】
同図において、以上の衛星による撮影画像は、観測画像４として、地上に送られる。かかる観測画像４は、衛星の軌道や姿勢の揺らぎ，センサの光学系の歪み，地球表面の形状などの要因により、その画像の形状に歪みが生じている。このため、衛星による撮影画像の利用のためには、この観測画像４の歪みを補正し、所定の地図投影法に従った補正画像５に変換する必要がある。ここで、かかる補正は、観測画像４の画素Ａ（ｐ，ｌ）からこれに対応する補正画像５の画素Ａ’（ｘ，ｙ）への写像を行なうものであるが、これを写像Ｆという。
【０００７】
なお、観測画像４は、衛星１でのセンサのセルの配列方向（即ち、地表観測線３の長手方向）をｐ座標、衛星１の進行方向Ｌをｌ座標とするｐ−ｌ座標系で表わされ、また、補正画像５は、観測画像４のｐ座標に対応する方向をｘ座標、観測画像４のｌ座標に対応する方向をｙ座標とするｘ−ｙ座標系で表わされる。
【０００８】
図８はかかる観測画像４の補正処理を示すフロ一チャートである。
【０００９】
図７に示す補正画像５では、画素が等間隔に配列・設定されているものであり、これに対し、観測画像４では、上記のような歪みのため、補正画像５の各画素に対応する画素位置は等間隔に配列されない。図８に示す補正処理は、補正画像５における各画素位置毎に、その画素位置に相当する観測画像４での画素位置の画素を見つけ出し、これを補正画像５での該当する画素位置に位置付けるものである。
【００１０】
図８において、かかる補正は、歪みのモデリング１００とリサンプリング１０１とからなっている。歪みのモデリング１００では、補正画像５における各画素位置毎に、その画素位置に相当する観測画像４での画素位置を求めるものである。図７についてみると、予め作成された歪モデル６を用いて、補正画像５での画素Ａ’(ｘ，ｙ）の位置に対応する観測画像４の画素Ａ（ｐ，ｌ）の位置を求めるものである。そのためには、写像Ｆの逆写像Ｆ-1を求める必要がある。
【００１１】
ここで、歪モデル６は、かかる観測画像４の撮影時での衛星の軌道や姿勢の揺らぎ，センサの光学系の歪みなどから予測されるされるもの、あるいは衛星の観測結果などから得られるものである。
【００１２】
また、リサンプリング１０１では、歪みのモデリング１００で補正画像５での画素Ａ’(ｘ，ｙ）の位置に対応する観測画像４の画素Ａ（ｐ，ｌ）の位置が求まると、この画素Ａ（ｐ，ｌ）の位置での画素の強度を求め、それを補正画像５の画素Ａ’（ｘ，ｙ）での画素の強度とするものである。
【００１３】
図９は図８におけるリサンプリング１０１の説明図であって、図７に示した観測画像４と補正画像５とでもって説明する。なお、７は観測画像４での画素の中心位置（以下、画素位置という）である。
【００１４】
同図において、図８で説明した歪みのモデリング１００により、補正画像５の画素Ａ’（ｘ，ｙ）の位置に対応する観測画像４の画素Ａ（ｐ，ｌ）の位置が求まるのであるが、画素Ａ（ｐ，ｌ）の位置は、図示するように、一般にセンサによる観測データの画素位置７に一致しない。そこで、この画素Ａ（ｐ，ｌ）の位置における観測データの強度を、この画素Ａ（ｐ，ｌ）の位置の周辺の実際に得られた画素位置７の画素強度から補間によって求める必要がある。このための補間方法としては、Nearest Neighboｒ，Bilinear，Cubic Convolutionなどの補間式が一般に用いられる。
【００１５】
ところで、従来、上記の原理に基づく衛星観測画像データの第１の補正方法としては、オルソ補正の方法が知られている（例えば、非特許文献１参照）。
【００１６】
これは、図８の歪みのモデリング１００を次の数１で表わすものである。
【００１７】
【数１】

この式はRational polynomial methodと呼ばれるものであって、撮影した画像毎にその都度、ｘ，ｙ，ｚの多項式の係数を、地上基準点GCP（Ground ControlPoint）を用いることにより、以下に説明するように計算して求め、得られた係数を用いた数１により、補正画像５の各画素位置に対する観測画像４の位置を求めるものである。
【００１８】
図１０はかかる地上基準点を用いた第１の補正方法の説明図である。
【００１９】
同図において、観測画像４において、複数の地上基準点Ｇ_i(ｉ＝１，２，……）、例えば、３個の地上基準点Ｇ₁，Ｇ₂，Ｇ₃を選定する。かかる地上基準点は、観測画像４で形状や位置が鮮明に識別できるものであり、かつこれら地上基準点の位置に対応する補正画像５での位置座標が、地図からの読み取りにより、あるいは測量により、正確に得られるものである。補正画像５では、これら地上基準点Ｇ₁，Ｇ₂，Ｇ₃に対応する位置を夫々Ｇ₁’，Ｇ₂’，Ｇ₃’として示している。
【００２０】
このようにして、観測画像４での複数の地上基準点Ｇ_i（ｐ_i，ｌ_i）とこれらに対応する補正画像５での位置Ｇ_i’（ｘ_i，ｙ_i，ｚ_i）とが求まると、これら地上基準点Ｇ_i（ｐ_i，ｌ_i）の画素座標値と位置Ｇ_i’(ｘ_i，ｙ_i，ｚ_i）の地図座標値とから、数１で示す多項式の未知の係数を最小自乗法で求めることができ、これにより、逆写像Ｆ^-1が可能となる。なお、多項式の次数は、地上基準点ＧＣＰの数や歪みの状況に応じて適宜設定する。
【００２１】
この第１の方法によると、観測対象地域の標高データDEM（Digital ElevationModel）と地上基準点とがあれば、衛星の軌道データ，姿勢データ或いは衛星のセンサの光学モデルがなくても、歪みのモデリング及びそれを用いた補正処理が可能である、という利点がある。
【００２２】
また、従来の衛星観測画像データの第２の補正方法としては、非オルソ補正、即ち、地表の標高による起伏を考慮しない方法として、次のような方法も知られている。これは、衛星の軌道データや姿勢データ，衛星のセンサの光学モデルを用いた幾何モデルによる方法である。この方法は、補正画像を等間隔の格子でもってブロックに分け、ブロック毎に歪みのモデル式を求めることに特徴がある。
【００２３】
図１１はこの第２の補正方法での格子によるブロック分けを説明するための図である。
【００２４】
同図において、補正画像５において、縦横に等間隔に格子８を設定し、かかる格子８により、補正画像５を複数のブロック９に区分する。補正画像５での画素位置は、かかるブロック９により、ブロック分けされている。そして、これら格子８が交差する各格子点１０に対応する観測画像４の位置を逆写像Ｆ^-1により求め、得られたこれらの位置を格子点１２とする格子１３を設定し、かかる格子１３により、観測画像４もブロック１１によって区分されることになる。一例として、観測画像４での格子点Ａが補正画像５での格子点Ａ’の逆写像Ｆ^-1として求められたものとし、これら格子点Ａ’，Ａが対応するものであることを示している。
【００２５】
そして、観測画像４のブロック１１は補正画像５のいずれかのブロック９と対応することになり、観測画像４と補正画像５とでブロック１１，９間が一対一に対応することになる。なお、ブロック９，１１は、これらの中で歪みが近似的に均一となるように、設定されるものであり、隣り合うブロック間では、歪みが異なっていてもよい。但し、隣り合うブロックの境界で歪みが連続するように、逆写像Ｆ^-1が設定されるものである。
【００２６】
ここで、逆写像Ｆ^-1として、次の数２で示すバイリニア関数を用いる。
【００２７】
【数２】

ここで、数２では、ｘ，ｙの０次，１次の項とｘｙの項とからなり、これにより、隣り合うブロックの境界で歪みが連続するようにしている。なお、ｘ²，ｙ²以上の高次の項を含ませると、隣り合うブロックの境界で歪みの不連続性（段差）が生ずることになり、これをなくすためには、多くの高次の項を付加して画面全体を１つのブロックとみなして計算する必要があり、計算が非常に複雑で時間がかかるものとなる。ここでは、歪みが均一とみてよい程度の大きさのブロックに画像４，５を分割し、このブロック毎にバイリニア関数を求め、これによって逆写像Ｆ^-1を行なうものであるから、数２のような簡単な式でもって逆写像Ｆ^-1を行なうことができるのである。この場合、上記のように、ｘｙの項を付加することにより、隣り合うブロックの境界で歪みが連続するようにしているのである。
【００２８】
図１２は上記のバイリニア関数での係数ａ₀，ａ₁，ａ₂，ａ₃，ｂ₀，ｂ₁，ｂ₂，ｂ₃の算出方法を示すフローチャートである。
【００２９】
同図において、この計算は、格子点計算２００と係数算出２０１とから構成されている。
【００３０】
図１１を例にして説明すると、格子点計算２００は、補正画像５の各格子点Ａ’（ｘ，ｙ）に対する観測画像４での座標位置Ａ（ｐ，ｌ）を求めるものである。これには、写像Ｆ^-1が必要となるが、これは、一般に、地上基準点なしには直接得られないので、写像Ｆの幾何学計算の繰り返し収束演算を用いて求める。
【００３１】
図１３は写像Ｆの幾何学計算を示すフローチャートであり、図１４は観測画像４の一部の画素位置の配列を示す図である。また、図１５は写像Ｆの幾何学計算の説明のための図である。写像Ｆは、座標値Ａ（ｐ，ｌ）に対する格子点Ａ’（ｘ，ｙ）を求めるものである。なお、図１３における軌道データ１４，姿勢データ１５，センサ光学モデル１６，画素位置データ１７が、図８における歪モデル６に相当するものである。
【００３２】
図１４では、４個の画素位置Ａ₁(ｐ₁，ｌ₁)、Ａ₂(ｐ₂，ｌ₁)、Ａ₃(ｐ₂，ｌ₂）、Ａ₄(ｐ₁，ｌ₂）を示している。まず、画素位置Ａ₁(ｐ₁，ｌ₁）を例にして、図１３に示す計算を説明する。
【００３３】
図１３及び図１５において、画素位置Ａ₁（ｐ₁，ｌ₁）の座標値ｌ₁からその画素位置Ａ₁の撮影時刻ｔが得られる（ステップ３００）。衛星１の軌道データ１４と得られた撮影時刻ｔとから座標値ｌ₁の撮影時の衛星１の慣性空間における位置(衛星位置）が得られ（ステップ３０１）、この衛星位置と衛星１でのセンサの光学モデル１６とから座標値ｌ₁の撮影時のセンサの光学中心の慣性空間での位置（センサ光学中心）１９が得られる（ステップ３０２）。
【００３４】
また、ステップ３００で得られた撮影時刻ｔと衛星１の姿勢データ１５とからこの撮影時刻ｔでの衛星１の慣性空間における姿勢（衛星姿勢）が得られ（ステップ３０３）、この衛星姿勢とセンサの光学モデル１６と画素位置データ１７（ここでは、画素位置Ａ₁（ｐ₁，ｌ₁）での座標値ｐ₁）とから、センサ面２０上のこの座標値ｐ₁に該当するセンサ素子（即ち、セル）の位置２１を通る、画素位置Ａ₁（ｐ₁，ｌ₁）に対する視線ベクトル２２が慣性空間で求められる（ステップ３０４）。そして、ステップ３０２で求めたセンサ光学中心１９から延びる視線ベクトル２２と地表モデル２３との交点２４が得られる（ステップ３０５）。この交点２４が観測画像４の画素位置Ａ₁（ｐ₁，ｌ₁）の地表モデル２３上の緯度，経度で示される位置を表わするものであり、これと所定の地図投影法１８による計算により、補正画像５での画素位置Ａ₁（ｐ₁，ｌ₁）に対応する座標位置Ａ’（ｘ，ｙ）を求める。なお、地表モデル２３としては、通常、回転楕円体を用いる。
【００３５】
ところで、このようにして求めた画素位置Ａ₁（ｐ₁，ｌ₁）に対応する座標位置Ａ’（ｘ，ｙ）は、補正画像５で設定された格子点１０（図１１）とは一致しない。そこで、図１４での画素位置Ａ₁（ｐ₁，ｌ₁）に隣接する画素位置、例えば、画素位置Ａ₂（ｐ₂，ｌ₁）を選択し、これについて同様の計算を行ない、補正画像５での画素位置Ａ₂（ｐ₂，ｌ₁）に対応する座標位置Ａ’（ｘ，ｙ）を求める。そして、これでも補正画像５で設定された格子点１０と一致しない場合には、再び他の画素位置を選択して同様の計算を行ない、というようにして繰り返し計算を行なって計算結果が補正画像５の格子点１０に収束するようにしていき、例えば、４ないし５回程度の計算により、補正画像５の格子点Ａ’（ｘ，ｙ）に対応する観測画像４の画素位置Ａ（ｐ，ｌ）が近似的に求める。
【００３６】
そして、かかる計算を補正画像５の格子点１０毎に行ない、補正画像５の各格子点１０に対応する観測画像４の画素位置Ａ（ｐ，ｌ）が得られることになる。以上が図１２における格子点計算２００である。
【００３７】
図１２における係数算出２０１では、図１１に示す補正画像５の各ブロック９毎に上記数２で示すバイリニア関数の係数ａ₀，ａ₁，ａ₂，ａ₃，ｂ₁，ｂ₂，ｂ₃を求めるものである。つまり、ブロック９毎にバイリニア関数を求めるものである。
【００３８】
例えば、図１６に示すように、補正画像５でのブロック９に対して、観測画像４では、ブロック１１が対応するものとすると、ブロック１１内の画素をブロック９内に逆写像Ｆ^-1するためのバイリニア関数を求めるのが、この係数算出２０１である。このために、これらブロック９，１１の四隅の格子点の座標値を用いて数２のａ₀，ａ₁，ａ₂，ａ₃，ｂ₁，ｂ₂，ｂ₃を求める。
【００３９】
ここで、ブロック９での四隅の格子点Ａ₁’(ｘ₁，ｙ₁），Ａ₂’(ｘ₂，ｙ₁），Ａ₃’(ｘ₂，ｙ₂），Ａ₄’(ｘ₁，ｙ₂）はブロック１１の四隅の格子点（画素位置）Ａ₁（ｐ₁，ｌ₁），Ａ₂（ｐ₂，ｌ₂），Ａ₃（ｐ₃，ｌ₃），Ａ₄（ｐ₄，ｌ₄）に対応するものであり、従って、ブロック９，１１の四隅の格子点の座標対（(ｘ₁，ｙ₁），（ｐ₁，ｌ₁)）、（(ｘ₂，ｙ₁），（ｐ₂，ｌ₂)）、（(ｘ₂，ｙ₂），（ｐ₃，ｌ₃)）、（(ｘ₁，ｙ₂），（ｐ₄，ｌ₄)）から、次の数３で示す連立方程式、即ち、
【数３】

を解くことにより、所望の係数ａ₀，ａ₁，ａ₂，ａ₃，ｂ₁，ｂ₂，ｂ₃が得られる。このようにして得られたこれら係数を数２に代入することにより、これらブロック９，１１に対するバイリニア関数が得られ、これがブロック９内の画素位置に対するブロック１１内の画素位置を示すことになる。なお、かかるバイリニア関数により、ブロック９内の画素位置がブロック１１内の画素位置に一致しない場合には、図９で説明したような補間処理を行なうことはいうまでもない。
【００４０】
そして、以上のような係数算出２０１（図１２）を補正画像５のブロック９毎に行なうことにより、これらブロック９毎のバイリニア関数が得られることになり、観測画像４から補正画像５への逆写像Ｆ^-1が可能となる。
【００４１】
この第２の方法には、（１）ブロック毎のバイリニアモデルなので、第１の方法のごとき最小自乗法に起因する誤差の発散がない、（２）ブロックの境界で歪みモデルが連続である、（３）数３を用いるので、補正画像の１点１点の位置座標の計算の演算量が数１の多項式に比べて小さい、（４）地上基準点がなくてもモデリングができる、（５）なお、地上基準点がある場合にそれを用いて衛星の軌道データ、姿勢データに含まれる誤差を推定し、それを補正することにより、よりモデリングの精度が高められるという利点がある。
【００４２】
【非特許文献１】
Demystification of IKONOS by Thierry Toutin and Philip Cheng Earth Observation Magazine vol.9,No.7,pp.17-21,2000
【００４３】
【発明が解決しようとする課題】
ところで、上記第１の補正方法では、（１）モデリングのために地上基準点が必須であり、地図などが整備されていない地域などで地上基準点が得られない場合に適用ができない、（２）モデリングの精度が地上基準点の数，配置に依存し、かかる数，配置によっては、多項式の係数の推定精度が劣化し、特に、地上基準点以外での誤差が発散する恐れがある、（３）補正画像の一つ一つに対して、上記数１の計算が必要となり、補正処理の計算量が大きくなるという問題点があった。
【００４４】
また、上記第２の補正方法は、地表面が起伏しない回転楕円体モデルで近似しているので、観測対象地域の標高データDEMを考慮し、地面の起伏に起因する倒れ込み歪を考慮したオルソ補正には、そのまま適用できない、という問題があった。
【００４５】
ここで、図１７により、地面の起伏に起因する倒れ込み歪について説明する。
【００４６】
同図において、上記第２の補正方法によると、回転楕円体の地表モデル２３に対し、実際には、起伏した地表面２５がある。そこで、衛星１０の撮影により、上記のようにして観測画像４を得る場合、センサ面２０でのセルＳについてみると、このセルＳは、上記のようにして求めたセンサ光学中心１９と視線ベクトル２２とにより、地表モデル２３との交点Ｐ₁を所定の地図投影法で表わした座標を補正画像の画素位置とするものであるが、実際には、起伏する地表面２５との交点Ｑを通る鉛直線と地表モデル２３との交点Ｐ₂を所定の地図投影法で表わした座標が正しい補正画像の画素位置となるものである。これら交点Ｐ₁と交点Ｐ₂との位置ずれに相当する補正画像上の位置ずれを倒れ込み歪δという。
【００４７】
本発明の目的は、かかる問題を解消し、地上基準点を必要とせず、画像全体で誤差の発散がなくて計算量が過大となることなく、しかも、地表面の起伏に起因する倒れ込み歪を補正可能とした衛星画像のオルソ補正方法を提供することにある。
【００４８】
【課題を解決するための手段】
上記目的を達成するために、本発明は、基準回転楕円体地球モデルの海抜標高ゼロの基準面に対して、複数の相異なる標高の基準面を設定し、補正画像を等間隔の格子でブロックに分割し、ブロック毎に、かつ上記複数の標高の基準面毎に、バイリニア関数による歪みのモデル式を求めるものである。
【００４９】
また、本発明は、補正画像の各画素位置について、該当するブロックでの画素位置の標高の値を挟む上下の標高の２つの基準面のモデル式を用いて、夫々の標高値に対応する観測画像の画素位置を計算し、標高データの値により、観測画像のこれら２つの画素位置を線形補間して、観測画像の画素位置を求めるものである。
【００５０】
【発明の実施の形態】
以下、本発明の実施形態を図面を用いて説明する。
【００５１】
図１は本発明による衛星撮影画像のオルソ補正処理方法の一実施形態を示す説明図であって、２６ａ，２６ｂ，２６ｃ，２６ｄは基準面である。
【００５２】
同図において、地表面２５を断面で示しており、補正画像５上で複数の相異なる標高の基準面２６ａ，２６ｂ，２６ｃ，２６ｄを設定する。ここで、基準面２６ａは基準回転楕円体の地表モデル２３の海抜標高ゼロ（Ｈ_O＝０）の基準面とし、基準面２６ｂは標高Ｈ₁の基準面、基準面２６ｃは標高Ｈ₂の基準面、基準面２６ｄは標高Ｈ₃の基準面であり、Ｈ₀＜Ｈ₁＜Ｈ₂＜Ｈ₃である。なお、標高Ｈ₁，Ｈ₂，Ｈ₃は夫々、所望とするように適宜設定される。
【００５３】
図１では、ブロック１，２，３を示している。そして、ブロック毎に各基準面でのバイリニア関数を求め、これを用いることにより、各基準面毎にその画素位置に対する観測画像での画素位置を求める。このようにして求めた基準面での画素位置と観測画像４での画素位置との対応関係を用いて、地表面２５上での画素位置と観測画像４での画素位置との対応関係を求めるものであり、これにより、地表面２５の起伏に起因する歪をなくすようにするものである。
【００５４】
以下、かかる処理について説明する。
【００５５】
図２はこの実施形態の補正処理を示すフローチャートであって、図８に示した補正処理での歪みのモデリング処理１００に相当するものである。
【００５６】
補正画像５には、図１に示したように、複数の基準面２６ａ〜２６ｄが設定されており、これら基準面に、互いに対応するようにして、先の従来技術と同様、等間隔の格子によってブロックが設定されている。補正画像を等間隔の格子によりブロックに区分する。
【００５７】
いま、補正画像４での着目画素Ａ’の座標を、この着目画像Ａ’が起伏する地表面２５上にあって標高データDEMで与えられる標高ｚを有していることから、（ｘ，ｙ，ｚ）と表わすことにする。
【００５８】
そこで、図２において、まず、着目画素Ａ’の座標値（ｘ，ｙ）をもとに該当するブロックを求める（ステップ４００）。図１の場合、ブロック１が求められる。
【００５９】
次に、この着目画素Ａ’の標高ｚをもとに、この着目画素Ａ’の位置を直接上下に挟む２つの基準面を求める。この標高値ｚが標高値Ｈ_i（但し、ｉ＝０，１，２，……）の基準面と標高値Ｈ_i+1の基準面との間、即ち、
【数４】

であるとき、標高値Ｈ_i，Ｈ_i+1の基準面を求める。図１の場合、標高値Ｈ₀の基準面２６ａと標高値Ｈ₁の基準面２６ｂとが求まることになる（ステップ４０１）。なお、ブロック２の点Ｂでは、標高値Ｈ₀の基準面２６ａと標高値Ｈ₁の基準面２６ｂとが求まり、点Ｃでは、標高値Ｈ₁の基準面２６ｂと標高値Ｈ₂の基準面２６ｃとが求まる。
【００６０】
そして、求めた標高値Ｈ_iの基準面での上記数２で示すバイリニア関数の係数ａ₀，ａ₁，ａ₂，ａ₃，ｂ₀，ｂ₁，ｂ₂，ｂ₃を後述する方法で求め、これら係数を用いたバイリニア関数を用いて、標高値Ｈ_iの基準面での画素位置Ａ_i'(ｘ_i,ｙ_i）に対応する観測画像５での座標値Ａ_i（ｐ_i，ｌ_i）を求める（ステップ４０２）。同様にして、標高値Ｈ_i+1の基準面における画素位置Ａ_i+1’(ｘ_i+1,ｙ_i+1）に対応する観測画像４での座標値Ａ_i+1(ｐ_i+1，l_i+1）を求める。図１の場合には、標高値Ｈ₀の基準面２６ａでの画素位置Ａ₀'(ｘ₀，ｙ₀，ｚ₀）に対応する観測画像５での座標値Ａ₀（ｐ₀，ｌ₀）と、標高値Ｈ₁の基準面２６ｂにおける画素位置Ａ₁’(ｘ₁，ｙ₁，ｚ₁）に対応する観測画像４での座標値Ａ₁（ｐ₁，l₁）が求まることになる(ステップ１２３）。
【００６１】
次に、観測画像４において、座標値Ａ_i+1と座標値Ａ_iとから着目画素Ａ’の標高ｚに対する座標値位置Ａ（ｐ，ｌ）を、次の数５により、
【数５】

線形補間して求める（ステップ４０４）。
【００６２】
図３は図２におけるステップ４０２，４０３でのバイリニア関数の係数の算出方法を示すフローチャートである。
【００６３】
同図において、標高が異なる各基準面毎に、図１２で説明した計算を行なうものである。ここで、格子点計算５００は図１２での格子点計算２００と同様のものであり、係数算出５０１も図１２での係数算出２０１と同様のものである。格子点計算５００での写像Ｆの幾何学計算も、図１３に示したものと同様である。但し、図３に示す計算では、計算の基準となる基準面が標高値Ｈ_iに応じて変更されるものである。
【００６４】
これを図４でもって説明すると、補正画像５での同じブロックの標高値Ｈ₀ ，Ｈ_i，Ｈ_i+1の基準面２６ａ,２６_i,２６_i+1について、基準面２６ａでの位置Ｐ₀に対し、基準面２６_i，２６_i+1では、この位置Ｐ₀の真上の位置Ｐ_i，Ｐ_i+1（即ち、位置Ｐ₀と緯度経度が同じ位置）が対応する。そこで、基準面２６ａでの位置Ｐ₀に対するセンサ面２０でのセルをＳ₁とすると、基準面２６_iでの位置Ｐ_iに対するセンサ面２０でのセルはセルＳ₂であり、基準面２６_i+1での位置Ｐ_i+1に対するセンサ面２０でのセルはセルＳ₃である。
【００６５】
このように、同じ緯度，経度の位置であっても、標高が異なる基準面に応じてセルが異なることなり、従って、基準面に応じてセルを異ならせる処理、即ち、図１３での座標値ｐ及びｌを修正することが必要となる。
【００６６】
図５は複数の標高が異なる基準面に対する観測画像４での格子点の写像を示す図であり、２７ａ，２７ｂ，２７ｃ，２７ｄはブロックである。
【００６７】
同図において、ブロック２７ａは図１での基準面２６ａに対する領域でのブロック１に対応するものであり、ブロック２７ｂは同じく基準面２６ｂに対する領域でのブロック１に対応するものであり、ブロック２７ｃは同じく基準面２６ｃに対する領域でのブロック１に、ブロック２７ｄは同じく基準面２６ｄに対する領域でのブロック１に夫々対応するものである。このように、補正画像５の同じブロックであっても、基準面２６ａ，２６ｂ，２６ｃ，２６ｄ毎に観測画像４でのブロックが異なることになる。これは、図３での格子点計算５００において、図１３における座標値１７を標高に応じて修正したことによるものである。
【００６８】
このようにして、格子点計算５００により、補正画像５の基準面毎の格子点位置とこれに対する観測画像４での画素位置が求まると、係数算出５０１により、図１２に示した係数算出２０１と同様にして、各基準面でのブロック毎に上記数２で示すバイリニア関数の係数ａ₀，ａ₁，ａ₂，ａ₃，ｂ₀，ｂ₁，ｂ₂，ｂ₃が求められ、各基準面でのブロック毎のバイリニア関数が得られる。
【００６９】
そこで、例えば、図１での補正画像４でのブロック１の標高値ｚの画素位置Ａ’に対する観測画像４での画素位置Ａを求める場合には、得られたバイリニア関数により、この画素位置Ａ’を直接挟む基準面２６ａ，２６ｂのブロック１に対する観測画像４でのブロック２７ａ，２７ｂ（図５）において、基準面２６ａでの画素位置Ａ’に対応する画素位置Ａ_a’に対する観測画像４の画素位置Ａ_aをブロック２７ａから求め、基準面２６ｂでの画素位置Ａ’に対応する画素位置Ａ_b’に対する観測画像４の画素位置Ａ_bをブロック２７ｂから求め、これらの得られた画素位置Ａ_a，Ａ_bでもって画素位置Ａ’の標高値ｚに応じた補間処理を行なうことにより、補正画像５での画素位置Ａ’（ｘ，ｙ，ｚ）に対する観測画像４での画素位置Ａ（ｐ，ｌ）が得られることになる。
【００７０】
なお、基準面の標高Ｈ_iの選定には、次の条件、即ち、
（１）各ブロックについて、該ブロック内の標高値データDEMの最小値と最大値との範囲をカバーするに十分な数の相異なる標高を持つ基準面を設定すること
（２）隣接するブロックの境界で歪みモデルが連続となるように、隣接する２つのブロック夫々の標高値データDEMの変動範囲の共通部分を越える最小の基準標高値とこの共通部分を下まわる最大の基準標高値との間に合まれる全ての基準面の標高値は、これら隣接するブロックで共通の値を設定すること
を考慮する必要がある。
【００７１】
なお、この（２）について、図１のブロック２，３を例にして説明すると、ブロック２での標高値の変化はｈ₁〜ｈ₃であり、ブロック３での標高値の変化はｈ₂〜ｈ₄である。ここで、ｈ₁＜ｈ₂＜ｈ₃＜ｈ₄であると、標高の変化範囲ｈ₂〜ｈ₃は、その範囲でブロック２での標高の変動範囲とブロック３での標高の変動範囲とが重なっており、これがブロック２，３での「隣接する２つのブロック夫々の標高値データDEMの変動範囲の共通部分」である。また、上記の「共通部分を越える最小の基準標高値」は標高値Ｈ₂であり、また、「共通部分を下まわる最大の基準標高値」は標高値Ｈ₁である。そして、かかる標高値Ｈ₁，Ｈ₂間に基準面を設定するときには、かかる基準面の標高値はブロック２とブロック３とで同じでなければならないとするものである。
【００７２】
また、標高Ｈ_iの間隔は、ブロック内の標高値データDEMの変動の大きさ、ブロックのオフネデイア角（衛星１０の直下点から計った斜視角の大きさ）によって決定できる。
【００７３】
即ち、バイリニアモデルの近似誤差と倒れ込み量の標高値に対する線形近似誤差との合計が、補正後の誤差の許容値の範囲に収まればよい。定性的には、ブロック内の標高変動が大きいとき、斜視角が大きいとき、夫々倒れ込み量及びその線形近似誤差が大きくなるので、標高Ｈ_iの間隔を小さくする必要がある。予めこの近似誤差をモデル計算などにより推定して、間隔を決めることができる。
【００７４】
【発明の効果】
以上説明したように、本発明によれば、（１）与えられた標高データDEMを用いて、地表面の起伏に起因する倒れ込み歪を補正できる、（２）ブロック毎のバイリニアモデルなので、Rational polynomial methodの多項式を最小自乗法で求めることに起因する誤差の発散がない、（３）ブロックの境界で歪みモデルが連続である、（４）地上基準点がなくても、モデリングができる、（５）地上基準点がある場合には、それを用いて衛星の軌道データや姿勢データに含まれる誤差を推定し、それを補正することにより、モデリングの精度がより高められ、補正精度を向上できるという効果が得られる。これにより、精度と計算量の両面で実用的な衛星画像のオルソ補正方法を提供することができる。
【図面の簡単な説明】
【図１】本発明による衛星撮影画像のオルソ補正処理方法の一実施形態を示す説明図である。
【図２】図１に示す実施形態の補正処理を示すフローチャートである。
【図３】図２におけるステップ４０２，４０３でのバイリニア関数の係数の算出方法を示すフローチャートである。
【図４】図３における格子点計算５００の説明図である。
【図５】複数の標高が異なる基準面に対する観測画像の格子点の写像を示す図である。
【図６】衛星による画像撮影の説明図である。
【図７】衛星画像の補正の原理を示す図である。
【図８】観測画像の補正処理を示すフロ一チャートである。
【図９】図８におけるリサンプリングの説明図である。
【図１０】地上基準点を用いた観測画像の第１の補正方法の説明図である。
【図１１】観測画像の第２の補正方法での格子点，ブロックを説明するための図である。
【図１２】バイリニア関数での係数の算出方法を示すフローチャートである。
【図１３】図１２における格子点計算２００での写像Ｆの幾何学計算を示すフローチャートである。
【図１４】観測画像の一部の画素位置の配列を示す図である。
【図１５】図１３での写像Ｆの幾何学計算の説明のための図である。
【図１６】補正画像と観測画像とでのブロックの関係を説明するための図である。
【図１７】標高による倒れ込み歪を説明するための図である。
【符号の説明】
１衛星
２地表
３地表観測線
４観測画像
５補正画像
６歪モデル
７画素位置
８格子
９ブロック
１０格子点
１１ブロック
１２格子点
１３格子
１４軌道データ
１５姿勢データ
１６センサ光学モデル
１７座標値
１８地図投影法
１９センサ光学中心
２０センサ面
２１センサ素子位置
２２視線ベクトル
２３地表モデル
２４交点
２５地表面
２６ａ〜２６ｄ，２６_i，２６_i+1 基準面
２７ａ〜２７ｄブロック

Claims

補正画像の画素位置と衛星による観測画像の画素位置とを対応づける歪みのモデリング処理と、該歪のモデリング処理に従って該観測画像の画素強度を補間するリサンプリング処理とからなる衛星撮影画像のオルソ補正処理方法において、
該歪みのモデリング処理が、基準回転楕円体地球モデルの海抜標高ゼロの基準面に対して、複数の相異なる標高の基準面を設定するとともに、該補正画像を等間隔の格子でもって複数のブロックに分割し、
該補正画像の各画素位置に該当する該ブロックを求め、求めた該ブロックにおける該画素位置の標高の値を挟む上下の２つの基準面を求め、該ブロックの該基準面毎に求めたバイリニア関数による歪みのモデル式を用いて、夫々の標高値に対応する観測画像の画素位置を計算し、該補正画像の該画素位置の標高の値により、計算した２つの該画素位置を線形補間して、該補正画像での該画素位置に対する観測画像の画素位置を求めることを特徴とする衛星観測画像のオルソ補正処理方法。
請求項１において、
前記各ブロックでは、前記ブロック内の画素位置の標高の値の最小値と最大値とによる範囲を力バーするに十分な数の相異なる標高の基準面が設定され、
隣接する２つのブロックでの標高の値の変動範囲の共通部分を越える最小の基準標高値と該共通部分を下まわる最大の基準標高値との間に含まれる全ての基準面の標高の値がこれら隣接するブロックで共通の値であることを特徴とする衛星観測画像のオルソ補正処理方法。