JP6041126B2 - 反射パルス推定方法、及び反射パルス推定プログラム - Google Patents

Description

本願発明は、レーザー計測結果から、地盤など目的とする対象物の計測結果を抽出する技術に関するものであり、より具体的には、ウェーブフォームデータ（Ｗａｖｅ Ｆｏｒｍ Ｄａｔａ：ＷＦＤ）と呼ばれる大量の計測データの中から、モーメント値に基づいて所望の計測データを適切に抽出する反射パルス推定方法、及び反射パルス推定プログラムに関するものである。

地形図を作製するためなど広範囲に渡って地形計測する場合、従来では航空機から撮影した空中写真を利用するのが一般的であった。昨今では、航空レーザー計測や、衛星写真を利用した計測、あるいは合成開口レーダを利用した計測といった様々な計測手法が出現し、状況に応じて好適な手法を選択できるようになった。

このうち航空レーザー計測は、図１示すように、計測したい地形１の上空を航空機２で飛行し、地形１に対して照射したレーザーパルス３の反射信号を受けて計測するものである。航空機２には通常、ＧＰＳ（Ｇｌｏｂａｌ Ｐｏｓｉｔｉｏｎｉｎｇ Ｓｙｓｔｅｍ）などの測位計とＩＭＵ（Ｉｎｅｒｔｉａｌ Ｍｅａｓｕｒｅｍｅｎｔ Ｕｎｉｔ）などの慣性計測装置が搭載されているので、レーザーパルス３の照射位置（ｘ，ｙ，ｚ）と照射姿勢（ω，φ，κ）を把握することができ、その結果、照射時刻と受信時刻の時間差から計測点（レーザーパルス３が反射した地点）の３次元座標を得ることができる。

地形１で反射したレーザーパルス３は、航空機２に搭載されたセンサで受信される。このとき、反射してセンサまで戻ってきたレーザーパルス３の強度（以下、「反射強度」という。）を取得し、受信時刻と併せて記録される。この反射強度は、いわば受信したレーザーパルス３のエネルギーの大きさであり、直接的には電圧として計測され、電圧を換算することでエネルギーの大きさが得られる。図１に示すように、一回の計測（フライト）で多数のレーザーパルス３が照射されるので、そのレーザーパルス３の数に応じた数の照射強度が記録される。

ところで航空レーザー計測は、密に点群データを取得できるのが一つの特徴であり、例えば１ｍ^２当たり２〜３程度でレーザーパルス３の照射が可能であり、一回の計測で１００ｋｍ^２程度を計測することを考えれば、一度に数億のレーザーパルス３が照射されることとなる。また、照射した一つのレーザーパルス３から複数の信号（以下、対象物に反射した信号を「反射パルス」という。）を受信することもある。例えば、航空機２と地表面の間に障害物が存在しない場合は、照射したレーザーパルス３に対して一つの反射パルスを得るが、図６に示すように航空機２と地表面の間に樹木や草花がある場合は、枝や草などに反射するため複数の反射パルスが得られる。

このように一回のレーザー計測で極めて大量の反射信号を受信することになるが、従来では記憶領域（メモリの容量）の問題からこれらすべてを記録していなかった。つまり従来方式では、所定の条件を設定することで記録する反射パルスの数を制限していたのである。この条件とは、照射した一つのレーザーパルス３から記録する反射パルス数の最大値（３〜４以内）を設定することであり、記録する反射パルスの間隔の最小値（２．５〜３ｍ以上）を設定することである。

図６は、従来方式によるレーザー計測を示す説明図である。この図では、航空機２と地表面の間に樹木と草むらがあり、樹木の枝で４回（１ｓｔパルス〜４ｔｈパルス）、草むらで１回（５ｔｈパルス）、地表面で１回（６ｔｈパルス）の都合６回の反射パルスを受信し、本来なら６個の反射パルスデータを記録するはずである。しかしながら図６に示す従来方式では、４個の反射パルスデータを記録するにとどまっている。さらに、５ｔｈパルスと６ｔｈパルスの間隔が小さい（２．５ｍ未満）ため、最も重要な地表面の反射パルス（６ｔｈパルス）が記録されていない。

このように従来方式によるレーザー計測では、記録する反射パルスのデータを制限するあまり、目的とする対象物（例えば、地表面）の計測結果が記録されない場合があるという問題があった。そこで、このような問題を解決するため、新たにウェーブフォームデータ（以下、「ＷＦＤ」という。）方式によるレーザー計測が採用されるようになった。図７は、ＷＦＤ方式によるレーザー計測を示す説明図である。この図に示すように、ＷＦＤ方式によるレーザー計測では、全ての反射信号を波形（ウェーブフォーム）として記録することができる。したがって、従来方式のように地表面の計測結果の記録を漏らすことはなくなった。

ところがＷＦＤ方式によるレーザー計測は、従来方式とは別の問題を抱えている。ＷＦＤ方式が全ての反射信号を波形データとして記録するのは既に述べたとおりである。この全ての反射信号には、「ノイズ」と呼ばれる反射パルスではないものも含まれている。つまり、ＷＦＤ方式で得られた波形データの中からノイズを除去し、真の反射パルスのみを抽出する処理が必要となる。例えば図７では、反射パルスを６個（１ｓｔ〜６ｔｈパルス）だけ記録すべきところ、記録された波形データではノイズを含むため、この波形から真の反射パルスを抽出しなければならない。

