로마 프리에타 지진(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() 샘플링 레이트를 얻을 수 있습니다. 이렇게 하면 경과 시간을 표시하는 데 유용한 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')
데이터 스케일링하기
데이터를 중력 가속도로 스케일링하거나 테이블의 각 변수에 상수를 곱합니다. 변수 유형이 모두 같기 때문에(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
관심 있는 데이터 저장하기
이 구간의 데이터로 또 다른 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')
요약 통계량 계산하기
데이터에 대한 통계 정보를 표시하려면 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')
궤적 시각화하기
성분 값을 사용하여 궤적을 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)))
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
궤적의 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
참고 항목
timetable
| head
| summary
| varfun
| duration
| seconds
| timerange