基于hilbert变换的数字信号_基于Hilbert变换处理绝对重力仪测量数据
绝对重力仪一般采用激光干涉测距原理测量重力加速度,即每当落体下落激光半波长的距离时,接收器上便会相应产生一个干涉条纹的变化,由此得到落体的位移信息,结合条纹变化时间得到时间-位移数据,进而重建落体轨迹,解算重力加速度。重建落体轨迹的方法有基于硬件处理法与基于数字信号处理法。硬件方法中,使用电路过零点比较器检测干涉条纹信号中的过零点,在每个过零点处比较器产生一个脉冲信号,每N个零点脉冲信号生成一个脉冲,送入以外部铷原子钟为时间基准的时间间隔分析仪中,由此取得落体的位移-时间数据。数字信号处理法有双采样过零点检测法、复杂外差调制法、加窗二次微分法[以及零差干涉条纹正交法[。
本文基于Hilbert变换相位平移原理提取各采样点瞬时相位,利用瞬时相位与位移的关系取得各采样点位移,结合采样频率得到位移-时间数据,最后采用最小二乘法拟合落体轨迹,得到位移拟合残差、g值偏差。同时,本文采用数据加窗的方法抑制Hilbert变换过程中产生的吉伯斯现象,使用零相移滤波器与去线性趋势方法去除信号的噪声干扰,降低测量误差。
1 算法原理
1.1 Hilbert变换原理
Hilbert变换定义为:
$
\begin{array}{l}
\mathit{\tilde U}{\rm{(}}\mathit{t}{\rm{) =\mathit H}}\left\{ {\mathit{U}{\rm{(}}\mathit{t}{\rm{)}}} \right\} = h{\rm{(}}\mathit{t}{\rm{)}}\mathit{U}{\rm{(}}\mathit{t}{\rm{) = }}\\
\int_{ - \infty }^\infty U (\tau )h(t - \tau ){\rm{d}}\tau = \frac{1}{{\rm{ \mathit{ π} }}}\int_{ - \infty }^\infty {\frac{{U(\tau )}}{{t - \tau }}} {\rm{d}}\tau
\end{array}
$
(1)
式中,$\tilde U$(t)为U(t)与脉冲响应函数h(t)=1/πt的卷积, 其频率响应为:
$
H({\rm{j}}\omega ) = \left\{ \begin{array}{l}
{\rm{j}}, \omega < 0\\
0, \omega = 0\\
- {\rm{j}}, \omega > 0
\end{array} \right.
$
(2)
式(2)的频率响应如
图 1 Hilbert变换频率响应 Fig. 1 Frequency response of Hilbert-transform
此时,U(t)的解析信号为:
$
\mathit{q}{\rm{(}}\mathit{t}{\rm{) = }}\mathit{U}{\rm{(}}\mathit{t}{\rm{) + i}}\mathit{\tilde U}{\rm{(}}\mathit{t}{\rm{)}}
$
(3)
由此可得U(t)的瞬时相位:
$
\mathit{\Phi }(t) = {\tan ^{ - 1}}(H\{ U(t)\} /U(t))
$
(4)
在连续条件下对正弦信号进行Hilbert变换,可得H{sin(x)}=-cos(x)(x>0)。然而离散条件下进行Hilbert变换时([,相位误差两端产生纹波(g值的过程中,数据段长短不影响经过拟合与修正后最终所得的g0,因此在使用Hilbert变换对干涉条纹进行处理时,为降低误差,在数据拟合前对数据进行加窗选取,去除出现纹波误差的数据段。
图 2 正弦信号Hilbert变换的相移及误差 Fig. 2 Phase shift and error of sine signal's Hilbert-transfer
1.2 信号模型与时间修正
采用Murata[给出的条纹信号模型:
$
U(t) = {u_0} + \beta t + {A_0}\sin (\mathit{\Phi }(t))
$
(5)
式中,u0与β分别为直流偏移与线性漂移成分,A0为条纹信号幅值。根据激光干涉测距原理,当位移S(t)变化为激光半波长λ/2时,对应的相位变化为2π。此时,相位Φ(t)为:
$
\mathit{\Phi }(t) = 2{\rm{ \mathit{ π} }}\frac{{S(t)}}{{\lambda /2}} = \frac{{4{\rm{ \mathit{ π} }}}}{\lambda }S(t)
$
(6)
因此,在一个条纹周期内,条纹采样点位移为:
$
S(t) = \frac{{\lambda \mathit{\Phi }(t)}}{{4{\rm{ \mathit{ π} }}}}
$
(7)
整个自由落体过程中,位移变化为:
$
{S^\prime }(t) = \frac{{n\lambda }}{2} + S(t)
$
(8)
式中,n为条纹周期数,n=0,1,2,…。
在忽略重力梯度的情况下,自由落体规律为:
$
S(t) = {s_0} + {v_0}t + \frac{1}{2}{g_0}{t^2}
$
(9)
加入重力梯度影响后,自由落体模型为:
$
\begin{array}{l}
{S^\prime }(t) = {s_0} + {v_0}t + \frac{1}{2}{g_0}{t^2} + \\
\gamma \left( {\frac{1}{2}{s_0}{t^2} + \frac{1}{6}{v_0}{t^3} + \frac{1}{{24}}{g_0}{t^4}} \right)
\end{array}
$
(10)
式中,s0、v0、g0为所选坐标系统原点的初始条件,γ为重力梯度常数。当γ未知时,取γ=0,以二阶落体模型求解g值,并由式(11)修正得到g0:
$
g = {g_0} + \gamma {h_{{\rm{ef}}}}
$
(11)
式中,hef为有效测量高度。
在进行重力加速度测量时,当落体在高度H处下落Δx距离时,光束仍需要经历H-Δx距离才能到达接收器,而光的速度是有限的,因此光速的有限性会导致干涉条纹生成时间落后于落体到达该下落位置的时间,最终影响测量结果。因此,需要对时间进行修正[:
$
t_i^\prime = t - \frac{{\Delta x}}{c} = t - \frac{{{s_i} - {s_0}}}{c}
$
(12)
式中,c为光速,t′i为修正时间,si为第i个观测点的位移。
1.3 零相移滤波原理与去线性趋势处理
零相移滤波算法的核心分为滤波、时域翻转信号、再次滤波、再次翻转信号等4步。
首先是翻转信号与原信号的Z变换关系。假设一个有限长序列信号{x(n)},n=1,2,…,n0,将该信号延拓至整个x轴,则信号表达式为:
$
{x_1} = \left\{ {\begin{array}{*{20}{l}}
{x\left( {n\% {n_0}} \right), n\% {n_0} \ne 0}\\
{x\left( {{n_0}} \right), n\% {n_0} = 0}
\end{array}} \right.
$
(13)
其中,%表示取模。同样,假设信号序列{y(n)}为信号{x(n)}时域翻转后的序列,将信号{y(n)}延拓至整个x轴,则延拓信号y1(n)和延拓信号x1(n)关于y轴对称,即y1(n)=x1(-n)。
对以上2个信号进行双边Z变换:
$
X(z) = \sum\limits_{n = - \infty }^\infty {{x_1}} (n){z^{ - n}}
$
(14)
$
\begin{array}{l}
Y(z) = \sum\limits_{n = - \infty }^\infty {{y_1}} (n){z^{ - n}} = \sum\limits_{n = - \infty }^\infty {{x_1}} ( - n){z^n} = \\
\;\;\;\;\;\;\;\;\sum\limits_{n = - \infty }^\infty {{x_1}} (n){\left( {{z^{ - 1}}} \right)^n} = X\left( {{z^{ - 1}}} \right)
\end{array}
$
(15)
可知,对信号X(z)进行时域翻转后,该信号Z变换为X(z-1)。假设滤波器的传递函数为H(z),则零相移滤波器的输出Y1(z)为:
$
{Y_1}(z) = X(z)H(z)H\left( {\frac{1}{z}} \right)
$
(16)
式中,H(z)为实系数等式,故有H(ejω)=H(e-jω)。令z=ejω,则经过二次滤波和翻转的信号输出为:
$
\begin{array}{l}
{Y_1}\left( {{{\rm{e}}^{{\rm{j}}\omega }}} \right) = X\left( {{{\rm{e}}^{{\rm{j}}\omega }}} \right)H\left( {{{\rm{e}}^{{\rm{j}}\omega }}} \right)H\left( {{{\rm{e}}^{ - {\rm{j}}\omega }}} \right) = \\
\;\;\;\;\;\;\;\;\;\;X\left( {{{\rm{e}}^{{\rm{j}}\omega }}} \right){\left| {H\left( {{{\rm{e}}^{{\rm{j}}\omega }}} \right)} \right|^2}
\end{array}
$
(17)
可知,零相移滤波器的传递函数为|H(ejω)|2,输出信号Y1(ejω)相对于X(ejω)并没有发生相位上的偏移。
由信号模型式(5)可见,条纹信号中包含直流偏移与线性漂移,两者构成信号的线性趋势。采用最小二乘法对信号进行多项式拟合后,信号被分解为各阶数据。从中提取出信号的线性趋势数据,使用原信号减去该数据,便可得到去线性趋势的数据。
使用式(5),取A0=1,u0=0.01,β=0.01,Φ(t)=t,并引入5%幅值高斯白噪声作为外部噪声,得到信号(
图 3 信号高频噪声滤波以及线性噪声去除 Fig. 3 High frequency noise of signal filtering and linear trend removal
1.4 数据处理流程
数据处理流程如g值。
图 4 数据处理流程 Fig. 4 The flow chart of data processing
2 仿真与实验结果分析
2.1 仿真结果与分析
使用式(5)、式(9),取g=9.8 m/s2,落体的下落距离为H=25 cm,A0=1,u0=0.01,β=0.01,引入5%幅值高频高斯白噪声作为外部噪声,并使该信号采样频率为200 MHz,建立干涉条纹仿真信号模型。采用本文数据处理流程对数据进行处理,并以频率作为数据窗口依据,分别在去除噪声与保留噪声条件下取得仿真结果(
表 1
Tab. 1
表 1 g值仿真结果
Tab. 1 Simulation results of g values
频率范围
(去噪)/MHz
g值残差
/μGal
频率范围
(不去噪)/MHz
g值残差
/μGal
0~7
7.98×10-1
0~7
4.28×10-4
0.5~6.5
3.84×10-4
0.5~6.5
5.60×10-4
1~6
1.61×10-4
1~6
8.14×10-4
1.5~5.5
2.93×10-4
1.5~5.5
4.27×10-3
2~5
3.11×10-5
2~5
9.36×10-3
2.5~4.5
2.38×10-4
2.5~4.5
9.71×10-3
3~4
5.65×10-5
3~4
3.76×10-2
表 1 g值仿真结果
Tab. 1 Simulation results of g values
根据标准二阶线性模型最小二乘法[, g值不确定度为:
$
{\sigma _g} = \frac{{12\sqrt 5 }}{{{T^2}}} \cdot \frac{{{\sigma _S}}}{{\sqrt N }}
$
(18)
式中,N为加入计算的数据个数,T为整个数据段的时间,σS为最小二乘位移残差的标准差。
由σS。但是随着窗口的缩小,去除噪声时,σS的值极小,且随窗口缩小而变小,N与T的减小并不能使σg增大; 未去除噪声时,窗口的缩小并不能进一步改善σS,因此随着N与T的减小,σg不断增大。
2.2 实验结果与分析
在中国地震局地震研究所绝对重力仪上进行9次下落实验,对所得干涉条纹数据进行去噪声后,使用本文方法处理所得数据。数据窗口选取1~6 MHz,其中一次落体实验所得的最小二乘位移拟合残差如
图 5 最小二乘拟合位移残差 Fig. 5 Residual of the least squares fitting displacement
根据该残差可以直接判断数据处理方法的优劣[。使用各次实验最小二乘位移残差的标准差作为σS,参与拟合的数据段的时间为T,使用式(18)求解g值不确定度,结果如g值不确定度在0.2 μGal以内。
图 6 实验得到的g值偏差σg Fig. 6 σg of g values obtained from the experiment
分别使用双采样过零点检测法与本文Hilbert变换法进行数据处理,将结果与仪器真值比较,所得g值误差如μGal,优于双采样法。
图 7 双采样法与Hilbert变换法结果误差 Fig. 7 Errors of two-sample method and Hibert-transfer method
3 结语
本文基于Hilbert变换相位平移的原理,结合零相位低通滤波、信号去线性趋势处理,以及拟合数据加窗口选取,实现对重力加速度的精确解算。仿真结果与实验测量所得g值不确定度表明,该方法准确、可靠。实验结果所得平均测量误差为0.12 μGal,相比于双采样法,本文方法具有更高的测量精度。对比仿真结果,g值测量精度仍能通过研究得到进一步提高。
基于hilbert变换的数字信号_基于Hilbert变换处理绝对重力仪测量数据相关推荐
- python基于水色图像的水质评价_基于Python和遥感图像的膨胀与腐蚀操作
引言 膨胀与腐蚀是图像形态学中的基本操作,本文将从两个角度实现python的膨胀与腐蚀算法,分别是特征提取与分割后处理.目前,在RGB图像领域,这两种形态学算法经常用于分割结果的处理,例如连通区处理. ...
- 基于hsv的亮度调整算法_基于手绘工程图离线识别的预处理研究精品论文推荐
DOI: 10. 3969 / j. issn. 1009-9492. 2020. 11. 028 李春晓,田怀文,刘奇,等. 基于手绘工程图离线识别的预处理研究[J] . 机电工程技术,2020,4 ...
- arduino读取水位传感器的数据显示在基于i2c的1602a上_基于PC-CAN适配卡和总线技术实现潜水电机分布式监控系统的应用方案...
目前,国内外已有多家公司开发出针对个体潜水电机的智能电了保护装置,虽然可以达到保护电机的目的,但无法实现集中管理.因此,对潜水电机群的实时监控与统一管理的要求日益迫切.本文设计了基于CAN现场总线的潜 ...
- 基于python的房地产数据分析_基于Python的数据分析
转载 | CSDN 编辑 | 雷课小雷 下面来介绍一下基于Python的数据分析,主要介绍数据分析的概念.数据分析流程.Python优势.常用模块的用途以及使用Python进行数据分析的学习方法及步骤 ...
- 基于高分辨率影像城市绿地信息提取_基于ENVI的高分辨率遥感影像城市绿地信息提取研究_陈阳pdf.doc...
基于ENVI的高分辨率遥感影像城市绿地信息提取研究_陈阳pdf 第卷第期测绘工程,244 Vol.24No.4年月,20154Engineering of Surveying and MappingA ...
- 基于安卓的考试系统_基于安卓11定制!华为最新手机系统曝光:体验堪比苹果iOS!...
在最近的一场发布会上,华为正式宣布了自家的HMS和AppGallery服务,对标安卓Play商店和苹果Appstore商店,这一举措让华为再度登上风口浪尖.这种做法在业界人士眼里的目的只有一个,华为要 ...
- 基于数据库数据增量同步_基于canal实现分布式数据同步
应用场景 分布式架构中,数据同步常常是个大问题.例如,mysql中的数据,可能在ElasticSearch有一份索引,在redis有一份缓存,在Nginx有一份缓存,这时候只要你修改了mysql中的数 ...
- 基于c#的相关性分析_基于数字高程模型的城市地貌与地名相关性分析——以兰州市为例...
江西地名研究 基于数字高程模型的城市地貌与地名 相关性分析--以兰州市为例 文/张鹏丽,李育 提要:以兰州市为例,使用 ASTER 30米分辨率的数字高程模型DEM并提取了兰州市周边所有地名信息:通过 ...
- 基于ssm的用户管理系统_基于SSM的高校学生实习管理系统
前言:本文描述了基于SSM的高校学生实习管理系统的实战项目,分析了该项目的功能需求以及最终呈现的效果 基于SSM的高校学生实习管理系统应完成如下功能: 1. 用户角色划分 本系统主要用户有四类:实习企 ...
- 基于单片机自动升旗系统_基于视觉定位的机器人全自动冲击系统
产品图片 产品图片 01 产品简介 产品集成了力学冲击实验机.全自动冷却机.视觉定位.机器人.分拣收集装置及信息管理软件等系统,在人工进行批量试样的上料.组批及任务下达后,自动进行试样的降温保温.上料 ...
最新文章
- Android 动态修改参数配置
- 计算机模块考试注意事项,常考计算机模块考试操作注意事项.doc
- Socket之UDP服务器【Python】
- SAP CRM 和 Cloud for Customer 的 Document flow API 介绍
- 多选取值_机制砂如何控制MB值和石粉含量
- ntp symmetric_Python使用示例设置symmetric_difference()方法
- [daily] fedora用过光盘做dnf repo
- 一位清华学长的CTO成长之路 -- 不断的打破自己的边界
- 如何找到局域网内所有主机ip
- 济南铃木UU125无破线改装仪表盘详解
- win10计算机磁盘图标,主编教你win10系统无法正常显示硬盘图标的方法
- 互联网巨头“赢家通吃”,会产生哪些严重后果?
- nodeJS实现简易爬虫
- lsof的安装及使用
- WAP网站的推广方式(自整合篇)
- 基于vue+springboot的校园疫情健康打卡和离校审批系统的设计 (百度地图API对接)
- 剑指核心-TaoCloud全闪SDS助力构建高性能云服务
- 简单玩一下SM2258H+双贴BGA316固态U盘,还有开卡软件
- 短网址(Short URL)服务
- php pdo连接oracle乱码,php pdo oracle中文乱码的快速解决方法_PHP