CN102565541B - 一种用于电学层析成像系统的递推解调方法 - Google Patents

Abstract

本发明涉及一种用于电学层析成像系统的递推解调方法，其特征在于：其包括下列步骤：(1)根据系统的激励频率f和采样频率f_s建立初始方程和递推方程；(2)以激励信号零相位时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据按照步骤(1)中的初始方程计算出数据矩阵B和状态矩阵P的初始值；(3)将其初始值及新增采样点数据代入步骤(1)中的递推方程，以更新B和P；(4)判断解调结果是否满足系统精度要求，若不满足则返回步骤(3)，若满足则结束解调，并输出测量信号的幅值和相位数据。上述递推解调方法最少只需两个采样点数据，且增加采样点数可提高解调精度和抗噪声能力。

Description

一种用于电学层析成像系统的递推解调方法

技术领域

本发明涉及一种用于电学层析成像系统的递推解调方法，属于分布参数测量领域。

背景技术

电学层析成像的工作原理是根据不同的介质具有不同的电特性(电导率/介电常数/复导纳)，通过获取敏感区域边界测量数据并采用适当的图像重建算法，反演出敏感区域内的电参数分布，进而获得区域内介质的分布。电学层析成像系统的典型结构如图1所示。目前，为了获得较好的测量信噪比、线性度和分辨率，电学层析成像系统多采用正弦信号作为激励源，而对测量信号幅值和相位的检测一般采用相敏解调的方法。

由于传统的模拟相敏解调(乘法解调)方法受到低通滤波器建立时间的限制，不适合系统实时性要求较高的应用场合。相比之下，基于数字处理器的数字相敏解调方法以其良好的测量实时性而受到研究者越来越多的关注。数字相敏解调首先利用高速模数转换器对被测信号进行采样，之后利用高性能数字信号处理器件，如FPGA、DSP等，采用数值计算的方法提取被测信号的幅值和相位信息。

目前最为常用的数字相敏解调方法为正交序列解调法，但其要求采样序列长度必须为完整的信号周期，降低了灵活性，同时也限制了解调速度的进一步提高。

发明内容

本发明的目的在于提供一种用于电学层析成像系统的递推解调方法，可在少于一个信号周期的时间内得到精度较高解调结果，且随着代入递推过程的采样点数的增加，可提高解调结果的抗噪声性能。

本发明所提供的一种用于电学层析成像系统的递推解调方法，包括下列步骤：

步骤一、根据系统的激励频率f和采样频率f_s建立初始方程和递推方程：