ＷＦＤ方式で得られた波形データのうち真の反射パルスを抽出（推定）するには、人による目視等では限界があり、コンピュータ等による自動計算（あるいは半自動計算）で処理する必要がある。ところが、特許文献１のようにレーザー計測のフィルタリング（得られた点群データから地表面を推定する）技術についは数多く提案されているものの、ＷＦＤ方式で得られた波形データから反射パルスを的確に抽出（推定）する処理方法や処理プログラムに関しては、これまであまり提案されることはなかった。

特開２００４−２７２６８８号公報

ＷＦＤ方式で得られた波形データは、複数の極大値（上に凸の頂点）を持ち、規則性のない極めて乱れた形状を示す。このような波形データを分析する場合、近似する連続波形を推定することが考えられる。図８（ａ）は、実際に得られた反射強度の分布を示す波形であり、横軸が照射位置から計測点までの距離（時間）、縦軸が反射強度を表している。また、図８（ｂ）は、３つのガウス関数を推定して表したもので、横軸と縦軸は図８（ａ）と同じである。なお便宜上ここでは、実際に得られた反射強度の分布を示す波形を「実測波形（図８（ａ））」、この実測波形を近似するために推定する連続波形（例えば、ガウス関数）を「モデル波形（図８（ｂ））」という。

波形のうち極大値を示す位置が真の反射パルスの位置と推定できるが、図８（ａ）に示す実測波形では極大値の位置を明確に指定することができない。一方、図８（ｂ）に示すモデル波形では極大値の位置が明確であり、その結果、真の反射パルスの位置を容易に推定することができる。図８（ｂ）に示すモデル波形は３つのガウス関数で推定しているが、直接的にこの結果が得られるものではない。まずは、ガウス関数の数とおおよその中心を推定し、モデル波形を作成する。ガウス関数の数を指定したとき、最も実測波形を近似するモデル波形を得る手法の一つとしてＥＭアルゴリズムが知られている。ＥＭアルゴリズムは期待値最大化法とも呼ばれ、対数尤度の期待値を計算するＥステップと、尤度期待値を最大化するようなパラメータを設定し直すＭステップを繰り返し行い、最適解を導き出す解析手法である。しかしこの場合、ガウス関数の数や中心の初期値は人の判断が必要である。そこでまずは、一つのガウス関数として推定し、これと実測波形との誤差を算出し、この結果を踏まえて次に二つ目のガウス関数を推定し、再度、実測波形との誤差を算出する。このような処理を繰り返し行えば、最も実測波形と近似するケースをモデル波形として得ることができる。

実測波形を近似するモデル波形の推定手法は、ＥＭアルゴリズムに限らず種々の手法が知られている。また、用いる関数もガウス関数に限らず多項分布や、台形波、三角波等とすることもできる。しかしながら従来におけるモデル波形の推定手法は、いずれも実測波形との誤差にのみ着目するものであり、レーザー計測における地表面の反射パルスを抽出するには不向きな面があった。

図９（ａ）は、レーザー計測によって地形を計測した実測波形を示す図である。この実測波形では、図の左側に示すように最初に大きな反射パルス（１ｓｔパルス）が得られ、その後しばらくして地盤で反射したと思われる反射パルス（２ｎｄパルス）を取得している。一方の図９（ｂ）は、図９（ａ）の実測波形図に対し、ＥＭアルゴリズムによってモデル波形を推定した波形図である。この図からわかるように、推定したモデル波形では、地盤と思われるピークが消滅している。これは、実測波形との誤差にのみ着目して推定するが故に、大きなピークである１ｓｔパルスの影響を強く受けた結果である。すなわち、二つ目のガウス関数を推定する際、地盤と思われるピークは１ｓｔパルスから離れた位置にあってしかも反射強度が小さいため、その影響が極めて小さく評価されたわけである。

このように、従来からのモデル波形推定手法は、地表面の反射パルスの抽出漏れを生むことがあり、レーザー計測には不向きであった。また、既述のとおり、繰り返し解析して推定する手法では、一つ目のガウス関数の中心や追加するガウス関数の中心を人の判断や実測波形と推定波形との誤差によって指定することになるが、その指定によっては得られる収束解が異なり、その点では極めて不安定な解析方法であると指摘することができる。

本願発明の課題は従来手法が抱える問題を解決することであり、すなわちＷＦＤ方式によるレーザー計測で得られる波形データから、適切に反射パルスを推定する技術であって、しかもレーザー計測の特性を踏まえたうえで地表面からの反射パルスを推定する技術を提供することであり、具体的には反射パルス推定方法、及び反射パルス推定プログラムを提供することにある。

本願発明は、実測波形とモデル波形の反射強度差と、モデル波形の極大値位置（ピーク位置）と計測値（実測波形を構成する各点）との離れ（距離）に着目し、これらの積である「モーメント値」という新規な概念に着目して反射パルスを推定するという点に着目したものであり、従来にはなかった発想に基づいてなされた発明である。

