根据零极点文件画仪器响应频谱图

SeisPiano

已于 2023-03-17 09:29:16 修改

阅读量1.9k

点赞数 3

分类专栏： Matlab 地震文章标签： matlab

于 2021-09-09 12:06:10 首次发布

本文链接：https://blog.csdn.net/weixin_45693832/article/details/120197631

版权

Matlab 同时被 2 个专栏收录

11 篇文章 0 订阅

订阅专栏

地震

6 篇文章 3 订阅

订阅专栏

仪器响应简介

地震仪观测到的地面运动记录可以表示为
$u (t) = s (t) * g (t) * i (t)$
其中 $s (t)$ 代表震源项， $g (t)$ 代表路径效应， $i (t)$ 代表仪器响应， $*$ 代表卷积。

PZ（零极点）文件

一般来说，地震仪的仪器响应文件有RESP、PZ 和 FAP。
这里以PZ文件（即零极点文件）为例，相对于 RESP 文件，PZ 文件中仅包含仪器响应中的零极点和增益信息，在去仪器响应时更方便。
下面是某个台站的 Response PZs文件：

* **********************************
ZEROS	5
	+0.000000e+00	+0.000000e+00	
	+0.000000e+00	+0.000000e+00	
	+0.000000e+00	+0.000000e+00	
	-1.344330e-02	+0.000000e+00	
	-1.344330e-02	+0.000000e+00	
POLES	6
	-8.191920e-03	+0.000000e+00	
	-1.235380e-02	-1.234000e-02	
	-1.235380e-02	+1.234000e-02	
	-1.883430e-02	+0.000000e+00	
	-3.918000e+01	-4.912000e+01	
	-3.918000e+01	+4.912000e+01	
CONSTANT	1.752545e+13

根据零极点以及 CONSTANT，即可计算得到仪器响应函数：
$H(s)=C_{0}*\frac{(s-z_{1})(s-z_{2})...(s-z_{nz})}{(s-p_{1})(s-p_{2})...(s-p_{nz})}$
其中 $s = 2 πi f$ 。

利用Matlab绘制仪器响应曲线

clear;close all;

% ZEROS	5
z = [+0.000000e+00+0.000000e+00i	
	 +0.000000e+00+0.000000e+00i
	 -1.344330e-02+0.000000e+00i	
	 -1.344330e-02+0.000000e+00i];	
% POLES	6
p =	[-8.191920e-03+0.000000e+00i	
	 -1.235380e-02-1.234000e-02i	
	 -1.235380e-02+1.234000e-02i	
	 -1.883430e-02+0.000000e+00i	
	 -3.918000e+01-4.912000e+01i	
	 -3.918000e+01+4.912000e+01i];	
% CONSTANT	
k = 1.752545e+13;

sys = zpk(z,p,k);

P = bodeoptions;
P.FreqUnits = 'Hz';
P.Grid = 'on';
P.XLim = {[0.0001 1000]};

bode(sys,P);