CN103340625B - 一种电阻抗断层成像中快速最优的正则化方法 - Google Patents

Abstract

本发明公开了一种电阻抗断层成像中快速最优的正则化方法。该方法通过建立稀疏和密集两种剖分模型，采用稳健有效的正则参数最优化方法求得稀疏模型的最优正则化参数，然后将稀疏模型的最优正则化参数作为密集模型的正则化参数初值，将稀疏模型的残差作为密集模型下残差的估计，运用迭代类正则化方法求取密集模型的最优正则化参数。该方法实现了密集剖分模型的快速最优正则化，保证了成像的速度和精度。

Description

一种电阻抗断层成像中快速最优的正则化方法

技术领域

本发明属于电阻抗断层成像技术领域，涉及一种电阻抗断层成像中快速最优的正则化方法。

背景技术

电阻抗断层成像通过在被测目标某一断层表面施加一定的交流电流，并测量相应检测电极上的边界电压，然后根据一定的重建算法重构出目标内部电阻率分布图像——静态EIT图像，或电阻率变化的分布图像——动态EIT图像。中国专利申请（专利号：ZL99115885.5），公开了名称为一种电阻抗断层成像方法，对EIT成像技术方案进行了详细披露。该方法所述的奇异值分解法近年来已被其他正则化方法如阻尼最小二乘法所替代，但在阻尼最小二乘法中正则化参数的选取一直是一个难题，因为对于任何的正则参数α＞0，正则解x^α是数值稳定的，参数α越大正则解越稳定，但是从逼近的角度来讲，为了使残差||Ax_α-u_δ||很小，α应该越小越好。如何兼顾“好的逼近”和“好的数值稳定性”成为改进图像质量的关键。而目前已有的正则参数方法具有运算速度慢、对噪声敏感等缺点，所以我们迫切地需要设计一种快速、稳健、有效的正则化参数最优化方法。

发明内容

本发明解决的问题在于提供一种电阻抗断层成像中快速最优的正则化方法，该方法可以通过快速求取最优正则化参数，提高电阻抗断层成像的图像质量和成像速度。

本发明是通过以下技术方案来实现：

一种电阻抗断层成像中快速最优的正则化方法，包括以下操作：

1）根据待测体的基本信息，分别采用有限元方法剖分建立稀疏剖分模型和密集剖分模型的两种重构模型；

2）采用稳健有效的正则参数最优化方法求取稀疏剖分模型的最优正则化参数和残差；

3）将稀疏剖分模型的最优正则化参数作为密集剖分模型的正则化参数初值，将稀疏剖分模型的残差作为密集剖分模型下残差的估计；采用广义偏差原理求取密集剖分模型的最优正则化参数，在该最优正则化参数下用阻尼最小二乘法重构得到电阻抗断层成像的图像。

所述待测体的基本信息包括待测体边界信息、待测体内部结构信息和电极位置信息；

提取待测体的基本信息，将其转换成电阻抗成像所需的先验信息，然后结合所得的先验信息建立稀疏和密集两种剖分模型。

所述建立图像重建模型时的稀疏剖分模型和密集剖分模型，其依据的剖分形式是：电压变化越剧烈的区域的剖分单元规模越小。

所述所依据的剖分形式是：

区域中心的剖分单元最大，区域边缘的剖分单元较小，靠近电极位置的剖分单元最小。

所述重构稀疏模型的图像时采用的稳健有效的正则化参数最优化方法为无需估计误差的正则化参数最优化方法，包括先验选取正则化参数、L型曲线法、广义交叉检验法、最优分辨率法、拟最优准则法和CRESO准则法。

所述获取稀疏剖分模型的最优正则化参数和残差包括以下操作：

候选正则化参数是从小到大的等比数列，对不同的正则化参数α，采用阻尼最小二乘法求解稀疏剖分模型下的逆问题，求解公式为：

x=(J^TJ+αJ^TJ)^-1J^TΔU

其中J为稀疏剖分模型的雅克比矩阵，重建得到每个正则化参数对应的电阻率分布||x_α||和残差||Ax_α-U_δ||；

以所有候选正则化参数的电阻率分布的模||x_α||和残差||Ax_α-U_δ||分别为横坐标和纵坐标作图得到一条曲线，该曲线曲率最大值点对应的即为最优正则化参数；

在该最优参数下求解稀疏剖分下的逆问题，得到残差δ₀=||Ax_α-U_δ||²。

所述获取密集剖分模型的最优正则化参数是将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值α_n，然后运用广义偏差原理迭代公式迭代求解。

所述的密集剖分模型的最优正则化参数的求解包括以下操作：