本願発明の反射パルス推定方法は、レーザー計測によって得られる複数の計測値の中から、計測対象物に反射した反射パルスを推定する方法であり、実測波形生成工程、反射パルス設定工程、モデル波形作成工程、モーメント値算出工程、及び反射パルス推定工程を備えた方法である。実測波形生成工程は、複数の計測値に含まれる反射強度に基づいて実測波形を生成するものであり、反射パルス設定工程は、実測波形のうち極大値（ピーク）を示す位置付近に第１反射パルスを設定するものであり、モデル波形作成工程は、第１反射パルスに基づいて反射強度分布のモデルとなるモデル波形を作成するものである。また、モーメント値算出工程は、実測波形、第１反射パルス、及びモデル波形に基づいて計測値ごとにモーメント値を算出するもので、反射パルス推定工程は、モーメント値からモーメント分布を生成するとともに、このモーメント分布のうち極大値を示す位置に基づいて第２反射パルスを推定するものである。なお、モーメント値算出工程では、距離値（計測値の位置から第１反射パルスの位置までの距離の関数として求められる値）を算出するとともに、実測波形とモデル波形との差分により反射強度差を算出し、さらに距離値と反射強度差の積を求めることでモーメント値を算出する。

本願発明の反射パルス推定方法は、第３反射パルス、あるいは第４反射パルス以降を推定する方法とすることもできる。このとき、モデル波形作成工程、モーメント値算出工程、及び反射パルス推定工程が繰り返し行われ、２回目のモデル波形作成工程では、第１反射パルス及び第２反射パルスに基づいて反射強度分布のモデルとなるモデル波形を作成し、２回目のモーメント値算出工程では、実測波形、第１反射パルス、第２反射パルス、及びモデル波形に基づいて計測値ごとにモーメント値を算出する。そして２回目の反射パルス推定工程で、モーメント値からモーメント分布を生成するとともに、このモーメント分布のうち極大値を示す位置に基づいて第３反射パルスを推定する。さらに選択的に、モデル波形作成工程、モーメント値算出工程、及び反射パルス推定工程を繰り返し行うことで、第４反射パルス以降を推定することもできる。なおこの場合、モーメント値算出工程で算出される距離値は、複数の反射パルスのうち当該計測値の位置から最も近い反射パルスの位置に基づいて算出される。このとき、ガウス関数の標準偏差を用いた正規化距離を用いることもできる。

本願発明の反射パルス推定方法は、さらに、合計モーメント値算出工程と収束判定工程を備えた方法とすることもできる。合計モーメント値算出工程は、計測値ごとのモーメント値を総和することで合計モーメント値を求めるもので、収束判定工程は、合計モーメント値と収束閾値を比較し、合計モーメント値が収束閾値以下（又は収束閾値未満）となる場合に、反射パルス推定の終了を判定するものである。収束判定工程は、反射パルス推定工程の後に行う。なお、反射パルス推定の終了を判定しない場合は、さらに、モデル波形作成工程、モーメント値算出工程、及び反射パルス推定工程が繰り返し行われる。

本願発明の反射パルス推定方法は、正規化したモーメント値を用いた方法とすることもできる。この正規化したモーメント値は、距離値と反射強度差を乗算した値を、さらに複数の計測値の反射強度のうち最も大きい値の反射強度で除すことで得られる。または、実測波形の積分値で除することもできる。

本願発明の反射パルス推定プログラムは、レーザー計測によって得られる複数の計測値の中から、計測対象物に反射した反射パルスを推定する処理を、コンピュータに実行させるものであり、実測波形生成処理、反射パルス設定処理、モデル波形作成処理、モーメント値算出処理、及び反射パルス推定処理を備えたプログラムである。実測波形生成処理は、複数の計測値に含まれる反射強度に基づいて実測波形を生成するものであり、反射パルス設定処理は、実測波形のうち極大値（ピーク）を示す位置付近を指定することで、第１反射パルスを生成するものであり、モデル波形作成処理は、第１反射パルスに基づいて反射強度分布のモデルとなるモデル波形を作成するものである。また、モーメント値算出処理は、実測波形、第１反射パルス、及びモデル波形に基づいて計測値ごとにモーメント値を算出するもので、反射パルス推定処理は、モーメント値からモーメント分布を生成するとともに、このモーメント分布のうち極大値を示す位置に基づいて第２反射パルスを推定するものである。なお、モーメント値算出処理では、距離値（計測値の位置から第１反射パルスの位置までの距離の関数として求められる値）を算出するとともに、実測波形とモデル波形との差分により反射強度差を算出し、さらに距離値と反射強度差の積を求めることでモーメント値を算出する。

