SCILAB ~内部基盤面と解放基盤面~

引き続き地盤の周波数応答関数についてです。

f:id:S_E_Hyphen:20180304155155j:plain

第1層を対象層、第6層を基準層とする増幅スペクトルを計算します。すなわち地表面を分子、基盤上面を分母とする加速度 / 加速度応答関数です。

地表面に観測点があって、図の1番上の地震波形が記録されたとします。この観測波形をフーリェ変換した複素フーリェ級数xを、先に求めた増幅スペクトルz1で「割り算」します。そして更に逆フーリェ変換で時間領域に戻したら、基盤上面の揺れ方y1を推定することが可能となります(図の第2段)。

さらに、もしも第1層から第5層が無いものとして基盤上面が剝き出しの場合の揺れ方y2も推定可能で、このとき揺れは約2倍に大きくなります(図の第3段)。

図の第2段を内部基盤面、第3段を解放基盤面と呼んでいます。

 

 スクリプトについては多少の補足説明が必要で、9行目のexec('acc_resp.sce');は SCILABで周波数応答関数を計算する で示した関数acc_respを記述した外部ファイルを参照しています。

同様に10行目のexec('material_set.sce'); も外部ファイルの参照で、ここでは SCILAB ~ひずみ/加速度応答~ の51行目から56行目で示したエルセントロ地域の地盤情報(単位体積重量や剛性率など)の読み込みを行っています。

また 原著 では、内部基盤波を推定する際には観測点のある第1層を基準層、推定対象である第6層を対象層と考えて増幅スペクトルを掛け算しています。一方、解放基盤波の推定にあたっては基準層と対象層をひっくり返したうえで増幅スペクトルの逆数(1/z)を掛け算するという説明としています。

ここでは「逆増幅は割り算を使用する」ものとして、いずれも基盤を基準層、地表面を対象層に統一した説明としました。解放基盤面といった仮想の基盤面を例外と捉えるならば、原著のように「掛け算」を基本に考える方が分かり易いかも知れません。

#!/bin/bash
temp=`mktemp ./XXXX`
cat エルセントロ.txt |\
 sharp_cancel |\
 awk 'NR<=800{print $1,-$2}' > $temp

scilab -nw << _SCI_

exec('acc_resp.sce'); // function acc_resp の読み込み
exec('material_set.sce'); // el centro 地盤情報の読み込み

//エルセントロ波の読み込み
fd=mopen('$temp','r');
[n,t,a0]=mfscanf(-1,fd,'%f %f');
mclose(fd);
ns0=size(a0,'r'); // ns0=800
fs=100; // サンプリング周波数
// 後続のゼロ
a=[a0;zeros(1024-ns0,1)];
ns=size(a,'r'); // ns=1024(2の10乗)
nspe=ns/2+1; // スペクトルの要素数
// ここまででエルセントロ波の読み込み

freq0=[0:fs/ns:fs]; // 周波数軸の設定
x=fft(a,-1); // FFTの実行

obj=1; // 対象層(地表面)
ref=6; // 基準層(基盤上面) 
for i=1:nspe // DCからナイキストまで if i==1 then z1=1; z2=1; else f=freq0(i); // 周波数 z1=acc_resp(f,'inner',obj,ref); z2=acc_resp(f,'expose',obj,ref); end // if文の終わり st1(i)=x(i)/z1; // 内部基盤波の複素フーリェ級数 st2(i)=x(i)/z2; // 解放基盤波の複素フーリェ級数 end // 周波数ループの終わり st1=[st1;conj(flipdim(st1(2:nspe-1),1))]; // 内部基盤波の共役作成 st2=[st2;conj(flipdim(st2(2:nspe-1),1))]; // 解放基盤波の共役作成 y1=fft(st1,+1); // 時間領域に戻す y2=fft(st2,+1); // 時間領域に戻す fd=mopen('$temp','w'); for i=1:ns0 mfprintf(fd,'%f %f %f %f\n',t(i),a0(i),real(y1(i)),real(y2(i))); end mclose(fd); _SCI_ gmt gmtset FONT_TITLE 12p,GothicBBB-Medium-UniJIS-UTF8-H gmt gmtset FONT_LABEL 12p,GothicBBB-Medium-UniJIS-UTF8-H gmt gmtset FONT_ANNOT_PRIMARY 10p,GothicBBB-Medium-UniJIS-UTF8-H cat $temp |\ awk '{print $1,$4}' |\ gmt psxy -JX27c/6c -R0/8.1/-600/600 -BnWeS+t"解放基盤面" \ -Bxa1+l"時間(秒)" -Bya600f100g600 -W0.5p -X1.5c -Y1.5c -K > wave.ps cat $temp |\ awk '{print $1,$3}' |\ gmt psxy -JX27c/3c -R0/8.1/-300/300 -BnWeS+t"内部基盤面" \ -Bxaf -Bya300fg500 -W0.5p -Y8.2c -O -K>> wave.ps cat $temp |\ awk '{print $1,$2}' |\ gmt psxy -JX27c/4c -R0/8.1/-400/400 -BnWeS+t"地表面" \ -Bxaf -Bya400fg500 -W0.5p -Y5.2c -O >> wave.ps rm $temp