①将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值α_n；

②寻找候选正则化参数中与α_n最接近的数作为第i次迭代的初值α(i)，并得到α(i)左右的两个点α(i+1)和α(i-1)；

③应用阻尼最小二乘法求取α(i-1)、α(i)、α(i+1)3个正则化参数下的残差||Ax_α-u_δ||²为p(i-1)，p(i)，p(i+1);

④应用偏差原理迭代公式求取第n次迭代正则化参数：

${\mathrm{\α}}_{n+1}={\mathrm{\α}}_n-\frac{2p\left({\mathrm{\α}}_n\right)}{p^{\′}\left({\mathrm{\α}}_n\right)+\mathrm{sign}\left(p^{\′}\left({\mathrm{\α}}_n\right)\right){(p^{\′}{\left({\mathrm{\α}}_n\right)}^2-2p\left({\mathrm{\α}}_n\right)p^{\′\′}\left({\mathrm{\α}}_n\right))}^{\frac{1}{2}}}$

式中p表示求一阶导数，p″表示求二阶导数，sign(·)为符号函数；p(α_n)用稀疏模型在最优正则化参数下的残差作近似，p(α)关于正则化参数α(i)的一阶导数 $p^{\′}\left({\mathrm{\α}}_n\right)=\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}*\frac{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)}+\frac{p(i+1)-p\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}*\frac{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)},$ p(α)关于α的二阶导数 $p^{\′\′}\left({\mathrm{\α}}_n\right)=(p^{\′}\left(i\right)-\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)})*2/(\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)).$

⑤重复②～④过程3～5次，得到密集剖分模型的最优正则化参数α_best。

所述在最优正则化参数α_best下用阻尼最小二乘法重构并显示图像。

与现有技术相比，本发明具有以下有益的技术效果：

本发明提供的电阻抗断层成像中快速最优的正则化方法，其核心是通过提取成像目标的内部结构信息，将其转换成电阻抗成像所需的先验信息，结合所得的先验信息建立稀疏和密集两种剖分模型，将稀疏模型下的最优正则化参数作为密集模型的正则化参数初值，将稀疏模型下的误差作为密集模型下的误差估计，采用广义偏差原理求取密集模型的最优正则化参数，最终在该最优正则化参数下用阻尼最小二乘法重构得到图像。该方法实现了密集剖分模型的快速最优正则化，保证了成像的精度和速度。

本发明提供的电阻抗断层成像中快速最优的正则化方法，能够实时获取当前测量信号下的最优图像，提高了成像精度和成像速度。

附图说明

图1是实施例的仿真边界电压变化的剖分模型。

图2是实施例的内部阻抗变化的效果图。

图3是实施例的稀疏剖分模型。

图4是实施例的密集剖分模型。

图5是实施例的0.1%噪声下本方法的重建结果。

图6是实施例的0.1%噪声下L型曲线法的重建结果。

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

本发明给出一种电阻抗断层成像中快速最优的正则化方法，包括以下操作：

1）根据待测体的基本信息，分别采用有限元方法剖分建立稀疏剖分模型和密集剖分模型的两种重构模型；

2）采用稳健有效的正则参数最优化方法求取稀疏剖分模型的最优正则化参数和残差；

3）将稀疏剖分模型的最优正则化参数作为密集剖分模型的正则化参数初值，将稀疏剖分模型的残差作为密集剖分模型下残差的估计；采用广义偏差原理求取密集剖分模型的最优正则化参数，在该最优正则化参数下用阻尼最小二乘法重构得到电阻抗断层成像的图像。

该方法实现了密集剖分模型的快速最优正则化，保证了成像的精度和速度。

具体的，所述的待测体的基本信息包括待测体边界信息、待测体内部结构信息和电极位置信息；

提取待测体的基本信息，将其转换成电阻抗成像所需的先验信息，然后结合所得的先验信息建立稀疏和密集两种剖分模型。

建立图像重建模型时的稀疏剖分模型和密集剖分模型，其依据的剖分形式是：电压变化越剧烈的区域的剖分单元规模越小。具体的所依据的剖分形式是：

区域中心的剖分单元最大，区域边缘的剖分单元较小，靠近电极位置的剖分单元最小。

所述重构稀疏模型的图像时采用的稳健有效的正则化参数最优化方法为无需估计误差的正则化参数最优化方法，包括先验选取正则化参数、L型曲线法、广义交叉检验法、最优分辨率法、拟最优准则法和CRESO准则法。

在下面的实施例中，假定待测体为阻抗均匀分布的圆域，具体实施步骤如下：