本願発明の反射パルス推定プログラムは、第３反射パルス、あるいは第４反射パルス以降を推定するプログラムとすることもできる。このとき、モデル波形作成処理、モーメント値算出処理、及び反射パルス推定処理が繰り返し行われ、２回目のモデル波形作成処理では、第１反射パルス及び第２反射パルスに基づいて反射強度分布のモデルとなるモデル波形を作成し、２回目のモーメント値算出処理では、実測波形、第１反射パルス、第２反射パルス、及びモデル波形に基づいて計測値ごとにモーメント値を算出する。そして２回目の反射パルス推定処理で、モーメント値からモーメント分布を生成するとともに、このモーメント分布のうち極大値を示す位置に基づいて第３反射パルスを推定する。さらに選択的に、モデル波形作成処理、モーメント値算出処理、及び反射パルス推定処理を繰り返し行うことで、第４反射パルス以降を推定することもできる。なおこの場合、モーメント値算出処理で算出される距離値は、複数の反射パルスのうち当該計測値の位置から最も近い反射パルスの位置に基づいて算出される。

本願発明の反射パルス推定プログラムは、さらに、合計モーメント値算出処理と収束判定処理を備えたプログラムとすることもできる。合計モーメント値算出処理は、計測値ごとのモーメント値を総和することで合計モーメント値を求めるもので、収束判定処理は、合計モーメント値と収束閾値を比較し、合計モーメント値が収束閾値以下又は未満となる場合に、反射パルス推定の終了を判定するものである。収束判定処理は、反射パルス推定処理の後に行う。なお、反射パルス推定の終了を判定しない場合は、さらに、モデル波形作成処理、モーメント値算出処理、及び反射パルス推定処理が繰り返し行われる。

本願発明の反射パルス推定プログラムは、正規化したモーメント値を用いたプログラムとすることもできる。この正規化したモーメント値は、距離値と反射強度差を乗算した値を、さらに複数の計測値の反射強度のうち最も大きい値の反射強度で除すことで得られる。

本願発明の反射パルス推定方法、及び反射パルス推定プログラムには、次のような効果がある。
（１）ＷＦＤ方式で得られた波形データを利用することで、従来方式では抽出漏れの可能性がある地表面の反射パルスも、的確に抽出することができる。
（２）モーメント値を採用することで、レーザー計測の特性も踏まえたうえで、より的確に地表面の反射パルスを抽出することができる。
（３）合計モーメント値と収束閾値を比較することで、繰り返し計算の収束を自動的に判断することができる。
（４）モーメント値を正規化することで、反射強度の最大値と最小値が極端に異なる場合でも、適切に処理することができる。

航空レーザー計測の実施状況を示す説明図。 本願発明の主な流れを示すフロー図。 （ａ）は、計測値によって構成される実測波形を示すモデル図、（ｂ）は第１反射パルスに基づいて作成されたモデル波形を示すモデル図、（ｃ）は第１反射パルスと第２反射パルスに基づいて作成されたモデル波形を示すモデル図、（ｄ）は第１反射パルス〜第３反射パルスに基づいて作成されたモデル波形を示すモデル図。 繰り返し処理を行う場合の本願発明の主な流れを示すフロー図。 ３個のガウス関数がある場合の、ガウス分布中心と計測値との距離を表すモデル図。 従来方式によるレーザー計測を示す説明図。 ＷＦＤ方式によるレーザー計測を示す説明図。 （ａ）は実際に得られた計測値の分布を示す波形図、（ｂ）は推定した３つのガウス関数を示す波形図。 （ａ）はレーザー計測によって地形を計測した結果を表す実測波形図、（ｂ）は（ａ）の実測波形図に対して推定したモデル波形図。

本願発明の反射パルス推定方法、及び反射パルス推定プログラムの実施形態の一例を、図に基づいて説明する。

（全体概要）
本願発明は、レーザー計測によって得られたウェーブフォームデータから、反射パルス（対象物に反射した信号）を抽出するものであり、具体的には、プログラムによって必要な処理をコンピュータに実行させることで実施される。このコンピュータ装置は、ＣＰＵ等のプロセッサ、ＲＯＭやＲＡＭといったメモリを具備しており、さらにマウスやキーボード等の入力手段やディスプレイを含むこともあり、パーソナルコンピュータ（ＰＣ）や、ｉＰａｄ（登録商標）といったタブレットＰＣ、あるいはＰＤＡ（Ｐｅｒｓｏｎａｌ Ｄａｔａ Ａｓｓｉｓｔａｎｃｅ）などによって構成される。

図２は、本願発明の主な流れを示すフロー図である。この図に従って、本願発明の概要を説明する。まず、レーザー計測で取得した結果から実測波形を生成する。既述のとおりレーザー計測では、センサが受信した信号に基づいて計測点の３次元座標値が求められる。この受信信号には、反射強度、受信時刻などが含まれる。なお、レーザーの進行速度は一定と考えることができるので、照射時刻と受信時刻との時間差を、照射地点から計測点（つまりレーザーが反射した地点）までの距離に換算することができる。

実測波形は、照射地点からの距離と反射強度に基づいて生成されるもので、具体的には、横軸を照射地点からの距離とし、縦軸を反射強度として表した分布図であり、これがいわゆるウェーブフォームである。図３（ａ）に、実測波形の例を示す。レーザー計測では多数のレーザーパルスが照射されることは既に述べたとおりであるが、実測波形は一つのレーザーパルスに対して生成される。例えば、８ビットとして処理する場合、横軸が２５６に分割され、一つのレーザーパルスで得た信号を基にそれぞれの分割点に対して反射強度を与える。これによって、実測波形が作成される。つまり、実測波形は多数（８ビットの場合２５６個）の点によって構成されるのであり、ここでは便宜上、実測波形を構成する個々の点を「計測値」ということとする。言い換えれば、計測値は、一つのレーザーパルスに対して照射地点からの所定距離ごとに反射強度を与えたもので、距離（時間差）と反射強度の組み合わせと言える。

