TREND-POINT

地上レーザーのデータを可視化しています。

サーフェス化して断面を切っていますが、同時に航空LPのデータも切っています。大体、同じ形状になります。両者の精度、思ったより良いですね。

レーザーのデータから地表を抽出するには、木や草を除去する必要があります。また、座標の変換も必要な場合があるでしょう。これ、ReCap ではできません。

プロに聞いたら、ライセンスに空きのあるソフトがありました。
TREND-POINT
https://const.fukuicompu.co.jp/products/trendpoint/

Leica のデータも読めました。
が、樹木等を除去するのはかなり困難。プロもトライ＆エラーで解決するそうです。航測会社さんの処理はとても速いのですが、きっと独自のノウハウがあるのでしょう。
これはあきらめました。

座標変換は簡単。
図面から公共座標を拾って点群に与えるだけ。それを計算させると独立座標から変換してくれます。この機能だけでも使えそうです。

プロに聞いたところ、点群を平面図に変換できるソフトはないようです。が、いずれできるでしょうね。フルサーフェス化したり、ソリッド化したり。

新たなハードが出ると新たなデータが生まれ、それを効率よく処理するソフトが整備されます。昔からそうですが、近年は浸透が速く、多様化しているように思えます。
頑張って、ついて行きましょう。

プロの仕事？ その2

地上レーザースキャナーを初めて使いました。

機材所有部署から事前にレクチャーを受け準備していたのですが、事情により現場の前日に別機種に変更。再度1時間ほどレクチャーを受けました。
テスト中、手元では70m程の点群しか表示されていなかったのですが、射程は倍だそうで、映っていない点群も解析時には使用できるとのことでした。

で、現場。
今まで3次元地質モデルを作成する際には、地表面を作るのがネックになっていました。
測量図面から3次元に起こすのは手間。UAV ＋ SfM が楽なのですが、自分では積極的に操作したくはない。結果的に他部署にお願いする形になっていました。
ですが、このスキャナーであれば墜落の心配はなく、数分でデータを取得してくれます。何度か場所を変える必要は生じますが、GPSも付いており、簡単・楽です。

帰社後にデータを確認。
が、70m以遠のデータが出てきません。スキャナーを所有していた部署にも確認したのですが、「そのようなことはないはず」との回答。弱りました。

結論から言うと、データは取れていませんでした。いくつかのモードの中で、点群を最も密に取得できるモードを選択していたのですが、この場合は65mまでと’マニュアル’に記載がありました。そう、マニュアルに。
読まずに出かけました。反省。

で、現場はやり直し。
これがプロの仕事？と問われると恥ずかしくなるような内容です。ベテランですし、変化に追従する方法をよく知っていましたが、実施しておらず失敗。お粗末です。

次は失敗しません。十分に準備していきましょう。

点群から変位の抽出

LP の整備もあり、点群データは多くの方に認知されたようと感じています。

最近では、LP を Web GIS で扱いたい要望や、地上レーザーを購入したい要望などが聞こえてきます。CIM 推しにより、点群が "become  a thing" 状態になったということでしょう。

そうなると次のステップも普及しそうです。例えば、点群の差分表示。差分を取れば地表面の変化を表現できますし、体積変化を出せます。崩壊前後の時系列データで、時点の崩壊・堆積量を求めることが可能になります（崩壊面、すべり面は、これだけではダメですが）。
よく利用されている表現ですが、広域かつ維持管理にまで普及すれば容易に変化箇所を把握できるようになります。広域だとメッシュ処理後の対象より粗いデータになりやすいのでスケールの問題はありますが、PSInSAR や MMS 等の併用でカバーできるでしょう。

次に変位。
何らかのマッチング手法を用いて2時期の変位を表現します。6年前に optical flow を試した時期がありましたが、これはイマイチでした。
https://phreeqc.blogspot.com/2012/02/optical-flow-block-matching.html
近年では ICP による点群同士、あるいは点群とメッシュのマッチングを行い、その結果から変位を抽出する手法が主流でしょうか？これに関してはプラントや岩盤を扱う例の載っている商用利用可能なソフト（GPL）も出ています。
CloudCompare Ver.2.6.1 http://www.danielgm.net/cc/
関連する箇所はコチラ↓
Registration > Align (point pairs picking)
　photogrammetry clouds で有効
Registration > Fine registration (ICP)
Distances > Cloud/Cloud dist. (cloud-to-cloud distance)
Distances > Cloud/Mesh dist. (cloud-to-mesh distance)

