CN111352099B - 一种基于互相关信号相位分解的时延估计方法 - Google Patents

一种基于互相关信号相位分解的时延估计方法 Download PDF

Info

Publication number
CN111352099B
CN111352099B CN201811564624.0A CN201811564624A CN111352099B CN 111352099 B CN111352099 B CN 111352099B CN 201811564624 A CN201811564624 A CN 201811564624A CN 111352099 B CN111352099 B CN 111352099B
Authority
CN
China
Prior art keywords
signal
cross
correlation
monocycle
sequence
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811564624.0A
Other languages
English (en)
Other versions
CN111352099A (zh
Inventor
仇璐
吴振谦
曹士杰
赵辛东
蒋志迪
俞牡丹
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
College of Science and Technology of Ningbo University
Original Assignee
College of Science and Technology of Ningbo University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by College of Science and Technology of Ningbo University filed Critical College of Science and Technology of Ningbo University
Priority to CN201811564624.0A priority Critical patent/CN111352099B/zh
Publication of CN111352099A publication Critical patent/CN111352099A/zh
Application granted granted Critical
Publication of CN111352099B publication Critical patent/CN111352099B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/523—Details of pulse systems
    • G01S7/526—Receivers
    • G01S7/53—Means for transforming coordinates or for evaluating data, e.g. using computers
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01F—MEASURING VOLUME, VOLUME FLOW, MASS FLOW OR LIQUID LEVEL; METERING BY VOLUME
    • G01F1/00—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow
    • G01F1/66—Measuring the volume flow or mass flow of fluid or fluent solid material wherein the fluid passes through a meter in a continuous flow by measuring frequency, phase shift or propagation time of electromagnetic or other waves, e.g. using ultrasonic flowmeters

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Electromagnetism (AREA)
  • Fluid Mechanics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种基于互相关信号相位分解的时延估计方法，将接收信号传播延时分解为K个整周期和ψ的相位偏移；通过设置单周期函数信号，求得互相关信号与所述单周期函数信号间的相关系数的极值，再往0时刻方向进行整周期移相，求得互相关信号的第一个起振点，从而得到接收信号的时延。整个运算过程复杂度低，非常适合基于C语言在任何嵌入式芯片编程实现，且具有广泛的场景适应性，能够获取高精度的测量结果。

Description

一种基于互相关信号相位分解的时延估计方法
技术领域
本发明涉及时延估计领域，特别涉及一种基于互相关信号相位分解的时延估计方法。
背景技术
目前，对于信号时延的估计方法，有多种方式，如电平阈值比较法、小波分析法、互相关法、先估计整周期后进行精确相位分析法；以上方法各有其优缺点。如电平阈值比较法分为单阈值和双阈值法，信号到达阈值的时刻即为求解的时延，方法简单，但是对于干扰信号的抑制能力很差，同时阈值的选择需要针对特定的测量距离、测量温度和液体浓度进行标定，不具备应用的广泛适应性。
小波分析法，其计算结果较好，但是计算过于复杂，在嵌入式芯片中实现较难，运算时间过长，实时性差。
互相关将发射信号和接收信号进行互相关运算，求取每个互相关周期信号的局部极值，再将所有局部极值进行多项式拟合，求解全局极值时刻即为求解的时延。但是由于互相关运算的理论原理，决定所求解的结果与所用测量的发射信号的时长有强关联关系，同时互相关信号的全局极值出现时刻受超声换能器的接收信号起振速度的快慢影响，所以测量结果的标定必然需要在不同的测量条件如距离、温度和液体浓度进行标定，算法的应用适应性差。同时进行极值拟合的过程计算复杂，在嵌入式芯片中实现较难，实时性差。另外由于互相关运算基于离散数据进行，由于算法原理的原因，导致计算的结果一定比实际值大，偏差的数据受限于AD采样芯片的采样率，所以该方法的计算结果受限于AD采样率，要提高测量精度，比如提高系统成本。
先估计整周期后进行精确相位分析法，是先用一脉冲信号估计整周期，后用一连续同频率信号估计精确相位。该方法能得到较为精确的结果，但是需要分两次发射信号，对硬件要求高，测量过程复杂。
基于克服上述分析方法的缺点，提出一种运算过程复杂度低、测量精度高、适合基于C语言在任何嵌入式芯片编程的时延估计方法是亟待解决的问题。
发明内容
本发明就是为了解决运算过程复杂度高、测量精度低、不适合基于C语言在任何嵌入式芯片编程的问题，提出一种基于互相关信号相位分解的时延估计方法，先求解相位偏移，再往0时刻方向进行整周期移相，求得互相关信号的起振点，从而得到接收信号的时延。
为实现上述发明目的，本发明采取以下技术方案：
一种基于互相关信号相位分解的时延估计方法，将接收信号传播延时分解为K个整周期和ψ的相位偏移，通过设置标准单周期正弦波信号，求得互相关信号与所述标准单周期正弦波信号间的相关系数的极值，再往0时刻方向进行整周期移相，求得互相关信号的第一个起振点，从而得到接收信号的时延。
优选地，包括以下步骤：
S1、对接收信号与理想发射信号进行互相关运算得到互相关信号，构建单周期函数信号，并对单周期函数信号进行步进移相，求解互相关信号和单周期函数信号的相关系数的极值，极值处相位即为接收信号传播延时的精确相位偏移ψ；
S2、对求得的精确相位偏移为ψ的单周期函数信号进行整周期移相，确定具体的移相周期数K，从而求得互相关信号的第一个起振位置点，从而求得接收信号传播延时。
优选地，步骤S1中，求解互相关信号和标准单周期正弦波信号的相关系数极值处的精确相位偏移ψ，具体包括以下步骤：
A1、选取时长为Ts的理想发射信号x(t)；
A2、获取接收波信号y(t)；
A3、进行互相关运算，生成互相关信号f(τ)的序列f(mT)；
A4、寻找f(mT)为最大值时的Tmax；
A5、获取离Tmax最近的整周期点，记为Km；
A6、构成标准单周期正弦波信号g(τ)的步进移相的单周期函数序列g((mT)θ(i)(j))；
A7、求取序列g((mT)θ(i)(j))与步进移相的互相关信号序列f((mT)θ(i)(j))的相关系数ρgf(Km,θ(i))；
A8、求取ρgf(Km,θ(i))最小值时的θ(i)，即为精确相位偏移值ψ；
其中，步骤A1与步骤A2无先后顺序区别。
优选地，所述互相关信号，如下式所示：
Figure GDA0003537989730000041
式中，τ表示互相关信号时域自变量；
优选地，所述标准单周期正弦波信号g(τ)如下式所示：
Figure GDA0003537989730000042
式中，θ为待求解相位偏移自变量，ω为超声激励信号的角频率,Km为f(mT)序列最大值处零点方向最近的整周期点。
优选地，所述单周期正弦信号g(τ)与互相关信号f(τ)的相关系数，如下式所示：
Figure GDA0003537989730000043
式中，
Figure GDA0003537989730000044
优选地，所述互相关信号离散形式序列f(mT)计算方式如下：
Figure GDA0003537989730000051
其中，x(t)和y(t)由同一个脉冲触发信号P(nT)进行数字化采样，分别形成x(nT)和y(nT)：
Figure GDA0003537989730000052
式中，Nx为接收信号起振后第一个采样点序列值，Ns为发射激励信号最后一个采样点序列值，Nr为接收信号最后一个采样点序列值；x(t)表示激励信号；y(t)表示接收波信号。
优选地，选取步进量Δθ＝θ(i)-θ(i-1)，寻找f(mT)序列中满足
Figure GDA0003537989730000053
的所有点，记为(mT)θ(i)(j),代入标准单周期正弦波信号g(τ)，得步进移相的单周期正弦波信号序列g((mT)θ(i)(j))。
优选地，计算ψ的方法如下：
Figure GDA0003537989730000054
优选地，步骤S2中，对在相关系数极值的相位ψ处的单周期正弦波信号进行整周期移相，确定具体的移相周期数Kx，具体包括以下步骤：
B1、从Km开始，将
Figure GDA0003537989730000061
往0时刻方向进行整周期移相至k个周期；
B2、寻找互相关信号离散形式序列f(mT)中满足
Figure GDA0003537989730000062
的所有点位置序列；
B3、求取所述位置序列的整周期移相的单周期值序列g((mT)k(j))；
B4、获取所述位置序列的整周期移相测量值序列f((mT)k(j))；
B5、计算单周期值序列g((mT)k(j))与测量值序列f((mT)k(j))的相关系数
Figure GDA0003537989730000063
B6、相关系数
Figure GDA0003537989730000064
是否小于设定值ε，若是，转B8；若否，进入下一步；
B7、k值减1，转B1；
B8、判定k-1为互相关信号周期Kx，根据当前k-1计算互相关信号的时延估计结果Tx。
优选地，选择k＝Km,Km-1,Km-2…，将
Figure GDA0003537989730000065
往0时刻方向进行整周期移相，寻找f(mT)序列中满足
Figure GDA0003537989730000066
的所有点，记为(mT)k(j),代入
Figure GDA0003537989730000067
得单周期正弦波信号序列g((mT)k(j))。
优选地，求解
Figure GDA0003537989730000071
Kx＝Kmin,从而得到
Figure GDA0003537989730000072
式中Ks为Ts包含的整周期个数，Ts为参与互相关运算选取的理想发射信号时间长度。
与现有技术相比，本发明的有益效果为：
本发明通过先利用单周期函数进行精确步进移相，求取互相关信号序列与单周期函数序列相关系数的最大值处的步进相位，获取精确相位；再进行整周期移相，对互相关信号的第一个周期进行整周期匹配，以此判定互相关信号第一个周期的出现位置，从而求取时延估计结果。整个计算过程利用了互相关运算，运算过程复杂度低，非常适合基于C语言在任何嵌入式芯片编程实现。
进一步地，本发明通过获取第一个互相关信号的出现位置，与超声波接收信号的起振速度没有关系，算法有广泛的场景适应性，不需要考虑测量的温度、浓度和距离等因素，在实施过程中免去了标定过程。
进一步地，本发明通过第一步的精确移相，能够获取高精度的测量结果，测量精度高达10-8秒数量级。
附图说明
图1是现有技术中超声波测量系统模型图；
图2是现有技术中超声波换能器接收端等效电路；
图3是现有技术中互相关结果与接收信号时间轴差示意图；
图4是现有技术中互相关运算时起振点与相邻采样点相互关系示意图；
图5是本发明具体实施方法互相关精确相位步进和整周期匹配算法原理示意图；
图6是本发明具体实施例互相关精确相位步进和整周期匹配算法流程示意图；
图7是本发明具体实施例互相关精确相位步进流程示意图；
图8是本发明具体实施例整周期匹配算法原理流程示意图；
图9、10、11是本发明具体实施方法的测量起振波极值与最大值的比值变化曲线示意图；
图12是现有技术中互相关算法计算相对传播延时示意图；
图13是本发明具体实施方法的互相关精确相位步进和整周期匹配法计算相对传播延时示意图。
具体实施方式
采用互相关进行接收信号的时延计算的系统很多，本发明以超声波测量系统进行说明，但并不限于此，以此类推。
超声波测量系统模型如图1所示，换能器S为超声波源发生器，通过发射激励信号x(t)发射超声波，通过介质传到换能器R，换能器R将接收到的超声波转换成电信号形成接收波信号y(t)。x(t)和y(t)由同一个脉冲触发信号P(nT)进行数字化采样，分别形成x(nT)和y(nT)，即满足x(t)和y(t)在0时刻产生第一个脉冲触发信号的条件，其中T为脉冲触发周期，即为采样周期。
压电式超声波换能器接收端的等效电路图如图2所示，图中Xi是超声波接收端的动态阻抗，C0是超声波接收端的静态电容，R是动态电阻。根据动态电路时域分析，如果发射端超声波为一频率恒定的正弦信号，那么接收信号具有指数形式的起振过程，可以表示为：
Figure GDA0003537989730000091
其中ω为超声激励信号的角频率，A和B分别为发射和接收信号的幅度，Tx为超声波测量目标TOF(time of flight飞行时间)，c为换能器起振参数，Ts为发射激励信号时间长度，Tr为接收信号结束时间点。
应用经典互相关算法对上述发射激励信号和接收信号进行互相关运算，由于x(t)和y(t)均为有限时长脉冲信号，互相关运算公式定义为：
Figure GDA0003537989730000092
式中，τ表示互相关信号时域自变量。
假设接收信号时长大于发射激励信号时长，即满足Tr-Tx>Ts,得到分段积分结果：
Figure GDA0003537989730000101
其中Ts为选取的参与互相关运算的理想发射信号时长，Tx为待求的实际时延
当τ∈[Tx-Ts,Tr-Ts)时，求解得：
Figure GDA0003537989730000102
其中F1,F2,F3,F4,F5,F6为常系数。如果满足
Figure GDA0003537989730000103
即Ts是激励波形周期的整数倍，容易证明得，
Figure GDA0003537989730000104
说明f(τ)在τ＝Tx处为局部极值点。根据f(τ)的解析形式得出以下结论：
1)y(t)从Tx开始起振，f(τ)从Γx＝Tx-Ts处开始起振，刚好相差Ts，其中f(τ)的起振点记为Γx；
2)f(τ)与y(t)的频率相同；
3)f(τ)在τ＝Tx为局部极值点；
如图3所示，上部分为接收信号波形，下部分为互相关计算波形，计算得到的互相关波形比接收波形在时间轴上往前移了Ts的长度。
经典互相关算法，应用于带起振过程超声信号TOF即Tx的求解时，由寻求f(τ)的极值点转变为在f(τ)中求解偏移Ts处的局部极值点。
由于超声波压电换能器特性参数无法直接获取，同时在测试过程中由于测试距离、工作电压和震动介质等因素的变化，换能器的起振过程无法准确标定，因此上述超声波收发信号的连续数学模型无法在实际测量计算中直接应用，必须基于数字化的测试结果进行离散互相关极值算法进行TOF(time of flight飞行时间)求解，下面对离散互相关算法进行分析。
互相关算法的离散形式为：
Figure GDA0003537989730000111
x(t)和y(t)由同一个脉冲触发信号P(nT)数字化采样后，经过由式⑴得：
Figure GDA0003537989730000112
Nx为接收信号起振后第一个采样点序列值，Ns为发射激励信号最后一个采样点序列值，Nr为接收信号最后一个采样点序列值。
接收信号中的起振点不一定刚好落在采样点上，如图4所示，为起振点落在两个相邻采样点之间的情况：
分别计算起振点前后两个采样点的互相关值，并求他们的差值：
Figure GDA0003537989730000113
Figure GDA0003537989730000121
t1和t2分别为前后采样点离起振点的时间差的绝对值。由于NrT-NxT>NsT＞＞t1,t2,可以忽略Δf(MT)中sint1和sint2的乘积项，易得Δf(MT)>0。如果采样率很高，T很小，那么Δf(MT)越接近0。
以上推导说明，当采用离散互相关算法通过极值判断获取时间差估计的方法，极值一定靠近(M+1)T采样点，计算结果总会落在实际时延的外侧，即计算结果比实际值偏大，测量结果直接受限于系统的采样率。为了提高测量精度，只能提高采样率，从而增加了系统成本。
基于此，本发明采用以下方法进行延时估计。
对接收信号与理想发射信号进行互相关运算得到互相关信号，设置单周期函数信号，求得互相关信号与所述单周期函数信号间的相关系数的极值，再往0时刻方向进行整周期移相，求得互相关信号的第一个起振点，从而得到接收信号的时延。即把待测传播延时分解为K个整周期和ψ的相位偏移。
具体地，包括以下步骤：
S1、对接收信号与理想发射信号进行互相关运算得到互相关信号，构建单周期函数信号，并对单周期函数信号进行步进移相，求解互相关信号和单周期函数信号的相关系数的极值，极值处相位即为接收信号传播延时的精确相位偏移ψ；
S2、对求得的精确相位偏移为ψ的单周期函数信号进行整周期移相，确定具体的移相周期数K，从而求得互相关信号的第一个起振位置点，从而求得接收信号传播延时。
本发明的一个实施例中，步骤S1中，先选取时长为Ts的理想发射信号x(t)，使得Ts<Tx并且Ts<Tr-Tx，作为互相关运算的数据源获取接收波信号y(t)，对两者进行互相关运算，得到互相关信号f(τ):
Figure GDA0003537989730000131
其离散形式为：
Figure GDA0003537989730000132
其中，x(t)和y(t)由同一个脉冲触发信号P(nT)数字化采样后，得到公式:
Figure GDA0003537989730000133
式中，Nx为接收信号起振后第一个采样点序列值，Ns为发射激励信号最后一个采样点序列值，Nr为接收信号最后一个采样点序列值。
若互相关之后的波形f(τ)的起振点记为Γx，根据图3及相关分析可得Γx＝Tx-Ts，因此求解Tx可以转化为求解f(τ)的起振点，即Tx＝Γx+Ts。
设：
Figure GDA0003537989730000134
其中Kx为待求解周期整数，
Figure GDA0003537989730000135
表示从0时刻开始最接近Γx的整周期点，ψ为待求解自该整周期点的偏移相位差。记Ks为Ts包含的整周期个数，则：
Figure GDA0003537989730000141
构造单周期函数信号：
Figure GDA0003537989730000142
式中θ为待求解相位偏移自变量，ω为超声激励信号的角频率,Km为f(mT)序列最大值处零点方向最近的整周期点。
计算单周期函数信号g(τ)与互相关信号f(τ)的相关系数ρgf(θ)：
Figure GDA0003537989730000143
式中，
Figure GDA0003537989730000144
求解相关系数的极值，得到精确相位ψ，即：
ψ＝argmax0→2πρgf(Km,θ) (15)；
步骤S2中，进行整周期移相，确定具体的Kx。将求得的精确相位偏移为ψ的单周期函数信号在
Figure GDA0003537989730000145
区间从Km往0时刻方向进行整周期移相，通过选择合适的阈值ε判断是否为互相关后信号的第一个波形，从而求得第一个起振波形位置点。即求解下式中k的最小值：
Figure GDA0003537989730000146
Kx＝Kmin (17)；
把Kx代入公式(12),即可得到Tx。
因为，互相关信号起振点Γx后第一个半周期不是正弦波，从第一个周期的
Figure GDA0003537989730000151
到第二个周期的
Figure GDA0003537989730000152
有较好正弦波特征，为提高判别率，故整周期移相选择单周期函数在
Figure GDA0003537989730000153
区间信号，故上述整周期移相判别原则(16)获得的是互相关后波形信号第二个周期的整周期偏移值，故准确起振点为该值减1(式(17))。
本发明的一个实施例中，延时估计测量采用离散方式进行，步骤S1中，包括如下步骤：
1)在发射信号中选取时长为Ts的信号，使得Ts＜Tx并且Ts＜Tr-Tx，作为互相关运算的数据源；
2)按照式(9)进行互相关计算，求得f(mT)序列；
3)在f(mT)寻找最大值的位置时刻Tmax；
4)按照式(13)构造单周期函数；
5)选取步进量Δθ＝θ(i)-θ(i-1)，寻找f(mT)序列中满足
Figure GDA0003537989730000154
的所有点，记为(mT)θ(i)(j),代入公式(13)，得步进移相的单周期函数序列g((mT)θ(i)(j))；
6)计算g((mT)θ(i)(j))与f((mT)θ(i)(j))的相关系数ρgf(Km,θ(i))；
7)求取相关系数ρgf(Km,θ(i))最小值时的θ(i),即为精确相位值ψ。
步骤S2中，包括如下步骤：
1)选择k＝Km,Km-1,Km-2…,将
Figure GDA0003537989730000155
往0时刻方向进行整周期移相k个周期；
2)寻找f(mT)序列中满足
Figure GDA0003537989730000161
的所有点位置序列，记为(mT)k(j),代入式(13)，得整周期移相的单周期函数序列g((mT)k(j))；
3)计算g((mT)k(j))与f((mT)k(j))的相关系数
Figure GDA0003537989730000162
若
Figure GDA0003537989730000163
Figure GDA0003537989730000164
则该k-1即为互相关信号起振周期Kx；根据式(12)求得Tx。
本发明的一个具体实施例中，步骤S1中的流程图如图6所示，包括以下步骤：
A1、选取时长为Ts的理想发射信号x(t)；
A2、获取接收波信号y(t)；
A3、进行互相关运算，生成互相关信号f(τ)的序列f(mT)；
A4、寻找f(mT)为最大值时的Tmax；
A5、获取离Tmax最近的整周期点，记为Km；
A6、构成步进移相的单周期函数g(τ)的序列g((mT)θ(i)(j))；
A7、求取序列g((mT)θ(i)(j))与步进移相的互相关信号序列f((mT)θ(i)(j))的相关系数ρgf(Km,θ(i))；
A8、求取ρgf(Km,θ(i))最小值时的θ(i)，即为精确相位值ψ；
其中，步骤A1、A2没有先后的分别。
步骤S2中的流程图如图7所示，包括以下步骤：
B1、从Km开始将
Figure GDA0003537989730000171
往0时刻方向进行整周期移相至k个周期；
B2、寻找互相关信号序列f(mT)中满足
Figure GDA0003537989730000172
的所有点位置序列；
B3、求取所述所有点位置序列的整周期移相的单周期值序列g((mT)k(j))；B4、获取所述所有点位置序列的整周期移相的互相关函数序列f((mT)k(j))；
B5、计算整周期移相的单周期值序列g((mT)k(j))与整周期移相的互相关离散形式函数序列f((mT)k(j))的相关系数
Figure GDA0003537989730000173
B6、相关系数
Figure GDA0003537989730000174
是否小于设定值ε，若是，转B8；若否，进入下一步；
B7、k值减1，转B1；
B8、判定k-1为互相关信号周期Kx，根据当前k-1计算互相关信号的时延估计结果Tx。
如果找到一个整周期位置的相关系数，其相关系数小于所述设定阈值ε，就说明此位置点已经不是正弦波了，那么上一个就是第一个起振整周期。
本发明的一个具体实施例中，如图6所示，本发明的一种基于互相关信号相位分解的时延估计方法，包括以下步骤：
D1、根据理想发射信号与接收信号求得互相关信号；
D2、对单周期函数信号进行步进移相，得到步进移相单周期函数信号：
D3、求取各步进移相单周期函数信号与所述互相关信号的相关系数序列最大值；
D4、求解相关系数序列最大值处的精确相位偏移；
D5、将精确相位偏移单周期信号往零时刻方向进行整周期移相；
D6、求取互相关信号与整周期移相后的相关系数，如果所述相关系数小于设定阈值，上一个整周期即为互相关信号起振周期；此处的相关系数为每个整周期位置处的相关系数；
D7、根据互相关信号起振周期、理想发射信号时长、精确相位偏移，求得接收信号实际传播延时。
本发明的一个具体实施例中，对本估计方法进行测试验证，具体地，
为了仔细研究换能器在不同工作条件下起振情况对传播时延计算的影响，设计三种自变量条件的对比实验，分别为液体质量浓度、液体温度和测量距离。
该测试方案采用型号为DYW-1M-01EA的插入式管道流量计用换能器，属于超声波液体换能器，中心频率1MHz，电容量2800pF，最小阻抗25Ω。设定换能器驱动电路工作电压12V，换能器驱动信号由FPGA通过直接频率合成DDS方法产生1MHz的方波信号，经由驱动电路放大后输入发射端换能器；接收端换能器接收信号后，经过信号调理和采样率为20MHz的AD芯片高速采样，数据输入FPGA以及STM32内核芯片进行传播时延算法估计。
测量时，将换能器发射和接收端刚性固定后，置于体积为1000ml的圆柱形溶液中，容器直径6cm。
2、测试项目
分别将温度、质量浓度和距离作为自变量进行分组测试，具体测试条件如表1所示。
表1测试项目表
Figure GDA0003537989730000191
Figure GDA0003537989730000201
按照表1测试点，进行波形测试，在每个测试点测量5组数据，共获取165组数据，每组数据包含1400个采样点，波形时长为70us，内含一个完整的接收信号波形。
1、测试结果分析
在图3中起振波形起振点到最大值点之间，具备9-10个幅度为指数增长的正弦波。前5个起振波极值相对于最大值的比值变化趋势，能反应不同测量条件下超声波起振情况，如图8、9、10所示，图中标号为1-5所示曲线，分别为第1-5个起振波形极值相对该波形最大值的比值变化曲线，可以观察到在测量范围内液体质量浓度越高起振越快，探头距离越大起振越快，而温度对起振过程影响不大。
a)传播时延计算算法比较
以上测试分析结果说明，在超声波换能器型号、工作电压等因素确定的条件下，超声波的起振波随着液体质量浓度、测试距离的不同而发生变化，而起振波受测试温度的影响较小。根据前文所述，在不考虑起振差异的情况下，直接应用传统互相关算法进行传播时延的计算，会导致较大的误差，下面分别用传统互相关算法和本文提出的互相关精确相位步进和整周期匹配法对测试数据进行相对传播时延估计。
相对传播延时，不是直接计算每次测试的绝对传播时延，而是选取某次测试的接收波形作为参考波形，将所有测试接收波形的起振点移至该参考波形的起振点，计算相对传播时延，其本质是各次测试波形相对于参考波形的测试误差。在质量浓度测试组中选取0g/L的接收波形作为参考波形，在距离测试组中选取9cm的测试波形作为参考波形，在温度测试组中选取31℃的测试波形作为参考波形。
针对三种数据测试组，分别进行传统互相关算法进行相对传播时延估计，如图11所示。在变化质量浓度情况下，应用传统互相关运算进行传播时延估计，由于起振波形的差异导致估计误差达到10-6s量级，接近甚至超过了超声波整周期，这在实际测试中不能接受。在变化距离测试中也会导致10-7us量级误差。
利用本发明提出的互相关相位分析法对三组测试数据组进行传播延时误差估计，如图12所示，该方法对三种数据测试组呈现较为一致的测量结果，测量误差精确到10-8s数量级，在不同浓度、不同距离情况下比传统互相关算法提高了两个数量级精度，极大的提高了算法在不同测量条件下的测量适应性。由此得出，从超声波起振情况不同导致传统互相关算法传播延时估计误差问题入手，通过理论和实测，证实了实际测量中测量条件比如液体浓度和测量距离导致起振速度不同，应用传统互相关算法引起较大的测量误差。在此基础上，本发明的互相关精确相位步进和整周期匹配法，简化了测量设备在每次进行超声波传播延时测量前的测量标定，给出了较为精确的测量结果。同时本方法时间复杂度仅为O(n3)，方便于在嵌入式微处理器实现，极大降低了超声波测量设备的制造成本。
以上内容是结合具体的/优选的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，其还可以对这些已描述的实施例做出若干替代或变型，而这些替代或变型方式都应当视为属于本发明的保护范围。

