统计信号分析与处理报告

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)); %秒

end

xx = 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

t_3(a)=xx(i);

a=a+1;

elseif ML(i)>3.0&&ML(i)<=4.0

t_34(b)=xx(i);

b=b+1;

elseif ML(i)>4.0&&ML(i)<=5.0

t_45(c)=xx(i);

c=c+1;

elseif 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 i=1:(length(t_3)-1)

tt_3(i)=t_3(i)-t_3(i+1);

end

for i=1:(length(t_34)-1)

tt_34(i)=t_34(i)-t_34(i+1);

end

for i=1:(length(t_45)-1)

tt_45(i)=t_45(i)-t_45(i+1);

end

for i=1:(length(t_56)-1)

tt_56(i)=t_56(i)-t_56(i+1);

end

for i=1:(length(t_6)-1)

tt_6(i)=t_6(i)-t_6(i+1);

end

x=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:n

i

if ML(i)<=3.0

less_3=less_3+1;

elseif ML(i)>3.0&&ML(i)<=4.0

less_3_4=less_3_4+1;

elseif ML(i)>4.0&&ML(i)<=5.0

less_4_5=less_4_5+1;

elseif ML(i)>5.0&&ML(i)<=6.0

less_5_6=less_5_6+1;

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('中国地震台网(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 on

plot3(data(Idx==2,1),data(Idx==2,2),data(Idx==2,3),'r*','MarkerSize',5)

hold on

plot3(data(Idx==3,1),data(Idx==3,2),data(Idx==3,3),'gx','MarkerSize',5)

hold on

plot3(data(Idx==4,1),data(Idx==4,2),data(Idx==4,3),'mo','MarkerSize',5)

hold on

plot3(data(Idx==5,1),data(Idx==5,2),data(Idx==5,3),'kp','MarkerSize',5)

hold on

plot3(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的点。X(Idx==1,1),为第一类的样本的第一个坐标；X(Idx==1,2)为第一类的样本的第二个坐标

figure(2);

plot(X(Idx==1,1),X(Idx==1,2),'r.','MarkerSize',14)

hold on

plot(X(Idx==2,1),X(Idx==2,2),'b.','MarkerSize',14)

hold on

plot(X(Idx==3,1),X(Idx==3,2),'g.','MarkerSize',14)

xlabel('震级');ylabel('深度');

%绘出中心点,kx表示是圆形

plot(C(:,1),C(:,2),'kx','MarkerSize',14,'LineWidth',4)

plot(C(:,1),C(:,2),'kx','MarkerSize',14,'LineWidth',4)

plot(C(:,1),C(:,2),'kx','MarkerSize',14,'LineWidth',4)

legend('Cluster 1','Cluster 2','Cluster 3','Centroids','Location','NW')