Cross-correlation coefficient between 3D arrays

 x,y (2次元) + time (1次元)の相関性をどう定量化しているか？

思いつくのは相互相関係数。動画のマッチングなどで使われているかな？と思い調べてみましたが、見当たらず。需要はないのでしょうか。

3Dで引っかかったのはコレ↓

Correlating two 3D arrays in Python
https://stackoverflow.com/questions/28861995/correlating-two-3d-arrays-in-python

なるほど。１次元にすればよいのか。確かにこれで目的は達せます。偉いねえ。

scipy.ndimage.correlate だと３次元のまま関数を計算できましたが、１次元の方が良さそう。

相関性の良さを比較するには正規化が必要です。が、正規化相互相関は各ライブラリに実装されていないようです。OpenCV にはありましたが、画像なので int 扱いになります。

難しい式ではないので、Zero Mean Normalized Cross-Correlation を関数として記述。2つの3d array を渡せば係数を返してくれるようにあしました。
これでクリアーです。

LiCSBAS

LiCSBAS: An Open-Source InSAR Time Series Analysis Package Integrated with the LiCSAR Automated Sentinel-1 InSAR Processor
https://www.mdpi.com/2072-4292/12/3/424

GitHub の wikiに従い install。

調査範囲を設定し、あとはbash を叩くだけ。これだけで DL から複数の計算まで止まらずに走ります。簡単です。

LiCSAR で2時期の差分画像までは処理済みですが、６～12日毎で数年分のデータがあるため、DLだけで1時間以上かかりました（ラップトップ＋4Gルーターです）。２度目からは不足分のみ DL される仕様になっていたので速い。ありがたいです。

計算自体は速いですね。ラップトップでも問題なし。ラズパイでも動くとのことでしたので、かなり軽量化に配慮されているのでしょう。1度だけ step02 で停止していましたが、安定性も問題なさそうです。

計算が終わったら、コヒーレンス等のマップを表示。その中から目的の位置をクリックすると、対応する時系列グラフが描画されます。異なる場所をクリックすると、それに応じて時系列グラフも変更されます。これは素晴らしい！

時系列値をテキストに吐き出すスクリプトが付属。変位速度マップを geotiff にしてくれるスクリプトも入っていました。GIS に持ち込むと、位置を視認できます。これも必要でしょう。 

LiCSBASを利用して、国内８か所の変位速度マップを描いてみました。１平方㎞程度の広さがあれば、判定し易くなります。平野や埋立地などは特に判定が容易。山はダメですね。

時系列解析に関する一通りの機能がそろっています。しかも軽量、安定。
作者らに感謝です。

 

InSAR time series analysis その2

時系列解析は、PSI 法 と SBAS 法に大別されます。
後者は聞いたことがあるけれど、中身を知らないなーと思って調べると、おおきな勘違いでした。GPSの誤差補正情報（SBAS：Satellite-based Augmentation Systems）ではなく、Small Baseline Subset algorithm です。

干渉SAR 時系列解析プロトタイプシステムの開発, 国土地理院時報
https://www.gsi.go.jp/common/000207740.pdf

GSITSAR
本システムでは，ある1 つの画像を共通のマスターとして干渉画像ペアを作成する解析方法（以下「シングルマスター」という．）と複数の画像をマスターとして干渉画像ペアを作成する方法（以下「マルチマスター」という．）を選択可能である（図-2）．一般的に，前者はPSI 法で，後者はSBAS 法で用いられる．

国交省さんから GSITSAR が公開されるとありがたいです。

干渉SAR 時系列解析による地盤沈下の検出, 国土地理院時報
https://www.gsi.go.jp/common/000079938.pdf

干渉 SAR 時系列解析は，PSI（Persistent Scatterer Interferometry ）と SBAS（Small Baseline Subset algorithm）という二つのカテゴリーに分類される．
PSI は全解析期間で後方散乱特性が変化しない点（PS 点）を抽出し，それらのピクセルのみで変動を推定する（Ferretti et al.，2000，2001）．PS 点においては，空間及び時間干渉度低下は小さく，長期間かつ長基線長のペアでも干渉度低下によるノイズが小さい干渉結果を得ることができる．人工構造物等のように長期間変化しない物体がPS 点となり得る．
一方，SBAS は，短い垂直基線長及び短い撮像日間隔のSAR 干渉画像を多数作成し，各観測時の変動量を推定する（Berardino et al.，2002）．これにより，空間及び時間干渉度低下の影響を最小限に抑えて，時系列的な変動を抽出できる．

LiCSBAS は、NSBAS となっています。SBASとの差異を理解できていません。

Nationwide urban ground deformation monitoring in Japan using Sentinel-1 LiCSAR products and LiCSBAS
https://progearthplanetsci.springeropen.com/articles/10.1186/s40645-020-00402-7