初始方程： $\left\{\begin{array}{c}P\left(1\right)={[V{\left(1\right)}^H\·V\left(1\right)]}^{-1}\\ B_1=P{\left(1\right)}^{-1}\·{V\left(1\right)}^H\·X\left(1\right)\end{array}\right.---\left(1\right)$

递推方程： $\left\{\begin{array}{c}P(k+1)=P\left(k\right)-\frac{P\left(k\right){V_{k+1}}^H{V}_{k+1}P\left(k\right)}{1+V_{k+1}P\left(k\right){V_{k+1}}^H}\\ B_{k+1}=B_k+P(k+1){V_{k+1}}^H(x_{k+1}-V_{k+1}{B}_k)\end{array}\right.,k\≥1---\left(2\right)$

在方程(1)和方程(2)中，k为采样点序数且k≥0；x_k为第k个采样点数据；P(k)为状态矩阵；B_k为数据矩阵；V_k为由激励信号频率和系统采样频率唯一确定的常数向量，其具体形式为：

$V_k=\left[\begin{array}{cc}{e}^{-k\·j2\mathrm{\πf}/f_s}& e^{k\·j2\mathrm{\πf}/f_s}\end{array}\right]---\left(3\right)$

V(1)为由V₀和V₁构成的辅助矩阵，即

V(k)＝[V₀ V₁]^T                    (4)

X(1)为由前2个采样点数据构成的测量向量，即

X(1)＝[x₀ x₁]^T                    (5)

步骤二、以激励信号相位为零的时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程(1)计算出包含测量信号幅值与相位信息的数据矩阵B和状态矩阵P的初始值B₁和P(1)；

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递推方程(2)中，逐步更新数据矩阵B和状态矩阵P的值；

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果B计算出测量信号的幅值A和相位

$\left\{\begin{array}{c}B={\left[\begin{array}{cc}{b}_1& b_2\end{array}\right]}^T\\ A=2\left|b_1\right|\\ \widehat{\mathrm{\θ}}=\arctan [\mathrm{Im}\left(b_1\right)/\mathrm{Re}\left(b_1\right)]\end{array}\right.---\left(6\right)$

本发明与现有技术相比的优点在于：可先利用从紧邻起始时刻起的2个测量信号采样点数据得到初步的解调结果，之后根据系统的测量精度和实时性等要求适当增加采样点数，通过递推的方法来提高解调精度和抗噪声能力；在无噪声的情况下，该方法不需要采用完整信号周期的采样点即可得到准确的解调结果，在存在噪声的情况下，亦可获得较为理想的解调结果，具有很好的实时性和灵活性。

附图说明

图1为电学层析成像系统的结构框图；

图2为本发明所提供的递推解调方法实施过程的流程图；

图3(a)为仿真试验中测量信号无噪声时的幅值解调结果；

图3(b)为仿真试验中测量信号无噪声时的相位解调结果；

图3(c)为仿真试验中测量信号含噪声时的幅值解调结果；

图3(d)为仿真试验中测量信号含噪声时的相位解调结果。

具体实施方式

本发明，即一种用于电学层析成像系统的递推解调方法，包括下列步骤：

步骤一、根据系统的激励频率f和采样频率f_s建立递推解调方程。

假设测量信号的表示形式为：

x_k＝Acos(2πkf/f_s+θ)                (1)

其中，k为采样点序数，A为测量信号幅值，θ为测量信号与激励信号之间的相位差。根据欧拉公式，x_k可表示为：

$x_k=\frac{A}{2}{e}^{-\mathrm{j\θ}}\·e^{-2\mathrm{\πkf}/f_s}+\frac{A}{2}{e}^{\mathrm{j\θ}}\·e^{2\mathrm{\πkf}/f_s}=\left[\begin{array}{cc}{e}^{-2\mathrm{\πkf}/f_s}& e^{2\mathrm{\πkf}/f_s}\end{array}\right]\left[\begin{array}{c}\frac{A}{2}{e}^{-\mathrm{j\θ}}\\ \frac{A}{2}{e}^{\mathrm{j\θ}}\end{array}\right]---\left(2\right)$

因此，考虑构造解调辅助矩阵V和含有测量信号幅值与相位信息的数据矩阵B

$V={\left[\begin{array}{ccccc}{e}^{-0\·j2\mathrm{\πf}/f_s}& e^{-1\·j2\mathrm{\πf}/f_s}& ...& e^{-k\·j2\mathrm{\πf}/f_s}& ...\\ e^{0\·j2\mathrm{\πf}/f_s}& e^{1\·j2\mathrm{\πf}/f_s}& ...& e^{k\·j2\mathrm{\πf}/f_s}& ...\end{array}\right]}^T---\left(3\right)$

$B=\left[\begin{array}{c}\frac{A}{2}{e}^{-\mathrm{j\θ}}\\ \frac{A}{2}{e}^{\mathrm{j\θ}}\end{array}\right]---\left(4\right)$

则测量信号的采样序列与解调辅助矩阵之间的关系可以表示为：

$\left[\begin{array}{c}{x}_0\\ x_1\\ .\\ .\\ .\\ x_k\\ .\\ .\\ .\end{array}\right]=\left[\begin{array}{cc}{e}^{-0\·2\mathrm{\πjf}/f_s}& e^{0\·2\mathrm{\πjf}/f_s}\\ e^{-1\·2\mathrm{\πjf}/f_s}& e^{1\·2\mathrm{\πjf}/f_s}\\ .& .\\ .& .\\ .& .\\ e^{-k2\mathrm{\πjf}/f_s}& e^{k2\mathrm{\πjf}/f_s}\\ .& .\\ .& .\\ .& .\end{array}\right]\left[\begin{array}{c}\frac{A}{2}{e}^{-\mathrm{j\θ}}\\ \frac{A}{2}{e}^{\mathrm{j\θ}}\end{array}\right]=\mathrm{VB}---\left(5\right)$

令

$V_0=\left[\begin{array}{cc}{e}^{-0\·j2\mathrm{\πf}/f_s}& e^{0\·j2\mathrm{\πf}/f_s}\end{array}\right],$ $V_1=\left[\begin{array}{cc}{e}^{-1\·j2\mathrm{\πf}/f_s}& e^{1\·j2\mathrm{\πf}/f_s}\end{array}\right],$ …， $V_k=\left[\begin{array}{cc}{e}^{-k\·j2\mathrm{\πf}/f_s}& e^{k\·j2\mathrm{\πf}/f_s}\end{array}\right],...\left(6\right)$

V(0)＝V₀， $V\left(1\right)=\left[\begin{array}{c}{V}_0\\ V_1\end{array}\right]=\left[\begin{array}{c}V\left(0\right)\\ V_1\end{array}\right],$ …， $V\left(k\right)=\left[\begin{array}{c}V(k-1)\\ V_k\end{array}\right],...\left(7\right)$

X(0)＝[x₀]，X₁＝[x₀ x₁]^T＝[X(0)  x₁]^T，…，X_k＝[X(k-1)x_k]^T，…      (8)

那么方程(5)可以表示成另一种形式，即

X(k)＝V(k)B                (9)

如果V(k)是方阵且其逆存在，则可通过

B＝V(k)^-1X(k)            (10)

直接计算出数据矩阵B，进而解调出测量信号的幅值和相位。然而，V(k)是一个k+1行2列的矩阵，除k＝1的情况以外均不能用方程(10)求解B。

考虑到矩阵V(k)的特点，重新构造求解数据矩阵B的公式，当k≥1时其可表示为：

B_k＝[V(k)^HV(k)]^-1V(k)^HX(k)            (11)

令

P(k)＝[V(k)^HV(k)]^-1                (12)

则

$P(k+1)={\left[V{(k+1)}^{H}V(k+1)\right]}^{-1}$

$={\left({\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]}^H\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]\right)}^{-1}$

$={\left(\left[\begin{array}{cc}V{\left(k\right)}^H& {V_{k+1}}^H\end{array}\right]\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]\right)}^{-1}---\left(13\right)$