Optical Flow では 2D に修正する手間を必要としましたが、コチラは 3D の点群で、そのまま処理・表現できるところが長所でしょう。
ただ、過去にある会社が特許申請していましたので、土木分野では気を付けないといけません（結果はどうなったのでしょう）。
こういった海外のソフトの説明書に既に載っている例を国内で特許申請かけようとする姿勢は、消化し難いですね。誰でも使えるソフトや技術が、ある分野では国内で使えなくなるというのもイマイチです。

いずれにしても環境は既に整っていますので、すぐに次のステップに進むでしょう。しばらく、動向を見守りましょう。

1mメッシュ（ LP） の精度

LP の 1m グリッドデータを使用してサーフェスを作成し、縦横断を切っていました。

その後、なぜか実測してもらう機会があったのですが、結果を見て LP の高さの精度に驚かされました。誤差10cm未満でした。
木や草を除去する技術が向上したのでしょうか？5m でなく 1m だからでしょうか？GPS の精度が向上したのでしょうか？
いずれにしても、測量屋さんの技術と努力に脱帽です。

よく利用している国土地理院の 5DEM では、精度が50cm以内と感じています。こうなると、1m DEM を整備してほしいですね。
https://phreeqc.blogspot.jp/2012/11/5m_4.html

せっかく高品質の成果を提供いただいたので、最低限、それを損なわないようにしなくてはなりません。
LPだけでなく、UAVや3次元設計、その他のノベーションは続くと思います。それらの変化に遅れを取らず対応できるよう、こちらも努力して備えなくてはなりません。

UAV で土量計算

UAVで土量計算を実施している例がありました。

大成建設がダム工事の現場で効果を実証
http://www.autodesk.co.jp/adsk/servlet/item?siteID=1169823&id=24316438
http://news.mynavi.jp/news/2015/04/24/101/

大成さんのスライドを見る限り、DJI F550 + GoProHero3+ + ReCap360＋Civil3Dですので、手持ちのツールで十分再現できます。が、Autodesk の発表を見る限り、点群作成はReCap360 とは書かれていません。ReCap360にイリノイ大学の新しいアルゴリズムが取り込まれたと解釈すべきかどうかわかりませんが、実務上は大成さんのスライド通り、ReCapを使用しているのでしょうね。
詳細は掲載されていませんが、この技術の素晴らしい点は、イリノイ大学のアルゴリズムでしょう。
SfMはLPと異なり相対位置を使用しますので、画素の位置を見定めるアルゴリズムと写真を撮る高度（解像度）さえ工夫すれば、今回のような精度を出せるのでしょう。また、SfM の場合、UAV の位置情報はLPに比べ軽視できる=安価な GPS で OK となりますので

UAV を地形測量に利用する場合、3つのノウハウが必要かと思われます。①撮影のノウハウ、②モデルを作るノウハウ、③ヘリを操縦するノウハウです。

①撮影はモデルを作る上で、どのカメラセッティングで、何処の写真が、どの程度の頻度（精度、カバー率）でほしいか、こだわりが出てくると思います。そうすると、どのスピードで、どの高度で飛ばしたいかなどといった要望が出てきます。何度かモデルを作れば、自然と身につくノウハウです。

②モデルを作るノウハウは、ほとんどありません。点群までは自動でできます。ググれば得られる情報です。スマホと 123D で遊んでいるのと同レベルです。ハードルはほとんどありません。
成果を売るのであれば、必要なのは「お金」です。金をかければ良いものが早くできるのが現状です。PhotoScanのようにGPU対応のモノもありますので。
なお、点群＞サーフェス＞コンター＞横断の流れは機械作業ですので、ＣＡＤの使えるパートさんでも可能です。

③飛ばすノウハウは・・・素人ですのでわかりません。基本はプロに任せようと考えています。
http://phreeqc.blogspot.jp/2015/04/uav.html