次に、実測波形に基づいて第１反射パルスを設定し、さらにこの第１反射パルスに基づいてモデル波形を作成する。第１反射パルスの設定とは、実測波形に描かれる計測値のうち極大値（上に凸の頂点）の中心付近に第１反射パルスの位置を入力することである。通常、実測波形に描かれる極大値は複数存在するので、そのうち最も大きな値（反射強度）を示す極大値を選択して第１反射パルス位置を入力する。

第１反射パルスの位置が定まると、これに基づいてモデル波形を作成する。モデル波形とは、実測波形を近似するもので、いわば反射強度分布のモデルとなるものであり、所定の関数を利用して作成する。この関数としては、ガウス関数や正規分布、多項分布、あるいは台形波や三角波など所望の関数を採用することができるが、ここでは便宜上、ガウス関数を利用した場合で説明することとする。

実測波形とモデル波形が作成されると、実測波形を構成する計測値ごとにモーメント値を算出する。実測波形は、離散型データ（計測値）を結線したものであり、一方のモデル波形は連続的な線形である。したがって、モデル波形も実測波形と同様に横軸方向に分割し、それぞれで反射強度を求めることとする。これによって、計測値ごとに実測波形とモデル波形の対応が可能となり、計測値ごとにモーメント値を算出することができる。なお、モーメント値については後に詳しく説明する。

計測値ごとのモーメント値が求められると、モーメント値分布図を作成する。このモーメント値分布図は、横軸を照射地点からの距離、縦軸をモーメント値として表すもので、横軸における分割は当然ながら実測波形の計測値と同一である。モーメント値はいわば実測波形とモデル波形の相違の程度を表すものであり、相当の相違がある区間には、新たに別のガウス関数を設定してモデル波形を作成することを示している。つまり、２つのガウス関数によるモデル波形で、実測波形を近似することを示しているわけである。２つ目のガウス関数を設定する位置（すなわちガウス分布の中心）は、実測波形とモデル波形の相違が顕著な区間とすべきであり、したがって波形として描かれるモーメント値分布図のうち極大値を示す位置に設定される。ここで設定された２つ目のガウス関数の中心が、第２の反射強度のピーク（極大値）と推定され、これをもって第２の反射パルス（つまり２ｎｄパルス）と推定することができる。

以下、本願発明の反射パルス推定方法、及び反射パルス推定プログラムを、構成する要素ごとに詳述する。

（実測波形の生成）
実測波形を作成するには、図２のフロー図に示すように、まず計測値を読み出す（Ｓ０１）。ここで計測値について詳しく説明する。既述のとおり一つのレーザーパルスに対して、ノイズを含めると多数の反射信号（リターン記録）が受信され、記録される。これらの受信信号は、所定の時間間隔（例えば１ｎｓ）で反射強度とともに記録され、例えば８ビットとして処理する場合は２５６の信号が記録される。この記録された個々の信号が、すなわち計測値である。なお、ここでは一つのレーザーに対して得られる計測値の個数をｎ個として説明する。

読み出した計測値を用いて実測波形を生成する（Ｓ０２）。計測値が持つ受信時刻は、照射地点からの距離と考えることができるので、横軸を照射地点からの距離、縦軸を反射強度とするグラフを設定する。このグラフ上に、ｎ個の計測値を配点して結線すれば、実測波形を生成することができる。図３（ａ）は、計測値によって構成される実測波形を示すモデル図である。ここで行われる内容が、反射パルス推定方法における「実測波形生成工程」であり、反射パルス推定プログラムにおける「実測波形生成処理」である。

実測波形が生成されると、反射パルスの初期値、すなわち第１反射パルスが設定される（Ｓ０３）。既に述べたとおり、実測波形に描かれる計測値のうち最も大きな極大値の中心付近を選び、この位置を入力する。ここで入力される位置は、実測波形を目視することで人の判断によって指定される。ここで行われる内容が、反射パルス推定方法における「反射パルス設定工程」であり、反射パルス推定プログラムにおける「反射パルス設定処理」である。

第１反射パルスが設定されると、第１反射パルスの位置を分布中心として、実測波形を近似するガウス関数が作成される（Ｓ０４）。この時点では、第１反射パルスの位置を分布中心としたガウス関数が、そのままモデル波形となる。図３（ｂ）は、第１反射パルスに基づいて作成されたモデル波形を示すモデル図である。ここで行われる内容が、反射パルス推定方法における「モデル波形作成工程」であり、反射パルス推定プログラムにおける「モデル波形作成処理」である。

実測波形及びモデル波形が作成されると、モーメント値が算出される（Ｓ０５）。このモーメント値は計測値ごとに算出され、ｎ個の計測値のうちｉ番目の計測値におけるモーメント値はＭ_ｉとして表現できる。このモーメント値Ｍ_ｉは、ｉ番目の計測値における反射強度差Ｅ_ｉと、同じくｉ番目の計測値における距離値Ｄ_ｉとの積によって求められる。

