CN104867178A - 一种二维血管内光声图像的建模与仿真方法 - Google Patents

Abstract

一种二维血管内光声图像的建模与仿真方法，所述方法首先建立含有粥样硬化斑块的血管横截面模型，成像导管位于模型中心，将血管横截面等角度分为份；然后对成像导管在每个角度从模型中心沿径向发射激光脉冲的过程进行仿真，对每个角度对应的多层血管壁组织进行光的蒙特卡罗模拟，得到空间电磁吸收分布函数；再根据仿真出多层血管壁组织的光声信号；最后根据个光声信号得到血管横截面图像。本发明可在有限的硬件设备和较短的时间内得到大量的IVPA图像，为血管内光声成像算法和IVPA图像后处理算法的研究和性能测试等提供数据源，为医学训练提供图像库，为获取更高质量的IVPA图像提供有益的参考。

Description

一种二维血管内光声图像的建模与仿真方法

技术领域

本发明涉及一种对二维血管内光声(intravascularphotoacoustic,IVPA)图像进行建模与仿真的方法，属于医学成像技术领域。

背景技术

血管内光声(intravascular photoacoustic,IVPA)成像是一种成像方法，它根据不同成分的粥样硬化斑块的光吸收系数不同，被光照射热膨胀后产生的声学信号强度也不同这一效应对血管壁及粥样硬化斑块成像。它结合了光声信号激发阶段光吸收较高的对比度和光声信号发射阶段超声检测较高的分辨率，可得到展示斑块类型和尺寸等信息的实时图像。

IVPA图像内容的差异由目标血管结构、斑块类型和重构图像的参数等造成，此外成像仪器的校准和微小的参数差异都会影响成像效果。测试者通过图像判断血管的准确程度、提升成像导管的参数质量、优化图像处理算法等都建立在分析大量的测试图像数据的基础上。然而，IVPA成像技术目前尚未广泛应用，可用的测试数据严重不足。此外，测试血管变异实践是一种比较有效的医学训练途径，但需要较多训练时间和有经验专家的专门指导，训练条件较为繁琐和严苛。

发明内容

本发明的目的在于针对现有技术之弊端，提供一种二维血管内光声图像的建模与仿真方法，在有限的硬件设备、较短的时间内得到大量的IVPA图像，为血管内光声成像算法和IVPA图像后处理算法的研究和性能测试等提供数据源，为医学训练提供图像库，为获取更高质量的IVPA图像提供有益参考。

本发明所述问题是以下述技术方案实现的：

一种二维血管内光声图像的建模与仿真方法，所述方法首先建立含有粥样硬化斑块的血管横截面模型，成像导管位于模型中心，将血管横截面等角度分为m份；然后对成像导管在每个角度从模型中心沿径向发射激光脉冲的过程进行仿真，对每个角度对应的多层血管壁组织进行光的蒙特卡罗(Monte Carlo，MC)模拟，得到空间电磁吸收分布函数A(r)；再根据A(r)仿真出多层血管壁组织的光声信号p(t)；最后根据m个光声信号得到血管横截面图像。

上述二维血管内光声图像的建模与仿真方法，所述方法按以下步骤进行：

a.建立血管横截面模型：

成像导管位于模型的中央，周围依次为管腔(血液)、斑块(纤维、钙质、脂质)、内膜/中膜(肌肉组织)和外膜(结缔组织)，将超声换能器看作理想的点换能器，其扫描轨迹为平行于成像平面、半径趋近于零的圆形轨迹，模型所在的坐标系为平行于成像平面的X-Y平面直角坐标系，其中坐标原点是成像导管中心；然后以血管横截面模型的中心为起始点将模型等角度分为m份，每份的中心线与X轴正半轴的夹角分别为θ₁,θ₂,...,θ_m，将每一份在θ-ρ直角坐标系下近似为多层血管壁组织，其中θ-ρ坐标系的θ轴表示以X轴正半轴为基准逆时针旋转得到的角度，ρ轴表示多层血管壁组织的厚度；确定每个角度对应的多层血管壁组织中每层的参数(包括吸收系数、散射系数、平均折射率、散射各向异性因子、厚度)，形成多层血管壁组织的光学参数模型；

