统计信号分析与处理报告

1、小组成员：杜磊，欧蒙川，员朝龙中国地质大学（武汉）统计信号分析处理实验报告小组成员: 魏彦斌 马全林 陈飞班 号: 075132 _院 系：_机电学院 专 业：_通信工程 指导教师：_侯强老师 一实验内容实验一、地震时间间隔的密度估计；实验二、地震震级与频度回归分析；实验三、地震空间分布聚类分析；实验四、地震优势深度聚类分析；二.实验要求及结果。 实验一、地震时间间隔的密度估计；读入数据后，把时间列（第一二列）转换成数据格式，然后分：3级以下，3到4级，4到5级，5到6级，6级以上等6个部分分别计算地震之间的时间差t，这个t就是地震时间间隔，而且是一个随机变量，对t这个随机变量进行密度估计。代

2、码：% clear all;close all; clc;filename = 中国地震台网(CSN)地震目录(1970-01-01至2015-09-31).xls;sheet = 1;xRange = A5:A8462;% xRange = A3:A8462;x2Range = B5:B8462;% x2Range = B3:B8462;yRange = H3:H8462;% ,x= xlsread(filename, sheet, xRange);% ,x2 = xlsread(filename, sheet, x2Range);,x= xlsread(filename, sheet, x

3、Range);ttt,x2 = xlsread(filename, sheet, x2Range);ml= xlsread(filename, sheet, yRange); %读取数据% X = x(isnan(ml);X2 = x2(isnan(ml);n = length(X2); %去掉无数据的日期和时间 for i = 1:n %将日期时间转化为数值形式 Xyy(i) = str2double(Xi(1:4); %年 Xmm(i) = str2double(Xi(6:7); %月 Xdd(i) = str2double(Xi(9:10); %日 XHH(i) = str2double

4、(X2i(1:2); %时 XMM(i) = str2double(X2i(4:5); %分 XSS(i) = str2double(X2i(7:8); %秒 endxx = datenum(Xyy,Xmm,Xdd,XHH,XMM,XSS); %将时间转化为数值形式ML = ml(isnan(ml) ); %去掉无数据项a=1; b=1;c=1;d=1; e=1;for i=1:n if ML(i)3.0&ML(i)4.0&ML(i)5.0&ML(i)=6.0 t_56(d)=xx(i); d=d+1; else t_6(e)=xx(i); e=e+1; end;end;%求个部分时间差for

5、 i=1:(length(t_3)-1) tt_3(i)=t_3(i)-t_3(i+1); endfor i=1:(length(t_34)-1) tt_34(i)=t_34(i)-t_34(i+1); endfor i=1:(length(t_45)-1) tt_45(i)=t_45(i)-t_45(i+1); endfor i=1:(length(t_56)-1) tt_56(i)=t_56(i)-t_56(i+1); endfor i=1:(length(t_6)-1) tt_6(i)=t_6(i)-t_6(i+1); endx=linspace(min(tt_3)-1,max(tt_3

6、),1200);p=Parzen(tt_3,x,15,);plot(x,p);grid on;figure(2);x=linspace(min(tt_34),max(tt_34),50);p=Parzen(tt_34,x,3,);plot(x,p);grid on;figure(3);x=linspace(min(tt_45),max(tt_45),100);p=Parzen(tt_45,x,3,);plot(x,p);grid on;figure(4);x=linspace(min(tt_56),max(tt_56),700);p=Parzen(tt_56,x,15,);plot(x,p);

7、grid on;figure(5);x=linspace(min(tt_6),max(tt_6),1600);p=Parzen(tt_6,x,25,);plot(x,p);grid on;实验二、地震震级与频度回归分析；读入数据后，对震级ML列进行分级统计，也是分为3级以下，3到4级，4到5级，5到6级，6级以上等6个部分，每个部分统计一下个数，然后进行线性回归拟合。拟合的公式是：LnN=a-bM式中N表示相应震级部分的个数， 比如3到4级的地震个数，M表示相关震级，比如3到4级就是4级。代码：clear all;close all; clc;graph,time=xlsread(中国地震台网

8、(CSN)地震目录(1970-01-01至2015-09-31).xls);%weidu=graph(:,2);%shendu=graph(:,3);%jingdu=graph(:,1);ML=graph(:,6);%time_day=time(:,1);%time_hour=time(:,2);n = length(ML);%figure(1);less_3=0;less_3_4=0; less_4_5=0;less_5_6=0;less_6=0;for i=1:n i if ML(i)3.0&ML(i)4.0&ML(i)5.0&ML(i)=6.0 less_5_6=less_5_6+1;

9、else less_6=less_6+1; end; end;N=log(less_3),log(less_3_4),log(less_4_5),log(less_5_6),log(less_6);M=3,4,5,6,7;plot(M,N,mo);lsline;xlabel(M);ylabel(LogN);实验三、地震空间分布聚类分析读入数据后，对经度列、纬度列和ML列分别取出来，不分震级大小，画出以经度纬度为横纵坐标，以震级为点的散点图，然后可以用我们学过的k-means进行聚类啦，通常k=6.代码：clear all;close all; clc;graph,time=xlsread(中国

10、地震台网(CSN)地震目录(1970-01-01至2015-09-31).xls);weidu=graph(:,2);%shendu=graph(:,3);jingdu=graph(:,1);ML=graph(:,6);%time_day=time(:,1);%time_hour=time(:,2);%n = length(ML);figure(1);plot3(jingdu,weidu,ML,.);grid on;xlabel(纬度);ylabel(经度);zlabel(震级);data=jingdu,weidu,ML;Idx,C,sumD,D=kmeans(data,6);figure(2

11、);plot3(data(Idx=1,1),data(Idx=1,2),data(Idx=1,3),b.,MarkerSize,5)hold onplot3(data(Idx=2,1),data(Idx=2,2),data(Idx=2,3),r*,MarkerSize,5)hold onplot3(data(Idx=3,1),data(Idx=3,2),data(Idx=3,3),gx,MarkerSize,5)hold onplot3(data(Idx=4,1),data(Idx=4,2),data(Idx=4,3),mo,MarkerSize,5)hold onplot3(data(Idx

12、=5,1),data(Idx=5,2),data(Idx=5,3),kp,MarkerSize,5)hold onplot3(data(Idx=6,1),data(Idx=6,2),data(Idx=6,3),y+,MarkerSize,5)hold on实验四、地震优势深度聚类分析；读入数据后，对深度列和ML列分别取出来，不分震级大小，画出深度-震级的散点图，然后可以用我们学过的k-means进行聚类啦，通常k=3.代码：clear all;close all; clc;graph,time=xlsread(中国地震台网(CSN)地震目录(1970-01-01至2015-09-31).xls

13、);%weidu=graph(:,2);shendu=graph(:,3);%jingdu=graph(:,1);ML=graph(:,6);%time_day=time(:,1);%time_hour=time(:,2);%n = length(ML);figure(1);scatter(ML,shendu,3,r);xlabel(震级);ylabel(深度); X=ML,shendu;K=3;Idx,C,sumD,D=kmeans(X,K);%X N*P的数据矩阵%K 表示将X划分为几类，为整数%Idx N*1的向量，存储的是每个点的聚类标号%C K*P的矩阵，存储的是K个聚类质心位置%sumD 1*K的和向量，存储的是类间所有点与该类质心点距离之和%D N*K的矩阵，存储的是每个点与所有质心的距离%画出聚类为1的点。X(Idx=1,1),为第一类的样本的第一个坐标；X(Idx=1,2)为第一类的样本的第二个坐标figure(2);plot(X(Idx=1,1),X(Idx=1,2),r.,MarkerSize,14)hold onplot(X(Idx=2,1),X(Idx=2,2),b.,MarkerSize,14)hold onplot(X(Idx=3,1),X(Idx=3,2),g.,MarkerSize,14) xlabel(震级);ylabel(深