ロマプリータ地震の解析

この例では、タイムスタンプ付きの地震データを timetable に保存する方法と、timetable 関数を使用してデータを解析し可視化する方法を説明します。

地震データの読み込み

この例のファイル quake.mat には、1989 年 10 月 17 日、Santa Cruz 山脈で発生したロマプリータ (Loma Prieta) 地震の 200 Hz のデータが含まれています。このデータは、Santa Cruz の California 大学、Charles F. Richter 地震学研究所の Joel Yellin 氏から提供されています。

データの読み込みから始めます。

load quake 
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 つの方向成分を参照します。データは、機器の応答に対する補正はされていません。

他のベクトルと同じ長さをもつ、200Hz でサンプリングされたタイムスタンプを含む変数 Time を作成します。関数 seconds と乗算によって正しい単位を表し、Hz ( $s^{- 1}$ ) のサンプリングレートを得ます。この結果が、経過時間を表すのに役立つ変数durationになります。

Time = (1/200)*seconds(1:length(e))';
whos Time

  Name          Size            Bytes  Class       Attributes

  Time      10001x1             80010  duration

timetable でのデータの編成

個々の変数をtableやtimetableに整理するとさらに便利です。timetable には、タイムスタンプ付きデータを操作するための柔軟性と機能が備わっています。時間変数および 3 つの加速度変数を使用して timetable を作成し、より意味のある変数名を指定します。関数 head を使用して最初の 8 行を表示します。

varNames = {'EastWest', 'NorthSouth', 'Vertical'};
quakeData = timetable(Time, e, n, v, 'VariableNames', varNames);
head(quakeData)

      Time       EastWest    NorthSouth    Vertical
    _________    ________    __________    ________

    0.005 sec       5            3            0    
    0.01 sec        5            3            0    
    0.015 sec       5            2            0    
    0.02 sec        5            2            0    
    0.025 sec       5            2            0    
    0.03 sec        5            2            0    
    0.035 sec       5            1            0    
    0.04 sec        5            1            0

ドット添字を使用して timetable の変数にアクセスし、データを調べます。(ドット添字の詳細については、table 内のデータへのアクセスを参照してください。)"East-West" の振幅を選択し、それを持続時間の関数としてプロット (plot) します。

plot(quakeData.Time,quakeData.EastWest)
title('East-West Acceleration')

Figure contains an axes object. The axes object with title East-West Acceleration contains an object of type line.

データのスケーリング

重力加速度でデータをスケーリングします。つまり、テーブル内の各変数に定数を乗じます。変数はすべて同じ型 (double) であるため、次元名 Variables を使用してすべての変数にアクセスできます。quakeData.Variables を使用すると timetable 内の数値を直接変更できることに注目してください。

quakeData.Variables = 0.098*quakeData.Variables;

調査対象のデータのサブセットの選択

衝撃波の振幅がゼロ近くから最大レベルまで増加し始める時間領域を調べます。上記のプロットを見ると、8 秒から 15 秒までの時間区間が該当の時間であることがわかります。より見やすくなるよう、選択した時点に黒い線を引いてこの区間を目立たせます。後のすべての計算にこの区間が関わってきます。

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

Figure contains an axes object. The axes object with title East-West Acceleration contains 3 objects of type line.

対象データの保存

この区間のデータを使用して別の timetable を作成します。対象の行を選択するには、timerangeを使用します。

tr = timerange(seconds(8),seconds(15));
dataOfInterest = quakeData(tr,:);
head(dataOfInterest)

      Time       EastWest    NorthSouth    Vertical
    _________    ________    __________    ________

    8 sec         -0.098       2.254          5.88 
    8.005 sec          0       2.254         3.332 
    8.01 sec      -2.058       2.352        -0.392 
    8.015 sec     -4.018       2.352        -4.116 
    8.02 sec      -6.076        2.45        -7.742 
    8.025 sec     -8.036       2.548       -11.466 
    8.03 sec     -10.094       2.548          -9.8 
    8.035 sec     -8.232       2.646        -8.134

3 つの加速度変数を異なる 3 つの座標軸で可視化します。

figure
subplot(3,1,1)
plot(dataOfInterest.Time,dataOfInterest.EastWest)
ylabel('East-West')
title('Acceleration')
subplot(3,1,2)
plot(dataOfInterest.Time,dataOfInterest.NorthSouth)
ylabel('North-South')
subplot(3,1,3)
plot(dataOfInterest.Time,dataOfInterest.Vertical)
ylabel('Vertical')

Figure contains 3 axes objects. Axes object 1 with title Acceleration contains an object of type line. Axes object 2 contains an object of type line. Axes object 3 contains an object of type line.

要約統計量の計算

データに関する統計情報を表示するには、関数summaryを使用します。

summary(dataOfInterest)

RowTimes:

    Time: 1400x1 duration
        Values:
            Min           8 sec      
            Median        11.498 sec 
            Max           14.995 sec 
            TimeStep      0.005 sec  

Variables:

    EastWest: 1400x1 double

        Values:

            Min       -255.09 
            Median     -0.098 
            Max        244.51 

    NorthSouth: 1400x1 double

        Values:

            Min       -198.55 
            Median      1.078 
            Max        204.33 

    Vertical: 1400x1 double

        Values:

            Min       -157.88 
            Median       0.98 
            Max        134.46

