やり残し事項 2021

今年は人生で1，2番目の悪い年でした。

年末のダメ押しが決定的で、再起動まで少し時間がかかりました。

仕事のストレスも相まって、余暇にPCを触る気力も失せてしまい、どうしようかと思案。
胃痛のため食欲に意識を振ることができず、長く寝られるわけでもない。
結果、PCを閉じ、何もしないことにしました。

といっても本当に何もしないわけではなく、大掃除をしながら体を動かし、静かな曲を聴いたり小説や漫画を読んだりしました（技術者でなければ、情報収集や勉強をしないというのも普通の生活なのかもしれません）。

救いは思いがけないところにあるもので、 漫画の主人公が人の死に思い悩んで立ち止まりかけた時に他のキャラが吐いたセリフ。これにハッとさせられました（我ながら安いものです）。

それを契機とし、PCを開き計算結果を確認したり、新たな計算を始められるまでには回復。胃痛は治らないものの、心は明確に切り替わりました。まだ何かやりたいことがあるわけではありませんが、手を止めないということは決定事項。後悔や反省は忘れることなく、考えながら2022年を進もうと思います。

やり残し事項です。

優先度中：流体＋個体（不連続体＋連続体）＋振動
優先度低：PSInSAR
優先度低：Dtransu の GPU 対応
優先度低：地表流＋地下水＋移流拡散

中期目標5年目は失敗。ま、成功しても結果的には意味がない状況に陥っていました。こちらは当分目標ナシです。しばらく考えましょう。

個別要素法の接触力モデル

Catherine O'Sullivan 原著「粒子個別要素法」3章より、接触力を表現するモデルです。
これまで見たり接してきたりしたモデルが、どのような分類に属するのかを理解できました。この図書、良いかも。


法線方向の接触モデル
・線形弾性ばね
・簡易ヘルツ接触モデル
　材料特性（G, ν）の導入
・降伏を考慮した法線接触モデル
　・ウォルトン―ブラウン線形モデル（半ラチェットばねモデル、履歴減衰モデル）
　・ばね、ダッシュポット（線形ケルビンモデル）：塑性変形によるエネルギー散逸を表現

接線力計算
・ミンドリン―デレシビッツ接触力モデル
　・簡略版：ブクオックモデル
　・ソーントン―イェンモデル（荷重履歴依存性の導入）：非線形、再載荷の考慮

引張り力伝達モデル
・並列結合モデル

転がり抵抗
・岩下―小田の回転抵抗モデル
・ジアンらの回転抵抗モデル

時間依存性
レオロジー：ばねとダッシュポッの組み合わせで表現
・バーガーモデル（ケルビンとマックスウェルの直列）：クリープの表現
・マックスウェルモデル（直列）：経時変化による劣化

 

個別要素法での減衰

個別要素法の減衰の扱いを整理。テキストは以下の２つ。
Catherine O'Sullivan 原著「粒子個別要素法」
古川ほか「個別要素法を用いた構造物の動的解析における減衰のモデル化に関する基礎的検討」
https://www.jstage.jst.go.jp/article/jscejseee/70/4/70_I_89/_pdf

F=Mx''+Cx'+Kx

1.要素間に作用する減衰力
F=Cδ'+Kδ
1-1.ダッシュポット(衝突エネルギーを消散)
・反発係数e←→減衰定数h←→ダッシュポットの粘性減衰係数C（衝突する粒子の質量・ばねに依存）

1-2.剛性比例型減衰：C=βK=(2h/ω)K
・速度項に乗じて力になる→単位：力/速度
・β固定（周波数ωのh固定）→周波数が高いほど減衰（定数）大
・単独で主要なモードの減衰を再現でき，剛性低下に応じた減衰力の低下も表現し，かつ高振動数領域の減衰性能も確保できる
・計算ステップを細かくする必要あり

2.要素重心に作用する減衰力
F=Ma+Cv
2-1.質量比例減衰(主要なモードの減衰)：C=αM=2hω
・α固定（周波数ωのh固定）→周波数が高いほど減衰（定数）小→１-1の併用で高周波数帯をカバー
・相対的に速度の大きな領域で誤った物体力→破壊モードに影響

