山野莽夫

  • 归档
    • 随笔
    • 建站资源
    • 分享
    • 代码
  • 地球物理学
    • 专业课
    • 概念解释
  • 计算机
  • 互联网
  • 教程
  • 规划
  • 实验室
    • 珍藏的软件
    • 贴吧云签到
    • A1账号自助申请
山野莽夫
小学生的挣扎的点点滴滴
  1. 首页
  2. 地球物理学
  3. 正文

工程物探实验报告-地震记录合成以及频谱分析

2016年5月30日 3531点热度 1人点赞 2条评论

实验要求

实验一(四课时): A 、地震信号频谱分析; 给定一个地震记录(随机或计算得出),编 程分析该地震记录的振幅和相位情况,绘制其振 幅谱,相位谱。 注:可用 C 语言编程实现或软件实现。振幅 谱,相位谱用 Excel 绘制。 B 、褶积合成地震反射记录 ; 选择适合自己的编程语言(C、C++等)或 软件合成地震反射记录。结果用 Excel 绘制。

实验报告

 一、实验目的

运用卷积方法合成地震记录,并对地震记录进行频谱分析,获取地震记录的振幅谱和相位谱,并进行解释。

二、实验原理

A 、地震信号频谱分析

获取地震记录的振幅谱和相位谱,进行频谱分析。

因为信号是复数,所以把实部和虚部分开来计算,并加以合成,便可以得到振幅谱和相位谱。

B 、褶积合成地震反射记录 ;

卷积方法可以得到地震反射记录。原始数据为地层的反射系数系列,通过把反射系数序列和雷克子波卷积既可以合成地震记录。

雷克子波的计算公式为。

反射系数序列为自己制定,本次实验取四层地层。共取200个点,给四个点赋值,其他均为零。

把得到的反射系数序列和计算出来的雷克子波进行卷积即可以合成地震记录。、

三、实验结果分析

结果分析:地震记录为反射系数序列和雷克子波卷积而成,应该保证反射系数序列比雷克子波序列足够的长,否者就会出现地震记录不圆滑等很多问题,而且雷克子波应该把频率选择高一些,这样周期才比较短,在一定长度的雷克子波序列比较完整。这样卷积的到的结果更加真实。

四、源代码

A、频谱分析

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define gcwt "file2.dat"

#define gcwt2 "file1.dat"

float pinpu(float*x,float *Xr,float *Xi);

int main()

{

float Xr[200],Xi[200],x[200],X[200],xiangwei[200];int i,j,k;

FILE *fp,*fp1;

fp=fopen(gcwt,"r+");

fp1=fopen(gcwt2,"w+");

for(i=0;i<200;i++)

{

fscanf(fp,"%f",&x[i]);

}//读取数据

fclose(fp);

pinpu(x,Xr,Xi);//计算实部虚部

for (i=0;i<200;i++)

{

X[i]=sqrt(Xr[i]*Xr[i]+Xi[i]*Xi[i]);//振幅谱

xiangwei[i]=atan(Xi[i]/Xr[i]);//相位谱

printf("%d\t%f\t%f\n",i,X[i],xiangwei[i]);

fprintf(fp1,"%d\t%f\t%f\n",i,X[i],xiangwei[i]);

}

fclose(fp1);

}

float pinpu(float*x,float *Xr,float *Xi)

{

int i,j;

for(i=0;i<200;i++)

{

Xr[i]=0.0;

Xi[i]=0.0;

for(j=0;j<200;j++)

{

Xr[i]=x[j]*cos(2*3.14/200*j*i)+Xr[i];//实部

Xi[i]=-x[j]*sin(2*3.14/200*j*i)+Xi[i];//虚部

}

}

return 0;

}

B、地震记录合成

#include<stdio.h>

#include<math.h>

#include<stdlib.h>

#define dz1 "file1.dat"

#define dz2 "file2.dat"

main ()

{ FILE *fp1,*fp2;

double dt=0.002,t[71];int fm1=10,i,j;double zibo[71],xishu[200],dzjl[300];

fp1=fopen(dz1,"w+");

fp2=fopen(dz2,"w+");

fclose(fp1);

fclose(fp2);

fp1=fopen(dz1,"r+");

fp2=fopen(dz2,"r+");

for(i=0;i<=35;i++)

{

t[i+35]=i*dt;

zibo[i+35]=exp(-pow(3.141592653*fm1*t[i+35],2))*(1-2*pow(3.141592653*fm1*t[i+35],2));

}

for(i=0;i<71;i++)

// printf("%f\n",zibo[i]);

for(i=0;i<=34;i++)

{

t[i]=-t[70-i];

zibo[i]=zibo[70-i];

}

for(i=0;i<=70;i++)

{printf("%f\t%f\n",t[i],zibo[i]);

fprintf(fp1," %f\t%f\n",t[i],zibo[i]);

}

for(i=0;i<=199;i++)

{

xishu[i]=0;

}

xishu[10]=0.2;xishu[65]=0.4;xishu[110]=-0.5;xishu[155]=0.7;xishu[190]=-0.36;

for(i=0;i<200;i++)

{ dzjl[i]=0.0;

for(j=0;j<70;j++)

{ if((i-j)>0)

dzjl[i]=dzjl[i]+zibo[j]*xishu[i-j] ; }

}

for(i=0;i<200;i++)

{

// printf("%d\t%f\n",i,dzjl[i]);

fprintf(fp2,"%d\t%f\n",i,dzjl[i]);

}

fclose(fp1);

fclose(fp2);

}

备注:代码可能不能脱离源码包执行,如果不能执行,请在下方留言,我把代码包发给你。当然如果有其他问题,都可以给我留言。

标签: 地震记录合成 工程物探 频谱分析
最后更新:2018年4月6日

小菜菜

菜鸟

打赏 点赞
< 上一篇
下一篇 >

文章评论

  • 爬行者

    大佬,能把代码包发一下吗?运行出错了

    2021年10月15日
    回复
    • 小菜菜

      @爬行者 年代久远,找不到了,这都是用的vc++6.0写的代码

      2022年1月26日
      回复
  • razz evil exclaim smile redface biggrin eek confused idea lol mad twisted rolleyes wink cool arrow neutral cry mrgreen drooling persevering
    取消回复

    此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据。

    标签聚合
    模板 ppt 宝塔面板 wordpress 虚拟机 c语言 地震学程序 onedrive
    最新 热点 随机
    最新 热点 随机
    Azure Student 微软云 学生订阅 免费12个月用量避坑注意点集合 MP3音频文件格式详细解析 python按固定采样点个数分割wav格式音频 愉快使用谷歌免费人工智能平台colab,训练你的神经网络模型,为你的学术生活添砖加瓦 华为云版轻量应用服务器-云耀云服务器简单体验评测 Cloudflare 免费CDN自定义节点ip之自选cloudflare 高速节点ip工具分享
    求爱网站源码 excel工作表(表格xlsx)去除工作表视图保护和工作簿保护 网盘挂载(映射)工具(支持onedrive、google drive)RaiDrive疫情期间限时使用专业版以及一点使用google网盘的说明 C语言黑窗程序执行完就退出的解决方法。 计算方法实习----牛顿前插公式 又一个功能强大的cloudflare workers onedrive 索引工具(OneDrive-Index-Cloudflare-Worker)部署教程

    COPYRIGHT © 2021 shanyemangfu.com. ALL RIGHTS RESERVED.

    Theme Kratos Made By Seaton Jiang

    蜀ICP备15031791号-2