이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.
로마 프리에타 지진(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)
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
를 사용하면 타임테이블 내에서 숫자형 값을 직접 수정할 수 있습니다.
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)
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
을 사용하여 계산할 수 있습니다. 이 함수는 테이블 또는 타임테이블의 각 변수에 함수를 적용하는 데 유용합니다. 적용할 함수는 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')
궤적 시각화하기
구성요소 값을 사용하여 궤적을 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