b.仿真IVPA成像导管从模型中心沿径向发射激光脉冲穿过血管内腔中的血液照射血管壁的过程：

忽略光的波动特性，将光照射多层血管壁组织看成大量光子与生物组织相互作用的过程，采用蒙特卡罗(MC)方法模拟多层血管壁组织中激光的传播，得到关于空间的电磁吸收分布函数A(r)；

c.仿真多层血管壁组织产生的光声信号：

根据下式仿真得到成像导管在角度θ＝θ₁,θ₂,...,θ_m时接收到的多层血管壁组织产生的光声信号p₁(t),p₂(t),...,p_m(t)：

$\left\{\begin{array}{c}{p}^{n+1}(j,k)=p^n(j,k)-P\·c^2(j,k)\·\mathrm{\Δt}\·\left\{\frac{1}{\mathrm{\Δ\θ}}\right[v_{\mathrm{\θ}}^{n+1/2}(j+1/2,k)-v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)]\\ +\frac{1}{\mathrm{\Δ\ρ}}[v_{\mathrm{\ρ}}^{n+1/2}(j,k+1/2)-v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)]\}+\frac{\mathrm{\β}\·c^2(j,k)\·\mathrm{\Δt}}{C_P}A(j,k)\·I^n\\ v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)=v_{\mathrm{\θ}}^{n-1/2}(j-1/2,k)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\θ}}[p^n(j,k)-p^n(j-1,k)]\\ v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)=v_{\mathrm{\ρ}}^{n-1/2}(j,k-1/2)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\ρ}}[p^n(j,k)-p^n(j,k-1)]\end{array}\right.$

式中，(j,k)是位于θ-ρ平面内的多层血管壁组织上一点的坐标；Δθ和Δρ分别表示平行和垂直于多层生物组织表面方向的离散平面间距；Δt表示离散时间间距；n表示离散的时刻；pⁿ(j,k)为n时刻(j,k)位置的光声信号强度；分别为n时刻(j,k)位置的质点在θ方向和ρ方向的振动速度；c(j,k)是光声信号在(j,k)位置的速度；A(j,k)是(j,k)位置的电磁吸收分布函数；I_n是n时刻的激光脉冲强度；P为血管壁组织的密度；β为等压膨胀系数；C_P为比热；

d.重建极坐标系下的血管横截面IVPA图像：

位置r处的血管横截面极坐标视图的灰度值即为电磁吸收分布A′(r)，由下式求得：