2-2.局所比例減衰力(主要なモードの減衰)：-α|Fp|sign(vp)
・加速させようとする運動のみを減衰→誤った減衰力が定常運動から生じない
・速度に依存しないので、履歴減衰と類似
・無次元のため、周波数に依存しない
・係数αを調整することで主要な固有モードの減衰定数を表現→減衰定数は高々数％→１-1の併用で高周波数帯をカバー

コロナ対策のバランス

忘年会があったようです。

私の周りの方々は参加されませんでしたが、他のフロアの方々を中心にバス2台でホテルへ向かわれていました。

勤め先の支店では、緊急事態宣言中も、その前後のテレワークが推奨されている期間も、基本、テレワークなしの出社体制を貫かれてきました。この支店のみのようなので理由は不明ですが、おそらく地域的な業務特性？を優先された結果だと想像します。

と思いきや、今月、昼食のための会議室利用が禁止さました。これまで使われていた方々は文句を言われていましたが、これは理解できます。が、今頃？という意見にも同感。
リモートはダメで全員出社を強いるのに、昼の少人数の集合はダメ、忘年会のような大人数会食はOKで、少人数会食はダメというように見えてしまう方針は、批判されても仕方ないでしょう。

なんだか、バランスが悪いですね。

メンテナンス その２

 引き続き、週末は計算を回しながら細々したメンテナンスを継続。

・自転車のブレーキ交換
　シュー、ワイヤー交換。快調です。

・自動車のタイヤ交換
　冬タイヤへ。ハブの錆が少し進んでいたので清掃。

・自転車チューブ交換
　パンクしたので開けてみると、異なるサイズのチューブが入っていました。即交換。

・大型ごみ搬出
　ヒーター、ストーブ、スキー、スノボ等々。

・図書整理
　配置換えとともに、不要なものを資源ごみへ。

・室外機のホース交換
　配管テープが破れていたので新品に巻きなおし。ついでに割れていたホース交換。

・Switch Joy-Con ステック交換
　代用品の DUALSHOCK3 がダメになり、壊れていた純正品を修理。パーツ交換で復活。

あとはバッテリーが気になるくらい。まだもつかな？

今年は仕事が暇で、土日を完全に休めています。このような年末は初めて。おかげで今までサボっていた身の回りの整理がはかどっています。晴耕雨読ではありませんが、日中動いて、夜文献整理。既に大掃除にも半分終わり。年末年始はゆっくりできそうです。

 

個別要素法でのパッキング

初めて3次元で DEM を扱っています。

理論は連続体に比べて単純。小さなモデルは特に問題なく流せました。
が、数万～数十万の粒子を扱いだしてから、パッキングで躓きました。結果が初期状態に大きく依存するのです。パラメータよりも、影響が大きいのかもしれません。

パッキングについて具体的なノウハウが書かれている図書、文献は見当たりません。手元の資料で関連する記述があったのは以下の通り。

伯野元彦「破壊のシミュレーション」よりDEM のパッキング方法です。
a） 規則的に配置する方法
b) ランダムに配置する方法
1. 棄却法
2. 局所移動法
3. 割込法
4. 落下法
5. 成長法

Catherine O'Sullivan 原著「粒子個別要素法」より
1. ランダム生成（半径拡大法、段階圧縮法）
2. 構築法
3. 三角分割
4. 重力堆積法

高橋ほか「個別要素法によるアスファルト混合物解析の要素パッキングに関する一検討」
https://www.jstage.jst.go.jp/article/journalpe1996/2/0/2_0_229/_article/-char/ja/
規則的にパッキングする場合、配置条件によって変形に方向性が現われやすい。
規則的な配置では、単一粒子径が望ましい。
粒子分割の粗密、すなわち粒子半径も解析結果に影響する。その程度は要素モデルによって異なる。
必要とする精度を満足する範囲で最も大きい粒子径を選定するのがよい。
要素モデルの組成(バネ要素とダッシュポット要素の組み合わせ)と、その物理定数のみではなく、パッキング方法も含めて個別要素の集合体モデルを構築しなければならない。

 

パッキング状態によって結果が大きく異なる事象は、DEM 界では良く知られているようです。この、「何とでもなる」という自由度が辛い。地盤を DEM で適切に模擬するのは、連続体で初期応力状態をつくるよりもかなり難しいと感じます。選択するパッキング方法、その組み合わせ、ノウハウは技術者の経験で培われているといえるでしょう。
が、逆に以下のような考え方もあるようです。

