2011年12月31日土曜日

解析時の定数設定

「地盤の地震応答解析」を読み終わりました。

といっても、最後の方は頭に入りませんでした。実際に手を動かしながら読む必要がありそうです。

この本、なかなか良かったです。文章が上手なのでしょうか、語りかけてくるようでどんどん読み進めることができました。技術的にも細かいことが書かれており、他の解析に通じる箇所もあります。メッシュサイズの決め方などは地震応答解析特有でしたが、他の解析でも気をつけなければならない点です。
実際に読みながら手を動かしたわけではないので、詳細は把握していないと思いますが、調査時に気をつけるべき個所、地震応答解析の概要や限界、困難さは理解できました。解析やその調査に携わる心づもりをする上でも、良書だと思います。

1つ気になったのは、弾性係数の拘束圧依存性があるのかどうかについて。
依存する場合もあれば、影響の少ない場合もあるとのこと。結局はそのあたりも意識して調査方法の選択・モデル化しないといけないということなんでしょうね。本では以下のように記載されています。
「深いところでは拘束圧の変化に伴う弾性係数の変化はそれほどないが、浅いところではかなり変化するので技術者は、簡単な計算でその影響を確かめておくなどの注意が必要である。」
まあ、他の解析でも現況再現でのパラスタ作業は、上記に該当するのでしょう。定数設定として共通する重要な作業です。解析のための定数設定やそのための調査は「解析がわからない」では済まないのです。

古い技術者は自分の経験や流儀で解決しようとする傾向があるようです。土木の場合、良い点も多くあります。が、「解析のような難しいことはわからない。」と平然と言うようになれば、技術者卒業ということでしょう。

2011年12月29日木曜日

Dtransu と CUDA

DtransuのCUDA化です。

作業をお願いしていたプロより連絡がありました。
C++にまで移植していただいたのですが、結果は精度面でOUT。

惜しいことに、デバイスエミュレーションで動かすと正しい値が出るようですが、実際にGPUに載せると正しい値が出ないとのこと。プログラミング自体は問題ないのでしょうね。TESLAで駄目なようですから、プロが仰るようにコンパイラかライブラリ側に問題があるのでしょう。素人でも分かりやすい説明でした。

ここまでしていただいて、本当に感謝です。
今回は残念でしたが、CUDAやコンパイラのVer.UPで解決できそうですね。気長に待ちましょう。

2011年12月28日水曜日

ソリッドのメッシング

ソリッドを全選択し、オートメッシュ機能で一気に切ってみました。


メッシングのテストですので、高次要素は使っていません。節点数: 24404、 要素数: 118925、所要時間: 3時間13分でした。なかなか良い感じで切れていますね。


すべり面を下から覗くとこんな感じです。
ボクセルメッシュ(GEORAMA でメッシング、SoilPlus で読み込み)↓


今回のメッシュ(GEORAMA でソリッドモデル作成、GTS で Parasolid 読み込み後、オートメッシュ)↓


ぜんぜん違いますね。素晴らしい。
これからは地下水でも妥協の必要が無くなりました。

AutoCAD 2012 のプラグインにより、ソリッド読み込みが実現したのですが、今後はCivil3Dにも搭載されるのではないでしょうか?


2011年12月27日火曜日

ソリッドからSTEP, Parasolid (AutoCAD 2012)

midas GTS は Nastran ファイルを読めません。

読めるのはSTEP, IGES, Parasolid DXF(ワイヤーフレーム)などのジオメトリ系のみです。GEORAMA から吐き出したボクセルメッシュは読めません。また、SoilPlusと同じプレですので、境界面を書き出して取り込むのも問題があります。
http://phreeqc.blogspot.com/2010/12/blog-post_16.html

CAD のソリッドモデルを STEP などに変換できないかと調べていると、Civil3d 2011 で AutoCAD 2011 用の subscription tool を適用し、IGES へ書き出せることが分かりました。

早速、DL して使ってみました。
が、駄目でした。書きだせるのですが、GTS 側でソリッドとして認識しません。コンパウンドとして認識されるのですが、分解するとサーフェスになります。地表面がメッシュで分割されてしまうので、後で地質毎にサーフェスをまとめるのが大変です。2012では解決しているのでしょうか?

