残念なオジサマ

講習会に来られていた少し残念なオジサマ。

小人数であったためか、講師が話されている途中で質問を投げかけます。熱心に質問・議論をされる姿勢は良いことなのですが、そのレベルの低さ、頻度、口調には閉口。もっと勉強してから来られるべきでしたね。
度々、進行が止まるため周りの方も迷惑そうでしたが（受講者がしびれを切らして回答される場面もありましたが）、一番前の席だったためか、元からそのような性格なのか、最後までお気付きになりませんでした。

残念ながら、私にも少し身に覚えがあります。
社内会議でレベルの低い技術報告が出てきたときには指摘すべきかどうか迷います。executives が（知ってか知らずか）指摘しないのだからこのままで良いのでは？などとも思うのですが、意見することはあります。が、たいていの場合は話になりません（基礎を疎かにしていなければ、最初から誤った報告は出されないでしょう）。で、結果的には多くの人が共有すべきでない時間を作ってしまいます。さらに空気を読む必要があるのでしょう。

若手なら育てる義務がありますが、オジサマを育てる義務はありません。
無知は罪。自分はそうならないよう、毎日コツコツ積み上げましょう。

SAR 緊急観測オーダーと課題

だいち防災 web。
久しぶりに訪れると、公開内容が増えたように感じました。気のせいでしょうか？
https://jaxa-dis.maps.arcgis.com/home/index.html

危機管理のための SAR 緊急観測オーダーが普及しつつあるようです。このサイトを見る限り、お盆前後の台風10号では、紀伊半島や中四国等にオーダーが出されていたようです。
今回の佐賀県豪雨では、28日昼に確認した時点で強度画像の合成結果と手動判定結果が公開されていました。が、使っているのは28日０時撮影のデータ。判定された浸水域は既に過去のものとなっており、危機管理の視点では情報が古く使えないでしょう。図らずも、次の課題（リアルタイム把握）が見えてきました。

ま、せっかくの衛星ですので、このような目的でも利用頻度は上げて頂きたいですね。

pooling

配列を pooling する方法です。

CNN ではなく、ただ単に配列を最大値で pooling したいと考え、その方法を探しました。
最初は numpy や scipy にあるだろうと簡単に考えていたのですが、見つけることができず。

最終的には、tensorflow の max_pooling を利用しました。このような感じ↓
******************************************

In [1]:

import numpy as np
from scipy.ndimage.filters import maximum_filter
import tensorflow as tf

In [2]:

x1=np.array([[0,1,2,3],[8,9,10,11],[4,5,6,7],[15,14,13,12]])
x1

Out[2]:

array([[ 0,  1,  2,  3],
       [ 8,  9, 10, 11],
       [ 4,  5,  6,  7],
       [15, 14, 13, 12]])

In [3]:

x = x1.reshape([1, 4, 4, 1])
x

Out[3]:

array([[[[ 0],
         [ 1],
         [ 2],
         [ 3]],

        [[ 8],
         [ 9],
         [10],
         [11]],

        [[ 4],
         [ 5],
         [ 6],
         [ 7]],

        [[15],
         [14],
         [13],
         [12]]]])

In [4]:

x.shape

Out[4]:

(1, 4, 4, 1)

In [5]:

maximum_filter(x1,size=[2,2])

Out[5]:

array([[ 0,  1,  2,  3],
       [ 8,  9, 10, 11],
       [ 8,  9, 10, 11],
       [15, 15, 14, 13]])

In [7]:

pool=tf.nn.max_pool(x, ksize=[1, 2, 2,1], strides=[1, 2, 2, 1], padding='SAME')
se = tf.Session()
pool = se.run([pool])
np.array(pool).reshape([1,2,2])

Out[7]:

array([[[ 9, 11],
        [15, 13]]])

メモリオーバー

興味本位でいくつかを対象に FFT や PCA をかけています。

これが興味深い。想像してなかったことや、「考えると当前」といったような結果が浮かび上がります。先の 1/f もそうですね。

この過程で、頻繁に遭遇するのがメモリーオーバー。
48GB 積んでいますが、簡単にオーバーします。GPUならなおさら。前処理としてスライスしたり、pooling したりするのですが、ダメですね。ビッグデータどころか1時間のデータも扱えません。

簡単なのは物量作戦。メモリが足りなければ積めば良い、です。が、ある程度で頭打ちになるでしょう。何とか、効率よい処理方法を考えないと。FFT はできそうですが、PCAは思いつきません。
うーん、どうしましょう。

1/f ゆらぎ

動画を時間軸で FFT。

予想通り、スペクトルにすると低周波が高く、高周波が低くなる形状です。微動では見慣れない形状ですので面白いなあと思いつつ、当初の着目点に戻って整理を進めました。

昨夜、文献を眺めていますと、この形状についての説明がありました。

寺西ほか「動画像の時系列周波数特性によるゆらぎ解析」法政大学情報メディア教育研究センター研究報告 19, 149-153, 2006

