CN103412189B - 一种用于电学层析成像系统的信息滤波解调方法 - Google Patents

Abstract

本发明涉及一种用于电学层析成像系统的信息滤波解调方法，其特征在于包括下列步骤：（1）根据系统的激励频率f和采样频率f_s建立递归方程；（2）以激励信号零相位时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据计算出误差相关矩阵的逆矩阵P^-1(k)和递归中间矩阵K的初始值；（3）将其初始值及新增采样点数据代入递归方程，更新P^-1(k)和K，并同时计算数据矩阵B_k；（4）判断解调结果是否满足系统精度要求，若不满足则返回步骤（3），若满足则结束递归过程，并输出测量信号的幅值和相位数据。上述递归解调方法最少只需2个采样点数据，且增加采样点数可提高解调精度和抗噪声能力。基于信息滤波算法的递归解调方法具有可以在无初值情况下启动、迭代简单、利于硬件实现、占用资源少的特点，具有广阔的应用前景。

Description

一种用于电学层析成像系统的信息滤波解调方法

技术领域

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

背景技术

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

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

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

信息滤波算法是1967年Fraser提出的滤波算法，是卡尔曼滤波器的一种实现形式，这种方法强调了滤波过程的递归最小二乘特性。基于信息滤波算法的递归解调方法可以在无初始信息情况下启动递归过程，且递归形式简单。

发明内容

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

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

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

其中，k为采样点序数且k≥0；x_k为第k个采样点数据；P^-1(k)是误差相关矩阵的逆矩阵；K(k)＝P^-1(k)B_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(2\right)$

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

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

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

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

步骤二、以激励信号相位为零的时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程计算出误差相关矩阵的逆矩阵P^-1(k)和递归中间矩阵K的初始值P^-1(1)和K(1)；

步骤三、将步骤二中计算出的误差相关矩阵的逆矩阵P^-1(k)和递归中间变量K(k)的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递归方程中，逐步更新误差相关矩阵的逆矩阵P^-1(k)和递归中间变量K(k)的值；

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递归过程，并根据以下公式：

${\hat{B}}_k{\left[P^{-1}\left(k\right)\right]}^{-1}K\left(k\right)---\left(5\right)$

计算出数据矩阵B的最终结果进而得到测量信号的幅值和相位

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

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

附图说明

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

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

图3为本发明所提供解调方法的仿真实验结果。

具体实施方式

本发明，即一种用于电学层析成像系统的基于信息滤波算法的递归解调方法，包括下列步骤：

步骤一、根据系统的激励频率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\left(0\right)=V_0,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^{-1}(k+1)=V{(k+1)}^{H}V(k+1)$

$={\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]---\left(13\right)$

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

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

方程(13)为信息滤波算法的第一个递归表达式。

在方程(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(14\right)$

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

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

方程(14)可变形为：

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

设递归中间变量K(k)＝P^-1(k)B_k，代入方程(15)，有：

K(k+1)＝K(k)+V_k+1 ^HV_k+1                                   (16)

方程(16)即为信息滤波算法的第二个递归表达式，这里是对递归中间变量K(k)＝P^-1(k)B_k进行递归，而不是对估计本身。

由方程(13)和(16)，可知递归的公式为：

$\left\{\begin{array}{c}{P}^{-1}(k+1)=P^{-1}\left(k\right)-{V_{k+1}}^H{V}_{k+1}\\ K(k+1)=K\left(k\right)+{V_{k+1}}^H{x}_{k+1}\end{array}\right.,k\≥1---\left(17\right)$

而其初始值可由：

$\left\{\begin{array}{c}{P}^{-1}\left(1\right)=V{\left(1\right)}^H\·V\left(1\right)\\ K\left(1\right)=V{\left(1\right)}^H\·X\left(1\right)\end{array}\right.---\left(18\right)$

计算得到。

步骤二、以激励信号相位为零的时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程计算出误差相关矩阵的逆矩阵P^-1(k)和递归中间矩阵K的初始值P^-1(1)和K(1)；

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

$\left\{\begin{array}{c}{P}^{-1}\left(1\right)=V{\left(1\right)}^{H}V\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}{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]\\ K\left(1\right)=V{\left(1\right)}^{H}X\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]\end{array}\right.---\left(19\right)$

步骤三、将步骤二中计算出的误差相关矩阵的逆矩阵P-¹(k)和递归中间变量K(k)的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递归方程中，逐步更新误差相关矩阵的逆矩阵P^-1(k)和递归中间变量K(k)的值；

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递归过程，并根据以下公式：

${\hat{B}}_k={\left[P^{-1}\left(k\right)\right]}^{-1}K\left(k\right)---\left(20\right)$

计算出数据矩阵B的最终结果进而得到测量信号的幅值和相位

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

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

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

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

（2）测量信号为在幅值为1、相位为30°、频率为100kHz的理想正弦信号中加入幅值为0.01的随机噪声的非理想信号，系统采样频率为6MHz，每个信号周期内有60个采样点。

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

由仿真实验结果可以看出，在理想正弦信号的情况下，本发明所提供的递归解调方法利用初始的两个采样点即可计算出标准无偏差的幅值和相位解调结果，且随着采样点数的增多解调结果无变化；而在测量信号含有随机噪声的情况下，本发明所提供的递归解调方法利用初始的两个采样点得到的初始幅值和相位解调结果的相对误差在迭代10次左右后，均小于0.5%，且随着采样点数的增多解调相对误差呈现衰减的趋势。仿真实验验证了本发明所提供的递归解调方法的良好效果。

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

Claims (1)

1.一种用于电学层析成像系统的信息滤波解调方法，其递归过程可在无初始信息情况下启动，且递归形式简单、利于硬件实现、占用资源少，其特征在于包括下列步骤：

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

其中，k为采样点序数且k≥0；x_k为第k个采样点数据；P^-1(k)是误差相关矩阵的逆矩阵；K(k)＝P^-1(k)B_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(2\right)$

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

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

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

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

步骤二、以激励信号相位为零的时刻为起始时刻，根据测量信号在紧邻起始时刻之后的2个采样点数据，按照步骤一中的初始方程计算出误差相关矩阵的逆矩阵P^-1(k)和递归中间矩阵K的初始值P^-1(1)和K(1)；

步骤三、将步骤二中计算出的误差相关矩阵的逆矩阵P^-1(k)和递归中间变量K(k)的初始值以及新增的第k(k≥2)个采样点数据x_k代入到步骤一中建立的递归方程中，逐步更新误差相关矩阵的逆矩阵P^-1(k)和递归中间变量K(k)的值；

步骤四、判断解调结果是否满足系统精度要求，若不满足则返回步骤三；若满足则停止递归过程，并根据以下公式：

${\hat{B}}_k={\left[P^{-1}\left(k\right)\right]}^{-1}K\left(k\right)---\left(5\right)$

计算出数据矩阵B的最终结果进而得到测量信号的幅值和相位

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