$A^{\′}\left(r\right)=-\frac{1}{\mathrm{\π}{c}^4}\frac{4\mathrm{\π}{C}_{P}c}{\mathrm{\β}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{z_0}{|r_i-r|}\frac{1}{t}\frac{\∂p_i\left(t\right)}{\∂t}{|}_{t=\frac{|r_i-r|}{c}}$

式中，A′(r)为图像重建得到的关于空间的电磁吸收分布函数；i＝1,2,...,m，z₀为位置r处距多层血管壁组织表面即θ轴的距离；p_i(t)为仿真得到的成像导管在角度θ_i处接收到的光声信号；r_i为θ-ρ平面中与成像导管在X-Y平面中的成像角度θ_i相对应的位置；c为光声信号的波速。

e.坐标转换

设θ-ρ坐标系中的一点坐标为(j,k)，其灰度值为f(j,k)，该点在X-Y坐标系中的对应点的坐标为(j′,k′)，灰度值为g(j′,k′)，其中j∈[0,360]，k∈[0,d]，j′∈[-d,d]，k′∈[-d,d]，d为步骤d所求极坐标图像的纵坐标的最大值，则有：

$g(j^{\′},k^{\′})=f(\mathrm{\θ},\sqrt{j^{\′2}+k^{\′2}})$

其中

θ为以X轴正半轴为基准逆时针旋转得到的角度。

本发明可灵活地调整血管横截面模型，通过改变模型中所含粥样硬化斑块的类型和设置多层血管壁组织的光学和声学参数，可精确地仿真声速不均匀的多层血管壁组织的声场，从而在有限的硬件设备和较短的时间内得到大量的IVPA图像，为血管内光声成像算法和IVPA图像后处理算法的研究和性能测试等提供数据源，为医学训练提供图像库，为获取更高质量的IVPA图像提供有益的参考。

附图说明

下面结合附图对本发明作进一步详述。

图1是含有薄纤维帽和脂质斑块的粥样硬化血管横截面模型示例；

图2是IVPA成像导管在角度θ_i处发射激光脉冲并接收光声信号的示意图；

图3是将角度θ_i处对应的成像区域近似为多层血管壁组织的示意图；

图4是将图1中的血管横截面模型示例分为256份后，其中的一份近似为多层血管壁组织的示意图；

图5是光子在多层血管壁组织中的MC模拟的流程图；

图6是光声图像重建的直线扫描模式示意图。

文中各符号为：X、Y、X-Y平面直角坐标系的横轴和纵轴，其中，成像导管的中心为坐标原点，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向；m、横截面模型被等角度分割的总份数；θ、以X轴正半轴为基准逆时针旋转得到的角度；θ_i、成像导管的第i个成像角度；θ_ia、θ_ib、成像导管在角度θ_i处进行成像时对应的角度范围的上、下限，其中，i＝1,2,...,m；θ、ρ、θ-ρ坐标系的横轴和纵轴，其中，水平向右的方向为θ轴正方向，垂直于θ轴向上的方向为ρ轴正方向；A(r)、激光照射血管壁组织得到的关于空间的电磁吸收分布函数；A′(r)、图像重建得到的关于空间的电磁吸收分布函数；哈密顿算子；θ-ρ坐标系θ方向和ρ方向的单位向量；t、时间；r、θ-ρ坐标系上的点；p_i(r,t)、成像导管在角度θ_i处进行成像时对应的血管横截面模型的角度位置θ_ia≤θ≤θ_ib区域的多层生物组织的位置r处的光声信号，其中，i＝1,2,...,m；A_i(r)、MC仿真得到的成像导管在角度θ_i处进行成像时对应的血管横截面模型的角度位置θ_ia≤θ≤θ_ib区域的多层生物组织的电磁吸收分布函数；I(t)、激光脉冲函数；c、光声信号的波速；β、等压膨胀系数；C_P、比热；p、声压；v、质点的振动速度；P'、密度变化量；带方向的质点振速；P、血管壁组织的密度；v_θ、v_ρ、质点在θ方向和ρ方向的振速；(j,k)、位于θ-ρ平面内的多层血管壁组织上一点的坐标；Δθ、Δρ、平行和垂直于多层血管壁组织表面方向的离散平面间距；Δt、离散时间间距；n、离散的时刻；pⁿ(j,k)、n时刻(j,k)位置的光声信号强度； n时刻(j,k)位置的质点在θ方向和ρ方向的振动速度；c(j,k)、光声信号在(j,k)位置的速度；A(j,k)、(j,k)位置的电磁吸收分布函数；I_n、n时刻的激光脉冲强度；p_i(t)、仿真得到的成像导管在角度θ_i处接收到的光声信号，其中，i＝1,2,...,m；z₀、位置r处距多层生物组织表面即θ轴的距离；r_i、θ-ρ平面中与成像导管在X-Y平面中的成像角度θ_i相对应的位置；d、极坐标图像的纵坐标的最大值；(j,k)、(j′,k′)、θ-ρ坐标系中的一点坐标及其在X-Y坐标系中对应点的坐标；f(j,k)、g(j′,k′)、点(j,k)和点(j′,k′)的灰度值。

具体实施方式

本发明方法包括如下步骤：

(1)建立血管横截面模型：

如附图1所示，血管横截面模型包括成像导管(接收光声信号的超声换能器位于成像导管顶端)、管腔(血液)、内膜/中膜(肌肉组织)、外膜(结缔组织)和斑块(纤维、钙质、脂质)四部分。其中成像导管位于模型的中央，周围依次为管腔(血液)、斑块、内膜/中膜和外膜。本方法忽略超声换能器的孔径效应(即超声换能器接收信号时，受感面各处不能同时达到一样的相位分布)，将其看作理想的点换能器，其扫描轨迹为平行于成像平面，半径趋近于零的圆形轨迹。根据斑块的类型(钙化、纤维化或者脂质斑块)和大小，以及血管内腔、血管壁内膜/中膜、外膜的厚度建立不同的血管横截面模型。模型所在的坐标系为X-Y平面直角坐标系，其中坐标原点是成像导管中心，水平向右的方向为X轴正方向，垂直于X轴向上的方向为Y轴正方向。

然后，如附图2所示，以血管横截面模型的中心为起始点将模型等角度分为m份，在对每一份进行激光照射和接收光声信号的仿真过程中，成像导管所处的成像角度为

θ_i＝360(i-1)/m             (1)

其中，i＝1,2,...,m。成像区域的角度范围为

θ_ia≤θ≤θ_ib           (2)

其中θ是以X轴正半轴为基准逆时针旋转得到的角度；

θ_ia＝θ_i-180/m          (3)

和

θ_ib＝θ_i+180/m             (4)

分别是成像导管在角度θ_i处进行成像时对应的成像区域角度范围的上、下限。

如附图3所示，将每一份在θ-ρ坐标系下近似为多层血管壁组织，其中θ-ρ坐标系的θ轴正方向为水平向右的方向，ρ轴正方向为垂直于θ轴向上的方向，表示多层血管壁组织的厚度。多层血管壁组织的表面平行于θ轴且垂直于ρ轴。确定每个角度对应的多层血管壁组织中每层的参数(包括吸收系数、散射系数、平均折射率、散射各向异性因子、厚度)，形成多层血管壁组织的光学参数模型。

例如，当m＝256时，将附图1中的含有薄纤维帽和脂质斑块的血管横截面模型示例以所在坐标系原点为中心等分为256份，成像导管在对每一份进行激光照射和接收光声信号的过程中所处的角度分别为：θ₁＝0,θ₂＝360/256,...,θ_i＝360(i-1)/256,...,θ₂₅₆＝360×255/256。以成像导管角度位置为θ₁＝0时为例，如附图4所示，将对应的成像区域近似为多层血管壁组织，该多层血管壁组织的光学参数如表1所列。i＝2,3,...,256时，成像角度为θ_i时对应的多层血管壁组织的光学参数与成像角度为θ₁时相比，仅厚度参数不同，当成像角度为θ_i时对应的多层血管壁组织不包含斑块时，c、d层的厚度参数为0。

表1

(2)仿真IVPA成像导管从模型中心沿径向发射激光脉冲穿过血管内腔中的血液照射血管壁的过程：

忽略光的波动特性，将光照射多层血管壁组织看成大量光子与生物组织相互作用的过程，由于折射率不同，光子在两层生物组织的交界处会发生散射或透射。本发明方法模拟5×10⁵个光子以总能量1J、半径0.01cm的高斯光束形式照射血管壁，参照文献(Wang Lihong.Biomedical optics:principles and imaging[M]//Wang Lihong.Monte Carlo modeling of photon transport in biologicaltissue.New Jersey:John Wiley and Sons,2007:37-65)中的蒙特卡罗(Monte Carlo，MC)方法模拟多层血管壁组织中激光的传播，得到关于空间的电磁吸收分布函数A(r)。进行蒙特卡罗模拟的流程图如附图5所示。

(3)仿真多层血管壁组织产生的光声信号

描述光声信号在声学均匀介质中传播的物理模型为：

$\left\{\begin{array}{c}{\▿}^2{p}_i(r,t)-\frac{1}{c^2}\frac{{\∂}^2}{\∂t^2}{p}_i(r,t)=-\frac{\mathrm{\β}}{C_P}{A}_i\left(r\right)\frac{\∂}{\∂t}I\left(t\right)\\ \▿=\frac{\∂}{\∂\mathrm{\θ}}\stackrel{\→}{j}+\frac{\∂}{\∂\mathrm{\ρ}}\stackrel{\→}{k}\end{array}\right.---\left(5\right)$

其中，为哈密顿算子；和分别为θ-ρ坐标系的θ方向和ρ方向的单位向量；t为时间；r为θ-ρ坐标系上的点；如附图3所示，成像导管在角度θ_i(i＝1,2,...,m)处进行成像时，成像区域的角度范围为θ_ia≤θ≤θ_ib(θ_ia＝θ_i-180/m，θ_ib＝θ_i+180/m)，将该成像区域在θ-ρ坐标系下近似为多层血管壁组织，p_i(r,t)为t时刻在位置r处的光声信号；A_i(r)为位置r的电磁吸收分布函数；c为光声信号的波速；β为等压膨胀系数；C_P为比热；I(t)为激光脉冲函数，可近似为冲激函数。

光声信号的本质是超声波，分析声波的传输需引入声压p、质点振动速度v、和密度变化量P'这三个物理量。声振动是一种宏观的物理现象，满足牛顿第二定律、质量守恒定律和描述压强、温度与体积等状态参数关系的物态方程这三个基本的物理定律，根据这些基本定律导出p、v和P'三者的关系方程：

$P\frac{\∂v}{\∂t}=-\▿p---\left(6\right)$

$\frac{\∂P^{\′}}{\∂t}=-\▿\·\left(P\stackrel{\→}{v}\right)---\left(7\right)$

p＝v²P'               (8)

式(6)是运动方程，式(7)是连续性方程，式(8)是物态方程，其中，P为血管壁组织的密度，假定为一个常数；为带方向的质点振速。

根据式(6)、(7)、(8)将式(5)改写为如下形式：

$\left\{\begin{array}{c}\frac{\∂p}{\∂t}=-P\·c^2(\frac{\∂v_{\mathrm{\θ}}}{\∂\mathrm{\θ}}+\frac{\∂v_{\mathrm{\ρ}}}{\∂\mathrm{\ρ}})+\frac{\mathrm{\β}\·c^2}{C_P}A\left(r\right)\·I\left(t\right)\\ \frac{\∂v_{\mathrm{\θ}}}{\∂t}=-\frac{1}{P}\frac{\∂p}{\∂\mathrm{\θ}}\\ \frac{\∂v_{\mathrm{\ρ}}}{\∂t}=-\frac{1}{P}\frac{\∂p}{\∂\mathrm{\ρ}}\end{array}\right.---\left(9\right)$

式中，v_θ和v_ρ分别为质点在θ方向和ρ方向的振速。

将式(9)离散化为下列形式：

$\left\{\begin{array}{c}{p}^{n+1}(j,k)=p^n(j,k)-P\·c^2(j,k)\·\mathrm{\Δt}\·\left\{\frac{1}{\mathrm{\Δ\θ}}\right[v_{\mathrm{\θ}}^{n+1/2}(j+1/2,k)-v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)]\\ +\frac{1}{\mathrm{\Δ\ρ}}[v_{\mathrm{\ρ}}^{n+1/2}(j,k+1/2)-v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)]\}+\frac{\mathrm{\β}\·c^2(j,k)\·\mathrm{\Δt}}{C_P}A(j,k)\·I^n\\ v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)=v_{\mathrm{\θ}}^{n-1/2}(j-1/2,k)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\θ}}[p^n(j,k)-p^n(j-1,k)]\\ v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)=v_{\mathrm{\ρ}}^{n-1/2}(j,k-1/2)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\ρ}}[p^n(j,k)-p^n(j,k-1)]\end{array}\right.---\left(10\right)$