$={(V{\left(k\right)}^{H}V\left(k\right)+{V_{k+1}}^H{V}_{k+1})}^{-1}$

$={(P{\left(k\right)}^{-1}+V_{k+1}^H{V}_{k+1})}^{-1}$

根据矩阵求逆的定理，有

(A+BCD)^-1＝A^-1-A^-1B(C^-1+DA^-1B)^-1DA^-1        (14)

令A＝P(k)^-1，

C＝1，D＝V_k+1，可得

P(k+1)＝(A+B·1·B^H)^-1

＝A^-1-A^-1B(1+B^HA^-1B)^-1B^HA^-1                (15)

＝P(k)-P(k)V_k+1 ^H(1+V_k+1P(k)V_k+1 ^H)^-1V_k+1P(k)

因此在方程(11)中，若用第k次计算的结果更新第k+1次计算的相关参数，则

$B_{k+1}=P(k+1)V{(k+1)}^{H}X(k+1)$

$=P(k+1){\left[\begin{array}{c}V\left(k\right)\\ V_{k+1}\end{array}\right]}^H\left[\begin{array}{c}X\left(k\right)\\ x_{k+1}\end{array}\right]$

$=P(k+1)\left[\begin{array}{cc}{V\left(k\right)}^H& {V_{k+1}}^H\end{array}\right]\left[\begin{array}{c}X\left(k\right)\\ x_{k+1}\end{array}\right]---\left(16\right)$

$=P(k+1)[V{\left(k\right)}^{H}X\left(k\right)+{V_{k+1}}^H{x}_{k+1}]$

$=P(k+1)P{\left(k\right)}^{-1}{B}_k+P(k+1){V_{k+1}}^H{x}_{k+1}$

将式(13)变换为

P(k+1)^-1＝P(k)^-1+V_k+1 ^HV_k+1                (17)

从而，

P(k)^-1＝P(k+1)^-1-V_k+1 ^HV_k+1                (18)

将式(18)代入式(16)，可得

B_k+1＝P(k+1)[P(k+1)^-1-V_k+1 ^HV_k+1]B_k+P_k+1V_k+1 ^Hx_k+1

＝B_k-P(k+1)V_k+1 ^HV_k+1B_k+P(k+1)V_k+1 ^Hx_k+1            (19)

＝B_k+P(k+1)V_k+1 ^H(x_k+1-V_k+1B_k)

由式(15)和(19)可知，递推方程为