結局は③がネックになり、①②を分かっているプロに頼むのが安心ということになります。もしくは、自分で③を身に付けるかですね。何度も墜落させて。


作ったコンターを地形測量としてお金をもらうとなると公共測量作業規定が絡んできますので、測量屋さんににお願いするようになります。
地質屋が単独で可能なことは、見せるモデルを作るくらいでしょう。4年前くらいでしょうか、当初は3次元解析用の表層データの取得を一つの目的として始めましたが、現段階では見せるモデルまでに止めています。
ただ、崩壊地形を色々な角度で見ていると、見えていなかった構造が見えてくる場合がありますので、UAV+SfM に限らず、3D可視化は一連の作業として実施すれば良いと思います。

とりあえず手を動かすことが必要でしょう、

崩壊前の地形

時間があったので、崩壊後、土砂がどの程度移動するかを再現していました（今度は LSFLOW です）。

最初は簡単に考えていました。

崩壊前後の LP を使用すれば、モデルは簡単にできてしまう、と考えていたのです。が、これが誤り。多くの場合、崩壊後の LP はありますが、崩壊前はありません。写真から復元したとしても、その精度はLPに比べて落ちます。場合によっては10m以上の標高差が生じます。今回も部分的に5～10m、最大20m近いギャップが生じました。

http://phreeqc.blogspot.jp/2012/11/5m_4.html

http://phreeqc.blogspot.jp/2012/11/5m-2.html

崩壊後のLPにしても、そのままではすべり面として使えません。冠頭部はそのままで使えるのですが、末端部は土砂で隠れています。そのため、ボーリング結果を使って面を作ってやらないといけません。こんな当たり前のことも、手を動かすまで気づきませんでした。迂闊でしたね。

こういった前処理に手間をかけることが、精度を大きく向上させるための一因になります。

案外、パラスタより時間がかかるかもしれません。

LPを利用した標高データ（5mメッシュ）の公開

一昨日より、基盤地図情報サイトにて、LPから作成した標高データ（5mメッシュ）が追加されています。
http://www.mlit.go.jp/report/press/mizukokudo03_hh_000481.html

早速、DL！と思っていたのですが、アクセスが集中しているようでサーバーに繋がりません。昨夜、やっとつながり確認できました。標高(DEM)が大幅に増えています。選択もしやすくなっていますね。数か月前よりメインサーバーのメンテナンス中でしたが、大幅なデータ追加のためだったのかもしれません。

図化してみましたが、構造物の消し方や図郭のつなぎ目に多少の違和感が残るものの、山中は良いですね。以前、落石調査に行った山を見ると、崖の位置や小さな谷筋が分かりやすくなっていました。シミュレーションはもちろん、踏査の基図や地質断面作成用のサーフェスとして利用価値が高いと思います。

国交省管轄の河川沿いのLPデータは、すべて撮り終わっていると聞いていましたが、公開されるとは思っていませんでした。昨年の地震の影響なのでしょうか？平成24年度は津波シミュを控えていますので、広範囲の詳細データが必要になったのでしょう。

今後は1mメッシュの公開を期待しましょう。

LPデータによる地表面変位表示

Optical Flow による地表面変位のベクトル表示ですが、LPデータを先輩からお借りし、試してみました。

いえ、正確には試そうとしましたが、無理でした。
移動前が10mメッシュ、移動後が1mのグリッドデータですと、見た目が全く異なります。そうなると、Optical Flow では対応できません。10mの解像度に合わせても無意味です（移動量にもよりますが）。

手を動かしてようやく気付いたのですが（遅い！）、この手法は実務での適用性が狭そうです。地すべりに対し2回LPを撮るような現場は少ないでしょうし、崩壊には適用できません。

ただ、画像として表現できれば、写真であれ数値であれ、材料を選ばず比較できます。当然、動画も問題ありません。泥流の流速を自動取得したいとか、動画に映っている落石の個数を取得するとかといった利用の方が向いているかもしれません。
技術的にはレベルが低いですから、何ができるかはアイデア次第ですね。

OpenCV とベクトル表示