式中，(j,k)是位于θ-ρ平面内的多层血管壁组织上一点的坐标；Δθ和Δρ分别表示平行和垂直于多层生物组织表面方向的离散平面间距；Δt表示离散时间间距；n表示离散的时刻；pⁿ(j,k)为n时刻(j,k)位置的光声信号强度；分别为n时刻(j,k)位置的质点在θ方向和ρ方向的振动速度；c(j,k)是光声信号在(j,k)位置的速度；A(j,k)是(j,k)位置的电磁吸收分布函数；I_n是n时刻的激光脉冲强度。

根据式(10)即可仿真得到成像导管在角度θ＝θ₁,θ₂,...,θ_m时接收到的多层血管壁组织产生的光声信号p₁(t),p₂(t),...,p_m(t)。

(4)重建极坐标系下的血管横截面IVPA图像

从数学的角度来讲，光声成像的图像重建问题即由光声信号p(t)计算出生物组织电磁吸收分布函数A′(r)。具体步骤如下：

首先利用步骤(3)得到的一组m个光声信号p₁(t),p₂(t),...,p_m(t)求出其中i＝1,2,...,m。然后进行多个测量位置带权重的反投影累加，近似得到电磁吸收分布A′(r)：

$A^{\′}\left(r\right)=-\frac{1}{\mathrm{\π}{c}^4}\frac{4\mathrm{\π}{C}_{P}c}{\mathrm{\β}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{z_0}{|r_i-r|}\frac{1}{t}\frac{\∂p_i\left(t\right)}{\∂t}{|}_{t=\frac{|r_i-r|}{c}}---\left(11\right)$

式中，A′(r)为图像重建得到的关于空间的电磁吸收分布函数；z₀为位置r处距多层血管壁组织表面即θ轴的距离；p_i(t)为仿真得到的成像导管在角度θ_i处接收到的光声信号；如附图6所示，r_i为θ-ρ平面中与成像导管在X-Y平面中的成像角度θ_i相对应的位置。根据式(11)求出的A′(r)即为位置r处的血管横截面极坐标视图的灰度值。

(5)坐标转换

将步骤(4)中获取的θ-ρ坐标系中的极坐标视图转换为X-Y直角坐标系下的横截面视图。具体步骤如下：

设θ-ρ坐标系中的一点坐标为(j,k)，其灰度值为f(j,k)，该点在X-Y坐标系中的对应点的坐标为(j′,k′)，灰度值为g(j′,k′)，其中j∈[0,360]，k∈[0,d]，j′∈[-d,d]，k′∈[-d,d]，d为步骤(4)中所求极坐标图像的纵坐标的最大值。那么

$g(j^{\′},k^{\′})=f(\mathrm{\θ},\sqrt{j^{\′2}+k^{\′2}})---\left(12\right)$

其中

如附图2所示，θ为以X轴正半轴为基准逆时针旋转得到的角度。所求g(j′,k′)即为X-Y坐标系中的血管横截面图像的点(j′,k′)处的灰度值。

Claims (2)

1.一种二维血管内光声图像的建模与仿真方法，其特征是，所述方法首先建立含有粥样硬化斑块的血管横截面模型，成像导管位于模型中心，将血管横截面等角度分为m份；然后对成像导管在每个角度从模型中心沿径向发射激光脉冲的过程进行仿真，对每个角度对应的多层血管壁组织进行光的蒙特卡罗模拟，得到空间电磁吸收分布函数A(r)；再根据A(r)仿真出多层血管壁组织的光声信号p(t)；最后根据m个光声信号得到血管横截面图像。

2.根据权利要求1所述的二维血管内光声图像的建模与仿真方法，其特征是，所述方法按以下步骤进行：

a.建立血管横截面模型：

成像导管位于模型的中央，周围依次为管腔、斑块、内膜、膜和外膜，将超声换能器看做理想的点换能器，其扫描轨迹为平行于成像平面，半径趋近于零的圆形轨迹，模型所在的坐标系为平行于成像平面的X-Y平面直角坐标系，其中坐标原点是成像导管中心；然后以血管横截面模型的中心为起始点将模型等角度分为m份，每份的中心线与X轴正半轴的夹角分别为θ₁,θ₂,...,θ_m，将每一份在θ-ρ直角坐标系下近似为多层血管壁组织，其中θ-ρ坐标系的θ轴表示以X轴正半轴为基准逆时针旋转得到的角度，ρ轴表示多层血管壁组织的厚度；确定每个角度对应的多层血管壁组织中每层的参数：包括吸收系数、散射系数、平均折射率、散射各向异性因子、厚度，形成多层血管壁组织的光学参数模型；

b.仿真IVPA成像导管从模型中心沿径向发射激光脉冲穿过血管内腔中的血液照射血管壁的过程：

忽略光的波动特性，将光照射多层血管壁组织看成大量光子与生物组织相互作用的过程，采用蒙特卡罗方法模拟多层血管壁组织中激光的传播，得到关于空间的电磁吸收分布函数A(r)；

c.仿真多层血管壁组织产生的光声信号：

根据下式仿真得到成像导管在角度θ＝θ₁,θ₂,...,θ_m时接收到的多层血管壁组织产生的光声信号p₁(t),p₂(t),...,p_m(t)：

$\left\{\begin{array}{c}{p}^{n+1}(j,k)=p^n(j,k)-P\·c^2(j,k)\·\mathrm{\Δt}\·\left\{\frac{1}{\mathrm{\Δ\θ}}\right[v_{\mathrm{\θ}}^{n+1/2}(j+1/2,k)-v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)]\\ +\frac{1}{\mathrm{\Δ\ρ}}[v_{\mathrm{\ρ}}^{n+1/2}(j,k+1/2)-v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)]\}+\frac{\mathrm{\β}\·c^2(j,k)\·\mathrm{\Δt}}{C_P}A(j,k)\·I^n\\ v_{\mathrm{\θ}}^{n+1/2}(j-1/2,k)=v_{\mathrm{\θ}}^{n-1/2}(j-1/2,k)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\θ}}[p^n(j,k)-p^n(j-1,k)]\\ v_{\mathrm{\ρ}}^{n+1/2}(j,k-1/2)=v_{\mathrm{\ρ}}^{n-1/2}(j,k-1/2)-\frac{\mathrm{\Δt}}{P\·\mathrm{\Δ\ρ}}[p^n(j,k)-p^n(j,k-1)]\end{array}\right.;$

