로마 프리에타 지진(Loma Prieta Earthquake) 분석

이 예제에서는 지진 데이터를 분석하고 시각화하는 방법을 보여줍니다.

지진 데이터 불러오기

quake.mat 파일에는 1989년 10월 17일에 산타 크루즈 마운틴에서 발생한 로마 프리에타 지진(Loma Prieta Earthquake)의 200Hz 데이터가 포함되어 있습니다. 이 데이터는 산타 크루즈에 소재한 캘리포니아 대학교 찰스 F. 리히터 지진 연구소(Charles F. Richter Seismological Laboratory)에서 근무하는 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

캘리포니아 대학교 산타 크루즈 캠퍼스 자연 과학 건물에 있던 가속도계의 시간 추적 데이터가 작업 공간에 세 개의 변수로 나타납니다. 이 가속도계에는 지진파의 본진 진폭이 기록되어 있습니다. 변수 n, e, v는 단층(지각 변동으로 지층이 갈라진 지형)에 평행하게 정렬되었던 이 가속도계(N 방향은 새크라멘토(Sacramento) 방향을 가리킴)가 측정한 세 방향 각각의 구성요소를 나타냅니다. 이 데이터는 교정되지 않은 가속도계 반응 그대로입니다.

다른 벡터와 길이가 동일하고 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

데이터를 타임테이블로 정리하기

편의를 위해 개별 변수들을 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')

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

데이터 스케일링하기

데이터를 중력 가속도로 스케일링하거나 테이블의 각 변수에 상수를 곱합니다. 변수 유형이 모두 같기 때문에(double), 차원 이름 Variables를 사용하여 모든 변수에 액세스할 수 있습니다. 참고로, quakeData.Variables를 사용하면 타임테이블 내에서 숫자형 값을 직접 수정할 수 있습니다.

quakeData.Variables = 0.098*quakeData.Variables;

탐색할 데이터 서브셋 선택하기

충격파의 진폭이 거의 0에서 최대 수준으로 증가하기 시작하는 시간 영역에 대해 알아보려고 합니다. 위 플롯을 시각적으로 검토해보면 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. The axes with title East-West Acceleration contains 3 objects of type line.

관심 있는 데이터 저장하기

이 구간의 데이터로 또 다른 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')

Figure contains 3 axes. Axes 1 with title Acceleration contains an object of type line. Axes 2 contains an object of type line. Axes 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을 사용하여 계산할 수 있습니다. 이 함수는 테이블 또는 타임테이블의 각 변수에 함수를 적용하는 데 유용합니다. 적용할 함수는 varfun에 함수 핸들로 전달됩니다. 시간적 평균을 계산한 후 시간은 의미가 없기 때문에 아래에서는 mean 함수를 세 변수에 모두 적용하고 결과를 테이블 형식으로 출력합니다.

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);

아래에서는 세 변수 모두에 대해 적분을 수행하여 속도를 계산합니다. 함수를 생성한 후 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으로 생성한 타임테이블의 변수 이름에는 함수의 이름이 그대로 포함됩니다. 이러한 방식은 원래 데이터에 수행된 연산을 추적하는 데 유용합니다. 점 표기법을 사용하여 변수 이름을 원래 값으로 다시 조정할 수 있습니다.

pos.Properties.VariableNames = varNames;

아래에서는 관심 있는 시간 구간에 대해 세 구성요소의 속도와 위치를 플로팅합니다.

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. Axes 1 with title Velocity contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical. Axes 2 with title Position contains 3 objects of type line. These objects represent EastWest, NorthSouth, Vertical.

궤적 시각화하기

구성요소 값을 사용하여 궤적을 2차원 또는 3차원으로 플로팅할 수 있습니다. 다음 예제에서는 이 데이터를 시각화하는 여러 가지 방법을 보여줍니다.

2차원 투영부터 시작하겠습니다. 다음은 Figure에 일부 시간 값을 텍스트로 표시한 첫 번째 보기입니다.

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. The axes 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. The axes with title Position contains 2 objects of type line.

지원 함수

함수는 아래와 같이 정의됩니다.

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

참고 항목