周波数の低下とともにパワースペクトラムが増加 するような信号の中で、パワースペクトラムの振幅が周波数に対して反比例する信号を「１/ｆゆらぎ」と呼ぶ。

これでしたか！と目が覚めた気分。
昔、「１/ｆゆらぎ」が商品名についていたのを覚えています。が、意味も違いも分かっていませんでした（現在、そのような商品はありませんので、多くの方が違いを感じなかったのでしょう）。ここにきて、周波数との逆相関=「1/f」だと、ようやく理解しました。

これ、自然に多くあるとのことですが、本当でしょうか？
少なくとも、微動にはありません。
雨音は単位時間であれば高周波が卓越しているように思われます。が、音の場合、ランニングスペクトルで評価するのではなく、単位時間を代表する1つの数値を時間軸方向に並べ、そのスペクトルを取る必要があるのでしょう（時間軸での「ゆらぎ」ですから）。そうすると「1/f」に近づくのかもしれません。
画像だと、グレースケール変換、畳み込み、固有ベクトルなどが代表値として考えられます。風は風圧で良いでしょう。が、音は音圧ではダメでしょう（音色の方が重要な気がします）。知りたいですね。そのあたりは昔研究されているでしょう。探せばいくらでも出てきそうです。

まだまだ面白い会遇があるものです。
色々手を出してみるものですね。

numpy での処理

python で3次元配列を扱うことが立て続けにありました。

これらを空間方向、時間軸方向に扱う処理方法が求められます。が、最初は効率的な方法がわかりませんでした（そもそも、3次元配列に整形する段階で躓きました）。

意外と web 上に情報がなかったので、最終的に採用した numpy での処理方法を書き残しておきます。
******************************************

In [1]:

import numpy as np

In [2]:

x1=np.array([[1,1,1],[2,2,2],[3,3,3]])
x2=np.array([[11,11,11],[22,22,22],[33,33,33]])
x3=np.array([[111,111,111],[222,222,222],[333,333,333]])

In [3]:

x=np.array([x1,x2,x3])
np.shape(x)

Out[3]:

(3, 3, 3)

In [4]:

x

Out[4]:

array([[[  1,   1,   1],
        [  2,   2,   2],
        [  3,   3,   3]],

       [[ 11,  11,  11],
        [ 22,  22,  22],
        [ 33,  33,  33]],

       [[111, 111, 111],
        [222, 222, 222],
        [333, 333, 333]]])

In [5]:

x[:,0,0]

Out[5]:

array([  1,  11, 111])

In [6]:

xx=[]
xx.append(x1.copy())
xx.append(x2.copy())
xx.append(x3.copy())
xx=np.array(xx)
np.shape(xx)

Out[6]:

(3, 3, 3)

In [7]:

xx

Out[7]:

array([[[  1,   1,   1],
        [  2,   2,   2],
        [  3,   3,   3]],

       [[ 11,  11,  11],
        [ 22,  22,  22],
        [ 33,  33,  33]],

       [[111, 111, 111],
        [222, 222, 222],
        [333, 333, 333]]])

In [8]:

xx[:,0,0]

Out[8]:

array([  1,  11, 111])

In [9]:

np.mean(xx,axis=0)

Out[9]:

array([[ 41.,  41.,  41.],
       [ 82.,  82.,  82.],
       [123., 123., 123.]])

In [10]:

np.mean(xx,axis=1)

Out[10]:

array([[  2.,   2.,   2.],
       [ 22.,  22.,  22.],
       [222., 222., 222.]])

In [11]:

np.mean(xx,axis=2)

Out[11]:

array([[  1.,   2.,   3.],
       [ 11.,  22.,  33.],
       [111., 222., 333.]])

In [12]:

np.mean(xx,axis=(1,2))

Out[12]:

array([  2.,  22., 222.])

In [13]:

np.mean(xx,axis=(0,1))

Out[13]:

array([82., 82., 82.])

衛星データ

オープンな衛星データを扱えるサイト。AWS でも扱っているのは知りませんでした。

https://www.sentinel-hub.com/
https://earthengine.google.com/
https://registry.opendata.aws/

Tellus でも扱っていますが、日本の首都圏主体です。せめてALOS 程度は全国整備してほしいところです。

衛星データと防災

SAR分析チャレンジ
https://sardac.jp/

4月から5月にかけて開催されていた、情報通信研究機構 機構の啓蒙活動。基本的なコードを GitHub 上で配布し、それを組み合わせる＋新たに組むことで、SAR データと防災に関するアイデアをリンクさせるハンズオンだったようです。
公開された成果を見ると、皆さん面白いアイデアを出されていました。小さな自治体を対象としたアイデア主体で、これまで緊急観測の視点から広域利用を考えていた（そして宇宙村から外れている限り現実的でないと距離を取っていた）私にとって、新たな視点からの、驚かされる成果が多くありました。これからの利活用が期待されるところです。

