JP2003323611A

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

Info

Publication number: JP2003323611A
Application number: JP2003069813A
Authority: JP
Inventors: Fuminobu Furumura; 文伸古村; Hiroshi Ueda; 浩史上田; Takashi Hino; 隆司樋野; Takatoshi Kodaira; 高敏小平
Original assignee: Hitachi Software Global Technology Ltd
Current assignee: Hitachi Software Global Technology Ltd
Priority date: 2002-04-24
Filing date: 2003-03-14
Publication date: 2003-11-14
Anticipated expiration: 2023-03-14
Also published as: EP1357516A2; US20030179950A1; EP1357516A3; CN1453746A; CN100405398C; JP3869814B2; DE60330524D1; EP1357516B1; US6810153B2

Abstract

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

Description

【発明の詳細な説明】

【０００１】

【発明の属する技術分野】本発明は、衛星から撮影した
地上画像の歪みを補正する画像処理方法に係り、特に、
地表の起伏に起因する歪みを補正するオルソ補正処理方
法に関する。

【０００２】

【従来の技術】始めに、衛星から撮影した画像（衛星画
像）とその補正の原理について説明する。

【０００３】図６は衛星による画像撮影の説明図であっ
て、１は衛星、２は地表、３は地表観測線である。

【０００４】同図において、衛星１は、下方に向けて搭
載したセンサ（図示せず）で地表２からの電磁波を検知
することにより、地表２を撮影する。かかるセンサが衛
星用センサとして一般的なＣＣＤなどのラインセンサで
ある場合、かかるセンサは、それを構成する１列のセル
（センサ素子）の配列方向が衛星１の進行方向Ｌに直交
するようにして、配置されており、かかるセンサによっ
て同時に観測される地表２上の範囲が地表観測線３であ
る。かかる地表観測線３は衛星１とともにその進行方向
Ｌに移動し、この移動とともに各セルの検知データが所
定のサンプリング周期で取り込まれる。かかる検知デー
タの取込みタイミングは全セル同時である。この地表観
測線３全体のデータが衛星画像の１ライン分に相当する
ものであり、従って、１回のデータ取込み（１サンプリ
ング）で撮影画像の１ラインに相当するデータを取得す
ることになる。このようにして、衛星１の軌道上の進行
に伴い、地表観測線３による一定の幅（Swath）で地表
２が連続観測される。

【０００５】図７はこのようにして得られた衛星画像の
補正の原理を示す図であり、４が観測画像、５が補正画
像である。

【０００６】同図において、以上の衛星による撮影画像
は、観測画像４として、地上に送られる。かかる観測画
像４は、衛星の軌道や姿勢の揺らぎ，センサの光学系の
歪み，地球表面の形状などの要因により、その画像の形
状に歪みが生じている。このため、衛星による撮影画像
の利用のためには、この観測画像４の歪みを補正し、所
定の地図投影法に従った補正画像５に変換する必要があ
る。ここで、かかる補正は、観測画像４の画素Ａ（ｐ，
ｌ）からこれに対応する補正画像５の画素Ａ’（ｘ，
ｙ）への写像を行なうものであるが、これを写像Ｆとい
う。

【０００７】なお、観測画像４は、衛星１でのセンサの
セルの配列方向（即ち、地表観測線３の長手方向）をｐ
座標、衛星１の進行方向Ｌをｌ座標とするｐ−ｌ座標系
で表わされ、また、補正画像５は、観測画像４のｐ座標
に対応する方向をｘ座標、観測画像４のｌ座標に対応す
る方向をｙ座標とするｘ−ｙ座標系で表わされる。

【０００８】図８はかかる観測画像４の補正処理を示す
フロ一チャートである。

【０００９】図７に示す補正画像５では、画素が等間隔
に配列・設定されているものであり、これに対し、観測
画像４では、上記のような歪みのため、補正画像５の各
画素に対応する画素位置は等間隔に配列されない。図８
に示す補正処理は、補正画像５における各画素位置毎
に、その画素位置に相当する観測画像４での画素位置の
画素を見つけ出し、これを補正画像５での該当する画素
位置に位置付けるものである。