このままでは先に進みませんので、もう少し調べてみました。
結果、AutoCAD 2012 (Cicil3Dではありません)で Inventor Fusion plug-in を適用すると STEP, IGES, Parasolid で書き出せるとのこと。早速、インストールしようかと思いましたが、既に使われている方が変換して下さいました。

変換していただいたデータは STEP, Parasolid の2種。
GTS で読み込むと、なぜかスケールが狂っていましたが、きちんと solid として認識されました。もう少し使い込まないと分かりませんが、パっと見は OK です。
GEORAMA +Civil3Dで作成したソリッドモデルを変換しますので、上記のような不具合がなくなります。(作成法はコチラ>>>http://phreeqc.blogspot.com/2010/11/civil3d2010-georama2010.html

あとはメッシングにどれだけ時間がかかるかですね。
オートメッシュで試してみましょう。ボクセルよりは綺麗に切れると思いますし、要素が少なくて済むかもしれません。期待大です。

2011年12月26日月曜日

地形判読

「深いすべりが調査前に分かるか?」という疑問に対し、いろいろ意見を頂きました。

個人的な経験上、多くは10~30m、深くても50m程度にすべり層が出てきます。
が、それよりも深い70m、80mなどのすべりの存在が調査前に分かるか?それらに対し調査計画を立てることができるか?という点についての意見です。

「現在の技術水準で分からないという方が困難」など、多くの思想が聞けました。

個人的には、深いものに対し調査計画を提案する自信がありません(変状が明瞭であれば提案しますが)。対策も思い浮かびません。しかし、今までに深いすべりを経験して来られた方は、自信を持って提案されるようですし、調査のノウハウもあるようです。やはり、地すべりは①地形判読を行って規模を推定し、②そこを実際に歩いて変状や範囲を確かめ、③調査計画を立てる、という順序と、今までの経験が重要なのでしょう。

この中で、再度思ったのが地形判読の重要性。
一番最初に行うわけで、地すべりの規模が大きいほど、ミスリードの影響が後工程に大きくのしかかります。いくら高価なLPの取得や地形判読のための様々な図面を作っても、技術者の判読力が低ければ価値がありません。

時代は変わったようで、変わっていないと思います。
地形判読に関するツールは増え、新たな資格も作るようですが、本質はそこにないと思われます。

もう一度、基本に戻って、地すべりも含めた地形判読の復習をしておく必要があります。

EXIFの位置情報をCADへ

写真の EXIF には、位置情報(ジオタグ)が含まれています。

携帯やスマホにも GPS が付いていますので、ネットに UP された写真がどこで撮られたものか判別できる場合があります。また、位置情報を使って撮影位置を地図に示すことができるアプリも多くみかけます。

GPS がついてないデジカメでも、track データと撮影時間を比較して、EXIF に位置情報を埋め込むことが可能です。先日もカシミール3Dで位置情報を埋め込んだところまでUPしました。

問題はこの先で、これをCADへはめ込む方法が分かりませんでした。位置情報を読み取って、撮影場所をCAD上で図示したかったのです。

調べたところ、頭脳RAPIDPRO Civil3 で可能なようでした。
http://www.photron.co.jp/products/cad/pro_civil/function02.html
が、持っていません。

仕方ないので、フリーの「F6 EXIF」で位置情報の一覧データ(csv)を作り、Civil3Dでポイントデータとして読み込みました。
が、ここで問題。csv が WGS84 の LL (緯度経度)データとなってしまいます。track データを UTM でなく LL で書き出せば、ばっちり合います。が、平面図は通常 LL で表記されていないため、そのまま載せると位置が大きくずれます。これは何とかならないのでしょうか?

平面図の測地系・・・JGD2000
GPSで選択した測地系・・WGS84
EXIFの位置情報・・・・・・・WGS84(LL)

WGS84(LL)をcsvからCAD読み込みの段階で変換したいですね。Civil3Dで可能なような気もしますが、やり方が分かりません。


考え方を変えて、位置情報から Garmin の Waypoint データを作成し、それをDXF変換することにしました。GPSBabel を使用します。これ、かなりの種類のファイル形式を扱えます。
http://www.gpsbabel.org/
この後、GPS の付属ソフトを使用すれば、WGS84 (UTM) で DXF を書き出せます。

結果、平面図・GPSデータ・写真位置がCADで重なりました。
ちょっと力技っぽいですが、変換、変換で下地ができました。

2011年12月24日土曜日

ハンディーGPS

今年最後の現場は落石調査になりそうです。

2つの山、5つのエリアを数人で手分けして担当することになりました。私はそのうち1つのエリアを担当しています。
設計者からの強い勧めと、先行していた先輩から「GPSを持っていけ」と言われたため、今回初めてハンディーGPSを山の踏査に携行しました。GARMIN GPSMAP 60CSx です。古い機種です。

ハンディーGPSは10年近く前から会社にあります。今まで、道路際の写真を多く撮る際に使ったことはありますが、山の中の踏査では使ったことがありません。山の中で衛星を補足できるのか?という疑問と、「GPSに頼らなくても山は歩ける」といった地質屋としての(ちょっとした)自負からです。

ところが今回、携行してみて驚きました。期待以上に追いかけてくれます。
普通に腰にぶら下げて gps track を拾っていたのですが、精度は±10mに入っていたと思います。冬で草木が枯れていたからかもしれません。藪こぎに近い状態や、倒木の下を這ってくぐり続ける個所もあったんですけど。

まあ、こうなると期待してしまいます。
落石源や崩壊箇所など重要な位置は現場で図面に落としますが、そこで撮った写真の番号や、その間の写真撮影位置は落とすことをやめました。今までは全ての写真番号と位置を現場で図面に落としていました。これで踏査スピードがあがりますし、重要な個所に時間をとることができます。

あとは帰社後にデータを取りこむだけ。
GPSの waypoints と track データはカシミール3D と CAD(付属ソフトで変換:UTM WGS84の座標系をもったDXF) へ取りこみます。写真はカシミールで  track データ を利用し撮影位置に落とします。これで、自分の書き込んだ図面をスキャンして CAD へ貼り付けると、下地は完成です。あとは自分の書いた落石位置で track データを(頭の中で)補正して整理すればOKですね。

今回は落石調査でしたが、地質調査時のルートマップにも利用しやすいと思います。便利な道具だったんですね。なぜ使わなかったのでしょう。
現行機種はさらに感度が良くなっているとのこと。室内でも補足することができると聞きましたが本当でしょうか。俄然、物欲が出てきます。


写真位置をCADに落とせば完璧です。が、やり方が分かりません、どうするのでしょう?
これは、また後日ですね。

2011年12月19日月曜日

孔内水平載荷試験の再現 その4

以前、孔内水平載荷試験の結果を数値計算で再現すれば、ダイレイタンシー角ψを同定できるのでは?と考えたことがあります。その時は、SoilPlusが関連流れ則しか扱えなかったため、同定作業自体ができませんでした。
http://phreeqc.blogspot.com/2011/08/3.html

今、MIDAS社からGTSをお借りしていますので、そちらで確かめてみることにしました。プリはSoilPlusもMIDAS社から提供を受けているので、操作性は同じです。ただ、こちらは軸対象で水平方向の荷重を扱えますので、2次元で計算できます。3次元より計算時間が早いので、すぐに結果がわかります。

計算結果はこんな感じです。ゾンデ60cm間に1000kN/m2を10分割でかけました。

変位のグラフの例です。


ゾンデの先、余掘り分15~20cmをモデル化し忘れたので、メッシュを切りなおそうとしました。が、同じ領域でメッシュを切りなおした場合、前の節点の情報をクリアしないようです。バグでしょうか?
仕方ないので、そのまま比較することに。

比較は弾完全塑性の関連流れ則、非関連流れ則、ひずみ軟化で行いました。
パラメータはダイレイタンシー角ψと単体以外、試験での実測値や算出値を使用しています。パラスタはしていません。変形係数は試験でポアソン比を仮定して得られた値、K0は土被り圧とP0より算出、φ'peakは深川の方法、せん断抵抗力(ピーク、残留)はσr-σθで求めた値を使用しました。なお、ψは仮にφ-30°として設定してみました。
ピーク、残留をφで指定できないのは、ソフト上の制限です。実務でGTSのひずみ軟化を使用する場合、砂質土をモデル化するには、深度毎に分割しないといけないようです。

で、結果は以下の通り。

うーん、弾性変形の傾きがやや立っていますね。試験時のポアソン比の設定が良くなかったのでしょう。逆に、試験でGがわかるので、パラスタで適切なポアソン比(=変形係数)が設定できそうです。

この結果では、ひずみ軟化が最も合っていますね。試験で求めた値を使うとばっちり合うのでしょうか?結論を得るにはもっとやってみる必要がありそうです。非関連も、まあ、許容範囲でしょうか?関連は駄目(つまり、SoilPlusでこの土の解析はダメ)ですね。

ちなみに、これは載荷区間の変位の平均値をグラフ化しています。実際の孔内水平載荷試験も、土砂の場合は水の注入量から平均的な変形量を算出していますので、こういった簡易な整理も良いと思います。

載荷区間の下の方は、余堀をモデル化していないミスで変位が抑制されています。試しに下15cmをオミットして平均をとると、以下の通り。


それほど変わらなかったですね。ひずみ軟化が良さそうです。


除荷、再載荷をしてみました。


感動!ちゃんと除荷時に変位の大きい側で低下し、負荷時に小さい側で上昇するループを作っています!実際の土もこういった挙動になるのですが、理由はちゃんとあるんですね。構成則を作った昔の人は偉い。すばらしい。


結果としては当ケースの場合、ひずみ軟化が良いでしょう。ψの同定を目的に始めましたが、思わぬ収穫でした。孔内水平載荷試験の結果より、適切な土質定数(の範囲)を設計側に助言することは可能ですね。

うーん、しかし、この結果はひずみ硬化ですよね?弾完全塑性やひずみ軟化としてモデル化しても、結果はひずみ硬化のような挙動を示しています。硬化と軟化、弾完全塑性の区別が良く分からなくなってきました。
まだ理解していない証拠です。もう一度、勉強し直す必要がありそうです。



不飽和透水試験で求める飽和透水係数

雨水浸透施設技術指針に、不飽和帯での透水試験(現地浸透試験)が掲載されています。
http://www.arsit.or.jp/book/tyousa.html

式は非常に単純です。

ダルシー則を変形するとq/k=iAとなりますが、iAが一定であれば、q/kも一定です。基準では、いろんな水深・形状(iA)の浸透試験結果を整理しておけば、比浸透量(q/k)が求まる、さらに試験で最終浸透量(q)を測っておけば、不飽和浸透試験から飽和(に近い)透水係数(k)が求まると考えたようです。
実際、これらの検証を数値解析で行ったとあります。私はここが面白いと思いました。実は、私も不飽和浸透試験を数値解析でパラスタすれば、不飽和浸透特性が求まるのではないかと思い、いつかやってみようと思っていました。既に基準に取り入れられていたんですね。

しかし、この検証の具体的内容がどこにも書かれていません。文献も示されていませんので、適用範囲がどの程度のものなのか、信頼性はどの程度のものなのかがわかりません。仕方ないので、発売元の雨水貯留技術浸透協会に問い合わせてみました。
結果、具体的な内容はわかりませんでした。出典は東大のドク論だそうですが、東大の論文データベースを検索しても出てきません。発表論文を検索しても、基準以上のものは出てきませんでした。残念ですね。
ただ、対応していただいた方の話では、不飽和層の実験を全国から集めて検証したそうなので、飽和層への浸透は範囲外ということでした。もともと、浸透施設が飽和層を対象としたものではないとも言われていました。

非常に惜しいのですが、適用範囲や信頼性はわからないままです。現段階では雨水浸透以外で使用することは控えておいた方が良さそうです。

2011年12月17日土曜日

雨水浸透施設

雨水浸透施設が国交省や東京都で推奨されているようです。
http://www.mlit.go.jp/common/000113727.pdf
http://www.tokyo-sougou-chisui.jp/shishin/

河川のプロに聞くと、やはり河川への流入量を減らす目的で、推奨されているとの事。
ただ、私は実物を見たことがありません。

先日も、建築屋さんがこの施設について話されているのを聞きました。
河川堤防の横で、帯水層まで鉛直のドレインを作り、その上の地表部に既製品を設置する計画のようです。確かに、ポンプ施設を作るよりは安上がりでしょう。

ただ、地質屋からみると、非常に危険だと思います。堤防天端より堤内側が5mも低い状態なのです。計画高水位は天端より1.5m程度低いとすると、3m以上のヘッド差ができてしまいます。
洪水時に堤内側の井戸から水が吹いたり、パイピングがあったという話を聞きますが、雨水浸透施設でも帯水層までドレインを入れると、同じような現象が起こるでしょう。

雨水を処理するために、安い工法、お客様に喜ばれる工法をと考えられたのでしょうが、一つのことに目を奪われると、他が見えなくなる例ですね。
何事も引いてみたり、近寄ったりを繰り返す必要があります。

2011年12月16日金曜日

ボスポラス海峡 CM

大成建設さんのCMです。



以前にも書きましたが、尊敬する技術者の一人が現地に赴任されており、個人的に関心の高い事業です。


秒速5センチメートルも綺麗でしたが、これも綺麗!豪華ですね。

圧密沈下

圧密沈下について相談がありました。

e-logp、mv、cv の評価で迷っているとのこと。
確かに、データを重ねたグラフを見ると、ばらついています。

この中で、pcが明らかに上載圧依存となっており、e-logpを使用するのであれば、ある深度毎に区切って代表的なラインを設定する必要がありそうでした。mvの過圧密側もばらついて見えるのは同じ理由でしょう。

いずれにせよ、一軸のみ上載圧依存を考えるのは片落ちで、pc、e-logp、過圧密側のmvなどもチェックしておく必要があります。基本的な事ですが、試験数が増え煩雑になると、手を抜きがちになる所です。港湾基準だと、性能設計に基づいた決定法がありますが(これもすっかり忘れてましたが)、今回は鋭敏粘土上の河川堤防を扱っていましたので、設定値は深度を細分した平均的な値になりそうです(これも変ですけど)。
http://www.web-gis.jp/e-Forum/2007/033.PDF

最近では数値計算で解く方法に目を奪われていましたが、その定数がどのように設定されたのか、地盤のどこで試験をしたのか、各層を代表しているのか、といった根本的なチェックが重要ですね。

グイっと引き戻された気分です。

2011年12月14日水曜日

LEM, SRM, SAM その2

LEM, SRM, SAM の比較です。

適当なモデルで計算してみました。まずは、非関連流れ則。


自動で書き出される比較表の一部です。計画安全率は自分できめることができますが、今回、特に意味はありません。便宜上表示しているだけです。

結果はFellenius より Bishop、SRM より SAM で安全率が高く評価されます。後者は円弧を仮定している影響ですね。FEM とLEM は単純比較できないでしょう。SRMでの結果が円弧に近いと、案外、Bishopと同じような安全率が出ています。


関連流れ則でも実施しました。
結果、このモデルでは SAM、SRM ともに安全率が0.01前後上がります。ひずみが集中しにくくなる影響でしょう。教科書に載っている注意点ですね。土の評価が重要になります。


モデルを変えるとコチラ。材料値は上のモデルと同じで、非関連流れ則です。


解析領域や結果に問題ありですが、傾向は同じです。
異なる点は、SRM で非円弧が顕著になると、安全率が円弧よりかなり低くなるということでしょうか?まあ、当たり前ですか。


通して考えると、Fellenius だけはやはり異常ですよね。あらためて実感しました。Felleniusの計画安全率1.2を考慮すると、、他の計算手法で1.2は駄目ですね。駄目なのはわかりますが、ではどうしたら?というと実務的には以下のような答えしかありません。

海外ではどう評価しているのでしょう?

2011年12月13日火曜日

LEM, SRM, SAM

SoilWorks を触り始めました。

最初は GTS を触っていたのですが、思うように動きません。モデルも私が持っているのは SSRFEM には大きすぎるので、後回しにしました。

とりあえず、SoilWorks の斜面安定モジュールを動かすことに。
マニュが手元になかったのですが、2次元ですし、メッシュ作成や解析ケースの設定はSoilPlusと同じインターフェースということも手伝って、簡単にできました。

ところが、いざ計算をしようとすると、手法が3つあります。LEM, SRM, SAM です。Limit Equilibrium Method, Strength Reduction Method はわかります。が、SAM ってなんでしょう?ヘルプには「応力場解析技法」とありますが Stress Field Analysis Method の訳でしょうか?検索しても引っかかりません。このソフト、微妙に日本語変なんですよね。「斜面安定分析」とか。どうせなら全部英語にしてくれたら良いんですが。まあ、本質には関係ないところです。

マニュアルを DL して確認しましたが、SAM について書かれていません。

が、よく見てみると、極限平衡法の説明がFEMを使った計算になっています。たぶん、これのことでしょう。つまり、与えた物性値によって応力場を計算し、仮定した円弧に沿う抵抗力の積分値をせん断応力の積分値で割って安全率を出す、それを繰り返し最小安全率の円弧を見つけるという手法です。

こんな手法があるんですね。
どういう使い道があるんでしょうか?

うーん。思いつきません。

2011年12月11日日曜日

自然由来重金属マニュアルのリスク評価

先日、国土交通省が公開している「建設工事における自然由来重金属等含有岩石・土壌への対応マニュアル (暫定版) 」のリスク評価について相談がありました。
http://www.mlit.go.jp/sogoseisaku/region/recycle/pdf/recyclehou/manual/sizenyuraimanyu_zantei_honbun.pdf

公開当時、このマニュアルを読みました。リスク評価手法に興味があったのです。しかし、その箇所の内容は安全側かつ簡易なもので、本質を見極める上で満足できるレベルではありませんでした。1~2次元主体であったり、100年間溶出量が一定のプログラムであったり、「レベル2でのリスク」=「地下水中の濃度が基準を超えるかどうか」=「通常の移流分散計算の簡易版」であったり。移流分散方程式は遅延の掛け方が違ってますし、濃度一定の1次元理論解も左辺が違います。すぐに放り投げました。

私は土壌汚染問題のプロではないので、法やマニュアルに従った調査法について詳しく知りません。しかし、解くべき問題がはっきりすれば、science を応用した解決策は提案できるでしょう。それが法的調査と一致すればBESTですが、異なることもあると思います。残念ながら、マニュアルのリスク評価は後者のようでした。

相談は以下の2点。
・お客様が要望されているマニュアルのリスク評価をするためにどんな調査が必要か?
・もっと良い評価手法がないか?それに必要な調査は?

両方、お客様への提案書としてまとめました。
2点目は国交省マニュアルの範囲で試験法を選択するケースと、それ以上の調査(費用)が必要になるケースがあります。当然、前者が過大設計になります。後者が望ましいのですが、目先の費用の面で両手をあげて喜んでは頂けないでしょう。御要望のマニュアルがなければ良かったのですが。

お客様の「国交省」「マニュアル」のとらえ方次第です。

皆既月食

皆既月食ですね。

欠け始めから皆既食まで、デジカメで撮影してみました。
数年前のバカチョンですが、18倍光学ズーム付いてます。

カメラは詳しくないのでいろんなモードを試してみました。夜景モード、露出補正、シャッター速度云々。マニュアルでも設定可能なようですが、ド素人ですので、全部カメラ任せの自動モードです。

初めて月や星を撮ったのですが、後でPCで見ると、想像以上にきれいに撮れました。もちろん、ボケてない程度ですので、自己満足です。が、天文ファンの心理がわかります。カメラ様々。

人間にとって「災」に映る自然現象がクローズアップされた年ですが、こういった自然現象は、誰でも、いつでも大歓迎ですね。

2011年12月9日金曜日

SAGA 2.0.8 + 基盤地図情報(DEM)

SAGAで基盤地図情報(DEM)を読み込んで Protection Index (≒地上開度)を作る手順です。

① 基盤地図情報DEMのダウンロード〔JPGIS(GML)形式〕と変換
http://phreeqc.blogspot.com/2010/11/givil3d2010-2500-25000.htmlの①~③を参考に、DLしたファイルをコンバーターで変換します。「標高メッシュをシェープファイルに変換」です。コンバーターの表示が綺麗ですね。



② SAGAで読み込み、クリギングでグリッドデータに変換
・File-Shape-Load

・krigingは何でもよいでしょう。地形ですから。


・バリオグラムを設定。地形ですが、案外、相関があるんですね。意外でした。密度が細かいので自動でも良いかもしれません。

・グリッドサイズはDLしたメッシュサイズに応じて入力すれば良いでしょう。


・終了するとビープ音が出ます。画面には何も表示されません。


③ Protection Index の演算
・左上の Workspace のタブから Dataを選びます。
・Workspace に読み込んだファイル名が出ますので、ダブルクリック。画像が表示されます。


・Module から Morphometric Protection Index を選択。


・ L を設定します。

・演算が終わればビープ音が鳴ります。Workspaceに追加されたProtection Index をダブルクリックし、表示します。


これで出来上がりですが、地上開度とは角度の取り方が異なります。が、理屈がわかっていればOKですね。



2011年12月8日木曜日

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

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

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

合成の手法がマニュアルに書かれていません。合成の演算は山ほどあります。計算手法を示すことができないのであれば、せめて何のソフトの、どの演算を使ったかを書いてないと、追いかけられません。ましてや開度図の階調を見やすいように修正してから合成だと、純粋な計算結果を演算するのではなく、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 のいろいろな演算が何を意味するのかは、興味が残るところです。何の分野でどのような目的に使用されているのか分かりませんが、教科書や文献はあるはずです。こちらは今後の「技術力」として調べてみましょう。

2011年12月7日水曜日

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

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

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

ウェーブレット関数とは何ぞや?の状態ですので、ネットで調べてもよく分かりませんでした。これだけで1冊の本になるほど多肢にわたっています。
しかし、マニュのメキシカンハット関数やウェーブレット係数の算出式は非常に簡素なので、これに限ってはできるような気がしました。

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


あれ?最大値が11ページのグラフと違います。おかしい。
ネットで調べてみると、ψ(x, y) = (2-x^2・・・ の箇所が(1-x^2・・・の形になっている式が多く見られます。y = 0のときに(1-x^2と同じ形にならないといけないはずですから誤植でしょうか?
とりあえず2を1に修正してみると、



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

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

続きは後日。

2011年12月6日火曜日

SAGA

SAGAの続きです。

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

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



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


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

2011年12月4日日曜日

地上開度

以前、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データ解析マニュアルの対象となる地すべりは表層地形に現れる現象ですので、データ加工が役に立つのかもしれません。


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

2011年12月2日金曜日

剛性の拘束圧依存

E=2G(1+ν)         EはGの1乗に比例
G0=ρVs^2          G0はVsの2乗に比例
Vs=80~100・N^(1/3)    VsはN値の1/3乗に比例
Nc=1.7N/(σ'v/100+0.7)  N値は上載圧の1乗に比例
Nc=√(98/σ'v)・N      N値は上載圧の1/2乗に比例

これらの式が意味するところは、Ncが同じ材料の場合、せん断剛性や弾性係数が拘束圧の1/3~2/3乗に比例するということでしょうか?

吉田望「地盤の地震応答解析」p96、97あたりに拘束圧依存の説明が書かれています。図7.4をみると、確かにそんなものなのかなあという気がします。べき乗までひずみレベルに依存するのであれば設定は難しいですね。どこかで割り切らないと。

実務では0.5を使われているのでしょうか?機会があればプロに聞いてみましょう。

動的緩和法

以前、講習会でひずみ軟化について講師と雑談した際に、「動的緩和法」を勧められました。

詳細は知らないのですが、時間がたてば、ある点で減衰により静止するといった考え方で、静的解析に運動方程式を利用する手法です。

その時は何か抵抗がありましたし、なぜそのような手法を勧められたのかが分かりませんでした。
が、最近地震応答解析の本を読んでいて、「あ、なるほど」と思うところがありました。
ひずみの進行に応じて急激に剛性が低下するモデル。そう聞くと、動的変形解析をされている方は地震応答解析を思い浮かべるでしょう。そちらのソフトは市販されているわけですから、少し工夫すればソースを変えることなくひずみ軟化が表現できると考えられたのだと思います。

いろいろやってみるべきですね。