统计信号分析与处理报告

XCXDFSEWRV中国地质大学（武汉）统计信号分析处理报告小组成员: 魏彦斌马全林陈飞班号: 075132 _院系：__机电学院专业：_通信工程指导教师：_侯强老师一实验内容实验一、地震时间间隔的密度估计；实验二、地震震级与频度回归分析；实验三、地震空间分布聚类分析；实验四、地震优势深度聚类分析；二.实验要求及结果。

实验一、地震时间间隔的密度估计；读入数据后，把时间列（第一二列）转换成数据格式，然后分：3级以下，3到4级，4到5级，5到6级，6级以上等6个部分分别计算地震之间的时间差t，这个t就是地震时间间隔，而且是一个随机变量，对t这个随机变量进行密度估计。

代码：%%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, xRange);[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(X{i}(1:4)); %年Xmm(i) = str2double(X{i}(6:7)); %月Xdd(i) = str2double(X{i}(9:10)); %日XHH(i) = str2double(X2{i}(1:2)); %时XMM(i) = str2double(X2{i}(4:5)); %分XSS(i) = str2double(X2{i}(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:nif ML(i)<=3.0t_3(a)=xx(i);a=a+1;elseif ML(i)>3.0&&ML(i)<=4.0t_34(b)=xx(i);b=b+1;elseif ML(i)>4.0&&ML(i)<=5.0t_45(c)=xx(i);c=c+1;elseif ML(i)>5.0&&ML(i)<=6.0t_56(d)=xx(i);d=d+1;elset_6(e)=xx(i);e=e+1;end;end;%求个部分时间差for 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),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);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('中国地震台网(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:niif ML(i)<=3.0less_3=less_3+1;elseif ML(i)>3.0&&ML(i)<=4.0less_3_4=less_3_4+1;elseif ML(i)>4.0&&ML(i)<=5.0less_4_5=less_4_5+1;elseif ML(i)>5.0&&ML(i)<=6.0less_5_6=less_5_6+1;elseless_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('中国地震台网(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);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==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');%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的点。