地盤工学会誌にレーザー計測に関する講座が連載されています。

今月号では2回の測定データの差分による、土砂移動量や移動ベクトルの表示が掲載されていました。後者の「5.3.2地表面変位ベクトル解析（ｐ55）」については特許技術のようですね（特許技術を講座に載せられても勉強にならないのですが）。

いまや XBOX の kinect が USB で PC に接続できる時代です。画像追跡などはそれほど困難な技術のように思えません（kinect は驚愕！）。
CGのプロにその話をすると、「OpenCV で可能では？」と教えていただきました。調べてみると、映像処理や解析に特化したライブラリー集のようです。CVは computer vision でしょうね。

確かに、モーション解析と物体追跡は実装されています。
私が見たサンプルコードでは、2枚の画像の特徴点を検出し、それらの間の移動ベクトルを表示させるアルゴでした。地すべりに適用するには、LPデータの画像処理にひと手間必要かもしれません。

ま、ダメもとで、試してみましょう。

開度-ウェーブレット解析図

地上開度図とウェーブレット解析図の合成です。

と、ここまで来て、かなり落胆しました。

合成の手法がマニュアルに書かれていません。合成の演算は山ほどあります。計算手法を示すことができないのであれば、せめて何のソフトの、どの演算を使ったかを書いてないと、追いかけられません。ましてや開度図の階調を見やすいように修正してから合成だと、純粋な計算結果を演算するのではなく、2枚の画像（適度に補正後の数値）を適度に演算する世界でしょう（だからマニュに書かれていないのでしょうか？）。画像の合成法は見やすいようにヨキニハカラエということなのかもしれません。

仕方ないので、Photoshopで合成することにしました。

まず、使った素材です。

修正地上開度図

Wavelet 解析図

SAGA の場合、Protection Index は水平からの角度になります（地上開度は鉛直からの角度）。グレースケールを反転して、修正地上開度図としました。

以下、合成図です。

合成図1 : Over Lay

合成図2 : Hard Light

合成図3 : Linear Light

合成図4  : Vivid Light

合成図5  : Multiply
（コントラスト、明度補正済み）

Linear Light, Vivid Light, Multiply が見やすいでしょうか？


これらの合成手法が、正しいかどうか分かりません。しかし、確かめることもできません。判読しやすい画像を作成するのが目的で、それに適した合成手法を選択するのが正解であれば、前処理ルーチン（材料作成）の必然性のなさが落胆の最大の理由です。最初から画像に対してエッジ抽出などのフィルターをかけて合成し、判読しやすいように作っても同じことでしょう。まあ、マニュに示してある手法でも対応できるといった「営業力」がついたということで、良しとしましょう。


ところでマニュから離れますが、今回知った SAGA のいろいろな演算が何を意味するのかは、興味が残るところです。何の分野でどのような目的に使用されているのか分かりませんが、教科書や文献はあるはずです。こちらは今後の「技術力」として調べてみましょう。

Wavelet その2

メキシカンハット関数を使った、ウェーブレット解析（ただの演算ですが）の続きです。