反射強度差Ｅ_ｉとは、ｉ番目の計測値における実測波形の反射強度と、モデル波形の反射強度との差である。既述のとおり、計測値の横軸値（照射地点からの距離）に対応するようにモデル分布から反射強度を取得することで、計測値ごとに実測波形との比較が可能となる。ｎ個のうちｉ番目の計測値の反射強度をＥｎ_ｉとし、この計測値に相当するモデル波形の反射強度をＥｍ_ｉとすれば、反射強度差Ｅ_ｉは次式で求められる。

距離値Ｄ_ｉとは、ｉ番目の計測値におけるモデル波形を構成するガウス関数の分布中心までの距離に応じて求められる値である。すなわち、ガウス関数の分布中心の横軸値をｕ_ｊとし、ｉ番目の計測値の横軸値をｘ_ｉとすれば、距離値Ｄ_ｉは関数ｆ（ｕ_ｊ−ｘ_ｉ）として、次式のように表すことができる。なお、ｕ_ｊの添え字ｊは、複数のガウス関数によってモデル波形が形成される場合に、ｊ番目のガウス関数の分布中心横軸値を意味するためのものである。また式中に示すｐはべき乗を示すもので任意の正数（整数に限らない）を与えることが可能で、例えばｐ＝１、ｐ＝２、ｐ＝３などが与えられる。なお、式中の「ａｂｓ」は絶対値を求める意味である。

ｉ番目の計測値における、反射強度差Ｅ_ｉと距離値Ｄ_ｉが得られれば、モーメント値Ｍ_ｉは次式によって求めることができる。

ｎ個の反射強度のうち、その最大値と最小値が極端に異なる場合、適切に処理されない場合も考えられる。このような事態を回避するため、モーメント値Ｍ_ｉを正規化することもできる。この正規化は、次式のとおり反射強度の最大値Ｅｎ_ｍａｘで除すことによって行われる。

ｉ番目の計測値における、反射強度差Ｅ_ｉと距離値Ｄ_ｉを求め、モーメント値Ｍ_ｉを算出し、選択的にモーメント値Ｍ_ｉを正規化することが、反射パルス推定方法における「モーメント値算出工程」であり、反射パルス推定プログラムにおける「モーメント値算出処理」である。

計測値ごとにモーメント値Ｍ_ｉが算出できると、入力した初期値以外の反射パルス（ここでは第２反射パルス）が推定される（Ｓ０６）。第２反射パルスを推定するには、モーメント値分布図を作成する。モーメント値分布図は、横軸を照射地点からの距離、縦軸をモーメント値として表すもので、横軸には、実測波形の計測値の横軸値ｘ_ｉと同一の横軸値が設定される。この横軸値ｘ_ｉそれぞれに対してモーメント値を配点し、これらを結線したものがモーメント値分布図である。

このモーメント値分布図は波形として描かれ、このうち極大値を示す位置に第２反射パルスが設定される。第２反射パルスが設定されれば、この位置を分布中心とする第２のガウス関数を作成することができ、すなわち第１反射パルスに基づく第１のガウス関数とあわせて（２つのガウス関数で）モデル波形を作成することができる。図３（ｃ）は、第２反射パルスが設定され、第１反射パルスと第２反射パルスに基づいて作成されたモデル波形を示すモデル図である。こうして第２反射パルスの位置が特定できれば、その位置ｘ_ｉにおける計測値から第２反射パルスの反射強度、あるいはその計測点の座標値を求めることができるわけである。ここで行われる内容が、反射パルス推定方法における「反射パルス推定工程」であり、反射パルス推定プログラムにおける「反射パルス推定処理」である。

（繰り返し処理）
ここまで、第２反射パルスの位置を特定し、第１のガウス関数と第２のガウス関数で実測波形を最も近似できるとは限らない。そこで、既述したＥＭアルゴリズムのように、第３反射パルス、第４反射パルス、第５反射パルスというように、繰り返し処理を行って複数のモデル波形を作成し、その中から最も妥当なものをモデル波形として、反射パルスを推定することもできる。

図４は、繰り返し処理を行う場合の本願発明の主な流れを示すフロー図である。この図に示すように、反射パルスの設定（Ｓ０３）〜反射パルスの推定（Ｓ０６）を繰り返し行う。この場合、反射パルスの設定（Ｓ０３）では、これまでに得られた反射パルスを全て入力として扱う。例えば、２回目（ｋ＝２）の処理であれば、それまでに第１反射パルス（初期値）と第２反射パルス（１回目の処理で推定）が得られているので、これらの位置を設定し、２つのガウス関数を作成してモデル波形とする。その結果、２回目の処理では第３反射パルスが得られ、第１反射パルス〜第３反射パルスに基づいてモデル波形が作成することができる。図３（ｄ）は第１反射パルス〜第３反射パルスに基づいて作成されたモデル波形を示すモデル図である。同様に３回目（ｋ＝３）の処理であれば、第１〜第３反射パルスの位置を設定し、３つのガウス関数を作成してモデル波形とし、第４反射パルスが得られる。このように複数のモデル波形を作成し、所定の閾値判定を行うことで、最も実測波形を近似しているものをモデル波形として確定し（Ｓ０９）、その中に含まれる複数の反射パルスを計測点（例えば枝や草を示す）として推定する（Ｓ１０）こともできる。

繰り返し処理を行う場合、距離値Ｄ_ｉを求める際に、複数のガウス関数を対象にすることになる。つまり、ｒ個のガウス関数があれば、ｕ_ｊ−ｘ_ｉ（ｊ＝１〜ｒ）はｒ個だけ求められることになる。図５は、３個のガウス関数がある場合の、ガウス分布中心と計測値との距離表すモデル図である。このとき、できるだけノイズを排除することを目的に、最も近いガウス関数を選択してｕ_ｊ−ｘ_ｉを求める。図５に示す計測値ｘ_ｉの場合では、第３のガウス関数の分布中心であるｕ_３が選択され、Ｄ_ｉは関数ｆ（ｕ_３−ｘ_ｉ）として求められる。

（収束判定）
繰り返し処理を行う場合、あらかじめ所定の収束条件を定めることで、自動的に収束判定できて便宜である。そこで、本願発明では合計モーメント値に基づく収束判定を行うことができる。

図４に示すように、繰り返し計算を行うたびに合計モーメント値ΣＭ_ｉを求める（Ｓ０７）。この合計モーメント値ΣＭ_ｉは、計測値ごとのモーメント値Ｍ_ｉを総和したもので、次式により求められる。

ここで行われる内容が、反射パルス推定方法における「合計モーメント値算出工程」であり、反射パルス推定プログラムにおける「合計モーメント値算出処理」である。

合計モーメント値ΣＭ_ｉが求められると、収束判定を行う（Ｓ０８）。ここでは、あらかじめ定められた収束閾値σ_ｅと合計モーメント値ΣＭ_ｉとを比較することで判定され、合計モーメント値ΣＭ_ｉが収束閾値σ_ｅ以下となる場合、あるいは合計モーメント値ΣＭ_ｉが収束閾値σ_ｅ未満となる場合に、収束したとして一連の反射パルス推定処理を終了すると判定する。終了判定された場合は、それ以上の繰り返し処理を行わず、その時点のモデル波形を確定し（Ｓ０９）、その中に含まれる反射パルスを抽出する（Ｓ１０）。

一方、合計モーメント値ΣＭ_ｉが収束閾値σ_ｅを超えている場合、あるいは合計モーメント値ΣＭ_ｉが収束閾値σ_ｅ以上となる場合は、まだ収束しないとして終了判定は行わず、反射パルス設定（Ｓ０３）以降の処理が再び繰り返される。

ここで行われる内容が、反射パルス推定方法における「収束判定工程」であり、反射パルス推定プログラムにおける「収束判定処理」である。

本願発明の反射パルス推定方法、及び反射パルス推定プログラムは、航空レーザー計測のほか、車載型のいわゆるモバイルマッピングシステム（Ｍｏｂｉｌｅ Ｍａｐｐｉｎｇ Ｓｙｓｔｅｍ：ＭＭＳ）でのレーザー計測、固定式のレーザー計測など、様々なレーザー計測に利用できる。また、本願発明を用いたレーザー計測を行い２時期の地形を比較することで、経年の地殻変動に伴う地表面変化が把握できるとともに、断層活動の活動状況や地すべりの活動状況が把握できるので、自然災害を未然に防ぎあるいは被害を軽減させることが可能であり、産業上利用できるとともに、社会的に大きな貢献を期待し得る発明である。

１ 地形
２ 航空機
３ レーザーパルス

Claims (8)

1. レーザー計測によって得られる複数の計測値の中から、計測対象物に反射した反射パルスを推定する反射パルス推定方法であって、
   前記複数の計測値に含まれる反射強度に基づいて、実測波形を生成する実測波形生成工程と、
   前記実測波形のうち極大値を示す位置付近に、第１反射パルスを設定する反射パルス設定工程と、
   前記第１反射パルスに基づいて、反射強度分布のモデルとなるモデル波形を作成するモデル波形作成工程と、
   前記実測波形、前記第１反射パルス、及び前記モデル波形に基づいて、実測波形を構成する計測値ごとにモーメント値を算出するモーメント値算出工程と、
   前記モーメント値からモーメント分布を生成するとともに、該モーメント分布のうち極大値を示す位置に基づいて、第２反射パルスを推定する反射パルス推定工程と、を備え、
   前記モーメント値算出工程では、計測値の位置から前記第１反射パルスの位置までの距離の関数として求められる距離値を算出するとともに、前記実測波形と前記モデル波形との差分により反射強度差を算出し、さらに該距離値と該反射強度差を乗算することで前記モーメント値を算出する、ことを特徴とする反射パルス推定方法。
