ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

ロマプリータ地震のデモ

この例では、実世界の地震データをどのように解析し、可視化するかを示します。

Authors:C. Denham, 1990, C. Moler, August, 1992.

ファイル QUAKE.MAT には、1989 年 10 月 17 日に起きた Santa Cruz 山脈のロマプリータ (Loma Prieta) 地震の 200Hz のデータが含まれています。このデータは、Santa Cruz の California 大学、Charles F. Richter 地震学研究所の Joel Yellin 氏から提供されています。データの読み込みから始めます。

load quake e n v
whos e n v
  Name          Size            Bytes  Class     Attributes

  e         10001x1             80008  double              
  n         10001x1             80008  double              
  v         10001x1             80008  double              

今、ワークスペースには、UC Santa Cruz の自然科学棟にある加速度計で記録された時間波形を含む 3 つの変数があります。加速度計は、地震の本震を記録しました。変数 n、e、v は、断層に平行に配置され、その N 方向が Sacramento の方向を示す機器で測定された 3 つの方向成分を参照します。データは、機器の応答に対する補正はされていません。

重力加速度によってデータをスケーリングします。さらに、時間基準を含む 4 つ目の変数 t を作成します。

g = 0.0980;
e = g*e;
n = g*n;
v = g*v;
delt = 1/200;
t = delt*(1:length(e))';

加速度をプロットします。

yrange = [-250 250];
limits = [0 50 yrange];
subplot(3,1,1), plot(t,e,'b'), axis(limits), title('East-West acceleration')
subplot(3,1,2), plot(t,n,'g'), axis(limits), title('North-South acceleration')
subplot(3,1,3), plot(t,v,'r'), axis(limits), title('Vertical acceleration')

t=8 秒から t=15 秒までの区間を見てみましょう。選択された場所は黒い線で描画します。後の計算はすべてこの区間を使用します。

t1 = 8*[1;1];
t2 = 15*[1;1];
subplot(3,1,1), hold on, plot([t1 t2],yrange,'k','LineWidth',2); hold off
subplot(3,1,2), hold on, plot([t1 t2],yrange,'k','LineWidth',2); hold off
subplot(3,1,3), hold on, plot([t1 t2],yrange,'k','LineWidth',2); hold off

選択した時間区間を拡大します。

trange = sort([t1(1) t2(1)]);
k = find((trange(1)<=t) & (t<=trange(2)));
e = e(k);
n = n(k);
v = v(k);
t = t(k);
ax = [trange yrange];

subplot(3,1,1), plot(t,e,'b'), axis(ax), title('East-West acceleration')
subplot(3,1,2), plot(t,n,'g'), axis(ax), title('North-South acceleration')
subplot(3,1,3), plot(t,v,'r'), axis(ax), title('Vertical acceleration')

この区間の真中の 1 秒にフォーカスすると、"North-South" に対する "East-West" のプロットは水平の加速度を表示します。

subplot(1,1,1)
k = length(t);
k = round(max(1,k/2-100):min(k,k/2+100));
plot(e(k),n(k),'.-')
xlabel('East'), ylabel('North');
title('Acceleration During a One Second Period');

3 次元空間内の点の速度と位置を計算するために加速度を 2 回積分します。

edot = cumsum(e)*delt;  edot = edot - mean(edot);
ndot = cumsum(n)*delt;  ndot = ndot - mean(ndot);
vdot = cumsum(v)*delt;  vdot = vdot - mean(vdot);

epos = cumsum(edot)*delt;  epos = epos - mean(epos);
npos = cumsum(ndot)*delt;  npos = npos - mean(npos);
vpos = cumsum(vdot)*delt;  vpos = vpos - mean(vpos);

subplot(2,1,1);
plot(t,[edot+25 ndot vdot-25]); axis([trange min(vdot-30) max(edot+30)])
xlabel('Time'), ylabel('V - N - E'), title('Velocity')

subplot(2,1,2);
plot(t,[epos+50 npos vpos-50]);
axis([trange min(vpos-55) max(epos+55)])
xlabel('Time'), ylabel('V - N - E'), title('Position')

位置データで定義された軌跡は、3 つの異なる 2 次元の写像で表示することができます。これが 1 つ目で、t のいくつかの値に注釈が付けられています。

subplot(1,1,1); cla; subplot(2,2,1)
plot(npos,vpos,'b');
na = max(abs(npos)); na = 1.05*[-na na];
ea = max(abs(epos)); ea = 1.05*[-ea ea];
va = max(abs(vpos)); va = 1.05*[-va va];
axis([na va]); xlabel('North'); ylabel('Vertical');

nt = ceil((max(t)-min(t))/6);
k = find(fix(t/nt)==(t/nt))';
for j = k, text(npos(j),vpos(j),['o ' int2str(t(j))]); end

同様のコードで、さらに 2 つの 2 次元表示を生成します。

subplot(2,2,2)
plot(epos,vpos,'g');
for j = k; text(epos(j),vpos(j),['o ' int2str(t(j))]); end
axis([ea va]); xlabel('East'); ylabel('Vertical');

subplot(2,2,3)
plot(npos,epos,'r');
for j = k; text(npos(j),epos(j),['o ' int2str(t(j))]); end
axis([na ea]); xlabel('North'); ylabel('East');

4 番目のサブプロットは、軌跡の 3 次元表示です。

subplot(2,2,4)
plot3(npos,epos,vpos,'k')
for j = k;text(npos(j),epos(j),vpos(j),['o ' int2str(t(j))]); end
axis([na ea va]); xlabel('North'); ylabel('East'), zlabel('Vertical');
box on

最後に、10 位置点ごとにドットをプロットします。ドット間の間隔は、速度を表します。

subplot(1,1,1)
plot3(npos,epos,vpos,'r')
hold on
step = 10;
plot3(npos(1:step:end),epos(1:step:end),vpos(1:step:end),'.')
hold off
box on
axis tight
xlabel('North-South')
ylabel('East-West')
zlabel('Vertical')
title('Position (cms)')

この情報は役に立ちましたか?