$\left\{\begin{array}{c}P(k+1)=P\left(k\right)-\frac{P\left(k\right){V_{k+1}}^H{V}_{k+1}P\left(k\right)}{1+V_{k+1}P\left(k\right){V_{k+1}}^H}\\ B_{k+1}=B_k+P(k+1){V_{k+1}}^H(x_{k+1}-V_{k+1}{B}_k)\end{array}\right.,k\≥1---\left(20\right)$

而其初始值可由初始方程

$\left\{\begin{array}{c}P\left(1\right)={[V{\left(1\right)}^H\·V\left(1\right)]}^{-1}\\ B_1=P{\left(1\right)}^{-1}\·{V\left(1\right)}^H\·X\left(1\right)\end{array}\right.---\left(21\right)$

计算得到。

步骤二、以激励信号相位为零的时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程(21)计算出包含测量信号幅值与相位信息的数据矩阵B和状态矩阵P的初始值B₁和P(1)。

令紧邻起始时刻之后的2个测量信号采样点的值为x₀和x₁根，据方程(21)计算得到：

$\left\{\begin{array}{c}P\left(1\right)={\left[V{\left(1\right)}^{H}V\left(1\right)\right]}^{-1}={\left({\left[\begin{array}{cc}{e}^{-0\·2\mathrm{\πjf}/f_s}& e^{0\·2\mathrm{\πjf}/f_s}\\ e^{-1\·2\mathrm{\πjf}/f_s}& e^{1\·2\mathrm{\πjf}/f_s}\end{array}\right]}^H\left[\begin{array}{cc}{e}^{-0\·2\mathrm{\πjf}/f_s}& e^{0\·2\mathrm{\πjf}/f_s}\\ e^{-1\·2\mathrm{\πjf}/f_s}& e^{1\·2\mathrm{\πjf}/f_s}\end{array}\right]\right)}^{-1}\\ B_1=P\left(1\right)V{\left(1\right)}^{H}X\left(1\right)=P\left(1\right){\left[\begin{array}{cc}{e}^{-0\·2\mathrm{\πjf}/f_s}& e^{0\·2\mathrm{\πjf}/f_s}\\ e^{-1\·2\mathrm{\πjf}/f_s}& e^{1\·2\mathrm{\πjf}/f_s}\end{array}\right]}^H\left[\begin{array}{c}{x}_0\\ x_1\end{array}\right]=\left[\begin{array}{c}{b}_{11}\\ b_{12}\end{array}\right]\end{array}\right.---\left(22\right)$

此时，可根据初始数据矩阵B₁计算出测量信号的幅值和相角：

$\left\{\begin{array}{c}{A}_1=2\left|b_{11}\right|\\ {\mathrm{\θ}}_1=\arctan [\mathrm{Im}\left(b_{11}\right)/\mathrm{Re}\left(b_{11}\right)]\end{array}\right.---\left(23\right)$

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递推方程(20)中，逐步更新数据矩阵B和状态矩阵P的值。

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果B计算出测量信号的幅值A和相位

$\left\{\begin{array}{c}B={\left[\begin{array}{cc}{b}_1& b_2\end{array}\right]}^T\\ A=2\left|b_1\right|\\ \widehat{\mathrm{\θ}}=\arctan [\mathrm{Im}\left(b_1\right)/\mathrm{Re}\left(b_1\right)]\end{array}\right.---\left(24\right)$

下面结合附图和实施例对本发明做进一步详细说明。

在Matlab计算软件中对本发明所提供的一种用于电学层析成像系统的递推解调方法进行仿真实验。实验条件如下：

(1)测量信号为理想正弦信号，其幅值为1、相位为30°、频率为100kHz，系统的采样率为2MHz，每个信号周期内有20个采样点；

(2)测量信号为在幅值为1、相位为30°、频率为100kHz的理想正弦信号中加入均值为0、标准差为0.01的随机噪声的非理想信号，系统采样频率为2MHz，每个信号周期内有20个采样点。

实验条件(1)下测量信号幅值和相位的递推解调结果分别如图3中(a)、(b)所示，实验条件(2)下测量信号幅值和相位的递推解调结果分别如图3中(c)、(d)所示。

由仿真实验结果可以看出，在理想正弦信号的情况下，本发明所提供的递推解调方法利用初始的2个采样点即可解调出准确的幅值和相位，且随着采样点数的增多解调出的幅值和相位仍然准确；而在测量信号含有均方值为随机噪声的情况下，本发明所提供的递推解调方法利用初始的2个采样点得到的初始幅值和相位解调结果的相对误差均小于1％，且随着采样点数的增多解调相对误差呈现震荡衰减的趋势。

仿真实验验证了本发明所提供的递推解调方法的良好效果。

以上对本发明及其实施方式的描述，并不局限于此，附图中所示仅是本发明的实施方式之一。在不脱离本发明创造宗旨的情况下，不经创造性地设计出与该技术方案类似的结构或实施例，均属本发明保护范围。

Claims (1)

1.一种用于电学层析成像系统的递推解调方法，其特征在于：其包括下列步骤： 

步骤一、根据系统的激励频率f和采样频率f_s建立初始方程和递推方程： 

初始方程：

递推方程：

在方程(1)和方程(2)中，k为采样点序数且k≥0；x_k为第k个采样点数据；P(k)为状态矩阵；B_k为数据矩阵；V_k为由激励信号频率和系统采样频率唯一确定的常数向量，其具体形式为： 

V(1)为由V₀和V₁构成的辅助矩阵，即 

V(1)=[V₀ V₁]^T   (4) 

X(1)为由前2个采样点数据构成的测量向量，即 

X(1)=[x₀ x₁]^T   (5) 

步骤二、以激励信号相位为零的时刻为起始时刻,根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程(1)计算出包含测量信号幅值与相位信息的数据矩阵B和状态矩阵P的初始值B₁和P(1)； 

步骤三、将步骤二中计算出的数据矩阵B和状态矩阵P的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递推方程(2)中，逐步更新数据矩阵B和状态矩阵P的值； 

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递推过程，并根据数据矩阵B的最终结果

计算出测量信号的幅值和相位

Images