やり残し事項 2022

昨年の後悔や反省を引きづりつつも、何とか1年が無事に終わりました。
達成感はありません。が、低空飛行ながら失敗は回避できた年だったように思います。
中期目標はまだありません。その日暮らしから早く脱却せねば。

やり残し事項はそのまま+1。もう消して良いかなとも思えますが、もう一年残しておきましょう。

優先度中：MEMSセンサー＋マイコン
優先度中：流体＋個体（不連続体＋連続体）＋振動
優先度低：PSInSAR
優先度低：Dtransu の GPU 対応
優先度低：地表流＋地下水＋移流拡散

地震時の岩盤斜面の安定性

地震時の斜面安定性に関する数値解析を扱った文献３つです。
要約と結論以外は飛ばし読み、機械翻訳ベースです。

いずれも中国の方、Itasca社の製品を使われています。流行りですかね？**************************************

Study on Dynamic Response of Rock Slope With Inverse Non-Persistent Joints Under Earthquake

受け盤斜面の地震応答に対し、節理の傾斜角、節理間隔、節理が途切れている距離が与える影響を数値実験で検討した。使用ソフトはUDEC（2D）。

検討の結果、節理が斜面の動的応答に大きな影響を与えることが示された。節理の傾斜角が大きくなるか節理間隔が減少すると、斜面上のPGA増幅係数が増加。影響範囲は主に斜面の中央から肩に集中。
水平方向では、節理の傾斜角の増加と節理が途切れている距離&節理間隔の減少に伴ってPGA増幅係数が増加。
垂直方向では、PGA増幅係数曲線は最初に増加し、次に減少し、その後、標高の変化とともに増加。
応答は節理が密に発達した斜面で大きくなる。

Dynamic Analysis of the Seismo-Dynamic Response of Anti-Dip Bedding Rock Slopes Using a Three-Dimensional Discrete-Element Method

受け盤斜面の地震時おける動的応答解析。使用ソフトは3DEC。

【増幅率】
一般的に2次元解析では、標高が増加するにつれて増幅率は単調に増加する。しかし、3次元動的解析では、その「標高効果」を示さなかった（最大値は斜面の中央部と上部に表示された）。
同様に、斜面の表面近くで大きな増幅率を示す「表面効果」は、3次元条件下では現れなかった。
全体として、プロファイル1、2、および3に沿った3方向加速度は単調に増加しない。代わりに、最初に値を増やす(減少させる)パターンがあり、その後に値を減少(増加)させるパターンが見られた。

【周波数】
斜面上部ほど低周波は増幅され、高周波は低減される。支配的な周波数は減少する傾向がある。