・Cバンドの Sentinel-1データがベースになっているため、植生の密な地域では十分なコヒーレンスを保つことができない。
・Lバンドであれば期待できるが、観測頻度が低く（10回/年以下）、詳細な変形の時系列を検出することは困難。
・ALOS-4、NISARなど、次世代のLバンドSAR衛星が登場すれば、森林地帯であっても詳細な変形モニタリングが可能になる。

LiCSBAS は Sentinel-1の利用です。
融雪のための組み上げ→冬季の沈下をいくつかの地方で検出できています。夏場の農業用の組み上げで、沈下を起こしている地域もあるようです。確かに、（粗くても）密に変化を追うことが可能であれば、水位変化を面的に把握することも可能でしょう。このようなシステムが公開されてきたので、水位のインバージョンなどを研究される方々が現れたのでしょう。

**************************************

最初の文献の謝辞
「解析に使用したSAR衛星「ALOS-2」のデータは，ALOS-2に関する国土地理院とJAXAの間の協定に基づき提供されました．ここで使用しただいち2号の原初データの所有権は，JAXAにあります．「ALOS」のデータは，国土地理院とJAXAの「陸域観測技術衛星を用いた地理空間情報の整備及び高度利用に関する協定書」に基づき，国土地理院がJAXA  から購入したものです．「だいち」データの著作権は宇宙航空研究開発機構（JAXA）及び経済産業省にあります．」

 国交省さんでも研究には購入が必要なのですね。

InSAR time series analysis

先日の水位の話を聞いてから、連休中に InSAR の 時系列ツールの状況を把握しておこうと考えていました。

ANATIS
JAXA さんが開発を進められていたツールです。2019年に有償で利用者を募集されていましたが、その後の話は聞こえてきません。検索しても出てこないので、検証中ということでしょう。
経産省さんの下、Tellus でも ALOS2 のデータは蓄積されつつあります。が、解析に自由に利用できる状況ではありませんし、量も十分ではありません。従来通り ALOS2 データも解析ソフトも有償で購入してください、となると利用促進は期待できないでしょうし、海外に勝てない状況は続くでしょう。

SARPROZ
MATLABを利用したGUIを備えたコードです。HPのきれいな絵と時系列グラフに惹かれて使用しました。手順が分かりやすく、unwrap も実装されています。が、私の環境では最後まで動きませんでした。2時期の DInSAR としては初めてきれいな絵を得られたので、残念。
クラウド上のデータをサーバー側で解析するのではなく、データをすべてDLしないといけない点は、やや古い仕様です。

GEE
これはまだ SLC、DInSAR に未対応でした。

LiCSAR: An Automatic InSAR Tool for Measuring and Monitoring Tectonic and Volcanic Activity
https://www.mdpi.com/2072-4292/12/15/2430/htm
これは素晴らしい。このようなシステムが稼働していることを知りませんでした。蓄積され続けるSentinel-1 のデータを自動処理し、DInSAR などの結果をDLできるよう公開されています。速度面でも可能なのですね。感服です。
https://comet.nerc.ac.uk/comet-lics-portal/

This interactive map provides links to EU Copernicus Sentinel-1 InSAR products available for download from COMET-LiCS. Interferograms and coherence maps have been produced automatically using the LiCSAR processor, which builds on the Gamma SAR and Interferometry software.

速度計算・時系列表示に関しては、LiCSBAS を利用とのこと。
LiCSBAS: An Open-Source InSAR Time Series Analysis Package Integrated with the LiCSAR Automated Sentinel-1 InSAR Processor
https://www.mdpi.com/2072-4292/12/3/424
これも画期的。python コードを bash でつなぎ、注目個所の時系列データを自動で作成してくれる内容です。SIP 1期の研究を軽く超えたのでは？
国土地理院の地理地殻活動研究センターの方が第一著者です。気付いたときは「国がSentinel で良いのか？」と思いましたが、データ整備に取り組まれているのは経産省さんなので問題ないのでしょう。精度によっては、宇宙開発利用大賞を受けるに値する成果だと思います（組織的にありえないでしょうが、開発としてはここ数年でダントツでしょう）。LiCSAR の著者でもあるので、将来的に地理院地図にHPと同様のシステムが実装されるかもしれません（マップ表示のみのデモは既に公開されています）。楽しみです。

令和3年度 地質CIMの対応 その3

Civil3D 2021 + GEORAMA で作成した3Dモデルに、Navisworks 2021 + Navis+ で属性を付与する手順です（ボーリング＋サーフェスモデル）。
ガイドラインの属性例（表30）の約6割を付与できます。

GEORAMA
・「境界テーブル」を開いてセルを全選択し、CTRL+C。
・EXCELにペースト。
・不要列（B~Q）削除。
・2行目に'0'等を挿入。
※文字列:0、実数:2、整数:3、日付・時刻:4、リンク(外部参照ファイルやURL):5
・タイトル行をつける（A列は'General:Layer Name'）
・Navisworks データの保存予定先と同階層に'ATTRIBUTE'フォルダを作成し、csv保存。
※電子納品基準に沿ったフォルダ名。
・「ボーリング」を csv 出力。'Borlist.csv'が出力される。'ATTRIBUTE'フォルダに移動。
・最後の列に1行目’柱状図’、2行目’5’、3行目以降’SOURCE¥○○.xml’として柱状図xmlを指定。
※Navis データ保存予定先に'SOURCE'フォルダを作成し、柱状図を集めておく。
・「３Dダイアグラム」で地表面、断面、境界面、地表面、3Dボーリング、3Dソリッドボーリングなどを選択し、表示。
・必要に応じ「オプション」-「ボーリング」で柱状図の幅、3Dソリッドボーリングの幅、旗揚げの高さ、文字の大きさ等を調整。
・問題なければ「３Dダイアグラム」で書き出し。

Civil３D
・書き出された dwg を開いて、地表面のコンター間隔を調整、保存。

Navisworks
・「ホーム」タブ-「プロジュクト」-「追加」より書き出された dwg を追加。
・保存。
・「ホーム」タブ-同タブ「選択と検索」-「選択ツリー」を表示。
・選択ツリーパネルをプロパティ表示に切り替え、「ボーリング」をクリック。
※ソリッドボーリングが選択される）

Navis＋
・「分析ツール」-「属性分析ツール」より、IDを追加。
・「分析ツール」-「属性のエクスポート」より、GRMA_ID¥ID、ボーリング¥ボーリング名、ボーリング¥地質名をエクスポート。
※'ATTRIBUTE'フォルダに保存。
・Navis＋パレット内の「Navis＋」を右クリック、「属性テーブル作成」、'Boring'属性テーブルを作る。
・「Boring」を右クリック、-「属性割り当て」。属性ファイルにエクスポートしたcsvを選択し、「属性割り当て」。
※相対パスで指定。以降も同様。
・同ダイアログの「更新」も同ファイルを設定。
※後に更新用ファイルを用意したら、それを指定。
・同ダイアログの「項目テーブル作成・更新」に GEORAMA で書き出した'Borlist.csv'を選択。「ボーリング名」で「更新」。
※見失う場合、'Borlist.csv'の A列タイトルを’ボーリング名’に変更
・Navis＋パレットの「Navis＋」を右クリック、「属性テーブル作成」、'Surface'属性テーブルを作る。
・「Surface」を右クリック、-「属性割り当て」。属性ファイルにEXCELから書き出したcsvを選択。
・'Navisworksのプロパティ名称'に'Civil3D\General:Layer name'を指定し、「属性割り当て」。
・同ダイアログの「更新」も同ファイルを設定。
・保存。

この手順において付与できない項目が、地質情報の記事（概要）、深度（層厚）、孔内水位、Ｎ値、土質定数（土質試験結果）。基本、柱状図xmlと重複する項目なので、属性付与はコピペで済みます。が、該当区間のみの掘進時水位の選択（ケースとの位置関係）、2層にまたがるN値や原位置試験結果の振り分けの判断など、機械的付与の難しい属性も含まれています。
土質定数は設計段階で付与すべきでしょう。土質試験結果は多岐にわたるので、柱状図同様に試験結果一覧表PDFへのボーリング毎のリンクの方が使いやすいと考えます。
深度は柱状図を参照する際に必要になるので、付与しておくべき情報でしょう。
付与できない項目のうち、深度、色、記事のコピペは、Python で組みました（その後、Navis＋で更新）。ま、柱状図と重複する内容ですから、スクリプトを組んでまで例を再現される方は少ないかもしれません。
Navis＋ がないと外部の属性ファイルの更新を反映できないのもイマイチ。属性の埋め込みが前提なら、V-nas Clair ＋ Geo_kit の方が簡単です。実務上は後者の属性付与・表示方法で十分ですが。
ガイドラインには属性付与「例」しか掲載されていませんので、何を含めるかはお客様と話し合うことになるのでしょう（そして BIM/CIM実施計画書「4.3 (2)属性情報の項目」へ反映）。

現状、属性の付与されたボーリングモデルをソフトが地質として認識してくれません。外部参照される属性の役目を柱状図（xmlデータ）が担っている状況です。そのため、属性csvデータは各調査ステップ間で利用されないでしょう。ここだけ見ると非効率です。が、維持管理段階まで考えると、１つの統合モデルに地質の情報も含まれていた方が良いのは確か。事業の流れの中間段階はともかく、最終段階のボーリングモデルだけは、きちっと属性まで整理したいものです（しかも永続的なフォーマットで）。

ま、属性付与だけなら1時間程度の作業です。中間段階では無駄になるかもしれませんが、「手間だから実施しない」とも言えない作業です。ガイドラインが変更されるまでは、淡々と処理せざるを得ない内容でしょう。

CIM というワードをよく聞きだしてから７年が経過しました。その間、準備してこられた技術者は淡々と作業しているでしょうし、準備してこなかった技術者は慌てていることでしょう（いまだに他人事、人任せなオジサンもいますが）。

２年後の原則 CIM に向けて、周囲の重い腰がいくらか上がり始めたようです。

**************************************
20210807追記
地形のJ-LandXML書き出しです。

Civil3D から LandXML書き出し
・「ツールスペース」-「設定」タブ、一番上を右クリックし、「LandXML 設定を編集。。」を選択。
・「LandXML 設定」ダイアログで、識別要素を書き出しを「オン」、作成者等の情報を入力、「OK」をクリック。
V-nas の i-ConCIM キットでインポート・エクスポート

令和3年度 地質CIMの対応 その2

今年度のボーリングモデル作成時のソフト比較です。

Civil3D 2021 + GEORAMA でボーリング取り込み可（調査結果モデル）。
・取り込み後に推定・解釈モデルへ書き換え可能であるが、すべて手作業になり面倒。
・柱状図の情報（属性）は埋め込まれる。閲覧には上記ソフトが必要。
・属性を別ファイルで管理するため、Navisworks + Navis+を利用。
・ただし、対応できる属性に限りあり。ガイドライン掲載例の6割はこれらで対応可。
・外部ファイルの更新を反映するには上記ソフト（Navis+）が必要。

V-nas Clair　2020.6 ＋ Geo_kit でボーリング取り込み可（調査結果モデル）。
・取り込み後に推定・解釈モデルへ書き換え不可。
・柱状図の情報（属性）は埋め込まれる。閲覧には上記ソフトが必要。
・i-ConCIM - IFC ツールの属性付与は、すべて手作業になるので現実的ではない。
・dwgで保存すると旗揚げ・N値などが倒れる場合あり。⇒Navisworks の利用不可。

応用さんの製品は手元にありません。設計者が使わないので、上記ソフトに比べシェアは小さいでしょう。同理由で購入予定もナシです。

地層研の Geo-Graphia は計算のためのプリポストとして他支店で購入しています。私は現行品を使ったことはないのですが、10年以上前に販売されていたソフトの機能やサポート対応により、あまり良い印象を持っていません（文句を聞かないので現行品では改善されているのかもしれません）。
HPの紹介では IFC 出力が開発中、属性管理は謳われていませんので、まだ CIM の土俵に上がっていないようです。

Navis は初めて利用しました。特筆すべき点はありませんが、AECコレクションの一つ（とアドオン）なので、データ互換に関してはストレスなく扱えました。

現状では、Civil3D 2021 + GEORAMA ＋ Navisworks + Navis+ を選択せざるを得ない状況でした。サーフェス、ソリッドモデルの作成を考えると、汎用ソフトでは一択になりそうです。

令和3年度 地質CIMの対応

設計者から 地質の CIM 対応が含まれた仕事について相談を受けました。

概要説明、特記仕様書の提示があったので、それに沿った返答をしました。近くの地質担当部署は3次元に対応していませんので、元地質の私のところに来られたようです。

調査業務でも特記に CIM が入り始めたようで、GEORAMA の絵が欲しいと仰る地質担当者も出てきました（使用経験も絵もないのに提案して、どう対応されるのでしょう？）。
GEORAMA の絵って見栄えの良いモノはないんですよね。単なるCADデータですから。見栄えよくするには、私は他のソフトを使用していました。
それに GEORAMA 自体は CIM 要領やガイドラインに対応していません。「GEORAMA でCIM！」と言われる方は半分正解ですが、おそらくガイドラインを読まず、手を動かされていない方でしょう。（自社開発を除き）対応している汎用地質ソフトはまだないと思われます。これからです。

まず、今年度の要領、ガイドラインで、押さえておくべき点をピックアップしておきます。運用上の問題はいろいろありますが、それでも要領等に従うことは基本です。

BIM/CIM モデル等電子納品要領（案）及び同解説 令和3年3月
https://www.mlit.go.jp/tec/content/001395718.pdf

1.3 成果品の作成範囲
受発注者協議により作成する BIM/CIM モデル等を決定する。
① BIM/CIM モデル照査時チェックシート、BIM/CIM モデル作成 事前協議・引継書シート、BIM/CIM 実施計画書、BIM/CIM 実施（変更）計画書、BIM/CIM 実施報告書等
② BIM/CIM モデル：構造物や地形等の各 BIM/CIM モデル
③ 統合モデル：各 BIM/CIM モデルを統合したモデル
④ 動画等：イメージ画像や動画等のファイル
⑤ リクワイヤメント（要求事項）として特別な検討のために作成した BIM/CIM モデル
⇒②、③、④「成果物モデル」、⑤「要求事項モデル」

詳細設計
「成果物モデル」3次元モデル成果物作成要領（案）
「要求事項モデル」BIM/CIM モデルの作成方法、ファイル形式等は規定しない
詳細設計以外
「成果物モデル」BIM/CIM活用ガイドライン（案）を参考に設定
「要求事項モデル」BIM/CIM モデルの作成方法、ファイル形式等は規定しない
⇒表1-3の納品例：調査業務ではボーリングモデルの納品が必須、パネルダイアグラム・サーフェス・ソリッドは任意。地形、構造物、統合モデルは作成した場合に納品。

3. ファイル形式（電子成果品のファイル形式）
・地形モデル、線形モデル、土工形状モデルのファイル形式は J-LandXML 形式及びオリジナルファイル形式。
・構造物モデルは、IFC 形式及びオリジナルファイル形式。
・地質・土質モデル及び統合モデルは、オリジナルファイル形式。

2.2.2 GEOLOGICAL（地質・土質モデル）
地質・土質モデルは、地質ボーリング柱状図、表層地質図、地質断面図等の地質・土質調査の成果を、3次元空間に CAD データとして配置したもの。
・ボーリングモデル（調査結果モデル・推定・解釈モデル）
・準3次元地盤モデル（テクスチャモデル（準3次元地質平面図）・準3次元地質断面図）
・3次元地盤モデル（サーフェスモデル・ソリッドモデル（B-Reps、ボクセルル、柱状体））

サブフォルダ 格納される成果品
GEOLOGICAL ・地質・土質モデル（オリジナルファイル）
　VIEW ・確認用ファイル又はビューア
　SOURCE ・ボーリング柱状図やボーリング交換用データの XML ファイル等
　TEXTURE ・テクスチャファイル（TIF、JPG 等）
　ATTRIBUTE ・3次元モデルから外部参照される属性情報ファイルや参照資料ファイル（PDF、CSV 等）
　DOCUMENT ・管理情報等、地質・土質モデルに関する各種ファイル

2.2 BIM/CIM モデル
発注者が BIM/CIM モデルのデータを操作できる環境にない場合でも確認することができるよう、必要に応じて確認用ファイル又はビューアを格納すること。格納するファイル形式やビューア等の選定にあたっては、発注者と協議の上、決定すること。
発注者が 3 次元モデルを確認できる代表的な方法
・3D PDF
・イメージ画像：必要な方向や位置は、発注者と協議の上、決定する。
・3次元モデルビューア：3 次元モデルビューアは、インストール等が必要な場合があるため、利用可能か発注者と十分な協議が必要である。

＜BIM/CIM モデルのデータファイル名について＞
「BIMCIM_MODEL」フォルダに格納する BIM/CIM モデル（特に構造物モデル）のデータファイル名は、当該 BIM/CIM モデルが何を表現しているか発注者や後工程の受注者が分かるように、受発注者協議により設定することを推奨する。

BIM/CIM活用ガイドライン（案）第１編　共通編 令和3年3月
https://www.mlit.go.jp/tec/content/001395762.pdf

第3章地質・土質モデル
1 地質・土質モデルの作成・活用に関する基本的な考え方
地質・土質モデルを作成することで、本体構造物と地質・土質構成等における位置関係を立体的な把握が可能となり、各段階の地質・土質上の課題や地質・地盤リスクを関係者間で共有することにより、追加すべき補足調査や計画立案に関する検討を円滑に進めることが期待できる。
地質・土質モデルの作成・活用にあたっては、不確実性の程度やその影響について、関係者間で共有・引き継ぎを行う必要がある。なお、このような不確実性の取り扱いについては『土木事業における地質・地盤リスクマネジメントのガイドライン』が参考となる。
地質・土質モデルの品質は、モデル作成時点における地質・土質調査の質と量に依存するものであり、事業の進捗に応じて構造物等のモデル詳細度がより詳細になったとしても、それに応じて地質・土質モデルの品質を必ずしも確保できないため、構造物等で適用する「詳細度」と同様の考え方を適用することに無理があることから、地質・土質モデルに対しては「詳細度」を適用しないこととする。

1.1 地質・土質モデル作成における基本方針
モデルの品質を明確にするために、作成で用いた地質・土質調査成果やこれらに基づく推定の考え方について「BIM/CIM モデル作成 事前協議・引継書シート」へ必ず記録し、継承するものとする。

2.3 地質・土質モデルの活用時の留意事項
・ボーリングモデルのうち、調査結果モデルと推定解釈モデルのどちらを納品するかは発注者と協議するものとする。
⇒調査結果モデルを選択する。推定解釈モデルは作業量大。
・単にボーリング結果を数学的に補間する方法もひとつの方策であるが、推定精度が低下する恐れがあるので留意する。

3.2 属性情報
地質・土質モデルは、ボーリング調査結果から得られた各地層に対して、物理特性や圧縮強度等の力学特性のような土質試験結果等の様々な属性情報を扱うことが可能である。そのため、地質・土質モデルは、形状情報（オブジェクト）と属性情報で構成され、各事業段階へモデルを更新していく場合は、形状情報と属性情報を一体化するよりも、形状情報と属性情報を分離し、「共通 ID」を使用して、各々を個別に管理するのが有効である。
⇒Civil3D＋GEORAMA でオブジェクトを作成した場合、Navisworks ＋ Navis＋で属性情報を管理します。現状では一方通行の作業になる（Navis のデータをCivil に取り込んだ段階で GEORAMA は地質を認識しない）ため、ボーリングが追加になるとIDを降りなおした方が効率的です。

