Main Content

このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

ロマプリータ地震の解析

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

地震データの読み込み

ファイル quake.mat には、1989 年 10 月 17 日に起きた Santa Cruz 山脈のロマプリータ (Loma Prieta) 地震の 200Hz のデータが含まれています。このデータは、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 つの変数があります。加速度計は、地震波の本震の振幅を記録しました。変数 nev は、断層に平行に配置され、その N 方向が Sacramento の方向を示す機器で測定された 3 つの方向成分を参照します。データは、機器の応答に対する補正はされていません。

他のベクトルと同じ長さをもつ、200Hz でサンプリングされたタイムスタンプを含む変数 Time を作成します。関数 seconds と乗算によって正しい単位を表し、Hz (s1) のサンプリング レートを得ます。この結果が、経過時間を表すのに役立つ変数 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)
ans=8×3 timetable
      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 の変数にアクセスし、データを調べます (ドット添字の詳細については、テーブル内のデータへのアクセスを参照してください)。"East-West" 振幅を選択し、それを期間の関数としてプロット (plot) します。

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

データのスケーリング

重力加速度でデータをスケーリングします。つまり、テーブル内の各変数に定数を乗じます。変数はすべて同じ型 (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

対象データの保存

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

tr = timerange(seconds(8),seconds(15));
dataOfInterest = quakeData(tr,:);
head(dataOfInterest)
ans=8×3 timetable
      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')

要約統計の計算

データに関する統計情報を表示するには、関数 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   

速度と位置の計算

衝撃波の伝播速度を特定するには、加速度を一度統合します。時間変数に沿って累積和を使用し、波面の速度を取得します。

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

以下では、3 つのすべての変数の統合を実行して速度を計算します。関数を作成し、varfun を使用してその関数を timetable 内の変数に適用すると便利です。この例では、このファイルの終わりに関数を含め、その関数に velFun という名前を付けています。

vel = varfun(@velFun,dataOfInterest);
head(vel)
ans=8×3 timetable
      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)
ans=8×3 timetable
      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')

軌跡の可視化

軌跡は成分値を使用して 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)))

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

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

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

軌跡の 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')

サポート関数

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

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

参考

| | | | | |

関連するトピック