Claims (10)

1.一种基于互相关信号相位分解的时延估计方法，其特征在于：选取时长为Ts的理想发射信号，获取接收信号，对接收信号与理想发射信号进行互相关运算得到互相关信号，将接收信号传播延时分解为K个整周期和ψ的相位偏移；设置单周期函数信号，求得互相关信号与所述单周期函数信号间的相关系数的极值，得到相位偏移ψ；将单周期函数信号再往0时刻方向进行整周期移相，求得互相关信号的第一个起振点，从而得到接收信号的时延，所述单周期函数与理想发射信号函数为相同类型的函数，且角频率相同；
设：
Figure FDA0003537989720000011
其中Kx为待求解周期整数，
Figure FDA0003537989720000012
表示从0时刻开始最接近Γx的整周期点，ψ为待求解自该整周期点的偏移相位差；记Ks为Ts包含的整周期个数，则：
Figure FDA0003537989720000013
构造单周期函数信号：
Figure FDA0003537989720000014
式中θ为待求解相位偏移自变量，ω为超声激励信号的角频率,Km为互相关信号f(τ)的序列f(mT)序列最大值处零点方向最近的整周期点。
2.根据权利要求1所述时延估计方法，其特征在于：包括以下步骤：
S1、对接收信号与理想发射信号进行互相关运算得到互相关信号，构建单周期函数信号，并对单周期函数信号进行步进移相，求解互相关信号和单周期函数信号的相关系数的极值，极值处相位即为接收信号传播延时的精确相位偏移ψ；
S2、对求得的精确相位偏移为ψ的单周期函数信号，往0时刻方向进行整周期移相，确定具体的移相周期数K，求得互相关信号的第一个起振位置点，得到接收信号传播延时。
3.根据权利要求2所述时延估计方法，其特征在于：包括以下步骤：所述接收信号传播延时等于互相关信号的第一个起振位置点与所述理想发射信号时长的差。
4.根据权利要求2所述时延估计方法，其特征在于：步骤S1中，包括以下步骤：
A1、选取时长为Ts的理想发射信号x(t)；
A2、获取接收波信号y(t)；
A3、进行互相关运算，生成互相关信号f(τ)的离散形式序列f(mT)；
A4、寻找f(mT)为最大值时的Tmax；
A5、获取离Tmax最近的整周期点，记为Km；
A6、构建步进移相的单周期函数信号g(τ)序列g((mT)θ(i)(j))；
A7、求取序列g((mT)θ(i)(j))与步进移相的互相关信号序列f((mT)θ(i)(j))的相关系数ρgf(Km，θ(i))；
A8、求取ρgf(Km，θ(i))最小值时的θ(i)，即为精确相位值ψ；
其中，步骤A1与步骤A2无先后顺序区别。
5.根据权利要求3所述时延估计方法，其特征在于：所述互相关信号，如下式所示：
Figure FDA0003537989720000021
式中，τ表示互相关信号时域自变量；
所述单周期函数信号g(τ)与互相关信号f(τ)的相关系数ρgf(θ)，如下式所示：
Figure FDA0003537989720000022
式中，
Figure FDA0003537989720000023
计算中的方法如下：
Figure FDA0003537989720000024
6.根据权利要求3所述时延估计方法，其特征在于：所述互相关信号离散形式序列f(mT)计算方式如下：
Figure FDA0003537989720000031
其中，x(t)和y(t)由同一个脉冲触发信号P(nT)进行数字化采样，分别形成x(nT)和y(nT)：
Figure FDA0003537989720000032
式中，Nx为接收信号起振后第一个采样点序列值，Ns为发射激励信号最后一个采样点序列值，Nr为接收信号最后一个采样点序列值；x(t)表示激励信号；y(t)表示接收波信号，Tx为待求解的飞行时间，c表示换能器起振参数。
7.根据权利要求5所述时延估计方法，其特征在于：
选取步进量Δθ＝θ(i)-θ(i-1)，寻找互相关信号离散形式序列f(mT)中满足下式的所有位置点，记为((mT)θ(i)(j))，
Figure FDA0003537989720000033
代入单周期函数信号g(τ)，得步进移相的单周期函数信号序列g((mT)θ(i)(j))。
8.根据权利要求2所述时延估计方法，其特征在于：步骤S2中，对在极值ψ处的单周期正弦波信号进行整周期移相，确定具体的移相周期数K，包括以下步骤：
B1、从Km开始，将
Figure FDA0003537989720000034
往0时刻方向进行整周期移相至k个周期；
B2、寻找互相关信号离散形式序列f(mT)中满足下式的所有点位置序列：
Figure FDA0003537989720000035
B3、求取所述所有点位置序列的整周期移相的单周期函数序列g((mT)k(j))；
B4、获取所述所有点位置序列的整周期移相的互相关信号离散形式序列f((mT)k(j))；
B5、计算整周期移相的单周期函数序列g((mT)k(j))与整周期移相的互相关信号离散形式序列f((mT)k(j))的相关系数
Figure FDA0003537989720000041
B6、相关系数
Figure FDA0003537989720000042
是否小于设定值ε，若是，转B8；若否，进入下一步；
B7、k值减1，转B1；
B8、判定k-1为互相关信号周期Kx，根据当前k-1计算互相关信号的时延估计结果Tx。
9.根据权利要求8所述时延估计方法，其特征在于：
选择k＝km，km-1，km-2…，将
Figure FDA0003537989720000043
往0时刻方向进行整周期移相，寻找f(mT)序列中满足
Figure FDA0003537989720000044
的所有点位置序列，记为(mT)k(j)，代入下式
Figure FDA0003537989720000045
得到整周期移相的单周期函数序列g((mT)k(j))；
求解：
Figure FDA0003537989720000046
Kx＝Kmin，从而得到y(t)的开始起振点：
Figure FDA0003537989720000047
式中，Ks为Ts包含的整周期个数，Ts为进行互相关计算的理想发射信号时间长度。
10.根据权利要求1所述时延估计方法，其特征在于：包括如下步骤：
D1、根据理想发射信号与接收信号求得互相关信号；
D2、对单周期函数信号进行步进移相，得到步进移相单周期函数信号：
D3、求取各步进移相单周期函数信号与所述互相关信号的相关系数序列最大值；
D4、求解相关系数序列最大值处的精确相位偏移；
D5、将精确相位偏移单周期信号往零时刻方向进行整周期移相；
D6、求取互相关信号与整周期移相后的相关系数，如果所述相关系数小于设定阈值，上一个整周期即为互相关信号起振周期；此处的相关系数为每个整周期位置处的相关系数；
D7、根据互相关信号起振周期、理想发射信号时长、精确相位偏移，求得接收信号实际传播延时。
CN201811564624.0A 2018-12-20 2018-12-20 一种基于互相关信号相位分解的时延估计方法 Active CN111352099B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811564624.0A CN111352099B (zh) 2018-12-20 2018-12-20 一种基于互相关信号相位分解的时延估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811564624.0A CN111352099B (zh) 2018-12-20 2018-12-20 一种基于互相关信号相位分解的时延估计方法