4.2　ボーリングモデル
表30 ボーリングモデル（調査結果モデル）の属性情報（例）
基本情報（ボーリング名、ボーリング柱状図の種類、緯度及び経度、角度・方位、孔口標高、総削孔長、孔径（ｍｍ単位）、オリジナルデータリンク、改訂履歴）
地質情報（土質・岩石名、記事（概要）、深度（層厚）、孔内水位、Ｎ値、土質定数（土質試験結果））
⇒属性情報は柱状図や土質試験結果から取り込むため、冗長的作業になります。

CIM の本質は、３次元化よりも経年的な情報管理にあると考えられます。そのためか「4.2ボーリングモデル」は設計同様に気合の入った例が示されています。
が、個人的には疑問。ボーリングが増えると、現状では属性のついていないCADデータに柱状図データを追加し、ボーリングモデルやサーフェスモデルを作り、再度属性を振りなおすことになります。前段階の調査会社と使用ソフトが変わると、全柱状図を取り込みなおしになるので、データの使いまわしの点でも属性を作りこむのはやや無駄に感じます（あくまで現状では、です）。7年前にも似たような考え方を書き残していますので、状況は何も変わっていません。https://phreeqc.blogspot.com/2014/11/cim.html

ガイドラインにおいて属性の例ではなく、何を付与すべきかが示され、データフォーマットが決まれば、ソフトが対応してくれると思われます。それによってソリッドからの属性読み込み、柱状図認識・再現機能など逆向きの流れができるようになれば、活用の場は広がります。期待を込めて、その時代を待ちましょう。

長くなったので、対応方法は後日。

InSAR と地下水位

 地下水位を InSAR で推定するという文献を見かけました。

観測孔の水位と InSAR の結果を関連付け、面で推定するという内容。意外と古くからあるんだなーと思いながらフォルダに片づけました。

そして今週、SimPEG の 1st Seminar が開催されていましたので Youtube で 後追いしました。驚きですが、内容は上記と同じ。InSAR で地下水の賦存量を推定するものでした。

異なるのは物理モデルを踏まえている点。
圧密沈下計算では水位低下によって地表変位を計算します。では、変位から逆に水位を推定することも可能でしょう。InSAR で時系列の変位を求め、クラスタリングで複数のパターンに区分し、それぞれの水位変化を推定するといった流れでした。

テクトニックな動きは分離できないのでは？と疑問に思っていると、「岩盤の動きを見たらよい」という旨の意見がありました。なるほど。

InSAR では、個人的に良好な結果を得られたことがありません。PSInSAR でないと欲しい精度が得られないこともあり寝かせていましたが、このような使い方が知られているのであれば、自由に使えるようにしておきたいですね。

発表内容は研究中のモノでしたが、インパクトは十分にありました。既に関連文献が出ていましたので、意外と早くコードが公開されるかもしれません。InSAR と共に注目しておきましょう。

Hertz’s theory of impact

先日読んだ文献に、落石衝突による振動発生（室内・フィールド実験）について書かれていました。

模型実験の最盛期を経験してこなかった私にとって、相似則は単なる知識でしかありません。特に時間（の逆数になる周波数）の相似則についてはイメージが追いつきません。
この2文献では、いわゆる相似則を用いない平均周波数の推定例が掲載されていました。模型実験に造詣が深くなくとも、周波数の当たりを付けることが可能です。

Farin, M. et al. (2015), Characterization of rockfalls from seismic signal: Insights from laboratory experiments, J. Geophys. Res. Solid Earth, 120, 7102–7137, doi:10.1002/2015JB012331.

室内実験

• 周波数に対する粘性減衰の影響を定量化するために、式25および式26の合成スペクトルに、係数 exp(−𝛾(𝜔)r)をかけている （1/γ(ω)はエネルギー減衰の特性距離を表す）。
• エネルギー Wel や周波数 fmean および Δf を計算する前に、伝播中の波の分散や粘性減衰を評価し、測定された振動をこれらの影響から補正することが重要。実験では、この補正をシステマティックに行っている。
• 実験では1MHz以上のサンプリング周波数で信号を記録。
• 球状のビーズが粗い表面に衝突したり、礫が平らな表面に衝突したりすると、等価接触半径がインパクターの半径よりも小さくなることがある。板材では接触半径Rが1.15だけ小さければ、理論上の放射弾性エネルギー Wel は2倍小さくなる。ブロックでは、有効接触半径 R が2.1倍小さければ、放射弾性エネルギー Wel は10倍小さくなる。礫が粗いブロックに衝突すると、接触半径はさらに小さくなり、放射弾性エネルギー Wel も小さくなる。
• 周波数 fmean と Δf は半径 R に反比例し、放射弾性エネルギー Wel に比べて半径の変化による影響が少ない。図9上で礫が放射する周波数が、球状のビーズのそれに近いことからもわかる。

