JP4226360B2 - レーザデータのフィルタリング方法及びプログラム - Google Patents
レーザデータのフィルタリング方法及びプログラム Download PDFInfo
- Publication number
- JP4226360B2 JP4226360B2 JP2003063841A JP2003063841A JP4226360B2 JP 4226360 B2 JP4226360 B2 JP 4226360B2 JP 2003063841 A JP2003063841 A JP 2003063841A JP 2003063841 A JP2003063841 A JP 2003063841A JP 4226360 B2 JP4226360 B2 JP 4226360B2
- Authority
- JP
- Japan
- Prior art keywords
- mesh
- laser data
- cluster
- ground surface
- laser
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Instructional Devices (AREA)
- Processing Or Creating Images (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、レーザデータのフィルタリング方法及びプログラムに関するものである。
【０００２】
【従来の技術】
航空機レーザスキャナ等により取得したレーザデータから地表面データ、あるいはそれ以外の地物データを抽出することは、例えば、ディジタル標高モデル（ＤＴＭ：Digital Terrain Model）を作成する場合に必要であり、このような方法としては、例えば、特許文献１に記載されたものが知られている。
【０００３】
この従来例において、取得されたレーザデータ群はクラスタリングにより類似した特徴を持つグループ（クラスタ）に分類される。クラスタリングのために設定される特徴空間には、各レーザデータの平面位置、標高を要素とする３次元空間がそのまま利用される。
【０００４】
また、クラスタリングに際して定義される近傍にはボロノイ近傍が使用され、レーザデータ群が形成する特徴空間は、先ず、ボロノイ分割により各レーザデータの近傍を定義し、探索される。
【０００５】
【特許文献１】
特開２００２-３６６９７７号公報
【０００６】
【発明が解決しようとする課題】
しかし、上述した従来例において、クラスタリングに際する近傍点の探索は、全レーザデータに対して行う必要があるために、膨大な作業量を必要とするという欠点がある。
【０００７】
本発明は、以上の欠点を解消すべくなされたものであって、処理の効率化を図ることのできるレーザデータのフィルタリング方法及びプログラムを提供することを目的とする。
【０００８】
また、本発明の他の目的は、このようにしてフィルタリングした結果の表示方法の提供にある。
【０００９】
【課題を解決するための手段】
本発明によれば、上記目的は、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法を提供することにより達成される。
【００１０】
航空機レーザ測量で上空からのレーザ掃射により取得されるレーザデータは、平面位置（Ｘ，Ｙ）情報及び標高（Ｚ）情報を有するポイントデータであり、各レーザデータ単位では地表面等の識別を行うことはできない。このため、地表面等の識別には、隣接するレーザデータとの位置関係を利用することが必要であり、レーザデータ群に対する種々の処理が必要になる。
【００１１】
そして、レーザデータ群に対する処理は、上述した従来例をはじめとして、一般に各レーザデータを平等に扱い、これらに対する統計的処理として行われていた。
【００１２】
一方、レーザデータを種々の目的によって分類、抽出するレーザデータのフィルタリングに当たって、先ず、地表面データの抽出を優先することによって、地表面であることの特殊性をレーザデータ群の処理対象量の軽減に利用する可能性が生じる。すなわち、レーザデータ群に設定された好適なサイズのメッシュ内における最低標高点は、地表面である確率が高いために、これら最低標高点に処理対象を限定して処理対象量を軽減した方が、地表同士を検索の対象とするので、地表抽出精度は高くなる。
【００１３】
本発明は、以上の認識の下になされたものであって、メッシュ化及びメッシュ内でのメッシュ代表値の設定を行うことによって処理対象量を大幅に軽減させるとともに、クラスタリングに伴う近接点の検索を隣接メッシュに限定することによって処理効率を高め、かつ、正確な地表面抽出を可能にするものである。
【００１４】
すなわち、本発明において、先ず、レーザデータ群に対して平面座標に基づくメッシュを設定するステップ（Ｓ２）が実行される。メッシュサイズは、レーザデータが含まれない空メッシュが発生することのないように十分な大きさで、かつ、過度に周囲の条件に平均化されてメッシュの個性が失われない程度の小ささであることが必要であり、レーザデータの密度、正確にはポイントデータ間の間隔を考慮して、例えば、ポイントデータの間隔の２倍程度に決定される。
【００１５】
次いで、上記各メッシュ内のレーザデータのうち、最低標高を有するものをメッシュ代表点１として抽出するステップ（Ｓ３）が実行された後、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングするステップ（Ｓ４）が実行される。メッシュ代表点１を最低標高点データとすることにより、車、樹木等を予め除去することも可能になる。
【００１６】
クラスタリング操作は、標高要素を判断基準として隣接するメッシュとの間で行われる。近傍は、注目メッシュを包囲するように隣接する８個であっても、注目メッシュに対して縦横に隣接する４個であってもよいが、メッシュ代表点１が最低標高点により定義されているために、縦横の４個をとれば十分な精度を維持可能であり、かつ、処理量も少なくてすむ。
【００１７】
また、標高要素による判断は、例えば、隣接するメッシュ代表点１間の比高差をとることも可能であるが、注目メッシュのメッシュ代表点１との間の傾斜角度が所定値以下であるメッシュ代表点１の数を基準に行うことができる。
【００１８】
メッシュ代表点１は、メッシュ内のレーザデータから抽出されたものであり、Ｘ，Ｙ値及びＺ値をそのまま保持しているために、クラスタリングの基準として傾斜角度を利用すると、メッシュ全体を同一標高として考えることとなる比高差を利用する場合に比して、より現実の状態に即した判断が可能になる。
【００１９】
閾値として利用される傾斜角度は、経験則上、水平距離５（ｍ）に対する標高差が０．５（ｍ）、すなわち、５．７°程度に設定するのが望ましい。
【００２０】
以上のクラスタリングステップの後に、構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するステップ、すなわち、クラスタ内のメッシュ数による地表判定ステップ（Ｓ５）が実行される。
【００２１】
最大建物想定面積は、通常想定される建物の最大面積であり、通常、５０００ｍ2程度に設定しておけば足りる。最大建物想定面積に相当するメッシュ数、すなわち、最大建物想定面積をＳ、１メッシュの面積をｓ、クラスタ内メッシュ数をＭとした場合、
Ｍ＞Ｓ／ｓ
の条件を満たすクラスタは、全体に比較的なだらかな傾斜が建物と想定される場合の面積を超えて連続していることとなるために、非常に高い精度で地表部分であることを推定できる。
【００２２】
また、以上のフィルタリング方法は、なだらかで広い平野部の抽出に適しているが、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数未満、最小建物想定面積以上のクラスタ内のメッシュ代表点１の平面位置及び標高値の統計的代表値を平面位置及び標高値とするクラスタ代表点２を設定し、
次いで、前記クラスタ代表点２との間の傾斜角度が所定値未満である前記クラスタの近傍における地表面反射レーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法により、クラスタの面積により直ちに地表面と判断できない場合の地表抽出方法が提供される。
【００２３】
この方法において、クラスタリングまでのステップ（Ｓ２〜Ｓ４）は上述した方法と全く同様に行われ、この後、構成メッシュ数が最大建物想定面積相当数未満、最小建物想定面積以上のクラスタを検索するステップ（Ｓ６１）が実行され、次いで、該当クラスタ内にクラスタ代表点２を設定するステップ（Ｓ６２）が実行される。
【００２４】
クラスタ代表点２は、該当クラスタ内の複数のメッシュ代表点１を代表する平面位置及び標高値を有する点であり、メッシュ代表点１に対して適宜の統計的処理を施した統計的代表値として与えられる。統計的代表値を求めるためには、平均、あるいはメディアン値等を使用できる。
【００２５】
以上のようにして求められたクラスタ代表点２は、該当クラスを代表する平面位置及び標高情報を有しており、次いで、これらクラスタ代表点２とクラスタ近傍における地表面反射レーザデータとの間の傾斜角度を求めるステップ（Ｓ６３）が実行される。
【００２６】
傾斜角度を求めるに当たって、先ず、クラスタ代表点２から、適宜の手段によって地表面と判断されたレーザデータとの間で形成されるＴＩＮ（Triagulated Irregular Network）ネットワーク上の各点を近接点として定義することができる。また、この方法は、上述したメッシュ数による地表面抽出方法に続いて行うことが可能であり、この場合、傾斜角度は、メッシュ代表点１を対象として算出される。
【００２７】
該当クラスタの地表面判定は、以上のようにしてクラスタ代表点２と近傍の地表面反射レーザデータの間の傾斜角度を求めた後、傾斜角度が所定値未満の地表面反射データが全測定対象の地表面反射データに占める割合を求めるステップ（Ｓ６４）を実行することにより行われる。
【００２８】
この方法における目的は、建物が地表面からの判断が曖昧なクラスタからの地表面の判定であるために、判断の閾値となる傾斜角度は、地表面の想定される最大傾斜角度を基準に設定することができ、例えば、３０°に設定される。
【００２９】
また、割合の閾値は、判定する地形等を考慮して判定精度が高くなるように適宜決定可能であるが、１／３程度にするのが経験則上望ましく、閾値未満の傾斜角度の割合が上記閾値以上である場合には、当該クラスタは地表面と判定される。
【００３０】
このように、判定に当たってクラスタ内の各点を対象とすることなく、クラスタ代表点２を基準にクラスタ全体を判定することにより、計算量を少なくすることができる。
【００３１】
さらに、以上のフィルタリング方法に加え、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最小建物想定面積相当数未満のクラスタ内で、かつ、セカンドパルス以降のレーザデータからなるメッシュ代表点１との間の傾斜角度が所定値未満である前記メッシュ代表点１の近傍における地表面反射レーザデータ及び他のセカンドパルス以降のレーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法を採用することにより、樹木等の多い山間部での地表面の抽出精度を高めることができる。
【００３２】
この方法は、上述した方法において、構成メッシュ数が最小建物想定面積相当数未満となったクラスタに対する有効な処理を提供するもので、レーザスキャンに際して掃射されたレーザは複数の反射点から反射され、各々がリターン情報とともに受領される点が利用される。すなわち、対象物に到達したレーザは、所定のスポット径を有しているために、対象物に孔、隙間等がある場合には、最高点において先ず反射し（ファーストリターン）、この後、木漏れ日のように徐々に標高方向に進行し、セカンド、サードリターンとして反射する。このようなセカンドリターン以降のレーザデータは、性質上樹木が植生している山間部等において取得されることが多く、深度の深い、すなわち、リターン数の多いものほど地表面からの反射レーザデータである確率が高いこととなる。
【００３３】
一方、ファーストパルスが樹木の頂部で反射したとも限らないし、ラストパルスが地表面で反射したとも限らない上に、ラストパルスは反射強度が弱いために誤差成分が多く含まれやすい。したがって、これらセカンドリターン以降のレーザデータを扱うには、統計的解析が必要となる。
【００３４】
そこで、本方法において、先ず、当該クラスタ内でのセカンドパルス以降のレーザデータの抽出ステップ（Ｓ７１）を実行した後、抽出された各レーザデータと、近傍の地表反射レーザデータ及び他のセカンドパルス以降のレーザデータとの間で、上述した方法と同様に、傾斜角度算出ステップ（Ｓ７２）と割合による判定ステップ（Ｓ７３）を実施し、所定の傾斜角度未満の割合が所定の閾値以上である場合に地表面と判定する。
【００３５】
また、以上のようにして行われたフィルタリングの結果を確認するためには、
コンピュータが、
地表面反射成分として抽出されたメッシュ代表点１の属するメッシュ間を補間して形成した地表面モデルから等高線図３と、適宜の断面取得線４で切断した断面図とを作成し、
次いで、前記断面図上に前記断面取得線４が通過するメッシュを中心として断面幅方向メッシュ数を変化させながら対象メッシュ内の全レーザポイントデータを表示するという各処理を実行するフィルタリング結果の表示方法を使用することができる。
【００３６】
【発明の実施の形態】
図１、２にレーザデータからＤＴＭを作成する方法に適用された本発明の実施の形態を示す。レーザデータは、上空からのレーザ掃射により得られた地上物（地表、地物）の平面位置及び標高値を有するポイントデータであり、地表面判定に当たり、予め建物部分のデータをポリゴンデータ等として有している場合には、前処理として建物部分を除去し（ステップ（Ｓ１））、次いで、平面視において矩形のメッシュを設定する（ステップ（Ｓ２））。レーザデータ間の平均間隔がほぼ２．５（ｍ）であることを考慮して、メッシュサイズは５（ｍ）とする。
【００３７】
次いで、ステップＳ３において、各メッシュ内で標高値が最低のレーザデータを検索し、メッシュ代表点１として指定する。図３（ａ）にレーザデータ群に設定された各メッシュ内に１個のメッシュ代表点１を指定した状態を示す。メッシュ代表点１は、生データから抽出されたデータであり、抽出前の平面位置および標高値をそのまま保持する。
【００３８】
この後、ステップ（Ｓ４）において、クラスタリングを行う。クラスタリングは、高さ方向傾斜量が所定の閾値内であれば同一クラスタとする処理を行うことによってなされ、先ず、図３（ｂ）に示すように、注目するメッシュ代表点１と、このメッシュ代表点１が属するメッシュに対して縦方向及び横方向に隣接する４個のメッシュ内のメッシュ代表点１との間の高さ方向の傾斜角度を求め、所定の閾値と比較する。閾値には、５．７°が採用される。図３（ｃ）において縦横隣接するメッシュ内のメッシュ代表点１間の傾斜角度が閾値未満のもののみを線分で結んで示す。
【００３９】
また、クラスタリングに際して、各メッシュ代表点１からは、隣接するメッシュが空メッシュでない限り、隣接するメッシュ内のメッシュ代表点１に向けて４本の傾斜線が引かれるが、このうち、閾値未満の傾斜角度のものが所定の本数（この実施の形態においては１本）以上あるものが同一クラスタとして取り扱われる。
【００４０】
クラスタの境界は、直接隣接しているメッシュ間が閾値以上の傾斜であっても、回り込んで隣接する場合には同一のグループとして分類される。すなわち、図４（ａ）において、点ＡとＤとの間には閾値以上の傾斜があるが、点Ｂ、Ｃ、Ｄを結ぶ閾値未満の傾斜ルートが存在するために、Ａ〜Ｄが属するメッシュは同一クラスタに分類される。これに対し、点Ｅ、Ｆが属するメッシュは、Ａ〜Ｄが属するメッシュとは別のクラスタに分類される。なお、図４（ａ）においてメッシュ代表点１を黒点で、メッシュ代表点１間の傾斜が閾値未満の場合には、メッシュ代表点１間を実線で、閾値以上の場合には、破線で示される。
【００４１】
また、図４（ｂ）に、図３（ｃ）の結果から、クラスタリングした状態を、同一クラスタを連続する模様を施して示す。
【００４２】
この後、ステップ（Ｓ５）において、メッシュ数を基準とした地表判定を行う。上述したように、地表判定基準は、メッシュ数＝面積を示すことに着目し、最大建物想定面積以上のクラスタを地表面として、最小建物想定面積未満のクラスタを樹木等の建物以外の地物として抽出する。図に最大建物想定面積を５０００ｍ2＝２５（単位メッシュ面積）×２００（メッシュ）、最小建物想定面積を２００ｍ2＝２５×８（メッシュ）と設定した場合の地表面、建物及び地物を示す。
【００４３】
図４（ｂ）の結果から、メッシュ数２００以上で、地表と判定されたクラスタを図５（ａ）において右上がりのハッチングを施して示す。
【００４４】
以上のステップ（Ｓ１）からステップ（Ｓ５）まではコンピュータにより実行可能であり、これらのステップを記述したプログラムとして独立して使用することも可能である。
【００４５】
この実施の形態において、上記ステップ（Ｓ５）に引き続き、ステップ（Ｓ６）において傾斜量とメッシュ数による地表判定が行われる。このステップは、上記ステップ（Ｓ５）において最大建物想定面積と最小建物想定面積の中間にあるクラスタから地表面を抽出するもので、このステップにおける処理対象を図５（ａ）において右下がりのハッチングを施して示す。
【００４６】
傾斜量とメッシュ数による地表判定を行うに際して、先ず、ステップ（Ｓ６１）において上記条件を満足するクラスタを検索し、選択し、次いで、対象クラスタ内のメッシュ代表点１を相加平均してクラスタ代表点２を設定した（ステップ（Ｓ６２））後、ステップ（Ｓ６３）でクラスタ代表点２と、ステップ（Ｓ５）において地表面と判定されたメッシュ内のメッシュ代表点１との間の傾斜角度を算出する。図５（ｂ）に対象クラスタ（図右側）、及び傾斜角度算出の対象である地表面クラスタ内のメッシュ代表点１（地表面メッシュ代表点１；図左側）を示す。また、図５（ｃ）には、対象クラスタのクラスタ代表点２を白丸で、またクラスタ代表点２と傾斜角度算出の対象である地表面クラスタ内のメッシュ代表点１の一部とを線分で結んで示している。
【００４７】
次に、クラスタ代表点２と地表面メッシュ代表点１との傾斜角度を求める。図６は標高の断面を示しており、縦線はメッシュ境界を表している。図６（ａ）は、クラスタ内メッシュ数が１の場合を示しており、左から２本目と３本目のメッシュ境界の間にクラスタが存在する。頂点部分は算術平均により求めたクラスタ代表点２である。同様に、図６（ｂ）、（ｃ）、（ｄ）は各々クラスタ内メッシュ数が４，９，１６の例である。クラスタ内のメッシュ数が１個の場合にクラスタ代表点２と隣接メッシュの中心とが成す傾斜角度をαとすると（図６（ａ）参照）、クラスタ内のメッシュ数が４個（幅方向に２個）の場合には、クラスタ代表点２と隣接メッシュの中心とが成す傾斜角度βは、
β＝ tan- １ ( １× tan α／３ ) 、
となる。同様に、クラスタ内のメッシュ数が９個（幅方向に３個）の場合には、
β＝ tan- １ ( ２× tan α／４ ) 、
クラスタ内のメッシュ数が１６個（幅方向に４個）の場合には、
β＝ tan- １ ( ２× tan α／５ ) というように、順次低下する。
【００４８】
このようにαとβとの間には、一般に
メッシュサイズをＭ、要素数をＮとして、
β＝tan-１(２×tanα／(sqrt(N)+１))・・・（１）
の関係があるため、式（１）において適宜の管理角度（α）を代入して角度判定の閾値とする。なお、図６において縦線はメッシュ境界を示す。
【００４９】
このようにして傾斜角度の算出が終了すると、次にステップ（Ｓ６４）で地表面の判定が行われる。地表判定は、クラスタ代表点２とのなす傾斜量が上記閾角度以内である傾斜角度算出対象点の全体に対する割合が所定値以上か否かにより行われる。
【００５０】
図５（ｃ）は、ステップ（Ｓ６４）における判定方法を示すもので、図５（ｃ）においてクラスタ代表点２から上方に線分が引かれていることからも明らかなように、図の上方にも地表面クラスタが存在する。また、クラスタ代表点２から６本の線分（実線及び点線）が延びていることから、角度算出対象点が全部で６点あることが分かる。
【００５１】
上述したように、クラスタ代表点２は、平面位置及び標高値を有しており、この点と角度測定対象点との標高傾斜が、例えば角度判定の閾値α＝３０°未満である角度算出対象点を抽出する。図５（ｃ）において実線で示すように、上記条件を満たす角度算出対象点が４点であったとすると、全角度算出対象点数に占める割合は、
４／６＝０．６６６となり、
地表判定の割合閾値を１／３（＝０．３３３）とすると、
０．６６６＞０．３３３
より、このクラスタは地表面であると判定される。
【００５２】
なお、このステップ（Ｓ６）は、ステップ（Ｓ１）において建物部分が完全に除去されている場合には、省略することができる。また、以上のステップ（Ｓ１〜Ｓ４及びＳ６）を独立にコンピュータに実行させるようなプログラムを作ることも可能であり、さらに、ステップ（Ｓ１）からステップ（Ｓ６）までを連続してコンピュータに実行させるプログラムも作ることができる。
【００５３】
以上のようにして一応の地表面抽出を行った後、ステップ（Ｓ７）でリターン情報によるパルス解析を行う。
【００５４】
パルス解析は、セカンドパルス未満のパルスを利用して上記ステップ（Ｓ５）で建物以外の地物と判定されたクラスタから地表成分を分離するために行われる。パルス解析に際して、先ず、該当クラス内でセカンドリターン以降のレーザデータを抽出する（ステップ（Ｓ７１））。図７（ａ）においてセカンドリターン以降のレーザデータを黒丸で、ステップ（Ｓ６）までで地表面と判定されたメッシュ代表点１を白丸で示す。
【００５５】
次いで、ステップ（Ｓ７２）でセカンドリターン以降のレーザデータを中心とする傾斜角度を個別に算出する。傾斜角度の算出に先立って、図７（ｂ）において実線で示すように、セカンドリターン以降のレーザデータ間及び上記地表面反射レーザデータ間でＴＩＮネットワークを形成し、傾斜角度の算出は、ネットワーク対象点間を対象にして行われる。
【００５６】
地表面か否かの判定は、ステップ（Ｓ６）と同様の基準で行われ、処理方法がステップ（Ｓ６）とほぼ同一なので、このステップ（Ｓ７）はステップ（Ｓ６）と同時に行うことができる。
【００５７】
図７（ｂ）に地表面判定の条件を充たさなかったレーザデータを二重丸で、この結果、地表面として新たに指定された領域を図７（ｃ）においてハッチングを施して示す。
【００５８】
なお、以上のステップ（Ｓ１〜Ｓ４及びＳ７）を独立にコンピュータに実行させるようなプログラムを作ることも可能であり、さらに、ステップ（Ｓ１）からステップ（Ｓ７）までを連続してコンピュータに実行させるプログラム、あるいは、ステップ（Ｓ１〜Ｓ５及びＳ７）をコンピュータに実行させるプログラムも作ることができる。
【００５９】
以上のようにして地表面の抽出を行った後、ステップ（Ｓ８）で補間を行って地表面を求める。補間には、スプライン、クリギング等の種々の手法が使用できるが、トポグリッド手法を使用することによって、尾根・谷にめりはりをつけることができる。
【００６０】
さらに、この実施の形態は、以上のようにして求めた地表面の確認（ステップ（Ｓ１０））を含む。確認は、上述したフィルタリング結果を表示し（ステップ（Ｓ９））、この結果を目視することにより不自然な点を発見、修正して行われる。
【００６１】
確認時の手助けとなる結果表示に際し、先ず、ステップ（Ｓ８）により形成された地表面データから等高線図３を作成し（ステップ（９１））、さらに、ユーザの指定した断面取得線４に沿った断面図を作成する（ステップ（９２））。
【００６２】
この後、この断面図上に断面取得線４が通過するメッシュを中心として断面幅方向の所定数のメッシュ内の全レーザデータを表示する（ステップ（９３））。断面図上に表示するレーザデータ数、すなわち、断面取得線４に直交する方向のメッシュ数は適宜変更可能であり、当該メッシュ数を増減させた状態の変化から、地表面データの確からしさを確認できる。
【００６３】
図８はこの状態を示すもので、図８（ａ）に断面取得線４を一点鎖線で示した等高線図を、図８（ｂ）に断面取得線４に沿う断面図と、断面取得線４が通過する１列のメッシュ列内のレーザデータを表示した状態、図８（ｃ）に断面取得線４を中心に上下計３列のメッシュ列内の全レーザデータを表示した状態、図８（ｄ）に上下計５列のメッシュ列内のレーザデータを表示した状態を示す。
【００６４】
このレーザフィルタリング結果の表示方法は、上記ステップ（Ｓ５、Ｓ６、Ｓ７）のいずれかに続いて、ステップ（Ｓ８）を実行し、さらに、与えられた断面取得線４に沿って、与えられた断面幅で、ステップ（Ｓ９）をコンピュータに実行させるフィルタリング結果表示用プログラムを走行させることにより達成できる。
【００６５】
以上のようにしてフィルタリング結果の表示によりＤＴＭに問題がある場合には、閾値を変えてステップ（Ｓ３）からやり直すか、あるいはステップ（Ｓ１１）で手動修正が行われる。
【００６６】
これに対し、確認結果が良好である場合には、必要に応じて図化編集（ステップ（Ｓ１２））、および図面間接合（ステップ（１３））が実行されて作業が終了する。
【００６７】
【発明の効果】
以上の説明から明らかなように、本発明によれば、処理の効率化を図り、かつ、正確に地表面の抽出等のフィルタリングを行うことができる。
【図面の簡単な説明】
【図１】 本発明を示すフローチャートである。
【図２】 本発明を示すフローチャートである。
【図３】 クラスタリング方法を示す説明図である。
【図４】 クラスタリング結果を示す説明図である。
【図５】 傾斜量とメッシュ数による地表判定手順を示す説明図である。
【図６】 傾斜角度の算出を示す説明図である。
【図７】 リターン情報によるパルス解析の手順を示す説明図である。
【図８】 フィルタリングの結果表示を説明する図である。
【符号の説明】
１ メッシュ代表点
２ クラスタ代表点
３ 等高線図
４ 断面取得線
【発明の属する技術分野】
本発明は、レーザデータのフィルタリング方法及びプログラムに関するものである。
【０００２】
【従来の技術】
航空機レーザスキャナ等により取得したレーザデータから地表面データ、あるいはそれ以外の地物データを抽出することは、例えば、ディジタル標高モデル（ＤＴＭ：Digital Terrain Model）を作成する場合に必要であり、このような方法としては、例えば、特許文献１に記載されたものが知られている。
【０００３】
この従来例において、取得されたレーザデータ群はクラスタリングにより類似した特徴を持つグループ（クラスタ）に分類される。クラスタリングのために設定される特徴空間には、各レーザデータの平面位置、標高を要素とする３次元空間がそのまま利用される。
【０００４】
また、クラスタリングに際して定義される近傍にはボロノイ近傍が使用され、レーザデータ群が形成する特徴空間は、先ず、ボロノイ分割により各レーザデータの近傍を定義し、探索される。
【０００５】
【特許文献１】
特開２００２-３６６９７７号公報
【０００６】
【発明が解決しようとする課題】
しかし、上述した従来例において、クラスタリングに際する近傍点の探索は、全レーザデータに対して行う必要があるために、膨大な作業量を必要とするという欠点がある。
【０００７】
本発明は、以上の欠点を解消すべくなされたものであって、処理の効率化を図ることのできるレーザデータのフィルタリング方法及びプログラムを提供することを目的とする。
【０００８】
また、本発明の他の目的は、このようにしてフィルタリングした結果の表示方法の提供にある。
【０００９】
【課題を解決するための手段】
本発明によれば、上記目的は、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法を提供することにより達成される。
【００１０】
航空機レーザ測量で上空からのレーザ掃射により取得されるレーザデータは、平面位置（Ｘ，Ｙ）情報及び標高（Ｚ）情報を有するポイントデータであり、各レーザデータ単位では地表面等の識別を行うことはできない。このため、地表面等の識別には、隣接するレーザデータとの位置関係を利用することが必要であり、レーザデータ群に対する種々の処理が必要になる。
【００１１】
そして、レーザデータ群に対する処理は、上述した従来例をはじめとして、一般に各レーザデータを平等に扱い、これらに対する統計的処理として行われていた。
【００１２】
一方、レーザデータを種々の目的によって分類、抽出するレーザデータのフィルタリングに当たって、先ず、地表面データの抽出を優先することによって、地表面であることの特殊性をレーザデータ群の処理対象量の軽減に利用する可能性が生じる。すなわち、レーザデータ群に設定された好適なサイズのメッシュ内における最低標高点は、地表面である確率が高いために、これら最低標高点に処理対象を限定して処理対象量を軽減した方が、地表同士を検索の対象とするので、地表抽出精度は高くなる。
【００１３】
本発明は、以上の認識の下になされたものであって、メッシュ化及びメッシュ内でのメッシュ代表値の設定を行うことによって処理対象量を大幅に軽減させるとともに、クラスタリングに伴う近接点の検索を隣接メッシュに限定することによって処理効率を高め、かつ、正確な地表面抽出を可能にするものである。
【００１４】
すなわち、本発明において、先ず、レーザデータ群に対して平面座標に基づくメッシュを設定するステップ（Ｓ２）が実行される。メッシュサイズは、レーザデータが含まれない空メッシュが発生することのないように十分な大きさで、かつ、過度に周囲の条件に平均化されてメッシュの個性が失われない程度の小ささであることが必要であり、レーザデータの密度、正確にはポイントデータ間の間隔を考慮して、例えば、ポイントデータの間隔の２倍程度に決定される。
【００１５】
次いで、上記各メッシュ内のレーザデータのうち、最低標高を有するものをメッシュ代表点１として抽出するステップ（Ｓ３）が実行された後、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングするステップ（Ｓ４）が実行される。メッシュ代表点１を最低標高点データとすることにより、車、樹木等を予め除去することも可能になる。
【００１６】
クラスタリング操作は、標高要素を判断基準として隣接するメッシュとの間で行われる。近傍は、注目メッシュを包囲するように隣接する８個であっても、注目メッシュに対して縦横に隣接する４個であってもよいが、メッシュ代表点１が最低標高点により定義されているために、縦横の４個をとれば十分な精度を維持可能であり、かつ、処理量も少なくてすむ。
【００１７】
また、標高要素による判断は、例えば、隣接するメッシュ代表点１間の比高差をとることも可能であるが、注目メッシュのメッシュ代表点１との間の傾斜角度が所定値以下であるメッシュ代表点１の数を基準に行うことができる。
【００１８】
メッシュ代表点１は、メッシュ内のレーザデータから抽出されたものであり、Ｘ，Ｙ値及びＺ値をそのまま保持しているために、クラスタリングの基準として傾斜角度を利用すると、メッシュ全体を同一標高として考えることとなる比高差を利用する場合に比して、より現実の状態に即した判断が可能になる。
【００１９】
閾値として利用される傾斜角度は、経験則上、水平距離５（ｍ）に対する標高差が０．５（ｍ）、すなわち、５．７°程度に設定するのが望ましい。
【００２０】
以上のクラスタリングステップの後に、構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するステップ、すなわち、クラスタ内のメッシュ数による地表判定ステップ（Ｓ５）が実行される。
【００２１】
最大建物想定面積は、通常想定される建物の最大面積であり、通常、５０００ｍ2程度に設定しておけば足りる。最大建物想定面積に相当するメッシュ数、すなわち、最大建物想定面積をＳ、１メッシュの面積をｓ、クラスタ内メッシュ数をＭとした場合、
Ｍ＞Ｓ／ｓ
の条件を満たすクラスタは、全体に比較的なだらかな傾斜が建物と想定される場合の面積を超えて連続していることとなるために、非常に高い精度で地表部分であることを推定できる。
【００２２】
また、以上のフィルタリング方法は、なだらかで広い平野部の抽出に適しているが、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数未満、最小建物想定面積以上のクラスタ内のメッシュ代表点１の平面位置及び標高値の統計的代表値を平面位置及び標高値とするクラスタ代表点２を設定し、
次いで、前記クラスタ代表点２との間の傾斜角度が所定値未満である前記クラスタの近傍における地表面反射レーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法により、クラスタの面積により直ちに地表面と判断できない場合の地表抽出方法が提供される。
【００２３】
この方法において、クラスタリングまでのステップ（Ｓ２〜Ｓ４）は上述した方法と全く同様に行われ、この後、構成メッシュ数が最大建物想定面積相当数未満、最小建物想定面積以上のクラスタを検索するステップ（Ｓ６１）が実行され、次いで、該当クラスタ内にクラスタ代表点２を設定するステップ（Ｓ６２）が実行される。
【００２４】
クラスタ代表点２は、該当クラスタ内の複数のメッシュ代表点１を代表する平面位置及び標高値を有する点であり、メッシュ代表点１に対して適宜の統計的処理を施した統計的代表値として与えられる。統計的代表値を求めるためには、平均、あるいはメディアン値等を使用できる。
【００２５】
以上のようにして求められたクラスタ代表点２は、該当クラスを代表する平面位置及び標高情報を有しており、次いで、これらクラスタ代表点２とクラスタ近傍における地表面反射レーザデータとの間の傾斜角度を求めるステップ（Ｓ６３）が実行される。
【００２６】
傾斜角度を求めるに当たって、先ず、クラスタ代表点２から、適宜の手段によって地表面と判断されたレーザデータとの間で形成されるＴＩＮ（Triagulated Irregular Network）ネットワーク上の各点を近接点として定義することができる。また、この方法は、上述したメッシュ数による地表面抽出方法に続いて行うことが可能であり、この場合、傾斜角度は、メッシュ代表点１を対象として算出される。
【００２７】
該当クラスタの地表面判定は、以上のようにしてクラスタ代表点２と近傍の地表面反射レーザデータの間の傾斜角度を求めた後、傾斜角度が所定値未満の地表面反射データが全測定対象の地表面反射データに占める割合を求めるステップ（Ｓ６４）を実行することにより行われる。
【００２８】
この方法における目的は、建物が地表面からの判断が曖昧なクラスタからの地表面の判定であるために、判断の閾値となる傾斜角度は、地表面の想定される最大傾斜角度を基準に設定することができ、例えば、３０°に設定される。
【００２９】
また、割合の閾値は、判定する地形等を考慮して判定精度が高くなるように適宜決定可能であるが、１／３程度にするのが経験則上望ましく、閾値未満の傾斜角度の割合が上記閾値以上である場合には、当該クラスタは地表面と判定される。
【００３０】
このように、判定に当たってクラスタ内の各点を対象とすることなく、クラスタ代表点２を基準にクラスタ全体を判定することにより、計算量を少なくすることができる。
【００３１】
さらに、以上のフィルタリング方法に加え、
コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点１として抽出し、
次いで、隣接するメッシュのメッシュ代表点１を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最小建物想定面積相当数未満のクラスタ内で、かつ、セカンドパルス以降のレーザデータからなるメッシュ代表点１との間の傾斜角度が所定値未満である前記メッシュ代表点１の近傍における地表面反射レーザデータ及び他のセカンドパルス以降のレーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法を採用することにより、樹木等の多い山間部での地表面の抽出精度を高めることができる。
【００３２】
この方法は、上述した方法において、構成メッシュ数が最小建物想定面積相当数未満となったクラスタに対する有効な処理を提供するもので、レーザスキャンに際して掃射されたレーザは複数の反射点から反射され、各々がリターン情報とともに受領される点が利用される。すなわち、対象物に到達したレーザは、所定のスポット径を有しているために、対象物に孔、隙間等がある場合には、最高点において先ず反射し（ファーストリターン）、この後、木漏れ日のように徐々に標高方向に進行し、セカンド、サードリターンとして反射する。このようなセカンドリターン以降のレーザデータは、性質上樹木が植生している山間部等において取得されることが多く、深度の深い、すなわち、リターン数の多いものほど地表面からの反射レーザデータである確率が高いこととなる。
【００３３】
一方、ファーストパルスが樹木の頂部で反射したとも限らないし、ラストパルスが地表面で反射したとも限らない上に、ラストパルスは反射強度が弱いために誤差成分が多く含まれやすい。したがって、これらセカンドリターン以降のレーザデータを扱うには、統計的解析が必要となる。
【００３４】
そこで、本方法において、先ず、当該クラスタ内でのセカンドパルス以降のレーザデータの抽出ステップ（Ｓ７１）を実行した後、抽出された各レーザデータと、近傍の地表反射レーザデータ及び他のセカンドパルス以降のレーザデータとの間で、上述した方法と同様に、傾斜角度算出ステップ（Ｓ７２）と割合による判定ステップ（Ｓ７３）を実施し、所定の傾斜角度未満の割合が所定の閾値以上である場合に地表面と判定する。
【００３５】
また、以上のようにして行われたフィルタリングの結果を確認するためには、
コンピュータが、
地表面反射成分として抽出されたメッシュ代表点１の属するメッシュ間を補間して形成した地表面モデルから等高線図３と、適宜の断面取得線４で切断した断面図とを作成し、
次いで、前記断面図上に前記断面取得線４が通過するメッシュを中心として断面幅方向メッシュ数を変化させながら対象メッシュ内の全レーザポイントデータを表示するという各処理を実行するフィルタリング結果の表示方法を使用することができる。
【００３６】
【発明の実施の形態】
図１、２にレーザデータからＤＴＭを作成する方法に適用された本発明の実施の形態を示す。レーザデータは、上空からのレーザ掃射により得られた地上物（地表、地物）の平面位置及び標高値を有するポイントデータであり、地表面判定に当たり、予め建物部分のデータをポリゴンデータ等として有している場合には、前処理として建物部分を除去し（ステップ（Ｓ１））、次いで、平面視において矩形のメッシュを設定する（ステップ（Ｓ２））。レーザデータ間の平均間隔がほぼ２．５（ｍ）であることを考慮して、メッシュサイズは５（ｍ）とする。
【００３７】
次いで、ステップＳ３において、各メッシュ内で標高値が最低のレーザデータを検索し、メッシュ代表点１として指定する。図３（ａ）にレーザデータ群に設定された各メッシュ内に１個のメッシュ代表点１を指定した状態を示す。メッシュ代表点１は、生データから抽出されたデータであり、抽出前の平面位置および標高値をそのまま保持する。
【００３８】
この後、ステップ（Ｓ４）において、クラスタリングを行う。クラスタリングは、高さ方向傾斜量が所定の閾値内であれば同一クラスタとする処理を行うことによってなされ、先ず、図３（ｂ）に示すように、注目するメッシュ代表点１と、このメッシュ代表点１が属するメッシュに対して縦方向及び横方向に隣接する４個のメッシュ内のメッシュ代表点１との間の高さ方向の傾斜角度を求め、所定の閾値と比較する。閾値には、５．７°が採用される。図３（ｃ）において縦横隣接するメッシュ内のメッシュ代表点１間の傾斜角度が閾値未満のもののみを線分で結んで示す。
【００３９】
また、クラスタリングに際して、各メッシュ代表点１からは、隣接するメッシュが空メッシュでない限り、隣接するメッシュ内のメッシュ代表点１に向けて４本の傾斜線が引かれるが、このうち、閾値未満の傾斜角度のものが所定の本数（この実施の形態においては１本）以上あるものが同一クラスタとして取り扱われる。
【００４０】
クラスタの境界は、直接隣接しているメッシュ間が閾値以上の傾斜であっても、回り込んで隣接する場合には同一のグループとして分類される。すなわち、図４（ａ）において、点ＡとＤとの間には閾値以上の傾斜があるが、点Ｂ、Ｃ、Ｄを結ぶ閾値未満の傾斜ルートが存在するために、Ａ〜Ｄが属するメッシュは同一クラスタに分類される。これに対し、点Ｅ、Ｆが属するメッシュは、Ａ〜Ｄが属するメッシュとは別のクラスタに分類される。なお、図４（ａ）においてメッシュ代表点１を黒点で、メッシュ代表点１間の傾斜が閾値未満の場合には、メッシュ代表点１間を実線で、閾値以上の場合には、破線で示される。
【００４１】
また、図４（ｂ）に、図３（ｃ）の結果から、クラスタリングした状態を、同一クラスタを連続する模様を施して示す。
【００４２】
この後、ステップ（Ｓ５）において、メッシュ数を基準とした地表判定を行う。上述したように、地表判定基準は、メッシュ数＝面積を示すことに着目し、最大建物想定面積以上のクラスタを地表面として、最小建物想定面積未満のクラスタを樹木等の建物以外の地物として抽出する。図に最大建物想定面積を５０００ｍ2＝２５（単位メッシュ面積）×２００（メッシュ）、最小建物想定面積を２００ｍ2＝２５×８（メッシュ）と設定した場合の地表面、建物及び地物を示す。
【００４３】
図４（ｂ）の結果から、メッシュ数２００以上で、地表と判定されたクラスタを図５（ａ）において右上がりのハッチングを施して示す。
【００４４】
以上のステップ（Ｓ１）からステップ（Ｓ５）まではコンピュータにより実行可能であり、これらのステップを記述したプログラムとして独立して使用することも可能である。
【００４５】
この実施の形態において、上記ステップ（Ｓ５）に引き続き、ステップ（Ｓ６）において傾斜量とメッシュ数による地表判定が行われる。このステップは、上記ステップ（Ｓ５）において最大建物想定面積と最小建物想定面積の中間にあるクラスタから地表面を抽出するもので、このステップにおける処理対象を図５（ａ）において右下がりのハッチングを施して示す。
【００４６】
傾斜量とメッシュ数による地表判定を行うに際して、先ず、ステップ（Ｓ６１）において上記条件を満足するクラスタを検索し、選択し、次いで、対象クラスタ内のメッシュ代表点１を相加平均してクラスタ代表点２を設定した（ステップ（Ｓ６２））後、ステップ（Ｓ６３）でクラスタ代表点２と、ステップ（Ｓ５）において地表面と判定されたメッシュ内のメッシュ代表点１との間の傾斜角度を算出する。図５（ｂ）に対象クラスタ（図右側）、及び傾斜角度算出の対象である地表面クラスタ内のメッシュ代表点１（地表面メッシュ代表点１；図左側）を示す。また、図５（ｃ）には、対象クラスタのクラスタ代表点２を白丸で、またクラスタ代表点２と傾斜角度算出の対象である地表面クラスタ内のメッシュ代表点１の一部とを線分で結んで示している。
【００４７】
次に、クラスタ代表点２と地表面メッシュ代表点１との傾斜角度を求める。図６は標高の断面を示しており、縦線はメッシュ境界を表している。図６（ａ）は、クラスタ内メッシュ数が１の場合を示しており、左から２本目と３本目のメッシュ境界の間にクラスタが存在する。頂点部分は算術平均により求めたクラスタ代表点２である。同様に、図６（ｂ）、（ｃ）、（ｄ）は各々クラスタ内メッシュ数が４，９，１６の例である。クラスタ内のメッシュ数が１個の場合にクラスタ代表点２と隣接メッシュの中心とが成す傾斜角度をαとすると（図６（ａ）参照）、クラスタ内のメッシュ数が４個（幅方向に２個）の場合には、クラスタ代表点２と隣接メッシュの中心とが成す傾斜角度βは、
β＝ tan- １ ( １× tan α／３ ) 、
となる。同様に、クラスタ内のメッシュ数が９個（幅方向に３個）の場合には、
β＝ tan- １ ( ２× tan α／４ ) 、
クラスタ内のメッシュ数が１６個（幅方向に４個）の場合には、
β＝ tan- １ ( ２× tan α／５ ) というように、順次低下する。
【００４８】
このようにαとβとの間には、一般に
メッシュサイズをＭ、要素数をＮとして、
β＝tan-１(２×tanα／(sqrt(N)+１))・・・（１）
の関係があるため、式（１）において適宜の管理角度（α）を代入して角度判定の閾値とする。なお、図６において縦線はメッシュ境界を示す。
【００４９】
このようにして傾斜角度の算出が終了すると、次にステップ（Ｓ６４）で地表面の判定が行われる。地表判定は、クラスタ代表点２とのなす傾斜量が上記閾角度以内である傾斜角度算出対象点の全体に対する割合が所定値以上か否かにより行われる。
【００５０】
図５（ｃ）は、ステップ（Ｓ６４）における判定方法を示すもので、図５（ｃ）においてクラスタ代表点２から上方に線分が引かれていることからも明らかなように、図の上方にも地表面クラスタが存在する。また、クラスタ代表点２から６本の線分（実線及び点線）が延びていることから、角度算出対象点が全部で６点あることが分かる。
【００５１】
上述したように、クラスタ代表点２は、平面位置及び標高値を有しており、この点と角度測定対象点との標高傾斜が、例えば角度判定の閾値α＝３０°未満である角度算出対象点を抽出する。図５（ｃ）において実線で示すように、上記条件を満たす角度算出対象点が４点であったとすると、全角度算出対象点数に占める割合は、
４／６＝０．６６６となり、
地表判定の割合閾値を１／３（＝０．３３３）とすると、
０．６６６＞０．３３３
より、このクラスタは地表面であると判定される。
【００５２】
なお、このステップ（Ｓ６）は、ステップ（Ｓ１）において建物部分が完全に除去されている場合には、省略することができる。また、以上のステップ（Ｓ１〜Ｓ４及びＳ６）を独立にコンピュータに実行させるようなプログラムを作ることも可能であり、さらに、ステップ（Ｓ１）からステップ（Ｓ６）までを連続してコンピュータに実行させるプログラムも作ることができる。
【００５３】
以上のようにして一応の地表面抽出を行った後、ステップ（Ｓ７）でリターン情報によるパルス解析を行う。
【００５４】
パルス解析は、セカンドパルス未満のパルスを利用して上記ステップ（Ｓ５）で建物以外の地物と判定されたクラスタから地表成分を分離するために行われる。パルス解析に際して、先ず、該当クラス内でセカンドリターン以降のレーザデータを抽出する（ステップ（Ｓ７１））。図７（ａ）においてセカンドリターン以降のレーザデータを黒丸で、ステップ（Ｓ６）までで地表面と判定されたメッシュ代表点１を白丸で示す。
【００５５】
次いで、ステップ（Ｓ７２）でセカンドリターン以降のレーザデータを中心とする傾斜角度を個別に算出する。傾斜角度の算出に先立って、図７（ｂ）において実線で示すように、セカンドリターン以降のレーザデータ間及び上記地表面反射レーザデータ間でＴＩＮネットワークを形成し、傾斜角度の算出は、ネットワーク対象点間を対象にして行われる。
【００５６】
地表面か否かの判定は、ステップ（Ｓ６）と同様の基準で行われ、処理方法がステップ（Ｓ６）とほぼ同一なので、このステップ（Ｓ７）はステップ（Ｓ６）と同時に行うことができる。
【００５７】
図７（ｂ）に地表面判定の条件を充たさなかったレーザデータを二重丸で、この結果、地表面として新たに指定された領域を図７（ｃ）においてハッチングを施して示す。
【００５８】
なお、以上のステップ（Ｓ１〜Ｓ４及びＳ７）を独立にコンピュータに実行させるようなプログラムを作ることも可能であり、さらに、ステップ（Ｓ１）からステップ（Ｓ７）までを連続してコンピュータに実行させるプログラム、あるいは、ステップ（Ｓ１〜Ｓ５及びＳ７）をコンピュータに実行させるプログラムも作ることができる。
【００５９】
以上のようにして地表面の抽出を行った後、ステップ（Ｓ８）で補間を行って地表面を求める。補間には、スプライン、クリギング等の種々の手法が使用できるが、トポグリッド手法を使用することによって、尾根・谷にめりはりをつけることができる。
【００６０】
さらに、この実施の形態は、以上のようにして求めた地表面の確認（ステップ（Ｓ１０））を含む。確認は、上述したフィルタリング結果を表示し（ステップ（Ｓ９））、この結果を目視することにより不自然な点を発見、修正して行われる。
【００６１】
確認時の手助けとなる結果表示に際し、先ず、ステップ（Ｓ８）により形成された地表面データから等高線図３を作成し（ステップ（９１））、さらに、ユーザの指定した断面取得線４に沿った断面図を作成する（ステップ（９２））。
【００６２】
この後、この断面図上に断面取得線４が通過するメッシュを中心として断面幅方向の所定数のメッシュ内の全レーザデータを表示する（ステップ（９３））。断面図上に表示するレーザデータ数、すなわち、断面取得線４に直交する方向のメッシュ数は適宜変更可能であり、当該メッシュ数を増減させた状態の変化から、地表面データの確からしさを確認できる。
【００６３】
図８はこの状態を示すもので、図８（ａ）に断面取得線４を一点鎖線で示した等高線図を、図８（ｂ）に断面取得線４に沿う断面図と、断面取得線４が通過する１列のメッシュ列内のレーザデータを表示した状態、図８（ｃ）に断面取得線４を中心に上下計３列のメッシュ列内の全レーザデータを表示した状態、図８（ｄ）に上下計５列のメッシュ列内のレーザデータを表示した状態を示す。
【００６４】
このレーザフィルタリング結果の表示方法は、上記ステップ（Ｓ５、Ｓ６、Ｓ７）のいずれかに続いて、ステップ（Ｓ８）を実行し、さらに、与えられた断面取得線４に沿って、与えられた断面幅で、ステップ（Ｓ９）をコンピュータに実行させるフィルタリング結果表示用プログラムを走行させることにより達成できる。
【００６５】
以上のようにしてフィルタリング結果の表示によりＤＴＭに問題がある場合には、閾値を変えてステップ（Ｓ３）からやり直すか、あるいはステップ（Ｓ１１）で手動修正が行われる。
【００６６】
これに対し、確認結果が良好である場合には、必要に応じて図化編集（ステップ（Ｓ１２））、および図面間接合（ステップ（１３））が実行されて作業が終了する。
【００６７】
【発明の効果】
以上の説明から明らかなように、本発明によれば、処理の効率化を図り、かつ、正確に地表面の抽出等のフィルタリングを行うことができる。
【図面の簡単な説明】
【図１】 本発明を示すフローチャートである。
【図２】 本発明を示すフローチャートである。
【図３】 クラスタリング方法を示す説明図である。
【図４】 クラスタリング結果を示す説明図である。
【図５】 傾斜量とメッシュ数による地表判定手順を示す説明図である。
【図６】 傾斜角度の算出を示す説明図である。
【図７】 リターン情報によるパルス解析の手順を示す説明図である。
【図８】 フィルタリングの結果表示を説明する図である。
【符号の説明】
１ メッシュ代表点
２ クラスタ代表点
３ 等高線図
４ 断面取得線
Claims (7)
- コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点として抽出し、
次いで、隣接するメッシュのメッシュ代表点を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法。 - コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点として抽出し、
次いで、隣接するメッシュのメッシュ代表点を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最大建物想定面積相当数未満、最小建物想定面積以上のクラスタ内のメッシュ代表点の平面位置及び標高値の統計的代表値を平面位置及び標高値とするクラスタ代表点を設定し、
次いで、前記クラスタ代表点との間の傾斜角度が所定値未満である前記クラスタの近傍における地表面反射レーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法。 - コンピュータが、
上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定した後、
各メッシュ内の最低標高点データをメッシュ代表点として抽出し、
次いで、隣接するメッシュのメッシュ代表点を近傍として標高要素に基づいてクラスタリングし、
この後、構成メッシュ数が最小建物想定面積相当数未満のクラスタ内で、かつ、セカンドパルス以降のレーザデータからなるメッシュ代表点との間の傾斜角度が所定値未満である前記メッシュ代表点の近傍における地表面反射レーザデータ及び他のセカンドパルス以降のレーザデータの割合を求め、
前記割合が所定値以上である場合に前記対象クラスタ内のレーザデータを地表面反射成分として抽出するという各処理を実行するレーザデータのフィルタリング方法。 - 前記傾斜角度算定に際して使用される地表面反射成分は、請求項１記載のレーザデータのフィルタリング方法により地表面反射成分として抽出されたクラスタ内のメッシュ代表点である請求項２または３記載のレーザデータのフィルタリング方法。
- 前記クラスタリングが、注目メッシュのメッシュ代表点との間の傾斜角度が所定値未満であるメッシュ代表点の数を基準に行われる請求項１ないし４のいずれかに記載のレーザデータのフィルタリング方法。
- コンピュータが、
請求項１ないし５のいずれかに記載のレーザデータのフィルタリング方法により地表面反射成分として抽出されたメッシュ代表点の属するメッシュ間を補間して形成した地表面モデルから等高線図と、適宜の断面取得線で切断した断面図とを作成し、
次いで、前記断面図上に前記断面取得線が通過するメッシュを中心として断面幅方向メッシュ数を変化させながら対象メッシュ内の全レーザデータを表示するという各処理を実行するフィルタリング結果の表示方法。 - 上空を飛行する飛行体から地上をレーザ掃射して取得したレーザデータ群に平面座標に基づくメッシュを設定するステップと、
各メッシュ内の最低標高点データをメッシュ代表点として抽出するステップと、
隣接するメッシュのメッシュ代表点を近傍として標高要素に基づいてクラスタリングするステップと、
構成メッシュ数が最大建物想定面積相当数以上のクラスタ内に属するレーザデータを地表面反射成分として抽出するステップとをコンピュータに実行させるためのレーザデータのフィルタリングプログラム。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003063841A JP4226360B2 (ja) | 2003-03-10 | 2003-03-10 | レーザデータのフィルタリング方法及びプログラム |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2003063841A JP4226360B2 (ja) | 2003-03-10 | 2003-03-10 | レーザデータのフィルタリング方法及びプログラム |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2004272688A JP2004272688A (ja) | 2004-09-30 |
| JP2004272688A5 JP2004272688A5 (ja) | 2006-05-11 |
| JP4226360B2 true JP4226360B2 (ja) | 2009-02-18 |
Family
ID=33125322
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2003063841A Expired - Fee Related JP4226360B2 (ja) | 2003-03-10 | 2003-03-10 | レーザデータのフィルタリング方法及びプログラム |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4226360B2 (ja) |
Families Citing this family (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5071679B2 (ja) * | 2008-06-17 | 2012-11-14 | 朝日航洋株式会社 | 傾斜分析装置、方法及びプログラム |
| JP5356269B2 (ja) * | 2010-01-29 | 2013-12-04 | 株式会社パスコ | レーザデータのフィルタリング方法及び装置 |
| JP5750309B2 (ja) * | 2011-05-27 | 2015-07-22 | 株式会社東京建設コンサルタント | レーザプロファイラデータから浸水予測地点を抽出する方法 |
| JP6097182B2 (ja) * | 2013-08-28 | 2017-03-15 | 株式会社パスコ | データ解析装置、データ解析方法、及びプログラム |
| US10127709B2 (en) | 2014-11-28 | 2018-11-13 | Panasonic Intellectual Property Management Co., Ltd. | Modeling device, three-dimensional model generating device, modeling method, and program |
| WO2016139819A1 (ja) * | 2015-03-03 | 2016-09-09 | 株式会社岩崎 | 地盤形状推定プログラム、地盤形状推定装置および地盤形状推定方法 |
| JP7389883B1 (ja) * | 2022-11-21 | 2023-11-30 | 株式会社パスコ | 地表点抽出装置、地表点抽出方法およびプログラム |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP4619504B2 (ja) * | 2000-09-19 | 2011-01-26 | アジア航測株式会社 | ３次元デジタル地図作成装置 |
| JP2002236019A (ja) * | 2001-02-08 | 2002-08-23 | Pasuko:Kk | 地表面抽出処理システム |
| JP2002366977A (ja) * | 2001-06-08 | 2002-12-20 | Pasuko:Kk | 地物形状作成システム及び３次元地図作成システム |
-
2003
- 2003-03-10 JP JP2003063841A patent/JP4226360B2/ja not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JP2004272688A (ja) | 2004-09-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106780524B (zh) | 一种三维点云道路边界自动提取方法 | |
| Zolanvari et al. | Three-dimensional building façade segmentation and opening area detection from point clouds | |
| CN108427124B (zh) | 一种多线激光雷达地面点分离方法及装置、车辆 | |
| JP5161936B2 (ja) | データ解析装置、データ解析方法、及びプログラム | |
| CN110363861B (zh) | 基于激光雷达点云的田地作物三维重构方法 | |
| CN110969624A (zh) | 一种激光雷达三维点云分割方法 | |
| JP5762251B2 (ja) | 建物輪郭抽出装置、建物輪郭抽出方法及び建物輪郭抽出プログラム | |
| CN111667574B (zh) | 从倾斜摄影模型自动重建建筑物规则立面三维模型的方法 | |
| Hofmann | Analysis of TIN-structure parameter spaces in airborne laser scanner data for 3-D building model generation | |
| JP4226360B2 (ja) | レーザデータのフィルタリング方法及びプログラム | |
| CN111325138A (zh) | 一种基于点云局部凹凸特征的道路边界实时检测方法 | |
| CN104751479A (zh) | 基于tin数据的建筑物提取方法和装置 | |
| CN111598780A (zh) | 一种适用于机载LiDAR点云的地形自适应插值滤波方法 | |
| JP5197205B2 (ja) | 輪郭抽出装置と輪郭抽出方法及び輪郭抽出プログラム | |
| CN114612632A (zh) | 一种基于三维激光点云数据的排序、插值处理方法 | |
| CN109410183A (zh) | 基于点云数据的平面提取方法、系统、装置及存储介质 | |
| US20220114810A1 (en) | Forest resource information generation structure and forest resource information generation method | |
| CN114898118A (zh) | 基于多源点云的输电线路房屋拆迁量自动统计方法及系统 | |
| Kada et al. | Sub-surface growing and boundary generalization for 3D building reconstruction | |
| CN110216664A (zh) | 基于点云数据的障碍物识别方法及装置 | |
| KR100726333B1 (ko) | 라이다 데이터로부터 셰도-그리드를 이용한 건물 외곽선자동추출방법 | |
| CN110415210B (zh) | 一种基于点云贪婪三角投影构建模型的孔洞检测和修补方法 | |
| JP2003296706A (ja) | 建物抽出方法及び装置並びにサンプルデータのベクトル化方法及び装置 | |
| CN112967384A (zh) | 识别建筑测绘构件的点云智能化切分方法 | |
| CN115937226A (zh) | 一种基于无人机Lidar点云数据的果树单木分割方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060309 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20060309 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20080826 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20081027 Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20081027 |
|
| RD02 | Notification of acceptance of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7422 Effective date: 20081027 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20081125 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20081126 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111205 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4226360 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111205 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121205 Year of fee payment: 4 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131205 Year of fee payment: 5 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313531 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| LAPS | Cancellation because of no payment of annual fees |