（1）仿真获得内部阻抗分布变化引起的测量电压变化：采用二维圆域作为重建模型，采用图1的剖分模型进行边界电压仿真，剖分后得到节点和单元信息，对单元格的电阻率初值设为默认值1Ω·m。计算图1所示剖分的刚度矩阵，根据初始阻抗分布及刚度矩阵求解正问题，得初始边界电压分布U₁；参见图2，提高圆域中一部分单元的电阻率为2Ω·m，求解正问题，得边界电压分布U₂，并加入-60dB的高斯白噪声于U₂得U₃，标准化差分得到标准化的测量电压变化

（2）得到稀疏剖分下的最优正则化参数和残差：图3为本实施例的稀疏剖分模型。候选正则化参数是以10^0.1为比从0.00001到100的等比数列，对不同的正则化参数α，采用阻尼最小二乘法求解较稀疏剖分下的逆问题，求解公式为：

x=(J^TJ+αJ^TJ)^-1J^TΔU

其中J为模型的雅克比矩阵。重建得到每个正则化参数对应的电阻率分布||x_α||和残差||Ax_α-U_δ||。

以所有候选正则化参数的电阻率分布的模||x_α||和残差||Ax_α-U_δ||分别为横坐标和纵坐标作图，可以得到一条曲线，一般来说，该曲线呈L型，该曲线曲率最大值点对应的即为最优正则化参数。在该最优参数下求解稀疏剖分下的逆问题，得到残差δ₀=||Ax_α-U_δ||²。

（3）求得密集剖分下的最优正则化参数：图4为本实施例的密集剖分模型。将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值，运用广义偏差原理迭代公式迭代求解：

①将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值α_n。

②寻找以10^0.1为比从0.00001到100的等比数列中与α_n最接近的数作为第i次迭代的初值α(i)，并得到与α(i)相差很小的左右的两个点α（i+1)和α(i-1)，在本实施例中，取α(i+1)=α(i)*10^0.1，α(i-1)=α(i)*10^-0.1。③应用阻尼最小二乘法求取α(i-1)、α(i)、α(i+1)3个正则化参数下的残差||Ax_α-u_δ||²为p(i-1)，p(i)，p(i+1)。

④应用偏差原理迭代公式求取第n次迭代正则化参数：

${\mathrm{\α}}_{n+1}={\mathrm{\α}}_n-\frac{2p\left({\mathrm{\α}}_n\right)}{p^{\′}\left({\mathrm{\α}}_n\right)+\mathrm{sign}\left(p^{\′}\left({\mathrm{\α}}_n\right)\right){(p^{\′}{\left({\mathrm{\α}}_n\right)}^2-2p\left({\mathrm{\α}}_n\right)p^{\′\′}\left({\mathrm{\α}}_n\right))}^{\frac{1}{2}}}$

式中p表示求一阶导数，p″表示求二阶导数，sign(·)为符号函数。p(α_n)用稀疏模型在最优正则化参数下的残差作近似，而残差的一、二阶导数可通过离散化后求差分获得，其中p(α)关于正则化参数α(i)的一阶导数 $p^{\′}\left({\mathrm{\α}}_n\right)=\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}*\frac{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)}+\frac{p(i+1)-p\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}*\frac{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)},$ p(α)关于α的二阶导数 $p^{\′\′}\left({\mathrm{\α}}_n\right)=(p^{\′}\left(i\right)-\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)})*2/(\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)).$

⑤重复②～④过程3～5次，得到最优正则化参数α_best。

（4）在最优正则化参数α_best下用阻尼最小二乘法重构并显示图像。

为了说明本发明的方法的准确性和快速性，在实施例中添加了本方法与普遍使用的L型曲线方法在加入0.1%高斯白噪声情况下的对比。

图5、图6分别是添加0.1%高斯白噪声后的本方法和L型曲线法的成像效果，可以看出，两种方法获得的最优图像成像效果相当。

表1展示了两种方法计算所耗费的时间。

表10.1%高斯白噪声下两种方法的对比

表1、图5、图6说明了本方法能够实现正则化参数最优化，与普遍使用的L型曲线法效果相当，并节省了大量的计算时间。

需要说明的是，以上的实施例仅用于本领域的技术人员进一步理解本发明，本发明并不限于该实施例，凡是由本领域技术人员根据发明的技术方案做出的等效替换和增加，同样属于本发明保护的范围。

Claims (9)

1.一种电阻抗断层成像中快速最优的正则化方法，包括以下操作：

1)根据待测体的基本信息，分别采用有限元方法剖分建立稀疏剖分模型和密集剖分模型的两种重构模型；

2)采用稳健有效的正则参数最优化方法求取稀疏剖分模型的最优正则化参数和残差；