データに関する追加の統計情報は、varfunを使用して計算できます。これは、table 内または timetable 内の各変数に関数を適用する際に役立ちます。適用する関数は varfun に関数ハンドルとして渡されます。関数 mean を 3 つの変数すべてに適用し、結果を table 形式で出力します。これは、時間平均を計算した後では時間に意味がなくなるためです。

mn = varfun(@mean,dataOfInterest,'OutputFormat','table')

mn=1×3 table
    mean_EastWest    mean_NorthSouth    mean_Vertical
    _____________    _______________    _____________

       0.9338           -0.10276          -0.52542

速度と位置の計算

衝撃波の伝播速度を特定するために、加速度を 1 回積分します。時間変数に沿った累積和を使用して、波面の速度を取得します。

edot = (1/200)*cumsum(dataOfInterest.EastWest);
edot = edot - mean(edot);

次に、3 つすべての変数の積分を実行して速度を計算します。関数を作成し、varfun を使用してその関数を timetable 内の変数に適用すると便利です。この例では、その関数が例の最後に velFun という名前で含められています。

vel = varfun(@velFun,dataOfInterest);
head(vel)

      Time       velFun_EastWest    velFun_NorthSouth    velFun_Vertical
    _________    _______________    _________________    _______________

    8 sec           -0.56831             0.44642             1.8173     
    8.005 sec       -0.56831             0.45769              1.834     
    8.01 sec         -0.5786             0.46945              1.832     
    8.015 sec       -0.59869             0.48121             1.8114     
    8.02 sec        -0.62907             0.49346             1.7727     
    8.025 sec       -0.66925              0.5062             1.7154     
    8.03 sec        -0.71972             0.51894             1.6664     
    8.035 sec       -0.76088             0.53217             1.6257

同じ関数 velFun を速度に適用して位置を特定します。

pos = varfun(@velFun,vel);
head(pos)

      Time       velFun_velFun_EastWest    velFun_velFun_NorthSouth    velFun_velFun_Vertical
    _________    ______________________    ________________________    ______________________

    8 sec                2.1189                    -2.1793                    -3.0821        
    8.005 sec            2.1161                     -2.177                    -3.0729        
    8.01 sec             2.1132                    -2.1746                    -3.0638        
    8.015 sec            2.1102                    -2.1722                    -3.0547        
    8.02 sec              2.107                    -2.1698                    -3.0458        
    8.025 sec            2.1037                    -2.1672                    -3.0373        
    8.03 sec             2.1001                    -2.1646                    -3.0289        
    8.035 sec            2.0963                     -2.162                    -3.0208

varfun で作成された timetable の変数名に、使用した関数の名前が含まれているのがわかります。これは元のデータで実行された操作を追跡するのに役立ちます。ドット表記を使用して変数名を元の値に合わせます。

pos.Properties.VariableNames = varNames;

対象とする時間区間における速度と位置の 3 成分をプロットします。

figure
subplot(2,1,1)
plot(vel.Time,vel.Variables)
legend(quakeData.Properties.VariableNames,'Location','NorthWest')
title('Velocity')
subplot(2,1,2)
plot(vel.Time,pos.Variables)
legend(quakeData.Properties.VariableNames,'Location','NorthWest')
title('Position')

Figure contains 2 axes objects. Axes object 1 with title Velocity contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical. Axes object 2 with title Position contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical.

軌跡の可視化

軌跡は成分値を使用して 2 次元または 3 次元でプロットできます。このプロットは、このデータを可視化する異なった方法を示しています。

最初は 2 次元の写像です。これは 1 つ目で、いくつかの時間値の注釈が付けられています。

figure
plot(pos.NorthSouth,pos.Vertical)
xlabel('North-South')
ylabel('Vertical')
% Select locations and label
nt = ceil((max(pos.Time) - min(pos.Time))/6);
idx = find(fix(pos.Time/nt) == (pos.Time/nt))';
text(pos.NorthSouth(idx),pos.Vertical(idx),char(pos.Time(idx)))

Figure contains an axes object. The axes object contains 5 objects of type line, text.

plotmatrixを使用して、すべての変数の散布図を並べ、各変数のヒストグラムを対角線上に配したグリッドを可視化します。出力変数 Ax はグリッドの各座標軸を表し、xlabel と ylabel を使用してラベル付けする座標軸を識別するために使用できます。

figure
[S,Ax] = plotmatrix(pos.Variables);

for ii = 1:length(varNames)
    xlabel(Ax(end,ii),varNames{ii})
    ylabel(Ax(ii,1),varNames{ii})
end

MATLAB figure

軌跡の 3 次元表示をプロットし、10 位置点ごとにドットをプロットします。ドット間の間隔は、速度を表します。

step = 10;
figure
plot3(pos.NorthSouth,pos.EastWest,pos.Vertical,'r')
hold on
plot3(pos.NorthSouth(1:step:end),pos.EastWest(1:step:end),pos.Vertical(1:step:end),'.')
hold off
box on
axis tight
xlabel('North-South')
ylabel('East-West')
zlabel('Vertical')
title('Position')

Figure contains an axes object. The axes object with title Position contains 2 objects of type line.

サポート関数

関数は以下のように定義されます。

function y = velFun(x)
    y = (1/200)*cumsum(x);
    y = y - mean(y);
end

参考