Publications (2)

Publication Number Publication Date
CN111352099A CN111352099A (zh) 2020-06-30
CN111352099B true CN111352099B (zh) 2022-05-10

Family

ID=71192087

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811564624.0A Active CN111352099B (zh) 2018-12-20 2018-12-20 一种基于互相关信号相位分解的时延估计方法

Country Status (1)

Country Link
CN (1) CN111352099B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112698333B (zh) * 2021-03-24 2021-07-06 成都千嘉科技有限公司 一种适用于气体和液体的超声波飞行时间测定方法及系统
CN113489663B (zh) * 2021-07-06 2022-06-07 武汉大学 甚低频人工源信号的时延提取方法
CN114283857B (zh) * 2021-12-16 2024-05-28 上海艾为电子技术股份有限公司 分频信号的延时补偿、分频方法、系统和分频器

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103197318A (zh) * 2013-03-18 2013-07-10 中国科学院声学研究所 一种基于Pattern时延编码水下声定位的时延估计方法
CN107219376A (zh) * 2017-05-27 2017-09-29 华北电力大学 一种适应对象运动特性的互相关测速方法

Family Cites Families (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6594595B2 (en) * 2001-04-03 2003-07-15 Advantest Corporation Apparatus for and method of measuring cross-correlation coefficient between signals
DE102005015456A1 (de) * 2005-04-04 2006-10-05 Viasys Healthcare Gmbh Verfahren zur Bestimmung der zeitlichen Lage eines Wellenpakets sowie Flussmessgerät
EP1949558B1 (en) * 2005-10-05 2015-01-14 Telecom Italia S.p.A. Method and system for multiple antenna communications, related apparatus and corresponding computer program product
CN101315397A (zh) * 2008-06-25 2008-12-03 中国海洋石油总公司 一种幅相测量方法
CN101388001A (zh) * 2008-06-25 2009-03-18 天津大学 基于全相位fft的高精度瞬间相位估计方法
US9689726B2 (en) * 2012-07-09 2017-06-27 Texas Instruments Incorporated Flow meter
DE102012212901A1 (de) * 2012-07-24 2014-01-30 Continental Automotive Gmbh Verfahren und Vorrichtung zum Betreiben einer akustischen Messvorrichtung
CN104459461B (zh) * 2014-11-20 2017-09-12 云南电网公司大理供电局 一种基于工频正弦拟合相关系数的故障数据自动筛选方法
CN104501889B (zh) * 2015-01-23 2018-05-01 中煤科工集团重庆研究院有限公司 基于互相关时差法超声波流量的检测方法
US10309813B2 (en) * 2015-05-15 2019-06-04 Reliance Worldwide Corporation Method and system for fluid flow rate measurement
CN105572650A (zh) * 2015-12-15 2016-05-11 宁波大学 一种宽带复相关流速测量方法
CN106772393B (zh) * 2016-12-14 2019-08-02 湖北工业大学 一种改进的基于飞行时间检测的超声波测距方法
CN107870034B (zh) * 2017-10-24 2019-12-24 宁波大学科学技术学院 一种基于相位差的水声声速测量方法
CN108020818A (zh) * 2017-12-07 2018-05-11 东北大学 一种在噪声下基于滑动dft的正弦脉冲信号测距方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103197318A (zh) * 2013-03-18 2013-07-10 中国科学院声学研究所 一种基于Pattern时延编码水下声定位的时延估计方法
CN107219376A (zh) * 2017-05-27 2017-09-29 华北电力大学 一种适应对象运动特性的互相关测速方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
广义二次相关时延估计算法改进;景源等;《辽宁大学学报(自然科学版)》;20180515(第02期);7-13 *

Also Published As

Publication number Publication date
CN111352099A (zh) 2020-06-30

Similar Documents

Publication Publication Date Title
CN111352099B (zh) 一种基于互相关信号相位分解的时延估计方法
US11747181B2 (en) Extended range ADC flow meter
Andria et al. Digital signal processing techniques for accurate ultrasonic sensor measurement
CN105102924A (zh) 对壁表面的变化的超声波检测
WO2012031302A1 (en) Multiphase fluid characterization system
CN110799810B (zh) 用于测量流体速度的方法
US20210003436A1 (en) Time-of-flight generating circuit and chip, flow meter and method of the same
US11340100B2 (en) Method for evaluating measurement data sequences of an ultrasonic flow measuring device and ultrasonic flow measuring device
JPS62169073A (ja) 超音波ドップラー診断装置
Zhou et al. High-accuracy ultrasonic temperature measurement based on MLS-modulated continuous wave
US6422094B1 (en) Method for determining the flow rate and/or the molecular mass of liquid or gaseous media
US11137276B1 (en) All digital travel time flow meter using time reversed acoustics
JP2014516162A (ja) 超音波伝搬時間法による流体の流量検出方法
CN112147366A (zh) 用于使用超声波来测量流体速度的方法
KR100739506B1 (ko) 정합필터의 간략한 계산을 사용한 초음파 거리 정밀측정방법
US10727787B2 (en) Flow detection with quadrature demodulation
Tu et al. Phase and frequency matching-based signal processing method for Coriolis mass flowmeters
CN109632951A (zh) 用于使用超声确定流体类型的方法和组件
RU2436050C1 (ru) Способ определения скорости звука в жидких средах
Yong-jun et al. Fractional Fourier transform of ultrasonic chirp signal for range measurement
WO2005119182A1 (ja) 流体の流量測定方法及び流量測定装置
RU2330298C2 (ru) Способ определения места повреждения линий электропередачи и связи и устройство для его осуществления
Andria et al. Digital Methods for very accurate Ultrasonic sensor measurements
EP4067833A1 (en) All digital travel time flow meter using time reversed acoustics
Battaglini et al. Frequency modulated continuous wave ultrasonic radar

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant