(57)【特許請求の範囲】
【請求項１】
コンピュータが、
飛行しながら所定のカメラ装置で地上を順次撮影して得た第１の空中写真Ｉ１及び第２の空中写真Ｉ２のそれぞれの画素の画像データ、地形モデルを記憶手段に記憶しておく工程と、
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２のそれぞれの各画素の画像データを読み込む写真データ読込工程と、
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２の指標の入力に伴って、この指標と前記カメラ装置の前記撮影時の各パラメータを用いて前記第１の空中写真Ｉ１の座標及び第２の空中写真Ｉ２の座標と画像座標の変換係数を求める内部標定工程と、
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２を撮影した際の前記カメラ装置の概略位置と姿勢Ｅ１及び姿勢Ｅ２を評価関数Ｅａｌｌ（Ｅａｌｌ＝（Ｅ１，Ｅ２））として読み込む工程と、
前記第１の空中写真および第２の空中写真に写っているエリアの地形モデルである地表面上の平面座標（ｘ、ｙ）に標高ｚを与える座標を読み込む地形モデル取得工程と、
前記地形モデルを用いて、地上座標（ｘ、ｙ）に対応する地表面の三次元座標（ｘ、ｙ、ｚ）を算出し、前記概略位置と姿勢Ｅ１及びＥ２並びに前記第１及び第２の空中写真の座標と画像座標の変換係数を用いて前記三次元座標に対応する前記第１の空中写真Ｉ１及び第２の空中写真座標Ｉ２の各輝度値を算出し、これらを用いて前記第１の空中写真Ｉ１及び第２の空中写真座標Ｉ２に対応するオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）及びオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）を生成する工程と、
前記オルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）及びオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）の輝度差ｅ（ｘ，ｙ，Ｅａｌｌ、Ｉ１、Ｉ２）を、
［数１］
ｅ（ｘ，ｙ，Ｅａｌｌ、Ｉ１、Ｉ２）＝
Ｉｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）―Ｉｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）
として算出する工程と、
前記算出された差を残差とし、この残差を最小自乗法で前記姿勢Ｅ１及び姿勢Ｅ２を求める工程と
を行うことを特徴とする空中写真の位置及び姿勢の計算方法。
【請求項２】
前記概略位置と姿勢Ｅ１および姿勢Ｅ２の相対関係が写真測量のステレオカメラの相互標定工程によって与えられるばあい、これを拘束条件として前記最小自乗法による最適化を行う請求項１記載の空中写真の位置及び姿勢の計算方法。
【請求項３】
前記概略位置と姿勢Ｅ１および姿勢Ｅ２に含まれる撮影高度Ｚは、撮影縮尺Ｓおよび撮影地区の平均標高Ｈおよび焦点距離Ｃによって
【数２】
0004480212.tif 000020
として与えることを特徴とする請求項１記載の空中写真の位置及び姿勢の計算方法。
【請求項４】
前記地形モデルの各座標に対応する写真座標のオルソフォトの輝度を
【数３】
0004480212.tif 000021
と定義する工程と、
を行うことを特徴とする請求項１記載の空中写真の位置及び姿勢の計算方法。
【請求項５】
前記評価関数の収束判定は、
【数６】
0004480212.tif 000022
を用いて判定することを特徴とする請求項１記載の空中写真の位置及び姿勢の計
算方法。
【請求項６】
前記収束判定は、最適化時の外部標定要素の修正量と初期値とを加算し、この加算結果が所定のしきい値に到達するまで繰り返すことを特徴とする請求項５記載の空中写真の位置及び姿勢の計算方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は空中撮影されたステレオ写真を、所定の地形モデルの座標系に容易に定義できる空中写真の位置及び姿勢の計算方法に関する。
【０００２】
【従来の技術】
一般に写真の像は、被写体から反射された光がレンズ中心を直進してフイルムの平面に投影されたものである。この投影は、中心投影とも呼ばれる。すなわち、写真は中心投影像と言える。
【０００３】
しかし、実際は光は空気の密度や温度差、レンズの特性等によってわずかに曲げられる。このため、精度の高い測量結果を得るには、前述の中心投影の原理を乱す諸原因を取り除く必要がある。
【０００４】
また、例えばピンホールカメラに理想的なレンズを用いた場合は、図７に示すように実空間に存在する任意の点Ｐから出た光は針穴Ｏを通過してフイルム面上の点ｐに像として記録される。つまり、点Ｐ及び針穴Ｏ並びに像ｐが同一直線上に存在しなければならないという条件が成立する。この条件を共線条件と称する。
【０００５】
また、理想的なレンズではレンズの光学的性質を代表する２つの浮動な点Ｓ及びＳ′が存在する。すなわち、実空間の点Ｐから出た光は必ずレンズのＳ及びＳ′を通過して像ｐとしてフィルム面に投影されるので、点Ｐと点Ｓを結ぶ直線ＰＳと、点Ｓ′と像ｐとを結ぶ直線Ｓ′ｐは平行になる。
【０００６】
従って、仮想的にＳの点とＳ′の点とを一致させたとすると、共線条件が成り立つ。針穴Ｏ、ＳとＳ′とを一致させたときのような点を投影中心と称する。写真の回転はこのＳ′を中心にして行われる。
【０００７】
また、写真測量は写真が１枚の場合（以下単写真測量という）と２枚以上の場合（以下立体写真測量）があるが、本説明では立体写真測量について説明する。
【０００８】
立体写真測量というのは、図８に示すように対象物Ｐを２個以上の異なる位置から撮影した２枚以上の写真Ａ、Ｂを用いるものであり、２つ以上の共線条件が得られる。
【０００９】
すなわち、２枚以上の写真の投影中心を通る光束が互いに交わる前方交会の結果として、対象物の三次元座標が求められることになる。
【００１０】
そして、写真の外部標定（基準点の地上座標（絶対座標）と、対応する写真上の像の平面座標とから、撮影時のカメラ位置と傾きを求めること）が完了すると、２枚以上の対応する写真の像と投影中心を結ぶ光束は、実空間においてことごとく交会し、立体像が得られる。この立体像は写真測量ではモデルと呼ばれる。
【００１１】
ここで、図８の各座標系について定義する。図８に示す写真座標系（ｘ、ｙ）とは、主点を原点とする２次元座標であり、カメラ座標系とは投影中心を原点としたものである。
【００１２】
また、モデル座標系は２枚一組の立体写真から形成される立体像を定義するための３次元座標系であり、一般には左写真の投影中心を原点にとる。
【００１３】
前述のモデルの形成は、写真が撮影されたときの状態の空間位置に厳密に等しくしなくとも、撮影されたときの写真の位置及び傾きの関係が相対的に同じであればよい、このような相対的な相似なモデルを作ることを相互標定と称している。
【００１４】
つまり、相互標定とはモデルを形成することであり、立体写真の対応する２本の光束がことごとく交会するように、左右それぞれの投影中心の位置及び傾きを定めることにある。
【００１５】
射影幾何学の原理から２枚の立体写真の対応する光束がことごとく交会するためには、少なくとも５本の光束が交会すればよいことが知られている。
【００１６】
前述の相互標定の方法に共面条件を用いるものもある。これは図９に示すように、２本の光束が交会するためには、２本の空間直線が同一平面に含まれるという条件が成立することを利用した方法（以下共面条件という）である。
【００１７】
このような各種条件、座標系を用いてステレオ撮影された２枚の写真の撮影位置（Ｘｏ、Ｙｏ、Ｚｏ）及び地上座標系（絶対座標系）に対する回転角（ω、φ、κ）を空中三角測量によって求める。
【００１８】
この空中三角測量にあたっては、両方の写真間の相対的な姿勢や位置関係は写真内に共通に写っている対応点（パスポイントともいう）を観測することにより求めていた。
【００１９】
前述の絶対座標系に対する撮影位置及び姿勢を求めるには、写真内に写っている基準点（絶対座標が既知の点；国土地理院の情報、又は実際に計測した点）を用いて計算する。
【００２０】
この撮影位置及び姿勢の算出に当たっては、一般に共線条件若しくは共面条件を用いて連立方程式を立て、これに最小自乗法等の最適計算を適用する。
【００２１】
【発明が解決しようとする課題】
しかしながら、従来のような空中写真を用いた測量は、写真を撮るエリアに予め複数の基準点を設ける必要（作業員が現地に行って基準とする点の緯度、経度を計測又は国土地理院の地理座標を用いる）があるので非常にコストと時間を要するという課題があった。
【００２２】
このため、近年は空中写真測量以外の方法であるレーザスキャナ測量や衛星写真測量等で地形情報を計測することが多くなってきている。
【００２３】
また、国土地理院発行数値地図は５０ｍメッシュのものがあり、このようなメッシュ間隔のものを用いた地形情報の利用も検討されてきている。
【００２４】
このような地形データ（地形モデルともいう）と空中写真とを利用してＣＧ画像等を作成する場合、
写真測量は写真内に写っている地上座標系を基準とする複数の基準点、両方の写真のパスポイント（地上座標系、写真座標系）を用いて計算するものであるから、地形モデル上の各点が対応する空中写真上の位置を求めるためには、地形モデルの座標系に対する空中写真の撮影位置（カメラ位置ともいう）及び姿勢を求めることが必要になる。
【００２５】
このため、レーザスキャナ測量や衛星写真測量等で得られた地形モデルを用いた写真測量は容易に用いることができないという課題があった。
【００２６】
本発明は以上の課題を解決するためになされたもので、地形モデルを用いて地形モデル座標系に対する空中写真の位置及び姿勢を容易に得ることができる空中写真の位置及び姿勢の計算方法を得ることを目的とする。
【００２７】
【課題を解決するための手段】
本発明の空中写真の位置及び姿勢の計算方法は、
コンピュータが、
飛行しながら所定のカメラ装置で地上を順次撮影して得た第１の空中写真Ｉ１及び第２の空中写真Ｉ２のそれぞれの画素の画像データ、地形モデルを記憶手段に記憶しておく工程と、
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２のそれぞれの各画素の画像データを読み込む写真データ読込工程と、
【００２８】
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２の指標の入力に伴って、この指標と前記カメラ装置の前記撮影時の各パラメータを用いて前記第１の空中写真Ｉ１の座標及び第２の空中写真Ｉ２の座標と画像座標の変換係数を求める内部標定工程と、
前記第１の空中写真Ｉ１及び第２の空中写真Ｉ２を撮影した際の前記カメラ装置の概略位置と姿勢Ｅ１及び姿勢Ｅ２を評価関数Ｅａｌｌ（Ｅａｌｌ＝（Ｅ１，Ｅ２））として読み込む工程と、
前記第１の空中写真および第２の空中写真に写っているエリアの地形モデルである地表面上の平面座標（ｘ、ｙ）に標高ｚを与える座標を読み込む地形モデル取得工程と、
前記地形モデルを用いて、地上座標（ｘ、ｙ）に対応する地表面の三次元座標（ｘ、ｙ、ｚ）を算出し、前記位置と姿勢Ｅ１及びＥ２並びに前記第１及び第２の空中写真の座標と画像座標の変換係数を用いて前記三次元座標に対応する前記第１の空中写真Ｉ１及び第２の空中写真座標Ｉ２の各輝度値を算出し、これらを用いて前記前記第１の空中写真Ｉ１及び第２の空中写真座標Ｉ２に対応するオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）及びオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）を生成する工程と、
【００２９】
前記オルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）及びオルソフォトＩｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）の輝度差ｅ（ｘ，ｙ，Ｅａｌｌ、Ｉ１、Ｉ２）を、
【数１】
ｅ（ｘ，ｙ，Ｅａｌｌ、Ｉ１、Ｉ２）＝
Ｉｏｒｔｈ（ｘ、ｙ、Ｅ１、Ｉ１）―Ｉｏｒｔｈ（ｘ、ｙ、Ｅ２、Ｉ２）
として算出する工程と、
前記算出された差を残差とし、この残差を最小自乗法で前記姿勢Ｅ１及び姿勢Ｅ２を求める工程と
を行うことを要旨とする。
【００３０】
【発明の実施の形態】
図１は本実施の形態に係わる空中写真の位置及び姿勢の計算方法を説明する概念図である。
【００３１】
図１に示すように、航空機等で得られたステレオ写真（空中写真Ａ、Ｂ）のピクセル座標をカメラ座標系に変換する内部標定をして、絶対座標に対する空中写真Ａ、Ｂの相対的な関係を求めた後に、地形モデル（例えばサーフェースモデル）上の点Ｐ（ｘ、ｙ、ｚ）からのステレオ写真（空中写真Ａ、Ｂ）上における画像輝度が一致する点ｐ１及び点ｐ２（共通エリア）を求め、この点ｐ１及び点ｐ２並びに地形モデルの点Ｐを用いて、地形モデルに対する空中写真Ａの位置及び姿勢Ｅ１並びに空中写真Ｂの位置及び姿勢Ｅ２を求める。
【００３２】
【数８】
0004480212.tif 000002
このようなことを実現するために本実施の形態では図２に示すように、Ａ／Ｄ変換工程１ａ及び１ｂと、内部標定工程２ａ及び２ｂと、相互標定工程３と、初期値決定工程４と、地表モデルに対する絶対標定工程５と、地表モデル取得工程６とを備える。
【００３３】
Ａ／Ｄ変換工程１ａ及び１ｂは、それぞれスキャナ（図示せず）で読み込まれた空中写真Ａ、Ｂを、デジタル変換する。但し、デジタル式のカメラを用いている場合は、このＡ／Ｄ変換は不要である。
【００３４】
内部標定工程２ａ及び２ｂは、それぞれの空中写真Ａ、Ｂの画像のピクセル座標を写真座標に直すためのキャリブレーションデータを得る。具体的には、空中写真Ａ、Ｂ上の既知の点のピクセル座標を計測し、ピクセル座標の変換式（射影変換又はアフィン変換）に代入して未知の係数（座標変換係数）を求める。
【００３５】
相互標定工程３は、この相互標定は空中写真Ａ及びＢを撮影したカメラ間の相対的な関係（位置、姿勢）を求めるものである。
【００３６】
例えば、空中写真Ａ及びＢにおける共通エリアを定め、この共通エリアにおける数個のパスポイントを入力し（手動又は自動）、何れか一方の空中写真のカメラを固定とする。つまり、一方のカメラの投影中心は（０、０、０）、姿勢は（０、０、０）と仮定すれば、標定要素は他方のカメラの投影中心（１、ｂｙ、ｂｚ）及びカメラの姿勢（ω、ψ、κ）の５つのパラメータ（標定要素）で表現できる。但し、普通は初期においては初期近似値として零が与えられる。
【００３７】
そして、近似値のまわりにテーラ展開して線形化したときの微分係数の値を求めて線形化した方程式を作り、最小自乗法を用いて近似値に対する微小補正量を求める。そして、この処理を所定回繰り返して収束させる。次に、空中写真Ａ、Ｂの光束の交会点のモデル座標を共面条件を用いて求める。
【００３８】
初期値決定工程４は、後述する接続標定における各写真の投影中心を求めるための初期値である。この初期値の内Ｘ及びＹは、空中写真Ａ、Ｂの概略撮影位置（例えば中心）を地図上で計測することにより求める。
【００３９】
また、概略撮影位置の撮影高度Ｚは、撮影縮尺Ｓと撮影地区の平均標高Ｈ及び焦点距離Ｃから次式で求める。
【００４０】
【数９】
0004480212.tif 000003
また、空中写真の場合は一般に垂直撮影であり、各カメラの姿勢成分の内ω及びφの初期値は「０」として差し支えない。また、κは概略の飛行方向から求めることができる。
【００４１】
さらに、空中写真Ａ及びＢのペアからステレオ処理によって得られる地表面データと地形モデルとを比較し、相関係数が高くなるようマッチングを行うことにより、位置及び姿勢の初期値を求めることが可能である。
【００４２】
地表モデル取得工程６は、国土地理院発行の数値地図（５０ｍメッシュ）等に基づく地表モデルを取得する。
【００４３】
地表モデルに対する絶対標定工程５は、地表モデル、各写真の投影中心を求めるための初期値を入力し、相互標定で求めたモデルを地表モデルの座標系に変換する絶対標定を行う。
【００４４】
このモデルのブロックに対して、各写真の外部標定要素（カメラ位置（Ｘｏ、Ｙｏ、Ｚｏ）及び姿勢（κ、φ、ω））を最小自乗法により同時に求めている。
【００４５】
この処理を以下に図及び数式を用いて詳細に説明する。例えば、図３に示すように、初期値決定工程４の地表モデルに基づく初期値であ地理座標標高データｚ＝ｚ（ｘ、ｙ）及び画像Ｉ１、画像Ｉ２、各々の画像の外部標定要素Ｅ１、Ｅ２が与えられたとき、図３に示すようにオルソフォトの重なる部分が同じになるように外部標定要素を最適化する。
【００４６】
この図３における最適化処理を数１〜  を用いて記述する。初めに外部標定要素を定義するための初期評価関数について説明する。
【００４７】
【数１０】
0004480212.tif 000004
【数１１】
0004480212.tif 000005
【数１２】
0004480212.tif 000006
【数１３】
0004480212.tif 000007
次に、最小自乗法による最適化について説明する。
【００４８】
【数１４】
0004480212.tif 000008
【数１５】
0004480212.tif 000009
【数１６】
0004480212.tif 000010
次に、偏微分係数の求め方を説明する。
【００４９】
【数１７】
0004480212.tif 000011
【数１８】
0004480212.tif 000012
一方、相互標定が済んでいる場合について以下に説明する。
【００５０】
【数１９】
0004480212.tif 000013
【数２０】
0004480212.tif 000014
【数２１】
0004480212.tif 000015
【数２２】
0004480212.tif 000016
【数２３】
0004480212.tif 000017
【数２４】
0004480212.tif 000018
【数２５】
0004480212.tif 000019
すなわち、最小自乗法の概略処理は、図５及び図６に示すようにステレオ画像（空中写真Ａ、Ｂ）、初期値決定工程でのＥａｌｌの初期値、地形モデル、式２ー９の初期評価関数χｏを読み込む（Ｓ５０１）。
【００５１】
次に、最小自乗法による最適化に当たって行列Ｈ及びベクトルＶの計算を行い（Ｓ５０３）、修正量ΔＥａｌｌの計算（式３ー６）と初期値の更新（式３ー７）を行う（Ｓ５０５）。
【００５２】
そして、式２−９に基づいて評価関数χｌを計算する（Ｓ５０７）。次に、収束の判定を行い（Ｓ５０９）、収束していない場合は評価関数χｌを更新して処理ステップＳ５０３に戻す（Ｓ５１３）。ステップＳ５０９において、収束したと判定したときは、ステップＳ５０７で求めたＥａｌｌを送出する（Ｓ５１１）。
【００５３】
次に、ステップＳ５０３の行列Ｈ及びベクトルＶの計算を図６のフローチャートを用いて説明する。
【００５４】
行列Ｈ及びベクトルＶの計算は、ステレオ画像、Ｅａｌｌの初期値、地形モデルを読み込み（Ｓ６０１）、その後に全て「０」にセットして行列Ｈ及びベクトルＶの初期化を行う（Ｓ６０３）。
【００５５】
次に全ての地形モデル上の点を定義して、投影される写真座標Ｐｐｈｔ１及びＰｐｈｔ２の計算を行う（Ｓ６０７）。
【００５６】
そして、写真座標Ｐｐｈｔ１及びＰｐｈｔ２がともに写真範囲内かどうかを判定する（Ｓ６０９）。
【００５７】
ステップＳ６０９において、共に写真範囲と判定したときは、微分計算を行う（Ｓ６１１）。この計算は、相互標定結果を用いる場合は、式５−９、式５−１０、式５−１１、式５−１２を用いる。但し、一般には式４−１、４−２で計算する。そして、行列Ｈ及びベクトルＶへの加算（式３−４及び式３−５）を行う（Ｓ６１３）。
【００５８】
【発明の効果】
以上のように本発明によれば、空中を飛行しながら順次撮影された第１の空中写真、第２の空中写真が相互標定された後に、これらの空中写真に対応する地形モデルを読み込み、この地形モデル上の入力された点Ｐ（ｘ、ｙ、ｚ）に基づいて第１の空中写真、第２の空中写真の輝度が一致する地形モデルにおける位置、姿勢Ｅ１、Ｅ２を地形モデルに対する絶対標定で求めるようにしたことにより、空中写真測量以外の方法であるレーザスキャナ測量や衛星写真測量等又は国土地理院発行数値地図によって生成した地形モデルであっても、その地形モデル上に対する位置、姿勢Ｅ１、Ｅ２を容易に決定できるという効果が得られている。
【００５９】
また、地形モデルを用いて写真の位置、姿勢を決定できるので、従来のように空中写真をとるための基準点を実際に出かけて行って測量しなくともよいので、コストを低減できる。
【図面の簡単な説明】
【図１】本発明の概念を説明する説明図である。
【図２】本発明の実施の形態の概略工程図である。
【図３】本発明の実施の形態の共通エリアにおける空中写真のパスポイントの定義を説明する説明図である。
【図４】評価関数の加算対象を説明する説明図である。
【図５】最適化の処理を説明するフローチャートである。
【図６】最適化の処理を説明するフローチャートである。
【図７】共線条件を説明する説明図である。
【図８】各座標系を説明する説明図である。
【図９】共面条件を説明する説明図である。
【符号の説明】
２ａ  内部標定工程
２ｂ  内部標定工程
３  相互標定工程
４  地形モデル取得工程
５  地形モデルに対する絶対標定工程
【図１】
0004480212.tif 000023
【図２】
0004480212.tif 000024
【図３】
0004480212.tif 000025
【図４】
0004480212.tif 000026
【図５】
0004480212.tif 000027
【図６】
0004480212.tif 000028
【図７】
0004480212.tif 000029
【図８】
0004480212.tif 000030
【図９】
0004480212.tif 000031