Main Content

로마 프리에타 지진(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(s1) 샘플링 레이트를 얻을 수 있습니다. 이렇게 하면 경과 시간을 표시하는 데 유용한 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)
      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 object. The axes object 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 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, ylabel East-West contains an object of type line. Axes object 2 with ylabel North-South contains an object of type line. Axes object 3 with ylabel Vertical 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)
      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으로 생성한 타임테이블의 변수 이름에는 함수의 이름이 그대로 포함됩니다. 이러한 방식은 원래 데이터에 수행된 연산을 추적하는 데 유용합니다. 점 표기법을 사용하여 변수 이름을 원래 값으로 다시 조정할 수 있습니다.

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 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차원 투영부터 시작하겠습니다. 다음은 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 object. The axes object with xlabel North-South, ylabel Vertical contains 5 objects of type line, text.

plotmatrix를 사용하여, 모든 변수를 서로 비교한 산점도 플롯을 격자(그리드)로 그리면서 대각선에는 각 변수의 히스토그램이 표시되도록 시각화합니다. 그리드의 각 좌표축을 나타내는 출력 변수 Axxlabelylabel을 사용하여 레이블을 지정할 좌표축을 식별하는 데 유용합니다.

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, xlabel North-South, ylabel East-West contains 2 objects of type line. One or more of the lines displays its values using only markers

지원 함수

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

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

참고 항목

| | | | | |

관련 항목