地盤工学のための個別要素法_7_パラメータの設定と土と地盤の作り方(その2)
「その作業ができるところがDEM解析の強みでもある。解析者はその地盤作成の深淵を感じながらも、速く、効率的に解析対象となる供試体や地盤をつくる必要がある。」

経験を積んで深淵を俯瞰できるようになりたいものです。

**************************************
20211217
粒子個別要素法の記述を追記。

PSInSAR ＆ Landslide その2

こちらは地すべり分布図と PSInSAR の組み合わせ。
前者はDEM＋RF。その精度を高めるのに後者を使用。

Andrea Ciampalini et al. (2016)
Landslide susceptibility map refinement using PSInSAR data
https://www.sciencedirect.com/science/article/pii/S0034425716302759

Landslide Susceptibility Map
The landslide susceptibility map (LSM) of Messina Province was produced by implementing a Random Forests (RF) algorithm in Matlab (Catani et al., 2013).
To avoid subjectivity in the choice of explanatory variables, several parameters were considered among the DEM-derived products: Aspect, Planar Curvature, Profile Curvature, Curvature s.s., Elevation, Flow Accumulation, Topographic Wetness Index (TWI), Log Flow Accumulation and Slope.
The choice of parameters depends on the map unit resolution (MUR) (Catani et al., 2013). Among the possible parameters, we used those suggested for a map unit resolution of 100 m.
Using this MUR, the expected Area Under the Curve (AUC) value is 0.81. A training set was created by randomly selecting 10% of the landslide database.
The DEM had a 20 by 20 m spatial resolution. In the LSM, the average value inside a 100 by 100 m cell was calculated.
Each pixel was classified using four susceptibility classes:
(i) low to null (0–0.3)
(ii) moderate (0.3–0.55)
(iii) high (0.55–0.75)
(iv) very high (0.75–1)

PSInSAR
SqueeSAR: TRE ALTAMIRA, detection of millimetre surface displacements, improving previous PSInSAR algorithm.
Xband COSMO-SkyMed (CSK) in Stripmap mode (40 × 40 km in the range and azimuth directions and 3 × 3 m ground resolution) in both ascending and descending geometry

Integration
The VSlope intervals were determined based on the standard deviation (σ = 7 mm/year) of VSlope for the whole PSI data set after merging the ascending and descending data sets. The matrix in Table 3 shows the correction factors for each considered case. The velocity intervals were determined to increase the susceptibility degree from level 1 to 4. For example, an area characterized by a susceptibility degree of 1 will increase by 1 degree if its VSlope falls within 1 σ to 2 σ of the stability threshold (7–14 in this case), 2 degrees if the VSlope is between 2 σ and 3 σ and so on. The higher the ground deformation velocity, the higher the susceptibility degree. A cell characterized by the maximum susceptibility degree cannot increase further.

Discussion
To improve the PS/DS detection, L-band SAR sensors can be used because their wavelength (30–15 cm) is more suitable to investigate densely vegetated areas.
Both the landslide susceptibility map and the ground deformation velocity map can be used in different ways to forecast landslide occurrence in a selected area. Combining them may improve the reliability associated with predicting these types of phenomena, particularly for slow moving landslides

PSInSAR ＆ Landslide

PSInSAR の活用のお話。
技術者が時間を割くのは PSInSAR 自体ではなく、その結果をどう問題に応用するかというところ。活用例として見本になるような文献でしたので要点を( ..)φメモメモ。

Kamila Pawluszek-Filipiak et al. (2021)
Multi-temporal landslide activity investigation by spaceborne SAR interferometry: The case study of the Polish Carpathians
https://www.sciencedirect.com/science/article/pii/S2352938521001658