其特征在于，还包括以下步骤：

3)将稀疏剖分模型的最优正则化参数作为密集剖分模型的正则化参数初值，将稀疏剖分模型的残差作为密集剖分模型下残差的估计；采用广义偏差原理求取密集剖分模型的最优正则化参数，在该最优正则化参数下用阻尼最小二乘法重构得到电阻抗断层成像的图像。

2.如权利要求1所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，待测体的基本信息包括待测体边界信息、待测体内部结构信息和电极位置信息；

提取待测体的基本信息，将其转换成电阻抗成像所需的先验信息，然后结合所得的先验信息建立稀疏和密集两种剖分模型。

3.如权利要求1所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，建立图像重建模型时的稀疏剖分模型和密集剖分模型，其依据的剖分形式是：电压变化越剧烈的区域的剖分单元规模越小。

4.如权利要求3所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，所依据的剖分形式是：

区域中心的剖分单元最大，区域边缘的剖分单元较小，靠近电极位置的剖分单元最小。

5.如权利要求1所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，重构稀疏模型的图像时采用的稳健有效的正则化参数最优化方法为无需估计误差的正则化参数最优化方法，包括先验选取正则化参数、L型曲线法、广义交叉检验法、最优分辨率法、拟最优准则法和CRESO准则法。

6.如权利要求1所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，获取稀疏剖分模型的最优正则化参数和残差包括以下操作：

候选正则化参数是从小到大的等比数列，对不同的正则化参数α，采用阻尼最小二乘法求解稀疏剖分模型下的逆问题，求解公式为：

x＝(J^TJ+αJ^TJ)^-1J^TΔU

其中J为稀疏剖分模型的雅克比矩阵，重建得到每个正则化参数对应的电阻率分布||x_α||和残差||Ax_α-U_δ||；

以所有候选正则化参数的电阻率分布的模||x_α||和残差||Ax_α-U_δ||分别为横坐标和纵坐标作图得到一条曲线，该曲线曲率最大值点对应的即为最优正则化参数；

在该最优参数下求解稀疏剖分下的逆问题，得到残差δ₀＝||Ax_α-U_δ||²。

7.如权利要求1所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，获取密集剖分模型的最优正则化参数是将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值α_n，然后运用广义偏差原理迭代公式迭代求解。

8.如权利要求7所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，所述的密集剖分模型的最优正则化参数的求解包括以下操作：

①将稀疏剖分模型的最优正则化参数α₀作为密集模型的正则化参数初值α_n；

②寻找候选正则化参数中与α_n最接近的数作为第i次迭代的初值α(i)，并得到α(i)左右的两个点α(i+1)和α(i-1)；

③应用阻尼最小二乘法求取α(i-1)、α(i)、α(i+1)3个正则化参数下的残差||Ax_α-u_δ||²为p(i-1)，p(i)，p(i+1)；

④应用偏差原理迭代公式求取第n次迭代正则化参数：

${\mathrm{\α}}_{n+1}={\mathrm{\α}}_n-\frac{2p\left({\mathrm{\α}}_n\right)}{p^{\′}\left({\mathrm{\α}}_n\right)+\mathrm{sign}\left(p^{\′}\left({\mathrm{\α}}_n\right)\right){(p^{\′}{\left({\mathrm{\α}}_n\right)}^2-2p\left({\mathrm{\α}}_n\right)p^{\′\′}\left({\mathrm{\α}}_n\right))}^{\frac{1}{2}}}$

式中p′表示求一阶导数，p″表示求二阶导数，sign(·)为符号函数；p(α_n)用稀疏模型在最优正则化参数下的残差作近似，p(α)关于正则化参数α(i)的一阶导数 $p^{\′}\left({\mathrm{\α}}_n\right)=\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}*\frac{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)}+\frac{p(i+1)-p\left(i\right)}{\mathrm{\α}(i+1)-\mathrm{\α}\left(i\right)}*\frac{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)}{\mathrm{\α}(i+1)-\mathrm{\α}(i-1)},$ p(α)关于α的二阶导数 $p^{\′\′}\left({\mathrm{\α}}_n\right)=(p^{\′}\left(i\right)-\frac{p\left(i\right)-p(i-1)}{\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1)})*2/(\mathrm{\α}\left(i\right)-\mathrm{\α}(i-1));$

⑤重复②～④过程3～5次，得到密集剖分模型的最优正则化参数α_best。

9.如权利要求8所述的电阻抗断层成像中快速最优的正则化方法，其特征在于，在最优正则化参数α_best下用阻尼最小二乘法重构并显示图像。