【破壊】
斜面の下部で斜面が崩壊し始め、上部に進行する。破壊面は階段状であり、トップリングの特性を示している。斜面と節理の角度が大きくなると、破壊面の深さが大きくなり、斜面の不安定部分の範囲が広くなり、変位が増加。
Newmark 法は、一方向の地震波の作用のみを考慮しており(つまり、この方法は平面問題にのみ役立つ）、動的解析の変位結果よりもはるかに小さくなる。

Numerical analysis of evaluation methods and influencing factors for dynamic stability of bedding rock slope

流れ盤斜面の安定性評価方法と影響要因について検討。FLAC3Dを使用。

【安定性】
水平震度を加算する擬似静的法で計算された安全率は、他の解析方法で得られた値よりも大。
斜面高さの増加と地層の厚さの減少に伴って斜面の安定性が低下。ただし、これら2つのパラメータの影響には臨界値がある。本論文の斜面の例では、斜面の高さが200 mに達し、地層の厚さが10 mに達すると、斜面の安定性の低下は明らかではなくなる。
地層傾斜は、堆積性斜面の動的安定性に直接的かつ重要な影響を及ぼす。理論的および数値計算結果によると、斜面の安定性は斜面角度の増加とともに継続的に低下する。ただし、降下速度は、斜面角度の増加とともに大幅に低下。この論文のモデル勾配では、角度が55°に達すると斜面安定性の低下は明瞭でなくなる。

【破壊】
主に斜面の堆積面に沿って破壊が発生。

【周波数】
入力周波数が斜面の固有振動数に近いと、動的安定性が悪化。
安全率は1Hzから3Hzへの周波数の増加とともに増加。増加率は比較的遅い。3Hz〜10Hzの範囲内では安全率は直線的に増加し、10Hz〜15Hz以内で上昇傾向は徐々に緩やかになる。斜面の動的安定性は、動的負荷の中低周波部分の影響を大きく受けるが、高周波部分の影響をあまり受けない。

PINN examples

NVIDIA Modulus の examples を読みながら動かしました。

以前、エラーを吐く example があったのですが、他にも動かないモノが混ざっていました。設定を変更すると動きますが、収束しないとかも。ユーザー側の GPU が様々ですので、細かい設定は環境に合わせて変えないといけないのかな？となると、GPU毎に答えが変わっているってこと？ありそうですね。

1Dはすぐに収束するのですが、2D や 3D は打ち切り回数がない場合に延々と回り続ける場合も。NVIDIA さんでも難しいのは変わらないようです。
書くべき内容は簡易なコードと大きく変わりません。込み入った部分は Modulus 側で対応してくれています。

困ったのは post。3次元の点に色を付けることはできますが、形状は困難。粒子法のように表面を抽出するアルゴリズムの追加が必要でしょうか？NVIDIAさん、どうやっているのでしょう。

まだまだ時間がかかりそうです。

機械学習コンペ

機械学習のコンペサイトは Kaggle が有名ですが、日本にもあります。

SIGNATE https://signate.jp/

土木がらみのコンペも開かれています。
昔は衛星データから崩壊個所を予測する内容や、海岸線を抽出するものが出ていました。最近では、広島県から河川水位と雨量から翌日の水位を予測するコンペ、山口県からドライブレコーダーの画像から要補修個所を拾うコンペが出ています。

仕事では、機械学習の性能を発揮できないようなダーティーなデータを扱うことが多いため、たまには綺麗なデータを扱ってみたいと思うことがあります。が、経験上、これらのコンペに出ているデータもそれほど綺麗なものではありません。ストレス解消にはならないのです。

が、発見がないというわけではありません。
山口県のコンペが成功すれば、今後、安価な代替え手法として流行るかもしれません。結果には注目です。

PINNや物理探査など、機械学習が様々な分野で身近になってきた１年だったように感じます。
来年はどこまで迫ってくるでしょうか。

PINN 2D は難しい

PINN 1D から 2D へ展開したのですが、うまくいきませんでした。

簡単なコードだからダメなのか、NNが難しいのか。

分かったこともあります。
商用ソフトでは境界条件を張る場所を指定し、経時変化を入力すれば自動で処理してくれるのですが、簡単なコードでは点数×時間分、境界条件を定めなくてはなりません。当たり前ですが、慣れ切っていましたね。同様に、初期条件はイメージ通りとしても、途中の観測点も時間分作らないといけません。良い復習になりました。

正しく解けたかどうかは分かり難いですね。loss が大きくてもそれなりの答えが出てしまうところが難点。十分小さな loss を目安にするしかないのでしょう。そのために、ニューロンを増やす、隠れ層を増やすなどの対応もあるのでしょうが、いつものように最適な構造を探すのは難しい。時間がかかりそうです。

それでも、試行錯誤したことで PINN の仕組みはわかりました。当初目的は達成したので、次のステップへ移りましょう。
NVIDIA さんの Example は動かなかったのですが、なんとかなると思います。いきなり分かり難くなるのですが、中間が見当たりません。当面、これを動かすことを目標に進めましょう。

PINN で 1D 移流拡散

続きです。

拡散項を追加。自動微分がありがたい。

良さそうです。

が、ここに至るまでにおかしな結果を何度か見ました。

拡散係数が小さいと、移流しなくなる答えを出してきます。
念の為、昨日の移流のみの計算を繰り返したところ、全く移流しないパターンが出てきました。NN的には loss が最小のものを出せば良いだけなので、それはそれで答えの一つなのかもしれません。が、怖い。かなり。

PINN で 1D 移流

GitHubから、 最も理解しやすかったコードをDL。
https://github.com/nanditadoloi/PINN

PDE を移流方程式に変えて計算してみました。

移流が綺麗に解けると伺っていましたが、本当に驚き。個々の値のブレが若干大きな気もしますが、全体形状を見るとオーバー、アンダーシュートが見られません。NNなので計算時間はかかりましたが、アリですね。

Pytorch ですので、一度計算をかけた後に観測点を増やして先の続きから計算、というのが手軽にできます。FEMでは計算途中でメッシュを切りなおして、なんていうと構えてしましますが、こちらはほぼ制限なし。メッシュレスの良い点でしょう。

数値解析の解法に関する知識がほぼ不要というのも受けそうです。機械学習をやっている方なら、問題なく扱え、すぐに計算結果を出すことができるでしょう。教師データが不要なのも入りやすいと思います。

残念ながら、解けない場合にどのように改善していくかというノウハウがありません。解けた場合でも、厳密解がない場合は、正しく解けているのかどうかわかりません。これから知見が蓄積されるのでしょうね。いえ、もう数年たっていますので、どこかにあるのかもしれません。

これ、面白いので、もう少し触ってみましょう。

Failed to Install linux-image

Ubuntu22.04 の１台で、apt upgrade が使えなくなって２ヶ月。
Teams を入れようとして支障が出たため、重い腰を上げることにしました。

「大量のエラーが発生したため、処理が停止しました。」

この場合の対応はweb上に色々と書かれています。手当り次第試しましたが駄目。そういえば、手がなくなって寝かせていたのでした。

最後にた取りついたのがコチラ。
you need to remove the post-installation script of the packages
https://askubuntu.com/questions/1371963/ubuntu-20-04-lts-failed-to-install-linux-image-5-4-0-89-generic

次は１発で済ませます。

 

PINN の動向

先日、PINN の話を聞く機会がありました。

初歩的な話でしたので計算に役に立つ、というようなレベルではなかったのですが、質問する機会は得られました。

PINNで気持ち悪いのが、時空間の関係があいまいになっている点。
FEMなら全体剛性マトリックス、SPHなら粒子間距離を使って近傍の関係性を明示できるのですが、PINN は NN の weight に押し込められています。時間方向も同様です。
これについて、時間方向は1タイムステップ毎に収束させてから次に進む研究があるとの回答を頂きました。ある時点で特異値を得ても全体で均されるとわからなくなるという状況を解決したいというところから始まっているようです。まだ論文は見つけていませんが、時間方向だけでも明確になれば空間方向は NN の構成と weight で、というのはありかなと思えてきます。

もう一つはGPU。
陰解法では倍精度主体の GPU が欲しいものの、機械学習や陽解法では単精度主体の GPU で良いと思っています。では、PINNはどうか？個人的には求める制度にもよりますが、基本は単精度で良さそうだと考えていました。
この点を伺うと、単精度で書かれたコードは多くみられるが、倍精度でないと結果が合わないケースもあるとのこと。ま、そうなんでしょうね。

今後はサロゲートモデルも選択肢の一つに入る時代になるでしょう。既存の手法を深めつつ、新しい手法も身に着ける必要があります。ま、新技術に対し真摯に取り組む必要があること自体は、時代が変わっても不変です。

Georisk

PINN の論文をたどりながら、思いかけず出会った雑誌。

Georisk: Assessment and Management of Risk for Engineered Systems and Geohazards | Taylor & Francis Online (tandfonline.com)

ジオリスクのみで雑誌ができるとは思ってもいませんでした。
興味を惹かれたので、Open になっている中から、一つ読んでみました。
Observational method as risk management tool: the Hvalfjörður tunnel project, Iceland 

フォールトツリー分析を用いたトンネル掘削時のリスクマネジメント例です。詳細までは書かれていませんが、スタンダードな手法でリスクマネジメントを行われたことは理解できます。

これを読んでいて思い出したのが、総監。「あなたは管理技術の何を使いましたか？」と聞かれたのを思い出しました。FTAのような技術を選択肢として複数持っていないと、管理は低レベルになり事故を引き起こしかねません。裸でリスク源に向かっていくようなものです（先のexecutivesです）。

適切な技術を選択し、観測データによるリスク対応をトンネル掘削前に決めておく。これが標準なのでしょう。管理技術を学んだ方々にとっては特に難しい話ではないのですが、訓練は必要でしょうね。

クラスター発生

職場のフロアーでコロナのクラスターが発生。

最終的には2週間で８人が陽性となりました。

私は2人目。3人目が出た時点でExecutivesに原因（ヒトの過密）と対策（在宅勤務の導入）を授けましたが、正常化バイアスが勝ったのかプライドが勝ったのか動けず、被害が拡大。1週間後、８人目が出た時点でようやく在宅に切り替わりました。遅い判断です。ご家族も入れるとまだ被害人数は多くなるでしょう。

私から見れば、リスクマネージメントも危機管理もできなかった、ただの人災です。が、当然のように責任を求めず、責任を取らせない組織なので人災という意識は芽生えないでしょう。このまま時間が過ぎると在宅勤務がなくなり、皆集まれ、出社しろ状態に戻ると思います。総監をお持ちでないので、リスクマネージメントを勉強されたことがないのかもしれません。無知は罪です。

基本的に管理能力がないなら Slack や Teams で情報を撒いて判断を個人に任せればよいのですが、それもできない。今まで通り情報は囲い込み、判断しない。それで何とかなれ～っと期待する。神頼みはマネージメントではありません。

ま、波風を立てない、特別なことや新たな改革に縁の遠い人ほど出世している会社です。有事に何もできないのは、仕方がないのかもしれません。

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

20221206追記

最終的に10人のクラスターとなりました。そして、「皆集まれ」で発生前の状態に戻りました。リスク源は保有らしい。
先日、同じフロア―に陽性者が出た場合に限り在宅勤務が認められるとお達しがありました。うーん。

PINN

先日、機械学習を利用した物理シミュレーションの話を伺いました。

近年多い話題なのですが、こちらは一味違いました。
PINN: Physics-Informed Neural Network と呼ばれているそうです。損失関数に支配方程式を組み込み、残差を０に近づけるという手法です。なので、一般的な機械学習のように学習データを大量に準備する必要がありません。
理解すると「そんな単純な話」と思ってしまいますが、まったく思いつきませんでした。凡人には思いつけない、素晴らしい発想ですね。

NVIDIAさんは2年前にソルバーを公開されています。SimNet と呼ばれています。今日、このUser Guide の1～３章までを読み、４章以降をざっと眺めました。理論部分を完全には理解できなかったのですが、使う部分に限っては難しい内容ではありません。むしろ、書かれていない部分（STL等形状の取り込み、計算点の作成、評価）の実装が想像できませんでした。こういう難しそうなところをサラッと流すのは、さすがNVIDIAさん。プロですね。

誰でも動かせるように作ってもらっているのはありがたいのですが、いきなりこれを使いだすと身につかないでしょう。何かしら簡単なネットワークを作って計算してみることから始めてみるべきでしょう。

メッシュレス、物理法則の取り込み、３次元への対応、いずれも有望です。これまでの数値解析手法に並ぶことができるのか、楽しみです。

論文作成

成り行きで、お客様ではないとある方の論文作成をお手伝いしています。

以前、結果的に先輩と2人でほぼ書き直したのですが、その査読対応でのフォローが始まりました。仕事ではないので、平日・休日問わずプライベートの時間で対応しているのですが、なかなか前に進んでくれません。

論文を書かれる研究者とお付き合いする際、逆に勉強させていただく場合があります。その場合、お相手からすると「なかなか前に進まない」と思われているのかもしれません。能力不足で御迷惑をお掛けするのは辛いですよね。

そうならないためには、身の丈に合ったお付き合いをするか、自分の能力を上げるしか選択肢がありません。いきなり能力を上げることはできませんので、長年の積み重ねが大事なのだと思います。

ま、研究職の兼務を外れましたので仕事で論文を書くことはないですし、プライベートの時間をつぶしてまで書きたいモノもないでしょう。余計な心配かもしれません。

送迎バスとセンサー

今夏、送迎バスに園児が取り残されて人命を失う、あるいは自力で脱出を余儀なくされるという人為的事故が続きました。居眠りや飲酒で児童の列に車ごと突っ込むヒトもいなくなりません。

ヒューマンエラーの一言で片付けるには重すぎる過誤です。
エラー回避策の立案と実行は簡単なのですが、それができなかったヒトがいるというのが根本的な原因ですので、ヒトに働きかけて事故を無くすという考え方では、残念ながら撲滅は難しいのでしょう。

私の自家用車には居眠りを検知するセンサーが付いているそうです。幸い、検知される機会がないのでどのような音なのか知らないのですが、おそらく色々な車種についているのでしょう。レーンキープ、衝突軽減ブレーキなど先進安全装置は既に多くの車種に装備されています。ロック忘れやハザードの消し忘れなども、アプリに通知される時代です。義務化が延期された社用車運転時の検知器による事前アルコールチェックも、いずれは車内環境をセンサーが常時測定しアプリ（サーバー）に記録、閾値を超えると責任者に通知される時代が来るでしょう。技術的には可能ですので。

園児の送迎バスにも座席毎に距離センサーを付けておけば、距離異常を検出可能になります。タイマーでスリープするようにしておけば、エンジン停止後のバッテリー上りを防ぐことが可能でしょう。20台で30万円くらい？
天井などにカメラを付けておけば、動体検知が可能です。これは5台で10万くらいでしょうか。いずれも検知したら アプリに通知やLINEで一斉配信すれば良いでしょう。これくらいのオプション、電子工作レベル（と最近わかりました）なので義務化してしまえば良いと思います。
プロならもっと安く、良い性能のオプションを作ることが可能でしょう。既に社外に向けた距離センサーは多くの車に付いていますので、それを車内に向けると良いわけですから。

性善説を前提としたヒューマンエラー対策、性悪説を前提とした予防、その両軸で進めないと、またどこかで繰り返されることになります。来夏までには対策を済ませ、皆が笑って過ごせるようにしてもらいたいものです。

*****************************
20221108 追記

今日もまた、小学生が置き去りに。
ドアを後ろに移動して、後ろからしか降りられないようにしたらダメなのでしょうか。昔は挟み込みの危険があったかもしれないけど、今は安価なセンサーがあるので解決できるのでは？

MEMS加速度センサーの +- 方向

MEMS加速度センサーの +- 方向についてよく理解できていなかったのですが、以下のサイトで理解できました。
https://lastminuteengineers.com/adxl335-accelerometer-arduino-tutorial/

測定器の中にボールが入っていて、それが測定器を押す方向をイメージすれば整合します。下向きに１gかかっていたら、下が＋。測定器を急に上に動かすとボールが測定器の下面を押すのでさらに＋の力がかかります。急に下に動かせば、ボールが浮いて上面を押すので-の値になります。
理解を難しくしていたのは、センサーの軸方向が google 先生の情報と逆だったこと。単純にセンサーを傾けてチェックすればOKでした。

ジャイロもあります。Z周りが相対的に弱いのも感覚的に掴めます。
https://lastminuteengineers.com/mpu6050-accel-gyro-arduino-tutorial/

あと、分からないのはノイズ。
https://www.analog.com/jp/education/landing-pages/003/td_accelerometer_specifications_definitions.html
https://www.analog.com/jp/analog-dialogue/articles/mems-vibration-monitoring-acceleration-to-velocity.html

ノイズ = ノイズ密度 * sqrt (BW * 1.6)

1.6は何でしょう？後者では考慮されていません。
いずれにしても、観測ではそこまで小さくなりません。なぜでしょう。

Scoops3D

Scoops3D
https://www.usgs.gov/software/scoops3d

試行円弧の3次元版です。FEMではなくて、LEMです。このようなコードが公開されるUSGS、さすがです。

現状、円弧（球）ですので実務では使用しません。これが楕円でも使用しないでしょう。山が円弧すべりをすることは稀ですし、仕様が2Dの間は利益に直結しません。しかし、技術者としては、こういった計算が海外では Open になっていることに一定の危機感を抱くべきでしょう。

マニュアルを読んでみましたが、難しいことは書かれていません。2D の試行円弧の知識があれば、問題なく読めますし結果も理解できるでしょう。2dで実施している一通りの機能は備わっていますし、ソースを少し改変すれば多くの基準に適合させることがますので、逆に加点要素として導入しやすいのかもしれません。何かしら円弧が適当な場で提案してみるのもアリでしょう。
難しくはないのですが、計算には時間がかかりそうです。単純に探索領域が3次元に拡大するためです。example の実行では比較的早く終わりましたが、ひと山となると大変でしょう。

残念ながら、こういった計算や解析に関し多くのサラリーマン技術者は他人事です。理論ではなく経験に基づいた仕事の処理で利益を生んでおり、新たなものを許容する技量と器量に欠けています。経験優先が悪いわけではないのですが、それによる思考停止は避けたいものです。

いずれにしても、海外ではコチラの内情に関係なく技術が発達します。止まっている日本の中でその恩恵を受けることができるならば、迷わずに学ばせてもらうべきでしょう。

ADXL312 & MPU6886

手元にあるMEMSセンサー（加速度）の比較です。

近所の公園に出かけ、ADXL312 と MPU6886 で振動を測定してみました。

測定条件は以下の通り。
・ADXL312: 1.5G, 2.9mG/LSB
・MPU6886: 2G, 0.06mG/LSB

・常時微動（10分程度）
・センサー周辺で連続ジャンプ

まずは常時微動。

ADXL312

MPU6886

ADXL312 の方が見かけのノイズが小さく見えます。Z（緑）がドリフトしています。
MPU6886 はその倍くらいのノイズになっています。ドリフトは小さい。

ちなみに、スペクトルは両者ともフラットなホワイトノイズのような形状でした。おそらく実際の微動はこれらの1/10以下でしょうから、微動計としては使えません。工業用センサーなら微動を測定できるのでしょうか？

次、ジャンプ。故意に振動を与えた場合です。

ADXL312

MPU6886

ADXL312は、鈍いですね。Zの揺れが小さい理由は何でしょう？

MPU6886は大きく反応しています。が、コチラは大きすぎるような。

絶対値として両者とも誤っている可能性があります。微動計とあわせて測定しないと評価できませんでした。

一概にMEMSセンサーといっても、モノによって特性が大きく異なるようです。使用する場合、その計測値を微動計の計測値と比較するなどして事前に特性を把握しておかないとダメですね。

M5StickC Plus + Blender

姿勢の算出方法に目途がつきましたので、M5StickC Plus をパッド代わりに Blender を操作します。いくつか、必要なライブラリ等をインストールして準備します。

Pyserial

https://github.com/pyserial/pyserial

https://shima2cgstudio.hatenablog.com/entry/2020/12/15/192418

Blender の modules フォルダに serial フォルダをコピー。 

C:\Program Files\Blender Foundation\Blender *.*\*.*\scripts\modules

Blender 側

Pyrhon スクリプト。

https://gist.github.com/BoKKeR/edb7cc967938d57c979d856607eaa658

M5StickC Plus 側

姿勢算出方法を書き込んだ M5 から、Bluetooth 経由で Win に Sertial 通信。ワイヤレスの方が動かしやすいかなと思いましたので、BTを選択しました。

XYZ の回転角をコンマ区切りで投げると、Blender に角度が入りました。M5 を90度傾けると Blender 内の矩形モデルが約45度傾斜します。Blender 側の Python スクリプト内の係数を2倍に設定すると M5 とモデルがシンクロするようになりました。OK、これで完成です。

Blender に地形やバイクの3Dモデルを取り込んで動かしてみました。が、テスト時の矩形と異なり安定性がイマイチでした。原点がずれている、あるいは大きなモデルだと、ビビりが目立ちます。マウスの方が安定して動かせるので、不自由ささえ感じます。制御が甘かったようです。

その点、ゲームはよくできていると改めて感心。さすがプロですね。そういえば、近年のドローンではジャイロも使って制御していたはずです。MEMSセンサーやGPSセンサー、小さなセンサーを組み合わせて姿勢や位置を制御しているのでしょう。

今回、M5 と Blender で実用的なモノは作れませんでしたが、制御技術の一端に触れることができました。また、一連の作業でセンサーの特性だけでなく、その利用法、接続の仕方、コードの書き方なども（浅いなりに）学ぶことができました。仕事で触っていた MEMS センサーに対しても理解が進みました。
やはり手は動かすべきです。今後も続けましょう。

Madgwick filter

M5StickC Plus の MPU6886 には、ジャイロがついています。

ジャイロといえばパッドで利用されています。これまで加速度しか需要がなかったのですが、世間ではジャイロがコントローラーの1機能として使われています。ふと、機能があるなら作ってみようかと思い立ちました。

Google先生に聞いてみると、同じようなことを考えられる方はいらっしゃるというか、ロボットの姿勢制御等に利用されているようでした。自動車にも入っているのでしょうね。制御工学の範疇になるようです。
現在時刻の姿勢を算出する方法として、加速度、角速度から直接回転角度を求める方法と、何らかのフィルターを通す方法があるようで、一般的には後者の方が安定しているとのこと。このあたり、制御工学を勉強しないと理解できないでしょう。

で、M5StickC Plus でどの程度できるのかを試してみることに。
角速度は時間積分しただけ、加速度からの算出式は web の中で分かりやすい説明があった以下から、フィルターは Arduino で配信されているモノ（MadgwickAHRS algorithm）のできあわせ比較です。
https://watako-lab.com/2019/02/15/3axis_acc/
https://github.com/arduino-libraries/MadgwickAHRS

結果はコチラ。

加速度からの算出値もそこそこ良い挙動ですが、細かな振動が現れています。その点、Madgwick filterは安定しています。そこそこ再現できそうです。どういう理屈なのか知りたいところですが、ひとまず恩恵を受けて進めましょう。

M5StickC Plus + LINE Notify

SGP30 を利用したタバコ検知の続きです。

窓を閉める季節になってしまい、音のデータ取得を中止しました。貯めたデータでは、低周波側でいくらか特徴のあるピークを得られたものの、ノイズと綺麗に分離できるほどでもない程度でした。

今期のコードは TVOC と eCO2 の両方が閾値を超えたら ライト点滅、ブザー発信、LINE に観測値を送る、という仕様にとどめました。eCO2 も急上昇することが分かったので、両方を組み合わせました。

LINE に関しては、異なる部屋にいたときになかなか気づかないことがあったので、その対応策です。ただ、これにハマりました。
API が用意されていますのでそれを叩くだけ & コードはネットに転がっている & どれも同じ内容、なので間違えようがないのです。でも、スマホに着信しない。Windows から叩いてみると着信しますので、環境的に使えないわけではありません。

３日ほど試行錯誤した結果、原因がわかりました。
常時データを送信している Ambient は WiFiClient でないとダメ、LINE Notify は WiFiClientSecure でないとダメ、でした。
もう一つ、void setup()  に setInsecure(); も必要でした。 

Ambient ambient;

WiFiClient client;

void ambientConnect(){

・・・

}

WiFiClientSecure clientS;

void lnotify(String msg){

・・・

}

void setup() {

　　・・・

    clientS.setInsecure();

これをクリアすれば、着信しました。

今期はココまで。また春に出番が来たら、起こしましょう。

M5StickC Plus 加速度センサー（MPU6886）その2

時間ステップはコレ↓で数10μsecの誤差に収まりました。CPUクロックは80MHzでもOKです。10msec=100Hzサンプリングの例です。

    //delay(9);

    while (millis()-10*sTime<10){

      delayMicroseconds(100);

    }

    sTime++;

EXCEL の代わりに MobaXterm でシリアル接続。最大表示＆他の作業の裏で動かしていると、欠測が生じました。また、M5側でLCDへの波形の作画自体は問題ないのですが、全消しタイミングで計測の遅れが生じました。MobaXterm 最小表示 ＆ 他の作業ナシ & M5 LCD 波形表示なしでは、欠測なく受信できました。うーん。

SDカードスロットを付けると良いのかな？
今回は、ここまでにしておきましょう。

********************************************************
20221008追記

コード全体です。Aボタンを押す毎に、波形表示・非表示を切り替えます。XYZ 軸の＋－方向と重力加速度の方向の扱いについては、イマイチ理解できませんでした。 ←解決。
https://phreeqc.blogspot.com/2022/10/mems.html

公式サンプルコードに加え、波形表示は以下を参考にさせていただきました。感謝。https://gist.github.com/TakehikoShimojima/d136e81e13eeea603a8e594a8f5ef90f#file-m5vibration-ino

/***********************************************************

Copyright (c) 2022 fal

Released under the MIT license

https://github.com/m5stack/M5StickC-Plus/blob/master/LICENSE

************************************************************/

#include <M5StickCPlus.h>

float accX = 0;

float accY = 0;

float accZ = 0;

float gyroX = 0;

float gyroY = 0;

float gyroZ = 0;

unsigned long sTime = 1;

int n = 0;

float ax[240];

float ay[240];

float az[240];

int amin = -1000;

int amax = 1000;

int btn = 0;

void header(const char *string, uint16_t color) {

    M5.Lcd.fillScreen(color);

    M5.Lcd.setTextSize(1); // フォントサイズ（文字倍率）

    M5.Lcd.setTextColor(TFT_WHITE, TFT_BLACK);//フォント色、背景色

    M5.Lcd.fillRect(0, 0, 120, 20, TFT_BLACK);// 四角形（線分）(左上x,左上y,幅,高さ,色)

    M5.Lcd.setTextDatum(TC_DATUM);

    M5.Lcd.drawString(string, 120, 2, 4);//(string, x, y, font)

    //font1:Adafruit 8ピクセルASCIIフォント

    //font2:16ピクセルASCIIフォント

    //font4:26ピクセルASCIIフォント

    //font6:26ピクセル数字フォント

    //font7:48ピクセル7セグ風フォント

}

void setup() {

    M5.begin();                // Init M5StickC Plus.

    Serial.begin(115200);

    M5.Imu.Init();             // Init IMU.

    M5.Lcd.setRotation(3);     // Rotate the screen.

    header("MPU6886",TFT_BLACK);

}

void loop() {

  M5.update();

  if (M5.BtnA.wasPressed()) {

    if (btn==0){

      btn=1;

    }else{

      btn=0;

    }

  }

  M5.Imu.getAccelData(&accX, &accY, &accZ);

  //M5.Imu.getGyroData(&gyroX, &gyroY, &gyroZ);

  ax[n]=(accX)*1000;

  ay[n]=(accY)*1000;

  az[n]=(accZ)*1000;

   

  Serial.print(millis());

  Serial.print(",");

  Serial.print(accX*1000, 2);

  Serial.print(",");

  Serial.print(accY*1000, 2);

  Serial.print(",");

  Serial.println(accZ*1000, 2);

  //Serial.print(",");

  //Serial.print(gyroX, 2);

  //Serial.print(",");

  //Serial.print(gyroY, 2);

  //Serial.print(",");

  //Serial.println(gyroZ*1000, 2);

  //M5.Lcd.printf("Temperature : %.2f C", temp);

  if (btn==0){

    if (n == 239){

      header("MPU6886",TFT_BLACK);

      n = -1;

    }else if (n != 0){

      int x0 = map((int)(ax[n - 1]), amin, amax, M5.Lcd.height(), 30);

      int x1 = map((int)(ax[n])    , amin, amax, M5.Lcd.height(), 30);

      int y0 = map((int)(ay[n - 1]), amin, amax, M5.Lcd.height(), 30);

      int y1 = map((int)(ay[n])    , amin, amax, M5.Lcd.height(), 30);

      int z0 = map((int)(az[n - 1]), amin, amax, M5.Lcd.height(), 30);

      int z1 = map((int)(az[n])    , amin, amax, M5.Lcd.height(), 30);

      M5.Lcd.drawLine(n - 1 + 5, x0, n + 5, x1, BLUE);

      M5.Lcd.drawLine(n - 1 + 5, y0, n + 5, y1, YELLOW);

      M5.Lcd.drawLine(n - 1 + 5, z0, n + 5, z1, MAGENTA);

    }

    n++;

  }

 

  while (millis()-10*sTime<10){

    delayMicroseconds(100);

  }

  sTime++;    

}

20240408追記
M5Tough + ADXL355 を用いた振動計のコードを以下に公開しました。FIRフィルタを実装していますので、周波数領域も扱えます。
https://github.com/T40O0/ADXL355_SPI_M5_SD_FIR.git

M5StickC Plus 加速度センサー（MPU6886）

M5StickC-Plus 内の加速度センサー（MPU6886）によるデータ取得。

加速度計の場合、msecオーダーの制御が必要です。受信側で時刻を付与すると、転送時間や受信側の処理速度により正しい計測時刻を設定できません。それが規則的にずれているなら補正できますが、ランダムであれば手の施しようがありません。センサー側でデータを保存するか、測定時刻を測定データとともに受信側に送るかでしょう。手元にはSDカードモジュールがないため、後者の選択になります。まずは EXCEL で取り込んでみましょう。

これまで使用してきた delay() を用いて、約10msec 間隔でデータを採取。1 loop の中で 10msec 止めるので、処理時間を加えるとそれ以上の間隔での測定になります。それが何m秒なのか知りたくて、RTC にて時刻を付与したデータをシリアル送信。

ダメでした。RTC の時刻付与は秒単位。msec単位で表示できないのでしょうか？
そうなると最初に RTC で計測開始時刻を定め、あとは一定間隔で値を取り出すしかないでしょう。仕事で使用する加速度計も同じ（計測開始時刻＋時間ステップ指定）かもしれません。

センサー側 CPU クロック 80MHｚ、delay(5)、EXCEL側 10msec で取ってみると、7ms間隔でデータが保存されていました。
センサー側 CPU クロック 240MHｚ、delay(5)、EXCEL側 10msec だと6ms間隔。
センサー側 CPU クロック 240MHｚ、delay(9)、EXCEL側 10msec だと10ms間隔。この程度でしょうか。
全てのケースで、たまに欠損があります。これはEXCEL側の問題。やはり表示がダメか。

delay()でもそれなりの結果を得られましたが、やはり正確なデータ間隔をプログラム中で設定したい。このままでは長時間の計測でズレが生じるでしょうし、EXCEL側でも他のソフト次第で欠測の多く発生することが予想されます。

時間刻み、データの保存方法について調べましょう。

Data Streamer

自宅にて測定を続けていると、窓から入る色々な匂い？にも反応していることがわかりました。全く匂わないときにもセンサーが反応することがあるので、ドリフトの影響を受けているのかもしれません。ひとまず対象を絞ることにしました。

まずは、タバコ。
隣家が換気扇を回している場合にタバコ臭が室内まで届くことがあるようで、この換気扇の音を利用する目論見です。

FactoryTest から音の FFT 部分を取り出します。整理しながらコードを見ていたのですが、サンプリング周波数がよくわかりません。FFTにかけるデータ数を指定できるのですが、それらが何秒ステップで採取するのか指定するところを把握できませんでした。include したライブラリの中に書かれているのかもしれませんが、ひとまず飛ばします。

抜き出したコードに FFT 結果のシリアル転送する部分を追加。これで値が PC にて見えるようになりました。が、転送結果をファイルに書き出す方法がわかりません。もしかしてないのか？

他の方法を調べてみると、簡易な方法がありました。EXCELです。Data Streamer という com アドインがありました。これで通信結果をターミナルのように取得できます。昔はVBAでしか通信できなかったのですが、いつの間にか実装されていました。Iot 時代の常識なのでしょうか。

使用してみると、少々難あり。１度目はデバイスの接続、解除に手間取り、接続できなくなりました。次はマウスとして認識されてしまい、計測中にカーソルが暴れて制御不能になりました。後者の現象と対処法は web 上でたくさん引っかかりました。Windows あるあるのようです。マウスとして認識されていたデバイスを削除し、COMポートの番号を変えてからは安定していますが、対処療法なのでいずれ復活するかもしれません。

EXCEL につなぐことができたら、データを容易に採取できます。サンプリングは0.01秒毎でも可能でした。200行分表示させていたのですが、その中身は時間とともに更新されます。その範囲で散布図を作っておくと、イコライザーとまではいきませんが、周波数領域のグラフが時間毎に更新されます（そういえば、振幅の単位もわからない）。EXCEL は表示が遅いので、データが飛んだりグラフに空白が生じたりしていましたが、手軽にここまでできるとは驚きでした。
FFT 結果は１度の処理で128行×３列のデータを流すようにしていたのですが、このままでは EXCEL 側で欠損します。そこで表示を 200 行から 5 行に変更。表示の負荷を小さくすることで欠損なく取得できました。これ、振動を図る際に都合が良いですね。次はRTCと振動かな。ま、時間方向にFFT結果が欠損したとしても、その数を集めてプロットすれば周波数領域での特性を把握できそうです。今回はこのまま進めます。

データを集めましょうと思いましたが、そう都合よく隣家が換気扇を回してくれるタイミングに出会えるわけではありません。自宅の換気扇でまずは試しました。
結果は上々。近いので当たり前なのですが、安いセンサーでも音が届けば分離できる可能性はあるようです。ということで、コードはひとまずこれでOK。

それでは、気長にデータを集めましょう。

換気

会社で環境測定です。

「在宅ワーク、なにそれ？」「換気？ああ、アナウンスが流れてるね」状態で人の集まるフロアーです。ひょっとするとセンサーが反応して音が鳴るかも？と思いながら計測してみました。閾値は昨夜と同じ TVOC 100ppb 。

室外でキャリブレーションを行い、 SGP30 を持って入った次の瞬間に鳴り出しました。あ、これは失敗。薬品庫のある部屋でした。感度が良い。

気を取り直し、再測定。
鳴らないですね。室内はタバコ臭よりマシでした。
と思った数分後に鳴り出しました。うーん。タバコと同じレベルか。数字にするとキツイ。オジサンたちが集まって何を出しているのか。うーんキツイ。

換気、大事です。

Ambient の制限

Ambient からデータをダウンロードしてみると、４～６秒毎の蓄積、最小時間単位が秒になっています。1Hzで発信しているので、Wi-Fiかサーバー側の問題ですね。次はMEMSで振動を計測しようと考えていましたが、これではダメかな。

バッテリー電圧は、Ambient側では小数点以下2桁まで保存されていました。LCD表示は整数です。これは CO2濃度と同様に drawNumber を使用していたのが原因。

M5.Lcd.drawNumber(vbat, 120, 110, 4);

下２桁まで表示するにするには、カーソルを合わせてから print。

M5.Lcd.setCursor(96, 110, 4);
M5.Lcd.print(vbat);

drawNumber のx座標は中央配置用の指定になっています。setCursorは左端指定でした。

テストをクリアする最小限の実装が終わったので、いよいよ実験開始。
風上の窓をほんの少し開けて、その隙間風を計測できるよう、窓の鍵にUSBケーブルをかけてセンサーをぶら下げておきました。USBケーブルは数年前に100均で購入していたモバイルバッテリーにつないでおきます。

22:46 にスタートして流していると、1時前に反応がありました。外から入ってきたタバコでした。匂うよりも先にM5StickC-Plus の発する音と光で気づきました。音が小さいので昼間に別の部屋にいると気づかないかも。これは要改善点です。

※窓全開で換気中です。

窓際に立つと、外から入ってきたタバコの匂いがします。が、部屋全体には広がっていません。屋内にセンサーを設置する関係上、これ以上は早く感知できないでしょうね。窓を全開にしている流入量が多くてと使い物にならないか？流速が小さくなって同程度？
ま、いきなりタバコに反応するということが確認できただけでも収穫です。いったい、タバコには何が入っているのでしょうね。

翌朝、バッテリーが切れてM5は停止していました。
Ambient側の記録を見ると 4:30 までのデータが保存されていました。バッテリー電圧は降下していましたが、3.7V でしたのでまだ動くはず。おかしいなあと思い確認してみると、データ数制限のようでした。

https://ambidata.io/refs/spec/

• 送信から次の送信まではチャネルごとに最低5秒空ける必要があります。それより短い間隔で送信したものは無視されます。
• 1チャネルあたり1日3,000件までデーターを登録できます。平均すると28.8秒に1回のペースです。
• 件数のカウントは0時に0クリアされます。
• チャネルデーターを削除しても1日の登録件数のカウントは0クリアされません。

0時以降、3363件のデータが保存されていました。４～６秒毎も同様に制限に引っかかっていたのですね。これ、有料版だと50Hzサンプリングの記録をミリ秒単位で保存できるのでしょうか？
https://ambidata.io/samples/vibration/vibration/

電圧降下開始が４：08でしたので、そこから40分稼働するとすると、モバイルバッテリー追加で6時間稼働する事がわかりました。

初回から良い結果を得られたのですが、閾値を低くしていたので誤検知も多くなっている状態かと思われます。Ambient の限界？もわかりました。代替えクラウドを探してみましょうか。

あとは、RTCを使ったセンサー側の時刻付与、観測時のディスプレーoff（Aボタン押下で切り替え）、Bボタン押下で音と点滅の停止、余計な電力消費のカット、マルチモーダル化（音の取り込み）、そして閾値の調整（誤検知軽減）ですかね。コツコツ実装してまいりましょう。

M5.begin()

続きです。

まず、ソースを確認。

https://github.com/m5stack/M5StickC-Plus/blob/master/src/M5StickCPlus.cpp
void M5StickCPlus::begin(bool LCDEnable, bool PowerEnable, bool SerialEnable) {
    //! UART
    if (SerialEnable) {
        Serial.begin(115200);//省略
}
    // Power
    if (PowerEnable) {
        Axp.begin();
    }
    // LCD INIT
    if (LCDEnable) {
        Lcd.begin();
    }

https://github.com/m5stack/M5StickC-Plus/blob/master/src/AXP192.cpp
void AXP192::begin(void) {
    // Set LDO2 & LDO3(TFT_LED & TFT) 3.0V

今は M5.begin(true,false,false)。AXP が begin されないと、各種電源関係が初期化されません。これでしょう。そもそも、false はどのようなケースで使用するのでしょうか？

M5.begin() = デフォ true*3 だと、電源を落とせます。USBケーブルでPCや電源につなぐと、自動で起動します。そういう仕様なのでしょう。
M5.begin(true,true,false)でも同様です。

調べてみると、M5.begin()の引数 はM5シリーズ内で異なるようです。純正のサンプルコードが誤っているのでしょうか？ま、意図したように電源が動作するようになったので解決です。

さらに続く。

M5StickC Plus + SGP30 + Ambient

続きです。

Ambient を使用するコードを転送。

書き込めましたが、Wi-Fiがダメ。

アクセス先の 5GHｚを 2.4GHz に切り替えるとつながりました。M5Stick 側で 5GHz を指定する方法はないのでしょうか？後回しにして、進めます。

サンプルコードを足しただけですので、動きました。TVOC, eCO2, 電圧が1秒毎に更新されます。バッテリー電圧の表示が一桁なのは御愛嬌。この修正も後回し。

Ambient を確認。
測定できています。グラフが書かれています。時間が経つと更新されます。素晴らしい。
この部分はPython並みにライブラリ依存ですので、間違いようがないのでしょう。

そのまま測定を継続すると、開始後40分で電圧が急降下し始め、電源が切れました。ほぼデフォで1時間持ちませんでした。残念。ま、粗々の流れは確認できましたので、良しとしましょう。

安心して再起動すると、LCDが暗いまま。

あれ？転送中の問題じゃなかった？

つづく。

M5StickC Plus + SGP30

今回購入したモノ

※イヤホンは購入していません。サイズ比較用です。小さい！

１．マイクロボード： M5StickC Plus
https://shop.m5stack.com/collections/m5-controllers/products/m5stickc-plus-esp32-pico-mini-iot-development-kit
このサイズで3軸加速度センサー、ブザー、マイク、Wi-Fi、LCDなどの欲しい機能が一通り実装されています。そして、この価格。さすが中国発。まさかの手作り？
https://wirelesswire.jp/2019/08/71960/

２．測定モジュール： TVOC/eCO2 Gas Sensor Unit (SGP30)
https://shop.m5stack.com/products/tvoc-eco2-gas-unit-sgp30
TVOCセンサーはいくつかありますが、SGP30 の感度が良さそうでした。https://www.jaredwolff.com/finding-the-best-tvoc-sensor-ccs811-vs-bme680-vs-sgp30/
価格が安く、ケーブル付きで繋ぐだけ。電圧調整も不要な M5 純正ユニットを選択。ケースに入っている分だけ反応は弱くなりそうですが。これらの組み合わせはネットでもいくつか紹介されています。
eCO2 への換算式は載っていません。テストがエタノールと水素ですので、測定された H2 と TVOC＝エタノールとみなして C2H6O＋3O2→3H2O＋2CO2 から推定しているのかもしれません。おまけみたいなモノと捉えておきましょう。
https://m5stack.oss-cn-shenzhen.aliyuncs.com/resource/docs/datasheet/unit/Sensirion_Gas_Sensors_SGP30_Datasheet.pdf

電源を入れるとFactoryTest が 始まるようです。ボタンを押すと、ジャイロ？による立方体の表示や音声のFFTが切り替わります。
用意していたHellowWorldの表示コード、加速度の表示コードでも正常動作を確認できました。初期不良はないようです。一安心。

SGP30を使用するコードをコンパイルし、転送。
動きました。SGP30も大丈夫みたいですね。

次に Ambient を使用するコードを転送。
と思いきや、転送途中でラップトップの電源が切れました。M5はどうなった？
電源をoff にした後、起動しなくなりました。ありゃ。

異なるラップトップにつなぐと認識します。ディスプレーがoffってるだけでしょうか？

再度、FactoryTestを転送すると、動きました。
USBケーブルを外して、もう一度電源 off, on。大丈夫です。

危なかった。

Arduino IDE その2

Arduino IDE で SGP30 のライブラリを検索すると、３つ引っかかりました。M5社のサンプルコードでは"Adafruit_SGP30.h"が使用されていますので、これをインストール。https://github.com/m5stack/M5StickC-Plus/blob/master/examples/Unit/TVOC_SGP30/TVOC_SGP30.ino

サンプルコードをコンパイルできるか？
通りました。問題ないようです。

次は Wifi に繋いで、Ambientにデータをアップする部分。
図書に記載してあるので、常套手段なのでしょう。Ambient さんのサンプルコードを利用します。https://github.com/AmbientDataInc/Ambient_ESP8266_lib/blob/master/examples/M5Stack_ESP32/Ambient_BME680/Ambient_BME680.ino
Arudino IDE では、Ambient_ESP32_ESP8266_lib をインストール。これでインクルードできました。

コンパイルできるか？
ダメ。wifiで引っかかります。
原因は私の記載ミス。WiFi.h を Wifi.h と記載していたので、どのライブラリを入れても認識してくれしなかっただけ。何も追加する必要はなく、素直にコピペしておけば最初から通っていました。これを修正すれば、通りました。

ひとまず、環境構築は終了。後追いだとサンプルコードが豊富にあるので楽ですね。

Arduino IDE

Arduino の開発環境です。
Arduino IDE
公式：https://docs.m5stack.com/en/quick_start/m5stickc_plus/arduino

まずは driver。これはリンクから exe 版を利用しました。

次に Arduino IDE。
インストール後に追加のURLを設定すると、依存関係を含むボード特有のライブラリを引っ張ってきてくれます。この点が支持される理由でしょう。

ボードマネージャーからボードを指定、必要なファイルをインストール。

ライブラリマネージャーも同様。GitHub から検索しているようです。引っかかったので、install を押しました。が、途中でエラー発生。

Error: 13 INTERNAL: Library install failed: moving extracted archive to destination

作成されたフォルダが読み取り専用になっていたため解除。IDE 再起動後にもう一度インストールで成功。https://github.com/arduino/arduino-cli/issues/723

この IDE では、sketchという単位でコードを管理するようです。inoファイルができていましたが、c++もインクルードされています。

'Hellow World' を sketch として保存、再読み込み、ボードを指定して検証・コンパイル。OKです。

つづく。

Autodesk Tinkercad

Arduino はイタリア発とのこと。

マイコンボードとしてはかなり売れているのでしょう。設計がOpen なので互換機もたくさん出ています。
今回は購入しなかったのですが、ちょっと試してみたい。そのような時に、思い出したのがコレ↓
Autodesk Tinkercad
https://www.autodesk.co.jp/solutions/circuit-design-software

回路作成シミュレーションが可能な web アプリです（以前は異なる名前でした）。これに Arudino が含まれています。Ardino やブレッドボード、センサー類を配置し配線することで、動きをシミュレーションできます。配置だけである程度のコードが作成されるのもありがたい。

センサーが豊富というわけではないので作成したいものをシミュレートできるわけではないのですが、使い方を把握するには手軽で便利です。

電子工作

MEMS センサーを扱ってから、電子工作に興味が出てきました。

RasPi や M5stack などにセンサーを付けるだけで（精度ははともかく）独立した計測器ができてしまします。機械学習済みのモデルを入れると、エッジAIカメラとしても使えるようですね。ネットや書店では数年前からよく見かけていましたが、これまで惹かれることがありませんでした。
が、この季節、窓を開けると時折風に乗ってくる匂いに迷惑を被っています。窓に臭気センサーかマイクを付けて、Arduino に接続してプログラムを組んでしまえば検知できるなあ、警報機はできるよね、できちゃうわ。ということで、ポチ。注文しました。ひとまずやってみようと。

パーツが届くまでにプログラムの組み方やクラウドへのデータ転送、表示の方法などを調べました。書籍にも記載されているくらい、メジャーなようです。AWS もサービスを供給していました。何とでもなりそうです。

言語は Python or C風の「Arduino言語」とのこと（by Wiki、私にはC++との違いが判りません）。後者の方が情報が多く、購入したセンサーの利用例もすぐに見つかりました。後追いは楽ですね。

ネットでは、実用化というよりもおもちゃとして作成されている方が多いように見えます。スマート農業におけるIoT機器の一躍を担っているので、すべてがおもちゃと片づけるわけにはいかないのですが、入り口が低くなってファンが増えているのは確かなようです。古参からすれば薄い「にわか」でしょうが。

MEMS を触って、C に触れた今がちょうど良いタイミングだったということでしょうね。
ま、おもちゃでも役に立つかもしれないなら試してみましょう。

3D display

3D display のデモがありました。

知らなかったのですが、既に2年以上前から裸眼での 3D display が販売されていたようです。いくつかの表示方式があるようですね。今日はそれらの中から最新版を体験。

驚きました。
画面から飛び出しているように見えます。VR ゴーグルを使わなくても、立体視ができます。面白い。

先輩が FBX のデータを持ち込まれました。
これは面白かった。トイカメラを覗いているようで、建物、橋、車が立体的に見えます。広い範囲を映したり、極手前にオブジェクトが来るような配置だと 3D 映像が崩れてしまいますが、そこそこの距離を保てば立体感抜群でした。

1枚の写真データからでも 3D 化できるということでしたの、地形の写真も試しました。が、これはダメ。専用ソフトのアルゴリズムがあっていなかったようで、高さがチグハグでした。2枚の写真で空中写真判読をする時のように並べると、正しく表示してくれるのかもしれません。

面白いのですが、仕事に有効かといわれると ？？？。
地形がダメでしたので、VR の代用くらいしか思いつきません。もう少し安ければ FPS 用に買えますが、まだまだお高い。

展示用に1台ですかね。 

open.intel

OSS の領域に Intel さんが参戦していたとは知りませんでした。

open.intel
https://github.com/intel

色々あります。
intel-extension-for-pytorch は気になりますね。当然ながら、CPU をうまく使うアプローチです。既存のGPU版とどちらが早いのでしょうか。

Intel Extension for Scikit-learn はありがたい。
https://github.com/intel/scikit-learn-intelex
https://medium.com/intel-analytics-software/save-time-and-money-with-intel-extension-for-scikit-learn-33627425ae4

(Note: The extension contains scikit-learn optimizations that were originally in the daal4py package, but it was decided to extract them into a separate Python package: scikit-learn-intelex. All future updates for scikit-learn will be available only in Intel Extension for Scikit-learn. We recommend using scikit-learn-intelex instead of daal4py.)

Intel oneAPI Data Analytics Library を使うとも書かれています。
https://github.com/oneapi-src/oneDAL

conda install scikit-learn-intelex を打ち込めば、 daal4py, dal も同時に install されました。daal4py が必要なのは説明と異なるようにも思えましたが、完全分離とまではいかなかったのでしょう。共通する部分は daal4py に残しているのかもしれません。

早速 SVC の example を動かしてみました。
が、import error を吐きます。
issue で検索すると、2021.6から修正されているようです。pip には上がっていますが、conda にはまだ上がってないとのこと。18日前の情報なので、これでしょう。issue 内の対応通り Scikit-learn を 1.0.2 に downgrade したら読めました。 conda に登録されたら、すべて update しましょう。

結果は歴然。学習時間は約1/9になっています。
extentionなし：877.６s
extentionあり：98.５s

予測が12倍。これは嬉しい。
extentionなし：13.4s
extentionあり：1.1s

xgboost, LightGBM も高速化されているようです。https://www.xlsoft.com/jp/blog/intel/2021/08/05/improve-the-performance-of-xgboost-and-lightgbm-inference/

これらは速いので、恩恵は Scikit-learn よりもかなり小さいでしょうね。いずれ、試してみましょう。

pybind11

手は動かしていませんが、忘れそうなのでメモ ( ..)φ

pybind11
https://github.com/pybind/pybind11

pybind11 is a lightweight header-only library that exposes C++ types in Python and vice versa, mainly to create Python bindings of existing C++ code. Its goals and syntax are similar to the excellent Boost.

C++ 側の書き込みは最低限で、私のような初心者でもわかりやすいのがgood!
ネットの書き込みを見ていますと、Python 側の list は C++ 側の vector や list に変換してもらえるとのこと。オーバーヘッドも Cython より小さくなりやすいと評価される方もいらっしゃいました。

本当ならありがたい。まずはこれらから確認してみましょう。。

カメラの位置推定

後輩君がマネージャーからの指示で点群を触っていました。

マネージャーがうまくビジョンを伝えられていないようで、指示を受けた内容に対し手を動かして少しづつ進んでいる状態のようでした。私も何度かマネージャーから話を聞きましたが、着地点から逆算できるような筋の通った説明は得られませんでした。これだとルートを考えられないので、指示をこなすしかありません。どちらかというと、後輩君の方が問題の本質と解決手段を理解しているような印象を受けましたが。

その中で、sfm の逆のように点群と画像からカメラ位置と角度を推定したいという要望を聞きました。これ、知らなかったのですが PnP 問題と呼んでいるそうです。
この辺とか https://daily-tech.hatenablog.com/entry/2018/01/21/185633 この辺が https://daily-tech.hatenablog.com/entry/2018/03/22/022218 わかりやすい。

ひずみを取る、既知点の座標を 2D & 3D から指定するという手作業が発生しますので、完全に機械化はできないのでしょう。簡単なモデルだと自動でできるようですが、限られているでしょうね。

OpenCVでは実装されているようです。種類が多い。https://docs.opencv.org/4.x/d5/d1f/calib3d_solvePnP.html

ま、後輩君なら大丈夫でしょう。

ArcGIS Proに入れない

朝、後輩君から「ArcGIS Pro に入れません」と報告を受けました。

画面を覗くと、サブスクリプション関連のエラーが吐かれています。
そういえば、全社でライセンスを統合すると聞いていました。

社内外の担当に調べてもらって理由がわかったのがその日の夕刻。
原因は米国ESRIの統合時の処理方法でした。ライセンス統合作業を進めるために 日本ESRI と ユーザー側に連絡なしにアカウントを使えない状態にしてしまったとのこと。日本ではありえないのですが、海外なら普通にあるかもと思わせる出来事でした。

新しくアカウントを作り直せば入れるように設定できるため、後輩君にはそれで対処してもらいました。

試行錯誤を１日続けていたため、特に担当の方々はお疲れだったと思います。
ま、今後も起こるであろう問題への対処法をつかめただけ、良しとしましょう。

QGIS のリレーション

QGIS からは PostgreSQL 内のリレーションが見えませんでした。

が、それらを全く扱えないわけではないようです。https://docs.qgis.org/3.22/ja/docs/user_manual/working_with_vector/attribute_table.html#creating-one-or-many-to-many-relations

QGIS 内であらためてキーを設定することで、リンク先の情報まで地物情報に表示されるようになりました。

簡単ですが、惜しいですね。もう少しです。

QGIS + PostgreSQL + 日本語

一息ついたので、QGIS + PostgreSQL を実践しています。

用意したのはガッツリ日本語の入った csv。

以前購入していた図書を参考に、PostgreSQL が日本語にどこまで対応しているのか知りたいというのと、QGIS から PostgreSQL へデータを取り込む方法のトレースが目的です。後者は図書を読むまで知らなかったんですよね。

まずはSQL。

ジオメトリ作成は日本語 out。””ナシだとエラー。”緯度”というように””で囲むとエラーはかからなかったのですが、ジオメトリとして認識されませんでした。日本語でカラムを指定する方法がわからず、Lat, Lon で csv を書き換え。

csvを取り込む際に小文字変換のチェックがあったのですが、大文字だと引っかかるみたいですね。ま、キホン英数小文字というのはわかるのですが、何か日本語を通す書き方はあるのだと思います。

ALTER TABLE s_aaa.tb_bbb ADD geom GEOMETRY(POINT,4612);

UPDATE s_aaa.tb_bbb SET geom=ST_GeometryFromText('SRID=4612;POINT('||lat||' '||lon||')');

※この後、空間インデックス作成

型変換は日本語カラムOK。

ALTER TABLE s_aaa.tb_bbb

ALTER COLUMN "あああ" TYPE int4

USING "あああ"::int4

ただし、Python と異なって、nan を含む場合は float にも変換できませんでした。これも何かコツがあるのでしょう。

キーの指定は日本語OK。テーブル名、カラム名ともに日本語でも問題ありませんでした。

次にQGIS

csv を shpにするとカラム名が半角８文字でぶつ切りにされます。が、csv を PostgreSQL に取り込むと、問題ナシ。SQL でジオメトリを作ってから QGIS に持っていくと、長い日本語カラムでも欠かすことなく取り込めました。

QIS からは geometry を持つテーブルしか見えません。オブジェクトをクリックするとその属性が pop up され、さらに外部キーから参照テーブルの情報を表示させる、といった機能はなさそうです。これに期待していたので、もう少し探してみようと思います。

一通り触ってみましたが、QGIS の機能の多さにあらためて驚き。
SQL は使わないかな。QGIS 側で地理情報のみの表示であれば、Python で作ってしまった方が柔軟かつ速いと思います。

Infraworks2023でPLATEAU その2

続きです。

1) のテキスト置換を Python で組んだのですが、配布しても使えない方が多い。

ということで、C++で作ってみました。

といっても、C++ は初心者。既存のコードを読んだりコンパイルしたりしたことはありますが、書いたことはありません。久しぶりにプログラミングの入門書を買ってきました。

読んでみると、独学で振り落としてきた重要なパーツがちらほら。半分ほど読んでみましたが、勉強になりました。

で、ネットでコードを拾いながら３日かかって完成（Python なら１時間かからないのに）。特にパスの型が難しかった。意外と、欲しい情報がなかったので、図書を買ったのは正解でした。

https://github.com/T40O0/CityGML-Replace

python も C++ も、バッチファイルを書いて ２) の投影変換と連携させました。ダブルクリックで１）２）を済ませてくれます。できれば C++ で Java まで取り込みたかったのですが、今は実力不足。いずれ整理しましょう。いえ、そのころには PLATAEU で 全都市 FBX が整備されているでしょうね。

Infraworks2023でPLATEAU

Plateau のデータを WebGIS に載せるところまでは確認していましたが、それ以降、触っていませんでした。

少し面倒だったのであらためて手順を整理しておきたいと思っていたところ、社内でも同様の考えを持った方と遭遇。これを機に、整理してみました。

まずは Infraworks 等での利用。
fbx が配布されている場合は、迷わずそちらを選びましょう。
問題は、CityGMLしか配布されていない場合。以下の3手が必要です。

GML⇒
①GML内のアドレスを置換
②citygml-tools で直角平面に投影変換。
③FMEフリー版でIMXへ変換（axis=2,1,3）
⇒Infraworks

1）GML内のアドレス置換

3D都市モデルのデータ変換マニュアル
https://github.com/Project-PLATEAU/Data-Conversion-Manual-for-3D-City-Model

ただし、i-UR1.4は2021年9月にi-UR1.5に改定されました。これに伴い、URLが変更されました。そのため、本ソフトウェアの利用にあたり、3D都市モデル（CityGML形式）及びソフトウェアに記述された旧URLを新しいURLに更新する必要があります。具体的には以下の手順に従い、更新してください。

これは置換のみです。大量の GML データに対しては Python で処理しました。これを飛ばすと3）でエラーを吐きます。

2）citygml-tools で直角平面に投影変換

日本のデータの変換について
https://github.com/tokyo-digitaltwin/citygml-tools

現在、CityJSON-Tools　が利用する空間参照系ライブラリでは、3D空間参照系間の参照系変換のための十分な環境整備ができていません。そのため、空間参照系の変換を行う際は、該当する2D空間参照系を利用して変換を行います。高さ方向は一般的に変換による差が小さいため、実用上変換は不要であるとして、垂直方向では元となるデータの値を変換後もそのまま利用します。以下で具体的なコマンドを上げて、変換方法を説明します。

投影変換にcitygml-tools を利用します。Javaの実行環境（JDKのインストール）が必要です。コマンドプロンプトから以下のコマンドを打つだけで変換してくれます。マニュアル通りにXYを入れ替えてもOK。

citygml-tools reproject --source-crs=6668　--target-crs=6677 --lenient-transform　*.gml

※6677はEPSGコード。東京⇒JGD2011-9系⇒EPSG:6677

3）FMEフリー版

CityGML Importer for Autodesk
https://www.safe.com/citygml-importer/

メールが送られてきます。

・インストーラーは有料版と同じ。

・ライセンスファイルをリンク先から落として記載の位置に入れる。

・ワークスペースをリンク先から落として記載の位置に入れる。

・FME Quick Translator を起動。

記載の手順でRun。

Axis は「2, 1, 3」。

2）citygml-toolsで投影変換と同時に軸を入れ替えた場合は「1, 2, 3」

これでIMXファイルの完成です。
あとは Infraworks に取り込むだけ、ですが時間がかかります。

配布されている都市単位で1）2）が10分程度、3）が20分程度、Infraworksの取り込みが1.5時間かかりました。再表示も30分以上かかります。

取り込んでも建物全体が見えない地域もありました。変換ファイル数を少なくしてもダメな場合があり、変換かInfraworks 側かわかりませんが、まだノウハウがあるのでしょう。
それは追々。

MEMS

先週、MEMSセンサーを購入。

今まで省電力、小型、現場で長く持つ、などと聞いて来ましたが、実物を見たのは初めてでした。

本当に、小さい！

加速度センサーだったのですが、5～6㎜四方。これ、はんだ付けする自信はありません。基盤についた商品を購入してよかった。手に負えないところでした。スマホを分解したことはありませんが、このような小型センサーが所狭しと並んでいるのでしょう。モノづくり技術者にリスペクトです。

制御ボード、インターフェースボードなどを加えて独立した観測機器とするのでしょうね。以前よりラズパイでの制御もネットで見ますが、それではセンサーの小型化といった特徴を損なっています。小型パーツを組み合わせ、工作されるのでしょう。ノウハウを学びたい。

初心者ですので、無理はしません。PCにつなげて初期不良のないことは確認できましたので、料理法をゆっくり考えましょう。

 

位相 unwrapping

位相の unwrapping について調べていました。

SAR ではなく、単なる信号の再現（逆フーリエ）でこの話題が出てきました。

調べてみると、ありますね。

位相アンラッピング(Phase unwrapping)
http://retrofocus28.blogspot.com/2013/12/phase-unwrapping_26.html

def mywrap(y):
    for i in range(1, len(y)):
        div = y[i] - y[i-1]
        if div < -pi:
            y[i:] += 2*pi
        elif div > pi:
            y[i:] -= 2*pi
        else:
            pass
    
    return y

numpy なら1行です。

np.unwrap() 

簡単ですね。
理想的な条件なら再現できるようですが、実際のデータでは難しいようです。

Cython

Pandas の公式に速度改善に関するチュートリアルがあります。
https://pandas.pydata.org/pandas-docs/stable/user_guide/enhancingperf.html

CPython ではなく、Cython の利用です。
C で書いたライブラリを Python で読むのではなく、Python で書いた関数を C に渡してコンパイルする形。渡す際に型を定義したり、型チェックを外すことで速くすることが可能になっているようです。これは、使い勝手が良さそう。

プロファイルは動かなかったのですが、それ以外は動きました。遅い処理はこのような簡易な処理方法があることも覚えておきましょう。

C++ & Python

C++ で Pythonのモジュールを作りたい。

初めは Python と C++ の exe をバッチファイルで繋ごうかと思いましたが、ライブラリにすればPython から直接読めます。

チュートリアルがあり、GitHub にテンプレートが公開されていました。それをトレースしましょう。https://docs.microsoft.com/ja-jp/visualstudio/python/working-with-c-cpp-python-in-visual-studio?view=vs-2022

まずは環境設定。
VS2019 は Python3.7 までのサポートのようです。これは古い。VS2019 を削除し、2022に入れなおし。これでデフォが3.9？でしょうか。デフォを3.10にしたかったのですが、入れ替え方がわかりません。
空のPythonプロジェクトを作り、新たに conda 環境を Python=3.10 で作成。これを新しいプロジェクトに対する既定の環境に設定。あらためて空の Python プロジェクトを作り、テンプレートのコードをコピぺしました。が、これではダメ。再度開くと3.9に戻ります。
Google先生に伺ったところ、[全般][インタープリター]で変更できました。
Python のみで動くように初期コードに書き換えた後、3.10で動作確認。OK でした。

次に C++。
既存のプロジェクトとして ’superfastcode’ を追加。VS2017 製でしたので[プロジェクト][ソリューションの再ターゲット]にてアップデート。’superfastcode’ というライブラリとして Python から呼び出せるとのこと。ソースファイル module.cpp は既に追加されている状態です。
プロジェクトのプロパティから、いくつかの設定を変更します。
・[全般][プロジェクトの既定値][構成の種類] dllを指定。
・[C/C++][全般][追加のインクルード ディレクトリ]新たに作成したenv下の include フォルダーを追加。
・[C/C++][コード生成][ランタイムライブラリ]マルチスレッドDLL (/MD
・[リンカー][全般][追加のライブラリディレクトリ] ]新たに作成したenv下の libs フォルダーを追加。

これで一度ビルド。OKです。
’superfastcode2’ も同様に取り込んで設定。採用していた conda 環境に pybind11 をインストールしておきます。念のため、Win11 を再起動してからビルド。
うーん、引っかかります。Python3.6が必要と出ます。どこでしょうか。

結局、原因がわからなかったので’superfastcode2’ はアンロード。
リリース環境でも同じ設定でビルド。で、試算。

Running benchmarks with COUNT = 500000
[tanh(x) for x in d] (Python implementation) took 0.512 seconds
[fast_tanh(x) for x in d] (CPython C++ extension) took 0.123 seconds

やはりCの方が速い。
チュートリアルを見ると pybind11 を使うと ’superfastcode’ より遅くなっていますので、ここまでにしておきましょう。

C++ でライブラリを作り、Python で読み込んで計算するチュートリアルでしたが、分かり難いですね。これ、gcc でも設定できるでしょうか。自信ないですね。

gcc と エンコード

C ++の初心者です。

C、C#とともに何度か入門しましたが、実務ではそれほど用がなかったため、覚える前に辞めてしまう状態が続いています。此度、実務（Windows10）で使うことになり、あらためて入門です。

Visual Studio 2019 を使ってコンパイルしていると、scanf でエラー。

error C4996: 'scanf': This function or variable may be unsafe. Consider using scanf_s instead. To disable deprecation, use _CRT_SECURE_NO_WARNINGS.

そうでした。https://phreeqc.blogspot.com/2021/05/c.html
scanf_s を使うと通りましたが、デバッグでエラー。

test.exe (プロセス *****) は、コード -****** で終了しました。

忘れていますね。scanf_sで文字列を入力する場合には第三引数で最大配列数を指定する必要がありました（引数不足のwarning が出ていましたが、無視していました）。 

引数を指定し、できたexe を実行！
動きました。

表示がすぐに消えてしまうので、「Enterkey を押して終了させたい。」をコピペ。
https://teratail.com/questions/140273
これでOK。

次はgcc
tdm64-gcc-10.3.0-2 を Win11 に入れて、ターミナル起動。gcc -v で動作確認。

>gcc -v
gﾝﾝ spec gpｵﾄ｢ﾜｷB
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/10.3.0/lto-wrapper.exe
・・・
T|[gｳﾄ｢ LTO ｳkASY€: zlib zstd
gcc o[W 10.3.0 (tdm64-1)

文字化けします。
コマンドプロンプト表示をUTF-8に変更

>chcp 65001
Active code page: 65001

>gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/10.3.0/lto-wrapper.exe
・・・
algorithms: zlib zstd
gcc version 10.3.0 (tdm64-1)

OKです。では、コンパイル。コードはUTF-8。

>g++ -o test test.c

できたexeを起動すると、文字化けしています。Win の実行ファイル側では shift-JIS 固定なのでしょう。
実行ファイルでのエンコードをCP932として指定します。

>g++ -o test test.c -fexec-charset=CP932

うまく変換できているように見えます。もう少しテストが必要かな。
ひとまず、環境構築は完了です。

PC Building Simulator

PC Building Simulator というゲーム。

購入したもののずっとインストールすらしておらず、昨晩、初めて触りました。
日々、PC のトラブルが持ち込まれ、それに対応していくモードが主でしょうか。ウイルス駆除だの、電源交換だの、メモリ増設だの。これ、やっていて辛くなってきます。実世界のシステム担当者の仕事を代わりにやっているようで、遊んでいる気になれません。ええ、あらためて担当者に感謝です。

フリーで作成できるモードは遊んでいる感じがします。高級な部品を思う存分使用できます。ケースが狭くてリザーバーとラジエーターが干渉したり、電源が足りなかったり。実機で組む前に試す、なんて使い方もできそうです。
3DMarkのベンチマークが動く、エラーで落ちるなんてところまでシミュレートしてくれます。CPU, GPU ともに水冷で組んでみましたが、GPUの温度はどこに出るのでしょう。クロック変更はどこかな？

これなら壊す心配がありません。不慣れな人の練習用に使っても良いかも。

テレワーク その2

第7波が全国的に猛威を振るっています。

第1波から現在まで、勤務先の私を含む周辺では一度もテレワークが許可されていません。executives はイロイロと対応できない？しない？まま、基本、全員出社が平常運転となっています。いえ、体調不良者や濃厚接触者、妊婦のいる家庭等一部の方は execcutives の許可によってテレワークが認められています。ただし在宅手当なし、と通達されています。
一部の支店ではテレワークを導入しているようですので、リクルート時には「在宅勤務導入」を表に出しているようです。ま、嘘ではないですね。

今日、知人から電話がありコロナ対応状況について話していました。
同業他社でも「小学生以下の子供がいる人はテレワーク」「普通にテレワーク」など。7波の特徴に応じた対応を取っている会社は、さすがと思います。就活では本質を見抜き、そういった会社を選んでいただきたいものです。

CIMへの対応？ その2

他支店で、地質CIMに対応せざるを得ない状況に陥った方々が現れました。

受注者希望だったそうですが、発注者要望を受けているという少し変わった状況。で、地質の executives 経由ではなく、CIM の executives 経由で相談が届きました。
地質調査業務であれば、詳細度はさて置き3次元可視化する程度しか要望されないと思います。であれば、特に難しい点はありません。が、皆さん3次元化を未だに避けていらっしゃるので、できる人を探している状なのでしょう。
https://phreeqc.blogspot.com/2021/08/cim.html

ちょうど GEORAMA の更新時期になりました。今年から地質部門は新しい executive に変わったのですが、保守更新については判断保留になりました。私が地質を離れてからどなたも使われていませんし、地質担当の方々は未だに2次元CAD 主体ですので、即決できなかったのでしょう。国が CIM を進めている中、何らかの理由でこれまで3次元化を避けてきた、今後も積極的にはかかわらないというのが、新しい executives が引き継いだ経営方針なのかもしれません。

設計部署は数年前から Autodesk 製品を使って人を育てていますし、CIM に対応しています。現実的に地質CIMでは GEORAMA が納品への近道の一つですので、私なら更新は即決ですが、この会社でその判断は誤っているのかもしれません。

さて、相談は来ても、ソフトはない、人は育っていない、誰も手を動かそうとしないでは、弱りました。いえ、executives は予測していた結果でしょうから、それほど困っていないのでしょう。
本当に困っているのは担当者のみ。あれ？昔から変わっていないですね、この悪習。
さて、どうしましょうね。

弾性論とFEM

FrontISTR の ハンズオンに参加しました。

FrontISTR のビルドやインストール部分は時間の都合で割愛。メッシング は商用ソフトからのインポートが主流になりつつあるとのことで、これも割愛。ジョブを投げて計算を回すのに慣れるといった5時間でした。うーん。東大のスパコン？を使わないのでそこのジョブ管理を覚える必要はなかったのですが。ま、スパコンの使い方は似たようなものだと確認はできました。

機能紹介や理論の部分に METIS による領域分割、FEMにおける弾性論など懐かしい話題が含まれていました。FEMでの展開を忘れかけており、帰宅してから一通り復習です。

**************************************
線形変位場を仮定
→要素変位ー節点変位関係・・・補間関数[N]を使用
{ue}=[N]{u}
→物理空間座標系と計算空間座標系（正規化座標系）の関係・・・写像関数
アイソパラメトリック要素:補間関数と写像関数が同一(形状関数)

ひずみ-変位関係（支配方程式１）
→FEMでは要素ひずみー節点変位関係
{ε}=[B]{u}・・・（要素内の）ひずみは（節点）変位の空間勾配
Bマトリックス（ひずみ-変位マトリックス）：[B]=∇[N]・・・補間関数の空間微分

応力-ひずみ関係（フックの法則、支配方程式2）
→FEMでは要素応力-節点変位関係
σ=E{ε}・・・1次元
{σ}=[D]{ε}＝[D][B]{u}
Dマトリックス（応力-ひずみマトリックス）:ばね。ラメ定数でも表現可。

仮想仕事の原理(支配方程式3)
外力の仮想仕事(力ｘ変位)と
δW=1/2{u}T{f}：節点変位ベクトルと節点外力ベクトルの積
仮想ひずみエネルギー（ひずみｘ応力）の
δU=1/2∫{ε}T{σ}dV=∫[B]T{u}T[D][B]{u}dV
つり合い(要素剛性方程式)
{u}T{f}={u}T∫[B]T[D][B]{u}dV
{f}=[k]{u}
kマトリックス:k=∫[B]T[D][B]dxdydz
→FEMではヤコビアン|J|+正規化領域でガウスの数値積分：k=∫[B]T[D][B]|J|dξdηdζ

全体剛性マトリクス
{F}=[K]{U}

変位の計算
{U}=[K]-1{Ｆ}

節点変位→要素ひずみ→要素応力
応力は積分点で計算。
節点変位は要素間で連続（適合）するが、ひずみと応力の連続性は保障されない。
応力・ひずみは変位の勾配に比例する物理量であるため、変位より精度が低くなる。
実用的には積分点の値を節点に外挿し平均した節点平均応力、節点平均ひずみが用いられる。

※アイソパラメトリック1次要素→せん断ロッキング→低減積分で回避→アワーグラスモード→非適合要素（要素間の変位の連続性が保たれないものの、実用的には問題ない程度）

※四辺形1次要素の完全積分:４点（2次積分）、2次要素9点（3次積分）→低減積分：1次要素1点、2次要素4点

Pyrhon ＋ GDAL

Raster を複数使用した安定計算。

解像度の違いでエラーを吐きます。
Arc で ポリゴンからラスターへ変換していたのですが、解像度が完全には一致しないようでした。仕方がないので、リサンプリングすることに。
以前、GDAL を使って Raster を読み込んでいたので、そこに追記。

# Open Raster
def open_src(file,band):
    src = gdal.Open(file, gdal.GA_ReadOnly)
    if src is None:
        print('file not found')
    else:
        print(src.GetDescription() )
        xy=src.GetGeoTransform()#
        print(xy) 
        print(src.GetProjection())
        print(src.GetMetadata() )
        arr_b1=src.GetRasterBand(band).ReadAsArray(0)
        print(arr_b1.shape)
        nodata=src.GetRasterBand(band).GetNoDataValue())
        print(nodata)
        return src,xy,arr_b1,nodata


#図化
def plot_raster(arr_b1, nodata):
    arr_b1 = np.where(arr_b1==nodata, np.nan, arr_b1)
    print('max:',np.nanmax(arr_b1))
    print('min:',np.nanmin(arr_b1))
    print(arr_b1.shape)    
    plt.imshow(arr_b1)
    plt.colorbar()
    return arr_b1

     
#測地系変換
#EPSG：4612 から 4326に変換して保存
ds = gdal.Warp(file_out, file_in, srcSRS="EPSG:4612", dstSRS="EPSG:4326")
ds=None


#2つのラスターの領域、解像度を（ほぼ）揃える
#切り出し領域を決める
minx=136
miny=35
maxx=137
maxy=36
#領域 をforium で確認
fmap = folium.Map(location=[(miny+maxy)/2,(minx+maxx)/2], zoom_start=8)
line = [[miny,minx],[miny,maxx],[maxy,maxx],[maxy,minx],[miny,minx]]
folium.PolyLine(locations=line).add_to(fmap)
fmap
#ピクセル数を決める
width=18000
height=11100
#resampling
cut_file0="./output/cut0.tif" 
cut_file1="./output/cut1.tif" 
ds=gdal.Translate(cut_file0, file_in0, projWin=[minx,maxy,maxx,miny], width=width, height=height, noData=nodata)
ds=gdal.Translate(cut_file1, file_in1, projWin=[minx,maxy,maxx,miny], width=width, height=height, noData=nodata)
ds=None


#保存
def save_dst(rasterfile, arr_new, nodata, src):
    driver = gdal.GetDriverByName("GTiff")  
    dst = driver.Create(rasterfile, 
                       xsize=arr_new.shape[1],
                       ysize=arr_new.shape[0],
                       bands=1,
                       eType = gdal.GDT_Float32)
    dst.SetGeoTransform(src.GetGeoTransform())
    dst.SetProjection(src.GetProjection())
    band1=dst_ds.GetRasterBand(1)
    band1.WriteArray(arr_new)
    band1.SetNoDataValue(nodata)
    dst.FlushCache() 
 

Arc でラスターを2枚作成・書き出していましたので、マルチバンドにしてから1枚出力すれば解像度が完全に一致したのですが。
https://pro.arcgis.com/ja/pro-app/2.8/tool-reference/data-management/composite-bands.htm

ま、今後使えるかもしれませんので残しておきましょう。

 

float, double

C#で作成されたソフトと、Pythonで書かれたスクリプト。
同じテキストファイルを読んで計算したのですが、答えが異なりました。

原因は変数の型。

C# は float（32bit） で読み込むように作られていました。Python の float は C の double (64bit)で作られているそうですので、そのまま読み込むと前者の桁が小さくなります。https://docs.python.org/ja/3/library/stdtypes.html

後者を np.float32 で変更すれば良いのでしょうが、そもそも、input ファイルの段階で有効桁数以上の数字を消していなかったのがマズイ。不要な桁を消して読み込めば、一致しました。

手を抜くとダメですね。すぐに帰ってきます。

 

斜面の向きと地震波

地形量と崩壊の関係を整理してみようかと考えていました。

DEM が １枚あれば、多くの地形量を作ることができます。それと相関性を見てみるだけです。GIS 仕事です。

文献を探してみると、地震時でトライされた内容がありました。誘因として雨の要因を考慮する必要がなくなるため、地震時の方が良いのでしょう。

安田ほか（2006）動的振動解析による地震時の加速度応答および斜面変位と地形効果に関する考察
https://doi.org/10.11475/sabo1973.59.4_3

統計的に整理するに数が少ないと思われますが、方針はこれと同様です。

地震時だと揺れの方向と斜面の方向が関連するので、数値実験の方が整理しやすいでしょう。斜面のモデルは少なくても、揺れの方向を自由に変えてその応答を見れば良いだけですので理屈は簡単です。時間はかかりそうですが。

実際の地震波では、斜面の法線方向から下向きの場合が最も危険になるでしょう。以下では、ある程度考慮されているようです。下向きの角度はよくわかりませんでしたが。

篠田ほか（2018）広域的な地震時斜面崩壊危険度図の作成方法
https://doi.org/10.2208/jscejge.74.177

どなたか整理されていないでしょうか。単純な話なのでありそうですが。

六甲式の改良

地震時崩壊を調べていると、国土地理院の方々が書かれた以下の文献を見つけました。

神谷ほか「地震による斜面崩壊危険度評価判別式「六甲式」の改良と実時間運用」写真測量とリモートセンシング 51 (6), 381-386, 2013

F＝0.075s－8.92c＋0.006a－0.3228 （六甲式)
G＝4.38log(s－119c)＋3.93loga－15.27（修正六甲式)

例えば曲率が0の場合（地表面の凹凸がない場合），傾斜が43度以上では  PGA が 0gal であっても崩壊と予想する。また，PGA が583gal（計測震度5.8に相当）以上では，傾斜が0度であっても崩壊すると予想する。本システムでは，このような場合でも，妥当な結果を得る必要がある。

言われてみると、納得。閾値を調整すれば六甲式でもカバーできそうですが。

六甲式，修正六甲式は，PGAに関係する項と，関係しない項の和である。PGAに関係しない項を，それぞれ「部分六甲式」，「部分修正六甲式」と呼ぶ。あらかじめ，小グリッドの全てのセルにおける部分六甲式の値を計算し，大グリッドの全てのセルにおいて，「大グリッドのセルに含まれる小グリッドの部分六甲式の値」のヒストグラムを作成する。地震発生時には，大グリッドの部分六甲式のヒストグラムと各階級の代表値を用いて，六甲式等のヒストグラムを計算する。この計算は，大グリッド単位で実行するため，小グリッド単位の計算と比較し，大幅に計算量を削減することができる。

あらためて10mメッシュの傾斜角や曲率を計算し表示してみましたが、2億メッシュでも数十秒～数分程度でした。うまい工夫だと思いますが、今では10mメッシュのままでも計算できるでしょう。日本全国の10mメッシュで s, c の項を計算しておけば、地震時に PGA や震度情報を取り込むことですぐに予測結果を得られます。

今後、高度化を目指すなら、まずは5mDEMの併用でしょうか。角度や曲率を計算する際に海岸線沿いの崖などが考慮できない点をカバーできると思われます。
シンプルかつある程度検証されてきた式ですので、地震のたびに検証・更新される方がいらっしゃっても良いように思います。一番手間の書かかるポリゴン制作も衛星画像等から精度よく自動作成できるようになれば、一気に実現しそうです。が、まだ先ですかね。

 

Pandas DataFrame の不一致

Pandas で EXCELから読み込んだ内容と、csv に変換してから読み込んだ内容が一致せず、どこが原因かを探っていました。

最初はエンコードを疑いましたが、違いました。
原因は数字の表記でした。

In [2]:

df_excel=pd.read_excel('Book1.xlsx')
df_excel

Out[2]:

	A	B
0	100	100

In [3]:

df_utf=pd.read_csv('Book1u.csv',encoding='utf-8')
df_utf

Out[3]:

	A	B
0	100.0	100.0

In [4]:

df_sjis=pd.read_csv('Book1s.csv',encoding='cp932')
df_sjis

Out[4]:

	A	B
0	100.0	100.0

In [5]:

df_excel.equals(df_utf)

Out[5]:

False

In [6]:

df_excel.equals(df_sjis)

Out[6]:

False

In [7]:

df_utf.equals(df_sjis)

Out[7]:

True

In [8]:

df_excel.dtypes

Out[8]:

A    int64
B    int64
dtype: object

 

EXCEL上で 100.0 と表記していると、csv 経由では float、excel からは int として読み込まれる仕様でした。