SARの活用は古くから検討されています。が、SIP 等で提案されているモニタリング関連の技術レベルは数年前から変化なし。ソフトを買えば容易に達成できるレベルですが、高額なためか流行っていません。
コンペ以外にも、上記のような啓蒙活動を通しユーザーを増やすことは、技術レベルを進める一歩かもしれません。
https://www.jst.go.jp/sip/k07_kadai_01.html
http://pixel.eri.u-tokyo.ac.jp/sarws2013/sarws2013.html

一方、ArcGIS のブログでは、可視画像を使った土砂災害発生個所の抽出例が載っています。これは航空写真だけでなく、衛星の光学画像でも使えるでしょう。LANDSAT での抽出は昔から知られていますので、データ配布だけでなく結果の見やすいサイトもできています。
https://blog.esrij.com/2018/09/20/post-31278/
https://www.jstage.jst.go.jp/article/journalofmmij/132/6/132_96/_article/-char/ja/
https://apps.sentinel-hub.com/sentinel-playground/　もしくは　https://www.sentinel-hub.com/explore/eobrowser

衛星データやDEMなど、身近に使えるデータはたくさんあります。'おいしい'データは制限がかかっていますが、それでもアイデア、プログラミングの腕、やる気があれば、まだまだビジネスチャンスの埋もれている可能性があります。「SAR分析チャレンジ」によって、あらためて認識させられました。

k-空間

東京電機大学出版会「リモートセンシングのための合成開口レーダーの基礎」p325より

g(x) のフーリエ変換 G(f) の f
1/距離：空間周波数（spatial frequency）
1/時間：周波数（frequency）

空間周波数 fの代わりに k=2πf がしばしば用いられる。
k-空間（またはフーリエ空間（Fourier domain/space)）：k=2πf
kは電磁波等の波数(k=2π/λ)と混同しないように。

先日の重力探査の文献では「kは波数」等と明記されていました。一方、Fourier domain の wave number と表記された文献もあったように思います。もう一度読みましょうか。

波の知識

微動の整理を進める中で、必要な知識も一緒に整理しています。

有限フーリエ近似、有限フーリエ変換とその複素表現、フーリエ積分、周波数領域・時間領域でのフィルタや積分、窓関数、畳み込みなど。電気回路や画像処理にもつながっており、必要に迫られている技術です。これらを実装しながら、理解していない箇所をつぶしています（全く分からないところも残っていますが）。

お盆ですがプロが出社してましたので、いくつか質問もできました。おかげで解決できた問題もあります。ありがたい環境です。

一つ山を越えると、次の山が見えてきます。
このまま進みましょう。

ObsPy その2

ObsPyで計算した変位が、何かおかしいと悩んでいました。

他のソフトで計算した結果と大きく異なります。EXCELで計算した値ともあいません。速度も微妙に違っています。
加速度はあってたはず・・・と思い再確認したところ、ch２とch3 が異なっていました。15分のうち、ある1分間のデータだけがおかしい。
Mergeの仕方がまずかったか？と、その部分だけを読み込みましたが、おかしな値が再現されます。では係数の乗算がまずかったか？いえ、読み込み直後の生データの段階でおかしくなっていました。
処理ではなくて、ライブラリの想定するWin データではなかったということでしょう。プロ曰く、Winデータには方言があるそうですので。

地震の世界も案外狭いようで、その折々、必要なソフトを必要な方が作られているようです。いくつかのソフトや頂いたコードでは、今回欲しい結果に至るまでに2種以上の組み合わせが必要になりました。ObsPy＋Python なら、一気通貫できると思ったのですが、ダメでした。他のソフトで処理した波形やスペクトル比を TXT で読み込み、図化するだけなら対応できますが。うーん、これでは今までと一緒です。
結局、既存ソフトの機能追加をプロに依頼しました。

結果は出せませんでしたが、手を動かすことで理解は深まりました。それだけでも価値はあったと考えます。
わかっていたつもりでしたが、まだまだ初心者。あらためて分かった「わかっていないこと」について、今後身に着けていきましょう。

ObsPy

微動の結果を整理しています。

専用のソフトはあるのですが、自分でも中身を追って理解を深めようと思い、この週末に組んでいました。使ったのは最近のお気に入り、python。既にfortran のコードがあるので、それを片手に組んでみました。

まずは win データの読み込みから。
ライブラリがあるでしょう、と探してみたら、ありました。

ObsPy
https://docs.obspy.org/

微動ではなく、地震波を処理するためのフレームワークのようです。公開サーバーからデータを入手し、分析・図化することが容易にできます。進んでいますね。

obspy を使って、1分の win データを読み込むと、データが1個の stream として認識されました。時間やサンプリング周波数などのヘッダ付きデータです。各チャンネルはその下位の trace に入り、data で読みだせました。
フォルダ内の全 win データを merge で連続データに集約し、plot で図化。integral 2回で加速度を変位に。あとは matplot で図化。particle motion の出来上がりです（この機能はなぜか見当たりませんでした）。簡単。

integral の有効数字が気になるのと、窓関数の使い方がわからなかったのですが、それは今後。
いずれにしても、便利な世の中になりました。

******************************************
20190807追記
winデータの読み込みに、問題がありました。残念。