【００１０】図８において、かかる補正は、歪みのモデ
リング１００とリサンプリング１０１とからなってい
る。歪みのモデリング１００では、補正画像５における
各画素位置毎に、その画素位置に相当する観測画像４で
の画素位置を求めるものである。図７についてみると、
予め作成された歪モデル６を用いて、補正画像５での画
素Ａ’(ｘ，ｙ）の位置に対応する観測画像４の画素Ａ
（ｐ，ｌ）の位置を求めるものである。そのためには、
写像Ｆの逆写像Ｆ-1を求める必要がある。

【００１１】ここで、歪モデル６は、かかる観測画像４
の撮影時での衛星の軌道や姿勢の揺らぎ，センサの光学
系の歪みなどから予測されるされるもの、あるいは衛星
の観測結果などから得られるものである。

【００１２】また、リサンプリング１０１では、歪みの
モデリング１００で補正画像５での画素Ａ’(ｘ，ｙ）
の位置に対応する観測画像４の画素Ａ（ｐ，ｌ）の位置
が求まると、この画素Ａ（ｐ，ｌ）の位置での画素の強
度を求め、それを補正画像５の画素Ａ’（ｘ，ｙ）での
画素の強度とするものである。

【００１３】図９は図８におけるリサンプリング１０１
の説明図であって、図７に示した観測画像４と補正画像
５とでもって説明する。なお、７は観測画像４での画素
の中心位置（以下、画素位置という）である。

【００１４】同図において、図８で説明した歪みのモデ
リング１００により、補正画像５の画素Ａ’（ｘ，ｙ）
の位置に対応する観測画像４の画素Ａ（ｐ，ｌ）の位置
が求まるのであるが、画素Ａ（ｐ，ｌ）の位置は、図示
するように、一般にセンサによる観測データの画素位置
７に一致しない。そこで、この画素Ａ（ｐ，ｌ）の位置
における観測データの強度を、この画素Ａ（ｐ，ｌ）の
位置の周辺の実際に得られた画素位置７の画素強度から
補間によって求める必要がある。このための補間方法と
しては、Nearest Neighboｒ，Bilinear，Cubic Convolu
tionなどの補間式が一般に用いられる。

【００１５】ところで、従来、上記の原理に基づく衛星
観測画像データの第１の補正方法としては、オルソ補正
の方法が知られている（例えば、非特許文献１参照）。

【００１６】これは、図８の歪みのモデリング１００を
次の数１で表わすものである。

【００１７】

【数１】この式は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 Thier
ry Toutin and Philip Cheng Earth Observation Magaz
ine 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 me
thodの多項式を最小自乗法で求めることに起因する誤差
の発散がない、（３）ブロックの境界で歪みモデルが連
続である、（４）地上基準点がなくても、モデリングが
できる、（５）地上基準点がある場合には、それを用い
て衛星の軌道データや姿勢データに含まれる誤差を推定
し、それを補正することにより、モデリングの精度がよ
り高められ、補正精度を向上できるという効果が得られ
る。これにより、精度と計算量の両面で実用的な衛星画
像のオルソ補正方法を提供することができる。

【図面の簡単な説明】

【図１】本発明による衛星撮影画像のオルソ補正処理方
法の一実施形態を示す説明図である。

【図２】図１に示す実施形態の補正処理を示すフローチ
ャートである。

【図３】図２におけるステップ４０２，４０３でのバイ
リニア関数の係数の算出方法を示すフローチャートであ
る。

【図４】図３における格子点計算５００の説明図であ
る。

【図５】複数の標高が異なる基準面に対する観測画像の
格子点の写像を示す図である。

【図６】衛星による画像撮影の説明図である。

【図７】衛星画像の補正の原理を示す図である。

【図８】観測画像の補正処理を示すフロ一チャートであ
る。

【図９】図８におけるリサンプリングの説明図である。

【図１０】地上基準点を用いた観測画像の第１の補正方
法の説明図である。

【図１１】観測画像の第２の補正方法での格子点，ブロ
ックを説明するための図である。

【図１２】バイリニア関数での係数の算出方法を示すフ
ローチャートである。

【図１３】図１２における格子点計算２００での写像Ｆ
の幾何学計算を示すフローチャートである。

【図１４】観測画像の一部の画素位置の配列を示す図で
ある。

【図１５】図１３での写像Ｆの幾何学計算の説明のため
の図である。

【図１６】補正画像と観測画像とでのブロックの関係を
説明するための図である。

【図１７】標高による倒れ込み歪を説明するための図で
ある。

【符号の説明】

１衛星２地表３地表観測線４観測画像５補正画像６歪モデル７画素位置８格子９ブロック１０格子点１１ブロック１２格子点１３格子１４軌道データ１５姿勢データ１６センサ光学モデル１７座標値１８地図投影法１９センサ光学中心２０センサ面２１センサ素子位置２２視線ベクトル２３地表モデル２４交点２５地表面２６ａ〜２６ｄ，２６_i，２６_i+1 基準面２７ａ〜２７ｄブロック

───────────────────────────────────────────────────── フロントページの続き (72)発明者古村文伸アメリカ合衆国コロラド州 80021 ウエストミンスターウエストムアドライブ 10355 スイート250 日立ソフトウエアグローバルテクノロジーリミテッド内 (72)発明者上田浩史アメリカ合衆国コロラド州 80021 ウエストミンスターウエストムアドライブ 10355 スイート250 日立ソフトウエアグローバルテクノロジーリミテッド内 (72)発明者樋野隆司東京都品川区東品川４丁目12番７日立ソフトウエアエンジニアリング株式会社内 (72)発明者小平高敏東京都品川区東品川４丁目12番７日立ソフトウエアエンジニアリング株式会社内Ｆターム(参考） 2C032 HA02 HB11 5B050 AA01 BA02 BA06 BA17 EA05 EA13 EA18 5B057 AA14 CA12 CA16 CB12 CB16 CD12 CD17 DA07 DB02 DB09 DC33

Claims

【特許請求の範囲】

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