However, due to the high temporal decorrelation resulting from vegetative cover and short wavelengths (X-band), it was only partially successful (Perski et al., 2009; Perski et al. (2011). This is mostly related to the low PS density due to the application of TerraSAR-X data. Therefore, the exploitation of C-band (Sentinel-1) and L-band (ALOS PALSAR) data can present more advantages, especially in rural areas.

The PSInSAR approach is based on the following steps:
(1)interferogram generation with respect to a common master image
(2) PS candidate selection based on amplitude dispersion index
(3) first estimation of atmospheric phase screen (atmospheric influence) and topographical and displacement components
(4) second estimation of interferometric components and final PS points selection
For a precise description of the algorithms applied within each step, see Ferretti et al. (2000; 2001).

The interferograms are then generated for each slave image always using the same master image.
Amplitude dispersion index of 0.25 was utilized as a threshold for PS candidate selection, and the final PSs were extracted with a coherence value higher than 0.75.

SARscape® software was used for SAR data processing (Sahraoui et al., 2006).

Persistent scattering (PS) postprocessing
R index: > 0.33
cosβ: >0.3
Moreover, we discarded PS points which showed Vslope>0 because positives values represent uphill movements and it is not representative for small landslide movements, even though positive values exist within landslides, especially in the toe area.

three categories based on PSI velocities :negligible, extremely slow, and very slow

It is worth emphasizing that the reality of the PSInSAR velocity estimation itself is not assessed due to the lack of external field measurements, but this measurement may be related to landslide activity. For this reason, field verification was performed.

In general, applying the PSI approach in rural areas is challenging. However, the results of this study prove that increasing the temporal sampling rate in view of launching the second Sentinel 1 satellite, provides a higher PS density and, therefore, this technique can deliver useful information for landslide activity assessment.
From another point of view, the application of two Sentinel datasets from ascending and descending orbits, using only one Sentinel satellite (C band), allows us to achieve a similar PS density in comparison to the one ascending dataset of ALOS PALSAR (L-band).

SAR の利用

勤め先でも、ようやくSAR の利用を考え 始めたようです。

といっても、担当はまだ本気ではないようで、
「使うための問題点や課題をこれから抽出します」
「社員の現状はこうでした」
といった報告をされていました。

SAR に関する私の要望はPSI一択です。実務なら数百万のソフトを買っていただく必要があります。手法の区別がつかない executives に承認権があるため、技術者も含めたそれらの方々の変化を待っている状態です。が、まだ時間がかかりそう。残念。

土木分野で見かける手法は以下の3種でしょうか。

強度画像（振幅のみ利用）
pros：RGB画像の合成のみなので、技術力や知識を必要としない。簡単。
cons：変位量が得られない。

DInSAR(differential interferometry SAR)
pros：広域の変位情報が得られる。
cons：土木構造物に対しては空間解像度が粗い。

MTInSAR (multi-temporal interferometric SAR)
SBAS (small baseline subset)
pros：広域の変位情報が得られる。変位の時間変化を把握できる。
cons：土木構造物に対しては空間解像度が粗い。
PSInSAR (Persistent Scatterer Interferometry)
pros：ピンポイントの変位情報を高精度で得られる。変位の時間変化を把握できる。
cons：欲しい位置の情報が必ずしもえられるとは限らない。

その他、TRE ALTAMIRA 社の SqueeSAR の利用を見かけます。精度が良いと謳われていますが、内容を知りません。パスコさんはじめ、国内でも提携する会社があるようですね。今後、報告が増えてくるでしょう。 楽しみです。

GMT

この年になると、新たなソフトの使い方を覚えるのが億劫になってきます。
できるだけシンプルに良いツールのみを使い続けられると楽なのですが、時間と人は止まってくれません。

今回は GMT。
OpenSWPCの結果が NetCDF で提供されており、図化に何を使用するか探っていました。

ArcGIS pro 2.8 から NetCDF を Basic でも扱えるようになりましたので、まずはこれで。
はい、読めませんでした。XYだと読めましたが、緯度経度ではダメ。後で中身を確認すると、後者の並びがダメでした。

で、Python で geotiff に変換。
これは Arc で読めましたが、全時間ステップを読み込んで動画を作成するのは手間でした。

Arc に持ち込まず、Folium 上にてプロットしようかと考えましたが、数が多いので却下。
https://phreeqc.blogspot.com/2020/02/folium.html

最終的には、Forium のかわりに PyGMT を使用。
地図、観測局、波形（NetCDF）、テキスト等を重ね合わせて図化。それを動画に変換して保存。計算後にスクリプトを動かすと動画が作成されるようにしました。
PyGMT は初めてだったのですが、Python っぽい書き方のためか違和感なく作業ができました。

今回は Ubuntu で環境を構築しましたが、Win 版も配布されていました。一昔前だと、Cygwin が必要だったように思います。ユーザーは増えているのかな？

得られた結果はイマイチ。
カラーマップを調整して何度かトライしましたがダメ。
OpenSWPC のソースを見ると、RGB の前2つに P,S が割り当てられていました。うーん、最初に見ておけば良かった。