フィールド実験

• 周波数fmean（およびΔf）については、mおよびVzに対する明確な依存性が認められない。
• 衝撃によって放出されるエネルギースペクトル全体をできるだけ多く測定する必要があり、そのためには、理想的には 3/Tc 以上の高いサンプリング周波数を使用する必要がある。例えば、今回の落石実験では、インパクタの質量が m=10kg から m=2000kg に増加すると、衝撃持続時間 Tc がそれぞれ0.01秒から0.06秒に増加するため、サンプリング周波数は少なくともそれぞれ 300Hz から 50Hz にする必要がある。

まとめ

• 平滑な板の上では、弾性波と粘弾性散逸がエネルギー損失の主なプロセスである。弾性波の放射は衝撃エネルギーの0.1％から0.3％に過ぎない。粘弾性散逸は、板厚の10％以下の直径のインパクターで主に見られる。直径がプレートの厚さよりも大きい場合、エネルギーのほとんどすべてが弾性波で放射される。
• 粗いブロック上では、弾性散逸は失われたエネルギーの0.03%から5%程度。一方、塑性変形などの他のプロセスで失われるエネルギーは、インパクターの質量が増えるにつれて、失われるエネルギーの50%から99%以上に増加する。また、粘弾性で失われるエネルギーは、ビーズの質量が大きくなるにつれて、失われるエネルギーの50％から2％に減少する。
• 粗い基板＋粗いインパクタを用いた衝撃実験でもヘルツのモデルが使用できることが実証された。
• 落石の際に失われるエネルギーのほとんどは、塑性変形、または衝撃物の並進モードや回転モードで散逸すると思われる。塑性変形や一般的に不可逆的な消散は、弾性波で放射されるエネルギーを減少させるもので、定量化は困難。しかし、衝突体の質量と速度にかかわらず、弾性波で放射されるエネルギーは、衝突エネルギーの数パーセント以下であると考えられる。衝突時に粘弾性の散逸で失われるエネルギーは、フィールド上の地震観測所で検出される質量の範囲では無視できると考えられる。
• 観測された振動から衝撃の特性を推定する際には、波動伝播において放射エネルギーの大部分が高周波で失われるという事実が大きな制約となる。落石の振動調査では、できるだけ衝撃に近い信号を記録し、少なくとも継続時間の3倍の高周波サンプリングを使用することが推奨される。

その他

• アクティブまたはパッシブ手法により、サイトのグリーン関数を評価することができる。このグリーン関数を式25、26に用いることで、正規化された周波数がレイリー波のグリーン関数を用いて計算された周波数とどの程度異なるかを推定することができる。
• 分散に加えて、伝搬中のエネルギーの粘性減衰は、特に高周波数の場合、フィールド上の測定周波数に大きな影響を与える可能性がある。Gimbertら[2014]は、河川の乱流によって発生する振幅スペクトルを調査し、音源からの距離 r が5mから600mに増加すると、その中心周波数が10分の1に減少することを示した。
  

Farin, M. et al. (2018).Link between the dynamics of granular flows and the generated seismicsignal: Insights from laboratory experiments. Journal of Geophysical Research:Earth Surface, 123, 1407–1429.https://doi.org/10.1029/2017JF004296

斜路の模型実験

• バルクの運動エネルギー Ec(t) は失われた位置エネルギー ΔEp(t) よりも約1桁低い。そのため、これらの実験で粒の流れによって失われた総エネルギー Etot(t)=ΔEp(t)+Ec(t) は、ΔEp(t) よりもわずかに高いだけである。
• Hibert, Ekström, et al. (2017) では、エンベロープ Env(t) の最大振幅と流れの最大運動量｜MV｜の間に定量的な相関関係があることを報告している。しかし、実験では Env(t) の最大振幅が X方向（斜面流下方向）の速度の最大値には対応せず、Z方向（垂直方向）の速度の最大値に対応することが観察された。
• 実験で得られた粒状流の周波数は，自然現象で記録されたものよりも100～1,000倍高い。この違いを理論的に説明するために、実験室と野外でヘルツのモデル（Hertz, 1882）で予測される衝撃の持続時間の逆数として定義される周波数を計算した。実験室で記録された振動の周波数を野外で観測されたものにアップスケールするには、ヘルツの衝突継続時間が適切な特性時間であるといえる。
• According to Hertz’s theory of elastic impact,　this duration can be written
  Tc ≃ 2.87 (m^2/（RE∗^2VZ）)^(1/5)
  where m, R, and VZ are the impactor’s mass, radius and impact speed in the Z direction, respectively, and E∗ is　an elastic modulus given by
  1/E∗ = (1−𝜈^2)/Ei+(1−𝜈^2)/Eg
  where 𝜈i and 𝜈g are Poisson’s ratios, and Ei and Eg　are Young’s moduli for the impactor and the ground, respectively (Johnson, 1985).
• This approach, however, can be imprecise because the frequencies of the seismic signal emitted by an impact also strongly depend on the roughness of the ground, which is not taken into account in Hertz’s expression of the impact time Tc. In the field, the measured mean frequency fmean should also decrease as the distance between the seismic source and the seismic station increases because high-frequency energy attenuates more rapidly than low-frequency energy in heterogeneous media (Aki & Richards, 1980).
• 地震効率（失われた位置エネルギーが放射された地震エネルギーに変換される割合）は，実験では傾斜角が0°から20°に増加するにつれて半分に減少した（0.033％から0.017％）。斜面の角度に依存するのは、斜面の角度が小さいほど粒子の衝突がより正常に斜面に向けられ、地震効率が高くなることによると考えられる。
• 粗い面に衝突した粒子の地震効率は、平滑面に衝突した場合よりも約10倍低く、状態（平滑，粗面，侵食可能）に大きく依存すると推論した。

E*の式は DEMでも見かけますね。ヘルツの接触応力として有名だったのは知りませんでした。Hzのイメージしかなかったもので。

電気探査と測線数

SiｍPEG で数値実験。

確認したのは、電気探査の2極法（pole-pole）と4極法（dipole-dipole）x 縦断（1測線）と縦横断（3測線）の組み合わせ。Tutorial の地形データ、モデルを利用しました。

まずは3次元地形のもとで測定データを作成します（forward modeling）。その後、各条件でinversion してみました。

縦断方向の導電率モデル。異常部は球になっています。

地形と測線位置

まずは、pole-pole。多少、差があるようです。
１測線

 ３測線

次は、dipole-dipole。 こちらは差がわかりやすい。
１測線

3測線

4極法の3測線が BEST なのは予想通りでした。
時間がかかるので2極法を選択した場合は、横断もいっしょに3次元で inversion をかける方が良い、という程度ですね。案外、いくつかのスパースモデリング結果をスタックした方が、実務上は判定しやすいかもしれません。 

*****************************************
20210714追記

pole-pole でノルムを変更してみましたが、ほぼ変化なし=スタックしても変化なしでした。

あと、seismic に3次元がないことに今頃気づきました。joint inversion (2D_seismic+) は pyGIMLi、sparse(3D_DC) は SiｍPEG のように使い分けましょう。

oneAPI

Intel Parallel Studio XE 2020 Composer Edition for Fortran のSSR更新に関する案内が届きました。

Parallel Studio が廃止され、後継の oneAPI に移行する案内は以前から届いていましたが、まだ前者を使い続けていました。で、更新のタイミングで入れ替えることに。
oneAPI へのライセンスの書き換えはボタンクリックのみ。無償アップグレードでした。

その後、SSR更新費用を確認すると、10万で約倍になっていました。VTune や Adviser が含まれた上位版の後継ですので妥当なのですが、毎年続くと思うと痛い。残念なことにアップする前だと今回に限り6万だったようです。うーん。ま、oneAPI 自体は無償なので、保守契約はやめましょうか。
https://software.intel.com/content/www/us/en/develop/articles/free-intel-software-developer-tools.html

この oneAPI の無償化ですが、大きな転機ですね。NVIDIA さんが Linux 向けに PGI の後継を無償で配布していますので、Intel さんも対抗されたのでしょうか。ユーザーにとってはありがたいことです。今後、Winのデファクトスタンダードになるでしょうね。

で、早速 Wiｎ 版をインストールしてみました。oneAPI 側の制限はなくなりましたが、Visual Studio 側の制限はそのまま。個人利用として継続します。
途中、C＋＋がはいっていませんと注意されます。C＋＋をいれて、Base、HPC の順でインストール。時間がかかります。更新がさらに面倒になりそうです。Linux 版だと apt や Docker で簡単に更新できそうですけど。

インストールが終わり、コンパイル。
exe を動かしてみると、今度は'dllがない'と叱られます。oneAPI からは インストールマシンでも path を通してくれないのでしょうか。うーん、面倒。
libiomp5md.dll cannot be found after upgrading to oneAPI for Windows
https://community.intel.com/t5/Intel-Fortran-Compiler/libiomp5md-dll-cannot-be-found-after-upgrading-to-oneAPI-for/td-p/1258331?profile.language=ja
...\OneAPI\compiler\latest\windows\redist\intel64_win\compiler

エラーが出るたびにdllを検索。都度、Path を追加すれば最終的に動きました。

この oneAPI ですが、OpenMP で Intel の GPU をサポートしているそうです。が、詳細な情報を拾えず。雑誌などで紹介されているかもしれませんので、ゆっくり調べてみましょう。

ま、ゴリゴリ使っているプロでもないので、当面はこの状態で良いかな。