式中，(j,k)是位于θ-ρ平面内的多层血管壁组织上一点的坐标；Δθ和Δρ分别表示平行和垂直于多层生物组织表面方向的离散平面间距；Δt表示离散时间间距；n表示离散的时刻；pⁿ(j,k)为n时刻(j,k)位置的光声信号强度；分别为n时刻(j,k)位置的质点在θ方向和ρ方向的振动速度；c(j,k)是光声信号在(j,k)位置的速度；A(j,k)是(j,k)位置的电磁吸收分布函数；Iⁿ是n时刻的激光脉冲强度；P为血管壁组织的密度；β为等压膨胀系数；C_P为比热；

d.重建极坐标系下的血管横截面IVPA图像：

位置r处的血管横截面极坐标视图的灰度值即电磁吸收分布A′(r)，由下式求得：

$A^{\′}\left(r\right)=-\frac{1}{{\mathrm{\πc}}^4}\frac{4\mathrm{\π}{C}_{P}c}{\mathrm{\β}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\frac{z_0}{|r_i-r|}\frac{1}{t}\frac{{\∂p}_i\left(t\right)}{\∂t}|}_{t=\frac{|r_i-r|}{c}},$

式中，A′(r)为图像重建得到的关于空间的电磁吸收分布函数；i＝1,2,...,m，z₀为位置r处距多层血管壁组织表面即θ轴的距离；p_i(t)为仿真得到的成像导管在角度θ_i处接收到的光声信号；r_i为θ-ρ平面中与成像导管在X-Y平面中的成像角度θ_i相对应的位置；c是光声信号的波速。

e.坐标转换

设θ-ρ坐标系中的一点坐标为(j,k)，其灰度值为f(j,k)，该点在X-Y坐标系中的对应点的坐标为(j′,k′)，灰度值为g(j′,k′)，其中j∈[0,360]，k∈[0,d]，j′∈[-d,d]，k′∈[-d,d]，d为步骤d所求极坐标图像的纵坐标的最大值，则有：

$g(j^{\′},k^{\′})=f(\mathrm{\θ},\sqrt{j^{\′2}+k^{\′2}});$

其中

θ为以X轴正半轴为基准逆时针旋转得到的角度。