昨日UPした通り、(1-x^2・・・から始めることにしました。
2次元関数はこのような感じです。

マニュアルではスケール S を1で固定し、DEM のグリッドサイズで波長が変わるといった内容です。λ=4グリッドといったところですね。実務的ですね。
x-aが影響するところまでを考えると、最低でも中心から左右に4グリッドずつ必要です。従って、フィルタは9*9のマトリックスになります。

早速、SAGA で User Defined Filter があるか調べてみました。結果、装備されていました。しかし、3*3までなので適用できません。仕方ないですね。
Surfer では9*9でも可能ですから、こちらでやってみましょう。ただし、Filter による演算結果は、最後にマトリックスの合計値で割る設定になっています。そのため、マニュアルの積分に合わせるには、後で Grid Math で合計値を掛け戻してやる必要があります。また、Filter の edge の処理設定もしないといけません。

演算結果を確認してみますと、あれ？といった感じの仕上がりです。全くマニュの図とは異なっています。係数のオーダーも大きく違います。チェックしましたが計算過程は問題ありません。
ということは、(1-x^2・・・ではなくて、（2-x^2が正しいということでしょうか？

(2-x^2・・・の2次元関数は以下のようになります。振幅の最大値は大きくなりますが、同じような形状です。

Surfer で演算し直し、できた GRD ファイルを SAGA で読み込んだ結果がこちら。

オーダーもあっています。（2-X^2・・・で問題ないようでした。
最終的に、地すべり地形を判読しやすいと実感したら、式の導出まで追いかけてみましょうか。このまま使うのは気持ち悪いですからね。


とりあえず、地上開度、ウェーブレット解析、両方ともクリアーです。
地上開度もよく読んでみますと、1グリッドで計算した場合も書かれていますので、最大値を取るアルゴを組まなくてもよくなっています。非常に実務的と言いますか、汎用ソフトの利用を前提に、適用しやすいよう単純化して書かれています。技術的に良いか悪いかは分かりませんが、裾野を広げることはできるでしょう。


あとは地上開度図とウェーブレット解析図の合成のみです。
こちらはまた後日。

Wavelet

地上開度がクリアできたので、俄然やる気が出てきました。

土研「地すべり地における航空レーザー測量解析マニュアル（案）」に従えば、あと視覚化で残るのはウェーブレット解析です。

ウェーブレット関数とは何ぞや？の状態ですので、ネットで調べてもよく分かりませんでした。これだけで１冊の本になるほど多肢にわたっています。

しかし、マニュのメキシカンハット関数やウェーブレット係数の算出式は非常に簡素なので、これに限ってはできるような気がしました。

で、とりあえず、マニュの式をEXCELでグラフ化。

あれ？最大値が11ページのグラフと違います。おかしい。

ネットで調べてみると、ψ(x, y) = （2-x^2・・・ の箇所が（1-x^2・・・の形になっている式が多く見られます。y = 0のときに（1-x^2と同じ形にならないといけないはずですから誤植でしょうか？

とりあえず2を1に修正してみると、

11ページと同じ図ができました。私の知識ではどちらが正解かわかりませんが、係数を出すときの重みの大きさが変わるだけなので、とりあえずこちらで行きましょう。今度、地形解析チームに出典の文献を持っていないか聞いてみましょう。

あとは積分の式をどう反映させるかです。これはソフト側でユーザー定義フィルターが使用できるかどうかも重要でしょう。SAGA で可能なのでしょうか？Surfer では可能ですし、ImageJでもMexican Hat の計算はできますから、最悪はそちらでの検討になるでしょう。

続きは後日。

SAGA

SAGAの続きです。

LPデータを読み込みました。
ナンバリングしてあれば、SurferでXYZのみに加工したデータをCSVとして保存しておけば読み込めます(EXCEL2010では行数に限度があります)。また、Surferのグリッドデータ読めます。

Terrain Analysis Module がたくさん装備されています。
http://saga-gis.org/geostat2011/05_saga_geomorphometry.pdf
が、まだその計算内容のまとまった説明を見つけることができません。先日のようにWikiで一つ一つ探さないといけないのでしょうか？また、これらの計算は何の目的で考案されたのでしょうか？何がわかるのでしょうか？これをツールとして使う分野は？地形学とも違うような気がします。知りたいですね。

とりあえず、正常に動きましたので、作業自体は問題ないでしょう。


これから計算内容を確認していきましょう。

地上開度

以前、LPデータに対し、開度の測定やウェーブレット変換ができるソフトは皆無と投稿したことがあります。http://phreeqc.blogspot.com/2011/04/arcgis.html

先日、「地上開度」について古い情報が引っかり、再度調べてみました。
検索で引っかかったのは、こちらの方のブログです。
http://kibouharuka.at.webry.info/201104/article_7.html

確認したところ、確かに現在も Dept. of Physical Geography, Hamburg で開発が続けられている SAGA （System for Automated Geoscientific Analyses) に地上開度のモジュールが装備されていました。
http://www.saga-gis.org/en/index.html

モジュールのソースはこちらです。positive openness と等価と書かれていますね。
http://saga-gis.cvs.sourceforge.net/viewvc/saga-gis/saga_2/src/modules_terrain_analysis/terrain_analysis/ta_morphometry/ProtectionIndex.cpp?revision=1.8&view=markup&pathrev=MAIN

ソースの内容が読めるのは、オープンソースの利点です。高さ/距離のarctanを求め、最大値を更新、8方向の平均を求めるといった内容のようです。非常に単純ですね。こんなことが全てのグリッドに対し順番にできるんですね。まあ、GISなら当たり前なのでしょうか？
このアルゴなら簡単に地下開度にも変更できますね。モジュールは自作できそうですのが、組み込めるのでしょうか？

しかし、このソフト、いろんなモジュールが組み込まれていますね。ArcGISだと地形解析オプションとVer.UP だけで100万程度、開度やウェーブレット解析の委託開発でさらに100万程度かかったような気がしますが。詳細は確認していませんが、ウェーブレット以外は大体同じような機能が入っているのではないでしょうか？無料で使えるなら、こっちを選ぶでしょう。

早速、基盤地図情報の5mメッシュを取り込み、地上開度を描いてみました。
正常に動作します。（が、こんなもん？といった感じです。）案外簡単にできますね。
手順は後日UPしましょう。


最近、傾斜区分図や地上開度図の作成を「地形解析」というようですが、できた「解析結果」を見ても、もうひとつピンときません。加工された図を使って人が判断するわけですから、判断しやすい材料に「データを加工する作業」と言ったところが妥当なのでしょう。解析と言いながら、問題に対する決定論的な結果を示せないところがピンとこない要因なのかもしれません。

個人的な経験では、地形から推定した地質構造は、本当に表層部でしか合わないと考えています。これは地形が地質構造を反映しないのではなく、人の推定能力が弱いためだと思います。地形から、構造の傾斜を考慮することが難しいのでしょう。経験上、深度100mとなると50%、つまり当たるも八卦・・・の世界です。そういう意味では、LPデータ解析マニュアルの対象となる地すべりは表層地形に現れる現象ですので、データ加工が役に立つのかもしれません。


ま、何事も加工する能力がないままでは、利用価値があるのかすら判断できませんので、とりあえずは手を動かしてみましょう。

地すべり監視のためのSAR干渉画像判読マニュアル

「地すべり監視のためのSAR干渉画像判読マニュアル」が国土地理院から出ています。

http://gisstar.gsi.go.jp/sar-jisuberi/manual-top.htm
http://vldb.gsi.go.jp/sokuchi/sar/index.html

合成開口レーダーは、約10年以上前に小池先生が地形解析（確かSTA）に使用されているのを見て、トライしようと考えたことがありました。ちょうど、LANDSATから脱却し始めた頃ですね。日本の場合、地表を植生が覆うため、LANDSAT画像の演算処理は適用しにくいことがあり、他の画像を模索していた時代でしょう。
いずれにしても、衛星画像は広域の解析が対象になるため、すぐに利用しなくなりました。通常業務で広域の構造物を扱うことはほとんどありませんのでしたし、当時、私がベクトル解析を深く理解していなかっこともあります。
ところが、国土地理院では広域から地すべり変動を抽出する目的で、SAR干渉画像を利用しています。やはり、広域を相手にせざるを得ないニーズもあるのでしょうね。

これを見て思いついたのですが、地すべりであれば、地上レーザースキャナかLPデータを1～2年毎に取得し、差分をとれば3次元の動きが分かりますね。
・・・・・いえ、何年も動きを追うくらいなら、対策してしまいますよね。逆に、何年も観測だけで良いなら、動きは少ないでしょうし。干渉SAR画像による地すべり監視とは、どのような地すべりを相手にするのか分からなくなってきました。

河川では、広域の監理（治水安全度評価）を目的として、LPデータが使用されています。
http://www.nilim.go.jp/lab/rcg/newhp/seika.files/lp/
ちょうど昨日の発表ですが、今回の震災前後のデータで比較した研究例もありますね。Civil3Dが使用されています。
http://www.mlit.go.jp/chosahokoku/giken/program/kadai/shitei.html
http://www.mlit.go.jp/chosahokoku/giken/program/kadai/pdf/shitei/shi2-01.pdf

LPで差分をとるというのも、もう新しい発想ではないのでしょう。