2. 前記モデル波形作成工程、前記モーメント値算出工程、及び前記反射パルス推定工程が繰り返し行われ、
   ２回目の前記モデル波形作成工程では、前記第１反射パルス及び前記第２反射パルスに基づいて、反射強度分布のモデルとなるモデル波形を作成し、
   ２回目の前記モーメント値算出工程では、前記実測波形、前記第１反射パルス、前記第２反射パルス、及び前記モデル波形に基づいて、実測波形を構成する計測値ごとにモーメント値を算出し、
   ２回目の前記反射パルス推定工程では、前記モーメント値からモーメント分布を生成するとともに、該モーメント分布のうち極大値を示す位置に基づいて、第３反射パルスを推定し、
   さらに選択的に、前記モデル波形作成工程、前記モーメント値算出工程、及び前記反射パルス推定工程を繰り返し行うことで、第４反射パルス以降を推定し、
   前記モーメント値算出工程で算出される前記距離値は、複数の反射パルスのうち当該計測値の位置から最も近い反射パルスの位置に基づいて算出される、ことを特徴とする請求項１記載の反射パルス推定方法。
3. 計測値ごとの前記モーメント値を総和することで合計モーメント値を求める合計モーメント値算出工程と、
   前記合計モーメント値と収束閾値を比較し、前記合計モーメント値が収束閾値以下又は未満となる場合に、反射パルス推定の終了を判定する収束判定工程と、を備え、
   前記収束判定工程は、前記反射パルス推定工程の後に行い、反射パルス推定の終了を判定しない場合は、さらに、前記モデル波形作成工程、前記モーメント値算出工程、及び前記反射パルス推定工程を繰り返し行う、ことを特徴とする請求項２記載の反射パルス推定方法。
4. 前記モーメント値算出工程では、前記距離値と前記反射強度差を乗算した値を、さらに複数の計測値の反射強度のうち最も大きい値の反射強度で除すことで前記モーメント値を算出する、ことを特徴とする請求項１乃至請求項３のいずれかに記載の反射パルス推定方法。
5. レーザー計測によって得られる複数の計測値の中から、計測対象物に反射した反射パルスを推定する処理を、コンピュータに実行させる反射パルス推定プログラムであって、
   前記複数の計測値に含まれる反射強度に基づいて、実測波形を生成する実測波形生成処理と、
   前記実測波形のうち極大値を示す位置付近を指定することで、第１反射パルスを生成する反射パルス設定処理と、
   前記第１反射パルスに基づいて、反射強度分布のモデルとなるモデル波形を作成するモデル波形作成処理と、
   前記実測波形、前記第１反射パルス、及び前記モデル波形に基づいて、実測波形を構成する計測値ごとにモーメント値を算出するモーメント値算出処理と、
   前記モーメント値からモーメント分布を生成するとともに、該モーメント分布のうち極大値を示す位置に基づいて、第２反射パルスを推定する反射パルス推定処理と、を前記コンピュータに実行させる機能を備え、
   前記モーメント値算出処理では、計測値の位置から前記第１反射パルスの位置までの距離の関数として求められる距離値が算出されるとともに、前記実測波形と前記モデル波形との差分により反射強度差が算出され、該距離値と該反射強度差の乗算により前記モーメント値が算出される、ことを特徴とする反射パルス推定プログラム。
6. 前記モデル波形作成処理、前記モーメント値算出処理、及び前記反射パルス推定処理が繰り返し行われ、
   ２回目の前記モデル波形作成処理では、前記第１反射パルス及び前記第２反射パルスに基づいて、反射強度分布のモデルとなるモデル波形が作成され、
   ２回目の前記モーメント値算出処理では、前記実測波形、前記第１反射パルス、前記第２反射パルス、及び前記モデル波形に基づいて、実測波形を構成する計測値ごとにモーメント値が算出され、
   ２回目の前記反射パルス推定処理では、前記モーメント値からモーメント分布が生成されるとともに、該モーメント分布のうち極大値を示す位置に基づいて、第３反射パルスが推定され、
   さらに選択的に、前記モデル波形作成処理、前記モーメント値算出処理、及び前記反射パルス推定処理を繰り返し行うことで、第４反射パルス以降を推定し、
   前記モーメント値算出処理で算出される前記距離値は、複数の反射パルスのうち当該計測値の位置から最も近い反射パルスの位置に基づいて算出される、ことを特徴とする請求項５記載の反射パルス推定プログラム。
7. 計測値ごとの前記モーメント値を総和することで合計モーメント値を求める合計モーメント値算出処理と、
   前記合計モーメント値と収束閾値を比較し、前記合計モーメント値が収束閾値以下又は未満となる場合に、反射パルス推定の終了を判定する収束判定処理と、を前記コンピュータに実行させる機能を備え、
   前記収束判定処理は、前記反射パルス推定処理の後に行われ、反射パルス推定の終了が判定されない場合は、さらに、前記モデル波形作成処理、前記モーメント値算出処理、及び前記反射パルス推定処理が繰り返し行われる、ことを特徴とする請求項６記載の反射パルス推定プログラム。
8. 前記モーメント値算出処理では、前記距離値と前記反射強度差を乗算した値を、さらに複数の計測値の反射強度うち最も大きい値の反射強度で除すことで前記モーメント値が算出される、ことを特徴とする請求項５乃至請求項７のいずれかに記載の反射パルス推定プログラム。

Images