层析成像

层析成像（精选十篇）

层析成像 篇1

我们都有这样的经历,在乌云密布的雷雨天,空中的闪电照亮了一大片云层,随之而来的是轰鸣的雷声。19世纪70年代,Few利用雷声信号到达不同雷声传感器的时间差以及雷声信号和雷电电磁信号的时间差首次定位了雷声声源的位置[1]。这个事实告诉我们,如果光源发光的同时产生声波,即使该光源处于光学散射介质中,例如生物组织,我们也可以通过观测声波获得光源的位置。本文要介绍的光声层析成像(PAT)就是利用这样的探测原理。

PAT成像的理论依据是光声效应(Photoacoustic effect,PA),该效应描述的是:当脉冲或经过调制的电磁波来照射物体时,有的物体会吸收电磁波能量并发热,伴随着的热膨胀产生向外传播的声波[2]。该效应最早由亚历山大·贝尔于1880年发现,到20世纪70年代开始广泛地应用于物理、化学、生物、医药等多个领域中(Rosencwaig 1980,Gusev et al 1993)。PAT正是利用这个效应,并结合声波在软组织中的低散射性,通过测量产生的声波获得体内光学吸收体的位置和形态等信息。

PAT最重要的优势就是突破了纯光学高分辨成像技术的成像深度壁垒:由于组织对光的强散射作用,光学成像分辨率随着深度的增加而急剧降低,使得纯光学技术的高分辨组织成像被限制在几个毫米深度,这在很大程度上限制了它的实际应用范围[3,4]。PAT正是通过光学与超声技术的结合,充分地利用了低散射的超声波,实现了几个厘米深的高分辨成像。以下我们先解释PAT的基本机制,然后介绍它的主要成像模式,在第三部分讨论了PAT在临床医学生的研究现状和潜在应用,最后总结并展望该技术未来可能的发展趋势。

2 光声层析成像原理及模式

在生物医学领域中,具有短脉宽(纳秒级)的脉冲激光常被用做PAT的光源[2]。就像云层中的闪电那样,激光射入组织后,光子被组织散射到很大的区域。如果在光子到达的区域中存在光学吸收体,那么这些部位吸收光能后产生的热量会引起吸收体升温膨胀挤压周围的组织从而产生超声波。就像雷声穿透厚厚的云层那样,产生的超声波在软组织中可以在软组织内自由传播很远。

在位置r处的物体吸收脉冲光能量后发生热膨胀会产生一个初始压强p0(r),其正比于吸收光的能量,数学表达式为:

其中A(r)是指位置r处单位时间和单位体积内吸收的光能量,Γ是格日尼森系数(Gruneisen coefficient),描述吸收的光能量转化为压力的关系,μa是吸收体的光吸收系数,F是到达r处的光子注量(正比于光子密度)。在一定的光照和环境条件下,初始压力波幅值正比于组织的光吸收系数。通过初始条件(公式1),根据声波的波动方程就可以得到超声波在软组织中的传播过程。由于对超声压力波的探测灵敏度足够高,PAT中引起的组织局部升温原小于1摄氏度,是一种非常安全的成像技术。,而不同组织的光吸收系数不同,就成为了PAT图像对比度的来源。图1列举了生物体内主要组织的消光系数(主要是吸收)随波长的变化。例如,在可见光和近红外范围内血管和黑色素相比于外周脂肪组织和水有更大的消光系数,这一特性使PAT在血管系统以及黑色素瘤等成像方面具有巨大的优势。

基于不同的重建算法,PAT可以分为两大类:光声显微(PA microscopy,PAM)和光声计算层析(PA computed tomography,PACT)。类似光学显微成像使用聚焦物镜那样,PAM使用聚焦的超声探头来探测在探头聚焦区域的原始信号,不同的是PAM可以根据超声信号的时域信息获得该方向上的一维深度信息。根据决定横向分辨率的方法不同,PAM又可以分为声分辨率的PAM (acoustic-resolution PAM,AR-PAM)和光分辨率的PAM (optical-resolution PAM,OR-PAM)。相比PAM的高分辨率,PACT更多的像临床用的B超那样多点采集信号,然后通过计算机的重建获得由深度散射光产生的光声信号,实现深度组织成像。图2展示了PAT不同模式的成像深度和分辨率的趋势。从图中可以看到,目前PAT的成像深度可以达到10-1～102mm量级,而分辨率跨越了10-1～103μm4个量级,同时成像深度和分辨率的比值高达200。多尺度高分辨率的成像优势使PAT在生物医学光学成像中成为了一颗耀眼的明星[4,5,6,7,8]。

2.1光声显微成像

我们已经知道,根据横向分辨率的决定因素不同,PAM可以分为OR-PAM和AR-PAM。下面将会一一介绍。

PR-PAM基于光学聚焦,在PAT成像的所有模态中提供了最高的成像分辨率。图3给出了一个OR-PAM模式的系统及成像结果[9,10]。在该式系统中(图3(a)和(b)),聚焦激光(532 nm)从样品底部照射,同时聚焦超声探头从样品上部接收。为了达到最优的信噪比,激发光和探测超声的要求“共聚焦”。由超声探头接收到的信号被放大器放大,经高速采集卡转换成数字信号存储并用于后续的图像处理。激发聚焦光路,超声探头以及放大器固定在二维平移台上,进行平面扫描。在PAM系统中,每一个光激发位置得到该深度方向的一维信息,随着平移台的二维移动,最终可以得到样品的三维信息。图3 (c)是用该系统对受精后3～4天斑马鱼进行成像以及对比照片结果。基于斑马鱼的自发对比度,我们可以清晰地看到斑马鱼全身的血管,黑色素,眼睛等精细结构。

由于光在组织中散射很强,依赖光学聚焦的OR-PAM的成像深度受到限制。AR-PAM的横向和纵向分辨率都是由超声探头决定。由于超声在组织中的散射很小,基于超声探头不同的带宽,可以实现从皮下血管到内脏的不同尺度成像。

图4给出了我们研发出的AR-PAM系统结构图以及对小鼠下肢血管的成像结果[11]。如图4(a),脉冲激发经多模光纤导出,经过锥透镜,聚焦透镜以及聚光器后在样品表面形成中空的多纳圈。被激发出的光声信号被聚焦探头接收,经过放大器后由计算机进行采集和处理。图4(b)是对(c)中黄框区域成像结果。在无任何标记的情况下,小鼠后腿部的大小血管都可以以高信噪比进行成像,并且部分血管是在皮下较深部位。

总体来说,光声显微成像利用纯光学聚焦或者超声聚焦,同时具备光学成像的高分辨率又利用超声的强穿透性弥补了光在组织中的强散射,从而提高了成像的深度。光声显微成像在细胞,表层以及皮下血管的成像中有不可替代的地位。

2.2 光声计算机断层成像

光声计算机层析成像(PACT)通过探测被成像物体周围不同位置处的光声信号p(r,t),经过重建算法得出组织内光学吸收体位置。最简单的重建过程与在不同地点记录雷声,从而得到闪电位置类似。

和PAM相比,PACT成像中收集信号的位置有更大的自由度,可以针对真实成像对象进行相应的设计,使其在临床前期动物研究和临床应用方面具有广阔的空间。早期的PAT成像利用单探头环形扫描的形式,成功获得了小动物脑皮层血管的高分辨无创成像。随着技术的发展,越来越多的定制PAT阵列用于不同的研究方向。例如,目前国际上已有多个实验室研发了PAT小动物全身成像系统,部分已经形成产品[12,13,14,15]。

图5是我们研发的一个PACT系统[16]。激发光经过棱镜、锥面镜以及反声透光膜水平照射于样品,由样品断层面发出的光声信号经过反声透光膜被聚焦超声探头接收,在由放大器放大后被计算机采集和处理。超声探头被固定在与入射光同轴心的圆形旋转控制台上进行环形扫描。利用该系统我们研究了成年斑马鱼。图5 (b)-(g)是(h)中相应黑色标记线的断层图像。眼睛、血管、鳃、鳍等器官结构在PACT上高对比度地成像。

除了介绍的这两种基本类型,根据不同的成像需要和条件,PAT己经发展了大量的成像系统,有的是这两类的综合。

2.3 光声功能与分子影像

作为光学成像的一种,PAT也可以通过多波长的进行光谱功能性成像。此外,和其他成像模态类似,当与靶向光声造影剂相结合,PAT可以针对性地对肿瘤、炎症等进行分子影像。

光声功能影像的代表是对血氧饱和度的成像。红细胞内的氧和血红蛋白(HbO2)和脱氧血红蛋白(Hb)对光的吸收有很大不同,所以血液的光学吸收系数为

其中[HbO2]表示氧合血红蛋白的摩尔浓度(mol L-1),[Hb]表示脱氧血红蛋白的摩尔浓度。若用εHb(λ)表示脱氧血红蛋白在波长λ下的吸收系数,εHbO2(λ)表示氧合血红蛋白在波长λ下的吸收系数。通过在多波长下进行PAT成像,就可以将[Hb](r)和[HbO2](r)计算出来,进而得到该位置的血氧饱和度

生物体局部血氧饱和度的变化是反映新陈代谢水平的重要标志,不少课题组在功能性光声成像中已取得了重要突破。Hu等(2011)利用高分辨率OR-PAM实现高精度小鼠耳部毛细血管血氧饱和度成像(如图6b)[7],Jiang等(2014)利用PAT实现乳腺血氧饱和度的定量检测,并监测新型辅助化疗后乳腺癌血氧饱和度的变化[17]。

光声分子影像把分子影像和PAT结合,通过对带有分子探针的光声对比剂成像,从影像学上研究在分子和细胞层次上的生理病理过程。为了使分子探针具有很好的光学吸收特性,利用生物化学的方法把分子探针(例如抗体蛋白)和光学吸收体结合(例如染料分子或者纳米粒子),成为了光声造影剂的主要合成方法。分子探针通过与靶分子的结合,将光学吸收体富集在靶分子的位置。通过对光声造影剂的高对比度成像,就可以得到靶分子的在体分布,实现分子影像的目的。目前,已有多种光学吸收体,包括碳纳米材料、金属纳米颗粒和染料等,被成功地用于合成光声造影剂,并用于研究包括肿瘤和炎症等生理病理过程[18,19,20,21,22,23,24]。以下以金纳米粒子做具体说明。

图7是对黑色素靶向的纳米金笼,随着观察时间的增加,在黑色素瘤处的PA信号逐渐增强,经过6 h后,达到最大并稳定[25]。通过在0h和6h的对比,我们可以看出对黑色素瘤具有靶向性的纳米金笼为PAT成像提供了极高的对比度,这也为肿瘤生理病理过程的研究提供了极大帮助。

3 光声层析成像在医学上的应用

PAT不仅突破了传统高分辨光学成像的深度壁垒,也可以通过光谱PAT实现功能成像。PAT成像在临床上有极其广泛的应用前景,是目前国内外重点研究的新型临床影像设备。根据PAT的成像特点,现有的重点研究方向包括:皮肤成像、乳腺成像等。

皮肤虽然平均厚度只有1 mm,但皮肤对光的强散射使得传统光学显微镜在皮肤病的诊断上作用非常有限,而纯超声对疾病组织的对比度又不够。因此皮肤科医生更多的是依靠经验而不是影像结果进行诊断。PAT可以对皮肤内和皮下的很多病变组织(如黑色素瘤、血管增生相关)进行高分辨率无创成像。目前,针对皮肤癌、烧伤、葡萄酒色斑等疾病的诊断和影像导引治疗已经开展前期研究[26,27,28]。除了皮肤成像外,由于血管增生是乳腺癌的重要特征,本身也是光学强吸收体,针对血管的光声乳腺成像一直是PAT研究的热点。国内外在这方面都做了大量研究,研发了不仅和超声检查结合的手持式光声/超声探测设备,还有定制的光声乳腺成像系统。目前,已经实现在无标记乳房4 cm深度的PAT成像[29]。多个实验室已经开展临床前期实验,并获得了重要进展[30,31,32]。

随着PAT的发展,越来越多的PAT系统面向临床问题需求进行研发,包括眼睛成像[10,33]、光声内窥镜[34]、脑成像[35]等。

4 展望与总结

本文系统性地介绍了PAT成像技术的原理,系统模式以及在基础研究和临床医学中的应用。PAT的重要优势是其利用软组织中超声弱散射的特点,将高分辨光学成像深度大大推进到几个厘米的范围,可以实现从细胞器到组织乃至器官的多尺度成像。另一方面它也保留了光学光谱成像的优势,可以通过多光谱实现重要的功能影像,如定量测量血氧饱和度。同时,光声造影剂的使用进一步提高了PAT的成像深度和对比度,并通过靶向造影剂实现了光声分子影像。作为一种新兴的成像手段,PAT成像技术在生命科学和临床医学上快速发展,并受到广泛关注。然而,PAT在临床的应用依然面临挑战,主要是存在对深层组织的图像的对比度不够的缺点。这一方面要求在激光(尤其是红外激光)技术和高灵敏超声探测技术的进一步发展,另一方面需要和临床研究更紧密结合,并结合图像分析和处理技术,对PAT的临床成像结果进行更加精确和全面的解读。

光声层析成像通过在体内制造“闪电”,并透过厚厚的组织探测产生的“雷声”,为生命医学成像开辟了一个新的道路。如何只让特定的组织(如癌症组织)产生“雷声”?如何更有效地探测“雷声”都是今后的研究重点。我们期望在今后的五到十年内,PAT可以实现在临床上的真正应用,成为像超声系统那样普遍采用的临床诊断影像设备。

致谢

电磁层析成像下的石油勘探论文 篇2

在石油勘探的过程中，需要对不同的岩层的电磁层析成像声波进行准确的分析，将分析结果作为石油勘探中图像重构的主要依据。重构结果的好坏直接影响到后续的工作:其详细原理如下:

1)针对所有的石油勘探中电磁层析成像信号进行信号分解，可以获取信号波变换值x(n，p)。

2)通过运算获取电磁层析成像信号分解尺度S(n，p)，根据下述公式可以计算上述成像信号分解结果的关联性:

3)针对所有的电磁层析成像信号分解结果进行归一化变换，能够得到下述结果PS(n，p)=S(n，p)Qx(n)Qs(n)，q=1，2，…，q(1)在上述公式中，Qx(n)=∑px(n，p)2QS(n)=∑pS(n，p)2

4)将PS(n，p)与x(n，p)的绝对值进行比较，假设|Ps(n，p)||x(n，p)|，则可以判断该位置的信号变换值是由初始的电磁层析成像信号中分离出来的，此时，需要将x(n，p)赋予xg(n，p)，同时将S(n，p)置零。否则，需要保留x(n，p)的初始值。

5)根据上面的阐述能够得知，x(n，p)都是由电磁层析成像信号中的声波引起的，则Qx(n)/(M-1)是信号变换尺度上对声波的均方差进行估计的结果。利用下述公式可以计算声波无偏差估计结果μ=Qx(n)/(M-1)ζn(2)将上述无偏差估计结果与阈值进行对比，假设大于阈值，则返回步骤(2)继续执行，否则，结束迭代处理。6)根据xg(n，p)进行电磁层析成像信号声波的信号逆变换处理。假设μ>1，能够得到sin(μ-1)>0，假设μ1，能够得到sin(μ-1)0。根据上述特性，可以得到最理想的估计结果如下所述d→d+sin[β(μ-1)](3)在上述公式中，β能够用来描述石油勘探中电磁层析成像信号步长调整因子。在第一次迭代处理的过程中，设置d的取值是1，假设μ1，则可以得知sin(μ-1)0，能够利用上述公式进行迭代处理，不断缩小d的取值，并且用该取值乘以对应的关联性系数，减少抽取样本的数目，直至μ的取值趋近于1。假设μ>1，则可以得知sin(μ-1)>0，根据上述公式进行迭代处理，d的取值将不断增大，将该取值乘以关联性系数，可以增加样本的数量，减少声波过滤的误差，直至μ的取值趋近于1。根据上面的阐述能够得知，假设在石油勘探电磁层析成像复杂结构声波过滤过程中，信号系数方差与声波的方差相等，则能够得到理想的声波过滤效果。在过滤过程中，为了保证声波过滤过程的稳定性和提高过滤时的收敛速度，需要将调整因子引入到过滤的过程中，假设该调整因子的取值比较大，则收敛速度得到大幅度提升，但是过滤的效果比较差，反之，收敛速度虽然会降低，但是过滤的效果更佳理想。根据上述内容，能够得到调整因子计算公式如下所述ζn2=ζ2×‖i0*i1*…*ik-2*hk-1‖2(4)根据上面阐述的方法，能够得到电磁层析成像声波过滤时的调整因子，完成声波过滤，得到高质量的成像信号，实现石油勘探中的电磁层析成像声波过滤。但是，通过声波对石油岩层进行勘探形成稳定图像一直存在一个难题，石油区域的岩层结构复杂，电磁层析成像技术在复杂的岩石结构中会产生声波图像变异，形成图像干扰。传统的石油勘探电磁层析成像技术在这种干扰下，正常电磁层析成像会造成干扰损失，影响成像效果。

2成像信号中声波过滤优化方法理论

利用传统算法进行石油勘探中电磁层析成像中声波过滤，假设勘探目标的岩层过于复杂，将导致成像信号中掺杂大量的声波，对成像信号造成干扰。为此，提出基于加权小波分析算法的石油勘探中电磁层析成像复杂结构声波过滤方法。

2.1计算成像信号的权重

针对石油勘探过程中采集的电磁层析成像信号，能够得到对应的非线性方程如下所述yl=gl(yl-1，xl-1)zl=il(yl，wl{)(5)在上述公式中，wl能够用来描述电磁层析成像信号概率密度函数。对上述成像信号中的声波进行过滤的详细流程如下所述

1)对成像信号进行初始化处理，在l=0时，需要使样本符合yj0～Q(y0)分布。

2)在l=1，2，…，U时，需要选择一组成像信号作为样本。

3)根据随机向量计算对应的.权重，随机向量是wj，j=1，2，…，P。4)利用下述公式，对电磁层析成像信号进行更新处理p(yl|Zl)≈∑Pj=1wjlε(yl-yjl)j=1，2，…，P(6)在上述公式中，ε是狄拉克函数。

2.2实现成像信号重构

现阶段，小波分析方法已经应用到各种不同的行业中，发挥着越来越重要的作用。将加权方法与小波分析方法相结合，能够完成成像信号声波的过滤。利用下述公式能够进行小波变换处理XUg(b，c)=1|b|∫-g(y)ζ*(y)dy=g，ζb，{}cζb，c(y)=1|b|ζ(y-cb)b，c∈S，b≠{0(7)在上述公式中，ζb，c(y)能够用来描述小波母函数，b能够用来描述尺度变换算子，c能够用来描述成像信号采样时间。石油勘探中电磁层析成像信号的时间与对应的位置关系密切，ζ*(u)能够用来描述ζ(u)的共轭函数。利用下述公式能够对石油勘探中电磁层析成像信号进行离散变换处理2k2ζ(2ky-l)(8)即:b=12kc=l2{k利用下述公式能够进行石油勘探中电磁层析成像信号分解处理dk，l=∑ni(n-2l)dk-1，nek，l=∑nh(n-2l)dk-1，{n(9)其中，n=0，1，2，…，P-1。利用下述公式，能够对所有的成像信号分解因子进行重构处理:dk，n=∑ldk+1，li(n-2l)+∑lek+1，lh(n-2l)(10)根据上面阐述的方法，能够进行成像信号声波过滤，获取清晰的石油勘探中电磁层析成像信号。

3实验结果分析

为了验证基于加权小波分析算法的石油勘探中电磁层析成像复杂结构声波过滤方法的有效性，需要进行一次实验。在实验的过程中，以岩层回波电磁信号为基础，采集的电磁波成像信号。利用传统算法进行电磁层析成像信息采集，得到的信号成像采集结果能够用图2表示。利用改进算法进行成像信号声波过滤，得到的成像信息采集结果能够用图3表示。根据上述两帧图像能够得知，利用改进算法进行成像信号声波过滤，获取的结果与实际情况更加接近，极大的提高了成像信号过滤结果的真实度。将上述实验数据进行整理分析，能够得到不同算法获取的成像信号真实度对比结果如下所述:根据上述两个表中的数据可以得知，利用改进算法进行石油勘探中电磁层析成像信号声波过滤，极大的提高了成像信号的真实度，降低了成像信号中的信噪比。

4结论

濮阳地震集中区双差层析成像研究 篇3

摘要：利用2000～2013年濮阳地震集中区的地震观测数据，采用双差层析成像方法，对地震进行精定位，并反演了濮阳地震集中区的三维速度结构。结果表明：小震震源主要位于聊城一兰考主干断裂上，优势深度集中在4～15km之间，与断层的产状吻合：濮阳地震集中区的速度结构与地质构造有着一定的相关性，地震大部分发生在高速异常体内，高速异常体内的脆性介质往往是应力最为集中的地方，这可能是导致濮阳地震集中区小震发生的主要原因。

关键词：双差层析成像；速度结构；重定位；濮阳地震集中区

中图分类号：P315.72 文献标识码：A 文章编号：1000-0666（2016）02-0255-06

0 引言

濮阳地震集中区处于山东鄄城、河南濮阳的交界地区，是华北地区地震活动较为集中的区域之一，主要受北北东走向的聊城一兰考断裂带和北东走向的长垣断裂带控制，同时该区还有多条次级断裂，构造特征比较复杂。历史上曾经发生过1830年河北磁县7.5级地震和1937年山东菏泽7.0级地震等多次7级以上强震。自1997年以来，该区连续发生多次4级以上中等地震，经历了两次小震活动的增强过程，目前小震活动仍比较活跃。因此构建濮阳地震集中区三维速度结构，有助于进一步了解该区域的地震活动特征。

地震层析成像是探测地下介质结构和深部断裂活动的重要方法。Zhang和Thurber（2003）在双差定位方法（Waldhauser，Ellsworth，2000）的基础上提出了双差层析成像方法。该方法考虑了介质速度结构在空间范围上的变化，采用层状网格模型，因此可以得到更加精确的定位结果。以往的研究结果表明，利用该方法可得到震源区精细速度结构，为揭示震源区介质分布及地震发生的深部构造研究提供了更加可靠的信息（zhang，Thurber，2005；王长在等，2011，2013；邓文泽等，2014）。本文采用双差层析成像方法，对濮阳地震集中区的地震序列进行重定位和速度结构的反演，为该区域的地震危险性判断和地震趋势研究提供一定的参考。

1 资料与方法

双差层析成像方法通过联合绝对走时数据和相对走时数据实现三维速度结构和震源参数的联合反演。假设有两个相邻的地震，对于第i次地震，台站k的观测到时与理论到时差rⁱ_k可近似表示为

对于符合组对条件的第i次地震，台站k同样可以得到上述的残差r^j_k，则第i和第j两次地震在同一台站的到时差之差可以表示为

通过联合反演就可以得到震源区精细的速度结构及地震重新定位结果。由于该方法联合使用绝对走时和相对走时，增加了双差方程，因此相比于传统方法可以揭示更多的精细结构信息。

本文选取濮阳地震集中区有数字记录资料以来的2000～2013年发生的315个地震。为了确保震相拾取的准确性，所选地震震级均大于1.0级，首先对震相数据进行筛选，依据P波、S波的时距曲线，剔除有明显错误或数据偏差较大的震相。图1给出了本文所使用的震相数据的时距曲线，同时保证每次地震至少被4个台站记录到，其中参加反演的P波绝对到时2539个，相对到时资料33972个，参与反演的台站有24个。台站和地震分布见图2。

本文反演采用的初始模型参考了该区域以往的研究结果及人工地震观测的结果（嘉世旭等，2001：莘海亮等，2011）。表1给出了P波速度初始模型。表1中的上界面深度表示为分层速度模型中各层的顶面深度，层厚度为下一个上界面深度与该上界面深度之差。海平面以上部分与第一层的速度相同。S波速度与P波速度的比值为1.72。

2 双差层析成像结果

濮阳地震集中区主要位于东濮凹陷地区，属于第三系裂陷盆地。东濮凹陷地处豫鲁边界，北起范县，南至兰考，西自濮阳，东到东明，为渤海裂谷系向西的分支。它是在结晶变质岩系基底上发育的地台型构造层系，受到新生代地壳水平拉张作用，逐渐裂解断陷而形成的断陷盆地（王婧韫等，2014）。图3为濮阳地震集中区主要地质构造图，东濮凹陷夹在鲁西隆起与内黄隆起之间，主要受NNE走向的聊城-兰考断裂和长垣断裂控制，其中聊城-兰考断裂带以东为鲁西隆起区，长垣断裂带以西为内黄隆起区。

本文利用双差层析成像方法反演濮阳地震集中区的三维速度结构，对研究区域内的315次地震进行重新定位，最终获得291次地震的精定位结果，重新定位后均方根残差平均值由原来的1.19s降为0.59s，通过对理论到时和观测到时残差的理论估计，震源位置的测定误差在E-w方向平均为O.36km，在N-S方向平均为0.40km，在垂直方向平均为0.49km。

对比重定位前后的震中分布（图4），可以看出双差定位后的震中位置更加集中，丛集性较为明显，优势分布方向为NNE向，与该区域断层的走向方向基本一致。重新定位后震中主要分布于聊城-兰考断裂带一侧，延展范围约20km，地震分布与断层的形态密切相关。胡长和和许坤福（1991）认为聊城一兰考断裂带走向23°～32°，倾角40°～70°，据钻孔资料揭露，该断裂东西两侧上第三系和第四系厚度相差660m，说明该断裂继承性差异运动十分强烈，这种差异运动可能与该地区小震活动有一定的相关性。

从反演结果看，绝大部分地震主要集中在深度15km以上，本文给出上界面深度为4km和10km的速度结构结果。深度为4km和10km速度结构结果基本反映了该地区上地壳速度的横向变化特征，与该地区的地表构造有着一定的关系。从深度为4km速度结构结果（图5a）可以看出，内黄隆起呈现为明显的高速异常，东濮凹陷呈现为低速异常，内黄隆起的高速异常有向东“侵入”的特征。从图5b同样可以看出，内黄隆起呈现为高速异常，高速异常体的形态有一定的收敛，而东濮凹陷呈现为相对更加明显的低速异常的形态，可以清楚地显示东濮凹陷的基底特征。该结果与莘海亮等（2011）所反演得到的4～12km深度层的速度结果比较吻合。聊城-兰考断裂带和长垣断裂带介于高低速异常的转换带上，表明断裂带两侧的介质速度结构存在差异性，这种差异性可能更有利于应力集中，从而产生地震。从整体的反演结果来看，该地区的浅层速度结构与地表构造基本一致。

3 讨论

濮阳地震集中区内的构造特征和断裂展布与速度结构的横向变化有着较好的对应，聊城一兰考断裂带和长垣断裂带的走向与高低速过渡带的形态比较一致，断裂带两侧的速度存在明显差异。深度为4km和10km的水平层析成像结果显示，地震主要发生在聊城一兰考主干断裂带上，并且处于高速异常体内。一般认为高速体与地下较脆、较强的岩体有关，这些部分能够积累更多的孕震能量，内黄隆起的高速异常体向东“侵入”特征可能与地震的孕育密切相关。

AA'剖面（图6a）的速度结构显示，内黄隆起的高速异常向东“侵入”聊城一兰考断裂带下方，可能导致聊城一兰考断裂带下方的应力比较集中，从而增强了断层的活动性。从BB'剖面（图6b）的速度结果可以看出，地震震源处的高低速异常过渡带比较明显，断裂带两侧介质速度结构具有明显的不均匀性，高低速异常体的过渡地带既是应力集中的地方，又是介质相对比较脆弱的地方，这样的环境具备了积累大量应变能的介质条件，容易发生破裂释放应力。

从图6可以清楚得看到，地震主要发生在高速异常体内，主要分布在25km以上，优势深度集中在4～15km之间，断层的破裂宽度约20km，地震分布比较密集，显示出断层的倾角比较陡立，这与于平等（2003）认为的聊城一兰考断面上部视倾角约为70°～80°的结果比较接近。对比速度结构与地震分布可以看出，速度结构的不均匀性可能是控制地震分布的主要因素。东濮凹陷发生过M≥6地震3次，即1502年10月27日濮城M61/2地震、1937年8月1日菏泽M7.0地震和M61/2强余震。濮阳地震集中区与1502年濮城M61/2震的震源破裂区的范围大致重合，而目前该区域内的地震活动的持续增强过程，可能表明历史地震的断层上又重新积累了一定水平的应变（郑建常等，2013）。

为了进一步检验研究区内成像结果的可靠性及分辨率，利用恢复性实验方法进行验证（zhao etal，1992）。采用实际地震资料的反演结果作为恢复性实验的输入模型，可以得到各震相的理论走时，然后将获得的理论走时作为观测值，基于初始模型，反演得到新的成像结果。通过对比新的成像结果和实际地震资料的反演结果，可以反映成像结果的可靠性。图7给出了研究区内的恢复性实验的检验结果。对比恢复性实验结果和实际地震资料的反演结果可以看出，二者速度异常形态非常相似，只有部分异常体的形态和幅值存在较小的差异，这表明成像结果中的主要构造特征是可靠的。

4 结论

本文应用双差层析成像方法反演得到了濮阳地震集中区的速度结构特征，反演结果与该区域的构造特征和断裂展布比较一致。濮阳地震集中区的小震活动主要发生在高速异常体内，地震分布比较密集。内黄隆起的高速异常体向东“侵入”现象导致了聊城一兰考断裂带处的速度结构的不均匀性，可能致使该区的应力比较集中，这表明速度结构的不均匀性可能是控制该区域地震分布的主要因素。

电容层析成像技术研究进展 篇4

ECT技术的基本原理是:根据被测非导电物质具有不同的介电常数,当混合物的组分或浓度发生变化时,引起混合流体的介电常数发生变化,从而引起传感器阵列的电容变化,经过数据采集和信号处理,通过电容敏感场的先验信息,利用电磁场中的反问题进行求解,建立流体内部的介质分布情况。ECT技术广泛应用于电力、石油、化工及医药等两相流及多相流检测领域。ECT技术具有非侵入及响应速度快等优点,被认为是解决连续相为非导电物质的两相流及多相流成像的最有效手段,同时ECT系统能获取复杂工业过程中管道或者容器截面的可视化信息,对于一些特性复杂多变,常规检测方法存在困难的参数检测提供了有效的手段。

国内外就ECT技术研究做了大量的工作,主要包括微小电容检测、传感器结构参数的优化设计、数据采集系统、传感器的敏感场分布、图像重构算法及成像质量评价等,取得了令人鼓舞的成果,笔者主要介绍微小电容检测及数据采集系统、图像重建和ECT系统的典型应用这三方面的研究进展。

1 电容层析成像系统的组成

电容层析成像系统主要由电容传感器阵列数据采集系统、成像计算机及相关的分析解释软件等部分组成,如图1所示。电容传感器阵列包括绝缘管道、管道外壁的检测电极、为减小杂散电容的径向电极及屏蔽电极等;传感器检测电极将管内流体的分布转换成传感器电容(电压)输出;数据采集系统进行A/D转换,将数据送入计算机,计算机通过相应的算法对传感器采集的信号进行图像重构,实现可视化;相关的分析解释软件对管道内介质的状态进行分析,提供相关的参数信息。

2 电容层析成像技术研究现状

2.1 微小电容检测及数据采集系统

微小电容检测是电容层析成像系统的重点环节,ECT系统的测量电容远小于电路中的杂散电容,因此测量电路应有较高的分辨率及较好的抗杂散电容能力;ECT系统的相邻电容和相对电容的值相差近两个数量级,因此要求电路的测量电容动态范围大;ECT系统成像不同于CT系统成像,ECT系统要求较高的成像速度,应满足流体的高速流动。

针对微小电容检测电路的这些特点,研究人员已经研制出多种微小电容检测方法,英国UMIST的Yang W Q和York T A最早研制了交流法电容/电压转换电路的ECT系统[4],Huang S M等利用直流充放电法微小电容检测进行了研究[5],Kühn F T和van Halderen P A设计了有源差分法电容检测电路[6],Fasching G E等用高压交流双边激励法研制了三维流态化成像的ECT系统[7]。在国内,律德财等采用高压370V交流激励电路,验证了电容采集电路具有良好的灵敏度及分辨率等静态特性[8],胡敏等利用集成芯片HT133设计了微小电容测量电路,验证了电路分辨率达0.5fF,且具有很强的抗干扰性[9],陈德运等采用前向补偿法测量微小电容,并以油水两相流为研究对象,验证了方法的可行性[10],周云龙等验证了交流法测量微小电容电路中减小反馈电容能够提高灵敏度[11]。此外,国内在传感器优化设计方面也进行了大量的研究。王化祥等研究了传感器的管道厚度、电极的张角和径向电极的插入深度对ECT敏感场的影响[12],Yan H等对ECT传感器的敏感场进行了三维分析[13],周云龙等分析了不同的电极覆盖率以及屏蔽电极与电极外表面距离对ECT系统电容值和敏感区域的影响[14]。

综上,微小电容的检测主要有高压交流双边激励法、直流充/放电法、有源差分法及交流法等。高压交流双边激励法的系统有很高的分辨率。利于其实现高速实时在线测量,能够通过频率或幅值的调整来提高系统的灵敏度。但该传感器结构设计中无屏蔽层,当有外界电磁场干扰时,干扰信号会通过电缆的屏蔽层作用于内保护极板,产生干扰信号,对微小电容检测不利。直流充放电法采用直流电源激励,信号处理采用直流放大器,放大器的零点偏移是影响电路稳定工作的最大问题,另外此种电路存在注入电荷问题;有源差分法与直流充放电法相比省去了CMOS开关,而这一元件的内部电容比被测电容大得多,虽然通过差分放大器能够部分抵消,但由于器件的离散性,对微小电容检测有很大的影响。交流法是目前较成熟的且应用最多的微小电容检测方法,其原理如图2所示。

图2中Cx为检测电容,Cs1、Cs2分别为激励极板,检测极板对地的杂散电容,Cs1与高频正弦激励信号源并联,它的存在并不产生通过Cx流向运算放大器的电流,对输出的影响可以忽略不计;s2与运算放大器的反相端相连,处于虚地状态,C它的存在也不会对输出产生影响,因而可以消除杂散电容的影响。依据电容的电学原理,正弦波信号的频率越高,对电容信号的检测越有利。此种电路的输入为高频正弦波,根据电路有:

式中w——角频率;

T——时间常数。

根据不同的应用场合选择式(2)或式(3),式(3)的输出与角频率有关,要求正弦波频率一定要稳定。

随着电子技术的发展,ECT系统数据采集的研究经历了从模拟电路到数字电路,串行采集到并行采集,串口通信到USB通信的过程。近几年,York T A等研制了可定制集成芯片的微型ECT系统,该系统由单片机和FPGA芯片控制,数据采集速度为750帧/s,测量精度达到300aF[15]。国内天津大学开发出基于FPGA的数字化电学层析成像系统,数据采集速度可达1 000幅/s,在线成像速度可达120帧/S[16];东南大学也开发了基于DSP、CPLD的数据采集系统。

2.2 图像重建算法的研究

ECT系统的图像重建就是利用测量的电容值求解介电常数分布,因为测量电容值的个数远远小于成像的像素数,即方程数小于未知数的个数,这是一个不定方程组,不定方程组的解不唯一,电容测量值的微小变化会引起图像灰度值很大的变化,方程组的解不稳定,因此图像重建技术是病态问题,需要特殊的算法。

自电容层析成像技术出现以来,国内外对图像重建算法进行了大量的研究。1996年,Isaksen O对ECT图像重建算法进行了综述[17]。近年来研究人员又提出和发展了许多图像重构算法,Yang W Q等给出了一种改进的Landweber迭代算法[18],Fang W和Cumberbatch E给出了Tikhonov正则化方法[19],Soleimani M和Lionheart W R B提出了正则化的Gauss-Newton方法实现ECT的非线性图像重建[20],Yan H等使用多元线性回归模型的正则化算法,用回归矩阵代替了灵敏度矩阵进行图像重建[21],Ortiz-Aleman C等研究了快速模拟退火全局优化算法进行图像重建[22],Marashdeh Q等用神经网络方法进行图像重建,取得了一定的成果[23]。在国内,陈德运等用遗传算法进行了ECT图像重建[24],王化祥等针对图像的不适定问题,研究了灵敏度矩阵的奇异值分解理论,提出正则化共轭梯度法[25],苏邦良等引入同步迭代图像重构技术(simultaneous iterative reconstruction technique,SIRT),针对该算法的收敛速度和图像模糊性,对算法进行了改进,取得了较好的效果[26]。

目前研究图像重建的算法可以分为三大类:非迭代算法、迭代算法和现代智能计算方法。非迭代算法包括线性反投影法LBP (Linear back Projection,LBP)、奇异值分解算法(Singular Value Decomposition,SVD)和Tikhonov正则化方法,这类算法的共同特点是计算量小、图像重构速度快,但成像效果不理想,主要用于在线成像;迭代算法包括迭代Tikhonov算法、Landweber迭代及同步迭代重建技术等,这类算法特点是计算量大,图像重构速度相对慢,成像效果逐步逼近真实,主要用于离线分析;现代智能计算方法包括人工神经网络、遗传算法及模拟退火算法等,这些算法在诸多的领域有应用,通常称为现代智能计算方法。

当前最常用的ECT图像重建算法有线性反投影法LBP、Tikhonov正则化法和Landweber迭代法,下面分别详细介绍这3种方法的原理和应用特点。

2.2.1 线性反投影法

最早用于电容层析成像的算法为LBP算法,在LBP中,电容值被认为是通过不同像素单元内的高电容率介质产生的不同电容值叠加而成的,可以写成关于传感器灵敏度矩阵:

式中C——归一化的电容值矩阵;

G——成像单元;

S——灵敏度矩阵。

由于是直接通过电容率分布,利用灵敏度矩阵线性投影到电容的归一化测量值,因此称为线性反投影算法。因为S矩阵通常不是方阵,不存在逆阵,用灵敏度的转置矩阵代替逆阵进行运算,式(5)可以写成:

为减少LBP算法带来的边缘模糊现象,还需对G设置门限滤波。这种重建算法是最早使用的算法,现在仍在实时成像中广泛使用。

2.2.2 正则化Tikhonov方法

正则化Tikhonov方法是处理病态逆问题的有效方法[27],病态逆问题是由于缺乏足够的信息而导致的,正则化方法是找到一个由先验信息约束的解集,然后从中选择一个解,本质上是最小化式(7)的目标泛函:

式中L——约束矩阵;

μ——正则化因子。

选取合适的μ,可求得正则化的G:

若L取单位矩阵,即为标准正则化:

I是单位矩阵,正则化的算法主要取决于μ,实际应用中根据实验方法来定。

2.2.3 Landweber迭代法

Landweber迭代法主要原理是在数据残差的负梯度方向上对解进行修正,数据残差的梯度表达式为:

迭代公式为:

G(0)是根据测量的归一化电容矩阵C通过LBP算法的成像,α为迭代步长,e(k)为图像与原图像的误差。G(0)的误差较大,为减小误差,再将重建的图像G通过式(11)换算成电容向量,电容向量与实测的电容向量相减得到一个误差向量。此误差向量再经LBP重建成一误差图像。此误差图像与原图像相加以进行图像修正,迭代后,看误差是否符合要求,符合要求停止迭代。Landweber迭代法的特点是该算法有着完备的正则化理论,它的缺点是收敛速度较慢且存在半收敛问题。这一算法在图像重构中有广泛的应用。

2.3 电容层析成像的应用

据国内外文献报道,ECT技术正在走出实验室,已在许多工业过程领域中得到初步应用,其应用涉及多相流动的流型识别、参数检测、气力输送系统和燃烧过程监测。

Xie C G等应用ECT系统对两相流的中心流、环流及层流等典型流型进行了识别,又将ECT系统成功应用于监测油气两相流,并初步给出了气液两相流相含率测量的思路和方法[28],IsaksenΦ等对环状流和层流进行了识别[29];随后,又将ECT技术成功应用于油气水三相分离器分界面的监测,可显示油、水、泡沫和气的界面[30]。

Zhu K等将ECT系统应用于气力输送系统,用来监视流态的发展[31]。Srivastava A等利用ECT技术研究了循环流化床在稳态和非稳态下气体与颗粒行为,指出在循环流化床建模时,需要考虑不同流相之间的关系[32];Dyakowski T等研究了ECT技术在气固两相流监测中的应用,利用ECT监视流化床的流态过渡,同时采用图像处理技术分析了用于流态控制的图像数据信息[33]。

Waterfall R C等用ECT技术对内燃机燃烧室中的火焰状态实现了可视化,重建图像可以确定火焰的面积及对火焰定位,可以监视各种空气燃料比下的燃烧效率[34];在国内,夏靖波等利用ECT系统重建火焰图像,从图像中提取特征值,采用模糊和神经网络控制策略对静态火焰进行闭环控制[35]。

目前已有部分商品化ECT系统,如英国的PTL公司、ITS公司。

3 ECT技术的发展与挑战

尽管ECT技术这些年取得了许多成果,但真正走出实验室,在工业现场得到普遍的应用,还面临着许多技术难题,主要表现在多模态成像、病态逆问题的求解、图像重建质量的评价及三维成像等方面。

3.1 多模态成像

ECT作为检测两相流和多相流的各种参数提供了有效的方法,但也不可避免地具有局限性,如水作为液相的两相流或多相流,水的体积流量就非常关键,水介于导电与不导电之间,不适合作为电容的电介质,需要和电阻层析成像(Electrical Impedance Tomography,EIT)技术结合起来,构成双模态成像,扩大两种模态成像技术的适用范围。目前,可以测量管道或容器截面的EIT和ECT双模态传感器,可以将两种模态传感器安装在管道不同截面上,能独立工作,但测量结果不能反映同一时间、同一截面上介质分布信息。ECT和ERT传感器安装在同一截面上激励模式耦合问题、多模态后的数据融合及图像融合等问题都是研究者努力的方向。

3.2 图像重建算法

在ECT系统的逆问题不适定问题上,多采用正则化的方法来解决,这涉及到正则参数的合理选择,合理有效地选择参数是值得深入研究的问题。在ECT的敏感场矩阵方面,目前所使用的逆问题求解方法都是出于实际经验的应用,而对ECT自身理论的研究还有待进一步的深入。ECT敏感场存在“软场特性”,即静电场电力线在介质变化时发生强度和方向的改变,如何从ECT基础理论中找出使敏感场矩阵变成良性的方法,或使用良态问题逐步逼近病态问题是一个需要攻克的难题。

3.3 重构图像的质量评价

当前对ECT重构图像质量的评价主要在相对图像误差、相对电容值残差及设定图像与重建图像间的相关系数等方面,所有的评价方法都是基于静态图像的评价。ECT系统是对被测断面的动态成像,静态的图像评估不一定满足动态成像的评估。成像质量的有效评估,可以验证系统硬件和软件的改进效果,对整个电容层析成像技术的发展都具有指导意义。

4 结束语

层析成像 篇5

天然地震层析成像相对于传统的反射地震方法而言是一种新的经济的勘探方法.这是由于天然地震层析成像所需的观测值直接来自于研究区下方发生的天然微地震,而反射地震却需要在研究区表面进行人工放炮.因此本工作是将天然地震层析成像方法应用于柴达木盆地西部某油田约100km2地区的深层构造的尝试性研究.626个地震事件的3289个P波到时的初步结果与研究区已知的.大的构造吻合较好.该模型中非常显著的特征就是观测到一个北西向的背斜.此外,微地震的分布也与研究区中活动断裂带的位置基本一致.

作 者：王小凤 冯梅 史大年 马寅生 陈宣华 区明益 霍光辉 王连庆 田晓娟 张西娟 李会军 李国岐 蒋荣宝 WANG Xiao-feng FENG Mei SHI Da-nian MA Yin-sheng CHEN Xuan-hua OU Ming-yi HUO Guang-hui WANG Lian-qing TIAN Xiao-juan ZHANG Xi-juan LI Hui-jun LI Guo-qi JIANG Rong-bao  作者单位：王小凤,冯梅,马寅生,陈宣华,区明益,霍光辉,王连庆,田晓娟,张西娟,李会军,李国岐,蒋荣宝,WANG Xiao-feng,FENG Mei,MA Yin-sheng,CHEN Xuan-hua,OU Ming-yi,HUO Guang-hui,WANG Lian-qing,TIAN Xiao-juan,ZHANG Xi-juan,LI Hui-jun,LI Guo-qi,JIANG Rong-bao(中国地质科学院地质力学研究所,北京,100081)

史大年,SHI Da-nian(中国地质科学院矿产资源研究所,北京,100037)

刊 名：地质通报  ISTIC PKU英文刊名：GEOLOGICAL BULLETIN OF CHINA 年，卷(期)：2006 25(9) 分类号：P61 关键词：天然微地震   层析成像   深层构造   地震事件   活动断裂带  

★ 均台辞,均台辞陈与义,均台辞的意思,均台辞赏析

★ 防台防汛事故应急准备与响应预案

★ 精河Ms5.1级地震前精河台的定点形变异常的初步研究

★ 接触网避雷器性能在线监测系统的研究与设计论文

电阻层析成像逆问题的仿真研究 篇6

1.1 电阻层析成像逆问题分析

ERT逆问题就是图像重建算法。关键在于用测量的边界电压重建和反映出物场内部电阻率分布, 成像算法的关键在于图像重建的速度和精度。目前ERT逆问题的求解过程中存在如下几个问题:场域内电阻率和边界电压都为非线性分布, 并且敏感场为软场;独立测量数据的有限性导致很多的不确定性;逆问题的不适定性。

ERT逆问题的非线性和不适定性给其求解带来很大困难, 主要变现在ERT逆问题固有的病态性、信息量小和计算量大等方面。考虑到这些因素, 为了提高重建图像的质量, 在进行ERT图像重建时可采取以下几个途径:降低病态程度;提高数据采集精度并消除数据中噪声;增大采集数据量;采用数值稳定性好的图像重建算法;对算法进行正则化改进。

1.2 电阻层析成像逆问题图像重建算法种类

ERT的图像重建算法可以有很多分类方法。根据所利用的测量数据分为动态式和静态式算法, 主要有以下三种方法:1) 反投影法。2) 灵敏度系数法。3) Newton-Raphson法。

M atlab实现重建算法步骤:本文主要实现修正的牛顿——拉夫逊类算法 (ModifiedNewtonRaphson算法, 简称MNR算法) 。它是基于解非线性最小二乘 (least-squared) 问题的著名的Gauss-Newton算法及其改进形式的最优化思想的静态电阻图像重建算法。它的算法流程如图1所示, 具体步骤如下:初始电阻率分布估计值设定为, 进入迭代;求解正问题得到边界电压计算值。假设第k步迭代的电阻率分布估计值为, 则计算得电压值

判断迭代次数是否大于N, 若大于N则终止迭代。若小于N, 将计算得到的边界电压值与实际边界电压测量值V0相比较进行迭代终止判别, 判别准则为:

当式1) 满足时, ρk!"即为所求电阻率分布, 迭代终止, 否则进入下一步对ρk!"进行修正;对ρk!"进行如下修正:

2 仿真实验数据分析

本次研究对Tikhonov正则化的牛顿-拉夫逊算法的图像重建效果进行了比较, 首先设定电导率的分布, 然后计算相应的测量电压, 再用仿真到的测量数据进行图像重建。设连续相介质的电导率设为0.5mS·cm-1, 离散相介质的电导率设为1mS·cm-1, 采用16电极相邻模式, 电极所对圆心角△=22.5°, 激励电流10mA。考察层流介质分布情况, 仿真结果如下。

层流介质图像重建效果, 如表1:

从图1、2中可以明显辨别流型的分布, 可以得出结论Tikhonov正则化的牛顿-拉夫逊算法对层流介质分布分辨较清楚, 即对层流介质分布的成像能力较好。

3 结论

电阻层析成像技术具有成本低、响应快、可实现非侵入测量的特性, 是一种优越的高新检测技术。本文的主要工作和结论如下:对ERT软场特性及逆问题的不适定性进行理论分析, 指出改进方法;针对逆问题, 利用牛顿-拉夫逊算法对其求解, 根据逆问题不适定性的改进方法, 根据Tikhonov正则化的牛顿-拉夫逊算法原理建立了ERT图像重建软件包, 并进行了大量的仿真分析。

参考文献

[1]魏颖.电阻层析成象技术 (ERT) 及其在两相流测量中的应用研究[D].沈阳:东北大学, 2001.

[2]胡海涛.电阻层析成像图像重建算法的研究[D].哈尔滨:哈尔滨理工大学, 2009.

[3]HAN Y H, YANG S H, JIN H B.Electrical resistance tomography for flowcharacterization ofvariousgas-liquid system[C].Book ofabstractsofthe 12th Asianpacific confederation ofchemical engineering congress, Japan, 2008.

[4]余金华.电阻层析成像技术应用研究[D].杭州:浙江大学, 2005.

光学相干层析成像旋转失真评价研究 篇7

关键词：医用内窥镜,光学相干层析成像,旋转失真,图像处理,质量控制

0 引言

光学相干层析成像(Optical Coherence Tomography,OCT)诞生于1991年[1],是生物医学光学领域较为成熟和常见的一大类可用于临床的技术,经过20多年的发展,由时域发展到了频域[2,3,4,5],成像方式由传统的体外检查(眼科)演变出了不同规格和用途的介入式光纤探头[4,5,6,7],成像速度普遍提高到了100帧/s量级以上,最快可达4000帧/s[8]。随着光源的发展和光路设计的进步,分辨率由10µm量级提高到了1µm量级[9],图像细节日益丰富,可提供更多的病理信息。在内窥成像领域,OCT的潜力巨大,国外学术界和工业界在心血管、消化道、呼吸系统、泌尿系统等方面[10,11,12,13,14]积极开展临床试验和推广,部分产品已经商业化,国内也有科研机构和企业在跟进。

目前,OCT内窥成像的空间扫描主要通过机电装置控制光纤探头进行螺旋式后退实现。在实际的系统设计和成像过程中,由于控制精度、稳定性、摩擦力、扭矩等方面的问题或故障,探头运动的周期性和连续性不理想,对应的图像失真叫做旋转失真,是OCT内窥成像的重要误差来源,也是评价机电装置运动和控制性能的重要参数。

目前,医疗器械质控领域缺乏评价OCT图像旋转失真的方法。OCT内窥成像领域当前并没有对应的国内外标准。OCT领域唯一的国际标准ISO 16971-2015[15]也仅针对眼科成像(体外扫描)。现有的对于OCT旋转失真的评价借鉴血管内超声的做法[16],通过在体模上每隔45°进行一个标记,评估实际图像中标记的偏移量来进行。由于不同OCT探头尺寸、形状、参数设定各有差别,设计通用体模检测旋转失真的难度较大,不适合推广。OCT学术界在开发运动伪影校正算法时关注相邻图像之间的相似性,提供了有益的启示[17,18,19,20]。

为了探讨通过图像分析直接评价旋转失真的可行性,本文通过实验观察了OCT成像中的旋转失真现象,针对相邻图像直接的相似性,计算了帧内相邻扫描线之间以及相邻两帧图像之间的差异,比对了旋转失真的不同量化表述方法,对于完善OCT内窥系统的质控和促进我国OCT内窥成像技术的发展提供了参考。

1 实验材料和方法

1.1 基本数学表述和实际应用案例

理论上,OCT探头应进行匀速螺旋运动。在每一帧OCT二维图像上,如果以OCT探头为柱坐标系原点,那么不同的OCT轴扫描线(A-line)对应不同的角度。以常见的扫频(Swept Source)OCT内窥系统为例,成像过程应满足以下关系:

其中,vA是光源的扫描频率(单位:Hz),N是每一帧图像的轴扫描线数,vr是探头旋转的速度(单位:转/s)。两幅图对应的螺距P满足:

其中,vp为水平后退的速度。

公式(1)、(2)是描述探头运动的基本公式。图像旋转失真在数学本质上意味着不符合公式(1)。旋转失真图像典型案例,见图1。两个例子的共同特点是旋转失真区域与周围正常图像之间有突变。本文的体模实验需要捕捉类似的突变,以分析图像特征。

注:a.冠状动脉OCT图像,虚线区域对应的组织是一个血管分支。由于旋转失真的发生,分支与冠状动脉之间的边界丢失,因而呈现出0点钟方向的跳变;b.食道OCT图像,虚线区域显示的是探头迟滞形成的模糊图像,也是旋转失真的一种。

1.2 实验材料

本实验使用聚乙烯和硅胶制作而成的软管充当体模,近似模拟人体内部的空腔结构(例如血管、气管等),其内径范围为1~2 mm,厚度为200~400µm,长度为50 cm。

1.3 仪器与方法

本实验使用商用扫频OCT内窥系统配合光纤探头进行测量,OCT光源的频率为40 k Hz,波长扫描范围为1250~1370 nm。每帧图像包含1024条扫描线,探头每秒钟旋转39次,因而系统成像速度为39帧/s。OCT探头远端采用球透镜设计,焦距为1 mm,横向空间分辨率为30µm,轴向空间分辨率为14µm。探头外面有透明塑料保护鞘,厚度约100µm,外径800µm。OCT体模成像实验示意图,见图2。

在实验中,保护鞘和体模位置固定,探头以10 mm/s的速度后撤,同时每秒旋转39转,从而对体模内部进行三维螺旋扫描,以模拟OCT在组织中的成像过程。根据实验参数设定,探头旋转的角分辨率为2pi/1024,对应体模上的圆弧长度为6~12µm,低于横向空间分辨率(30µm),这意味着相邻两条扫描线对应的组织存在重叠,扫描线本身的相关系数较高。为增加探头运动的阻力和需要的扭矩,体模软管经过一定的弯曲处理,曲率半径最小为1 cm。

1.4 数据处理方法

为规范数据格式,本文在柱坐标系下(原始OCT图像为矩形而不是圆形)提出两种评估旋转失真的思路,分别为帧内2-范数评估法和帧间2-范数评估法,其原理见图3。N为每一帧的扫描线数,M为每条扫描线上的采样点数,Aj指的是第j条扫描线。

帧内2-范数评估法指的是在每一帧图像里,计算相邻两条扫描线之差的2-范数,根据相邻扫描线的变化来跟踪旋转失真,由以下公式定义:

其中,Ij为第j条扫描线与第j+1条扫描线之差的2-范数,ai,j为第j条扫描线上第i个采样点,以此类推,最终Ij会组成一个矢量VI。帧内旋转失真通过两种不同途径进行评价:一是直接取VI的最大值VImax(帧内最大),作为每一帧内旋转失真的最大值;二是直接求VI的2-范数VIN2(帧内总体),作为每一帧内旋转失真的总体估计(也等效于平均值)。

帧间2-范数评估法指的是在相邻图像之间,计算相同角度对应两条扫描线之差的2-范数,根据帧间对应位置扫描线的变化来跟踪旋转失真,由以下公式定义:

其中,Bj为第j条扫描线与下一帧第j条扫描线之差的2-范数。Bj最终组成一个矢量VB。与上面类似,帧间旋转失真也通过VB的最大值VBmax(帧间最大)和2-范数VBN2(帧间总体)分别进行评价。

在进行计算之前,所有图像预先进行平滑滤波,以消除环境中的高频或脉冲噪声等对图像的干扰。

2 结果

2.1 体模中的旋转失真图像

在实验中,OCT旋转失真现象多次被观察到,见图4。

注:a.显示了来自聚乙烯软管的正常的OCT图像,其中最外侧的灰色环形为体模图像,内侧灰色环形为探头保护鞘的图像,两者均灰度均匀,过渡自然平滑;b.显示了具有旋转失真的OCT图像,在0点钟方向体模图像出现了明显的错位,以至于圆环不能闭合,与实际情况不符。说明本实验方法和平台可以产生旋转失真,从而提供数据比较旋转失真的评价方式。

2.2 旋转失真评估结果

本实验采集了硅胶软管体模的560帧连续图像,根据公式(3)~(6)分别计算了帧内最大、帧内总体、帧间最大、帧间总体四个矢量各自除以各自的最大值,得到归一化曲线的对比,见图5。

注:蓝色曲线代表帧内最大,主要分布在0.1~0.4之间,总体平稳;黑色曲线代表帧内总体,基线在0.9,整体平稳,对比度较小;绿色曲线代表帧间最大,基线在0.2附近,尖峰较多,说明在这些位置两帧之间在局部角度上差别变化大;红色曲线代表帧间总体,基线在0.6附近,同样伴随着较多的尖峰,说明从整个数据集合中帧与帧之间存在着明显的差别。

3 讨论

从图5的趋势来看,帧内最大和帧内总体两种方法反映的旋转失真情况截然不同,需要进行具体讨论,以61帧和178帧的对比为例,见图6。

图6(a)显示的是第61帧,其成因是成像过程中,OCT光线受到杂质碎屑的阻挡,在图像上形成了伪影(黄色星号标记)。图6(b)显示的是第178帧,在0点钟方向明显有旋转失真造成的图像突变(黄色星号标记)。肉眼观察下两帧均呈现类似的突变,但是帧内最大有明显的不同,前者为0.29,后者为1(帧内最大曲线的峰值)。究其原因,前者由于光线的绕射,在被阻挡区域和正常管壁之间依然存在过渡,后者直接出现突变,这说明帧内最大有助于区分旋转失真和由于光线阻挡造成的伪影。相反,61帧和178帧的帧内总体差别不大,分别为0.91和1,这说明由于平均效应,帧内总体难以捕捉旋转失真的变化。

另一方面,帧间总体和帧间最大两种评价结果曲线与帧内最大有较大差异,尖峰很多。经过检查具体图像,发现这两种评价结果反映的主要是组织不同位置之间的差异,而不是旋转失真本身,以第270和271帧为例(图7),对应的帧间总体和帧间最大分别为0.80和0.70。可以看出,两个图像自身并没有明显的旋转失真,组织本身差异较大,尤其是6点钟方向。本次实验中探头后退的螺距为256µm(10 mm除以39),意味着两幅图的间距已经远远超过了光斑直径(30µm),所以对应体模中的位置已有明显偏差,与成像结果相吻合。因此,帧间总体和帧间最大在这种设定下难以用于评价旋转失真。在实际的医学成像过程中,为缩短诊断时间,提升患者舒适度,探头后退速度和旋转速度都在不断提升,但缩小螺距本身依然面临很大困难,这意味着帧间评价方式同样很难应用。

4 结论

层析成像 篇8

关键词：医用光学,模拟眼,光学相干层析,三维分辨率

0 引 言

自从上世纪90年代首次用于离体视网膜检测以来[1],光学相干层析(OCT)成像技术已经获得了令人惊叹的进步和广泛的应用。尽管OCT在检测皮肤、牙齿、心血管、脑成像等医学领域都开展了研究,但至今为止,眼科仍然是OCT技术最匹配、应用最成熟的领域[2,3,4,5,6,7,8]。众所周知,OCT技术早已在美国、欧洲和亚洲等地获得临床诊断许可,但是与之相关的国际标准却迟迟没有发布。在OCT这个领域里,不管是制造商、临床使用者还是第三方检测机构都迫切需要一种有关OCT的标准测试方法与评价机制,从而确保注册检验、日常质控以及产品比对的顺利进行。

为了确保OCT设备的日常质控,可以定期检测一些关键的物理参数,如光源光谱和光束特性等。另一种测试方法就是采用模体检测,这种方法不依赖OCT的工作原理或某个特定内部模块,而且不用拆机属于非破坏性检测,其结果能更客观更全面的反映OCT设备的性能指标。国际上有不少实验室已经在从事OCT模体方面的研究,主要用于评价OCT的成像分辨率性能。英国国家物理实验室(NPL)的Tomlins和他的同事对OCT的横向与轴向分辨率进行了详细的理论论述,并采用透明树脂掺杂二氧化硅微球的办法制作了OCT点扩散函数模体用于评价OCT设备的三维成像分辨率性能[9]。美国食品药品监督管理局(FDA)的Agrawal等人在实验室研制了一款掺杂纳米尺寸颗粒的模体,用于评价OCT设备的三维点扩散函数[10]。美国国家标准与技术研究院(NIST)的Chang等人也尝试利用化学方法制作了层状光学组织模体用于测试OCT设备的轴向分辨率和成像对比度[11]。除了这些标准或监管机构的实验室之外,英国肯特大学的Avanaki等人也基于米氏散射理论研制了环氧树脂模体[12],并尝试掺杂了两种不同的散射颗粒:聚苯乙烯微球和金微球。针对眼科OCT设备,一些仿真视网膜的组织模体也成为了研究热点。Rowe和Zawadzki设计并制作了一种模体[13],包括5层透明组织,每层组织厚约60 ?m但具有不同的光学折射率,不仅如此,视网膜的中心凹也被仿真出来了。De Kinkelder等人也研制了一款层状模体用于仿真视网膜的生理结构[14],且具备很高的散射对比度。Baxi等人采用分层旋转涂覆掺有纳米颗粒的硅树脂,然后通过激光微刻蚀的方法去改变表面形貌,最终展示了一种视网膜模体[15],既仿真了视网膜厚度及近红外光学特性,也模拟了视网膜中央凹区的表面形貌。由于普通的平面加工工艺很难制作出微米尺度的三维图案去检测OCT设备的三维分辨率性能,一些新型的制造工艺如飞秒激光亚表面微刻技术[16]或改进后的刻蚀浇注方法[17,18]都被用于尝试制作OCT设备的三维分辨率测试工具。

本研究中,我们设计并制造了一种能够为眼科OCT设备提供三维分辨率测试图案的新型模拟眼,该测试图案的设计参考了USAF1951靶板并通过3D打印工艺制作模拟眼的眼底。不同于传统的USAF1951靶板的二维平面图案,新型模拟眼的测试图案都是三维的,不仅能够检测OCT设备的横向分辨率,而且也可用于评价其轴向分辨率。通过检测这种新型模拟眼,OCT设备的横向与轴向分辨率性能检测可以同时实现。为了验证这种新型模拟眼的有效性,使用一台科研级和一台临床用OCT设备分别对该模拟眼及三维分辨率测试图案进行了检测,研究结果表明该模拟眼能够对OCT设备的三维成像分辨率性能进行初步评估。

1 模拟眼的设计与制造

1.1 模拟眼设计

如图1所示,模拟眼的设计完全与真实人眼的光学结构相符,一些关键的光学元件如眼角膜和晶状体的结构都在该设计中实现。正常成年人眼的前后径平均为24 mm,垂直径平均为23 mm,最前端突出于眶外12～14 mm,为了简化模型,我们将包括视网膜的后半眼球设计成直径为24 mm的半球体。眼内腔包括前房、后房和玻璃体腔,眼内容物包括房水、晶状体和玻璃体,三者均透明,并与角膜一起共称为屈光介质。眼球最前面的光学结构是角膜,直径大约12 mm,其中垂直径略小于水平径。角膜的前后表面可以被近似的认为是球面,前表面的曲率半径约为7.8 mm,后表面的曲率半径约为6.8 mm,角膜中央区的厚度约为0.5～0.6 mm。

位于晶状体前,由虹膜构成的小圆孔,是人眼的孔径光阑,也叫瞳孔,它可以根据物体明暗调节进入人眼的光量大小,人眼瞳孔的直径可变动于1～9 mm之间。为了测试的方便,我们通过表面喷漆和一个橡胶圆环实现了7 mm直径的人工瞳孔。人工晶状体呈双凸透镜状,前表面曲率半径约12 mm,后表面曲率半径约为6 mm,中央厚约4 mm。真实人眼的晶状体富有弹性,当需要眼睛调节以看清远距离或近方物体时,主要通过改变晶状体前后表面的曲率半径来实现。

置于模拟眼眼底的三维分辨率测试图案参照了USAF1951靶板的设计,但与之不同的是,其所有的图案都是三维立体的。图案包括两部分:横向分辨率测试图案与轴向分辨率测试图案。其中横向分辨率测试图案包括6组测标,每组测标由垂直排列的三条纵向与三条横向的长条图案组成。所有的长条图案拥有相同的凸起高度100μm,长宽比5:1,但其宽度每组各不相同,从小至大依次为20μm,50μm,100 μm,200μm,300μm和500μm。纵向分辨率测试图案也包括6组测标,每组测标是一个矩形图案。所有的矩形图案拥有相同的长度(2 mm)和宽度(0.5 mm),但其凸起高度每组各不相同,从小至大依次为20μm,50μm,100 μm,200μm,300 μm和500μm。

1.2 模拟眼制造

为了方便模拟眼的制造工艺与后期装配,整个模拟眼被分成三个模块分别加工:眼前节部分、晶状体和视网膜半球体。眼前节部分和晶状体采用透明树脂通过注塑工艺成型,而包含三维分辨率图案的视网膜半球体加工则使用3D打印技术来实现。树脂注塑工艺提供了很好的光学透射性能,这一点对眼前节和晶状体模块至为重要,但这种工艺的缺点是需要使用模具。而3D打印技术不仅不需要任何模具,而且可以突破平面制造工艺的局限,能够制造复杂的三维结构。本研究中所设计的三维分辨率测试图案就是一种较复杂的三维结构,不仅长宽尺寸不同,而且高度各异,再加上又分布在一个半球面上,如果采用传统的平面制造工艺很难实现。因此,本研究采用了一种高精度的光敏树脂作为加工材料,一层一层的铺设,最小层厚可达16μm,使用紫外光对每层光敏树脂材料进行选择性的固化,从而最后实现整个视网膜半球体的加工。

三个模块都完成加工后,将眼前节与视网膜半球体模块都浸入到水中,在水面下将人工晶状体通过一个橡胶圆环固定在眼前节模块上,然后再与视网膜半球体模块通过提前设计好的卡槽进行衔接密封,如图2。通过这种水下装配的方法,确保了前房和玻璃体的中空结构都被水充满,从而模拟真实人眼的房水与玻璃体环境。不可否认,这种水下装配工艺稍显复杂,在后期的工艺改进上,我们在眼前节与视网膜半球体模块上分别都设计了一组进液与出液口,这样一来,眼球模型可以在空气环境中完成装配,模拟房水与玻璃体的液体可以在装配结束后通过进液口注入到前房与玻璃体结构中,既简化了装配工艺,增强了模型密封性,又增加了模拟体液的选择灵活性。

2 测试结果

2.1 OCT 系统

测试使用了一台科研级和一台临床用OCT设备,两台设备都是谱域OCT系统,其结构原理如图3所示。科研级OCT系统的工作中心波长是1 310 nm,通过结合两台超辐射发光二极管(SLD)光源得到近170 nm的带宽。此外,该OCT系统还配置了一套可见光探测的CCD系统,用于实时观察并对照OCT所检测的样品区域。该OCT系统的轴向分辨率和横向分辨率分别为6.5μm和13μm。临床OCT设备使用的是日本拓普康公司的3D OCT-1000,该谱域OCT系统的轴向分辨力可达5 μm,横向分辨力优于20μm。

2.2 OCT 测量结果与讨论

采用科研OCT系统对视网膜半球体上的轴向分辨率图案进行了检测,得到的B-scan图像如图4(a)所示,有四组高度顺序降低的图案沿着弧面清晰可见,但图像右侧高度最小的两组图案很难观察到,造成这种情况的原因很可能是所使用的3D打印工艺的加工精度不够。本研究所使用的3D打印工艺的精度在深度方向是16μm,在平面的X与Y方向上是42μm,而我们所设计的轴向分辨率图案中最小的两组的深度依次是20μm和50μm,比较接近该3D打印工艺的精度极限,正因为这个原因,在图4(a)中几乎观察不到高度最小的两组图案。值得注意的一点是,从图4(a)中能观察到的四组轴向分辨率图案的剖视图并不是很理想,上表面与两个侧面都不平整,与预先设计的沿弧面分布的矩形轮廓有差距。这样表面轮廓不规则的图案本身的几何尺度参数难以标定,如果用来作为测试或检验的标准,势必会造成困难。因此,未来的研究将会重点关注如何提高3D打印工艺的精度与稳定度,确保制造能得到足够精细的图案,并且形状规则,本身的尺度参数工艺可控且稳定。

图4(b)是采用临床OCT对装配好的模拟眼进行检测得到的视网膜三维图像,受限于临床OCT系统的观察视野,只能看到部分横向分辨率和轴向分辨率的图案。考虑到该模拟眼较复杂的整体结构,能得到清晰的OCT眼底图像证明了这种工艺的可行性,显示了采用3D打印技术制造OCT模拟眼的巨大潜力。

(a) 使用科研 OCT 系统得到的轴向分辨率测试图案的 B-scan 图像; (b) 使用临床 OCT 系统得到的眼底分辨率测试图案的三维重构图像

3 结 论

层析成像 篇9

电阻层析成像 (Electrical Resistance Tomography, 简称ERT) 技术是电学层析成像技术存在的四种基本形式之一, 是过程层析成像技术中的一个分支, 它与电阻抗层析成像 (Electrical Impedance Tomography, 简称EIT) 技术最重要的一点区别就是在重构场域内部的阻抗分布时, 电阻层析成像技术只保留了实部信息。因此从一定程度上讲, 电阻层析成像技术是电阻抗层析成像技术的一种简化形式。该技术中, 由于采用的是安全电流激励, 属非侵入式检测技术, 在研究人体生理功能与疾病诊断方面具有很重要的临床价值[1,2]。

为了提高正问题的计算效率, 目前国内外学者在有限元编号优化方面进行了大量的研究。本文首先通过直角坐标系对平面区域划分为四个象限, 在传统的电阻层析成像拓扑结构的节点编号方法上进行改动, 提供了一种新的节点编号方法, 该方法虽然在降低整体刚度矩阵带宽方面的效果不及传统方法, 但它提供了一种节点编号方法的思路, 并在此基础上研究出一种启发式全新编号方案, 有效降低了整体刚度矩阵带宽。

1 带宽的定义与优化方法

设有限元网络系统的节点总数为n, 有限元方程组的整体刚度矩阵K是一个对称方阵, 定义fi (K) 为K中的第i行最左边非零元素的下标, 那么有:

矩阵K的半带宽定义为:

半带宽β (K) 与矩阵带宽γ (K) 的关系为:

对于带状线性方程组AM×MX=B, 其计算时间受Mβ (K) 2影响[3,4,5,6]。

目前, 国内外学者对于该问题的研究普遍采用的一个方法就是从图形的几何中心节点出发, 逆时针方向按照层数由内而外进行编号。根据这种编号方法编写出的程序有一定的规律性, 编程思路也较为清晰, 但该方法在降低带宽方面的效果并没有从理论上证明是最优的。本文提供两种编号思路: (1) 将平面网格图形通过平面直角坐标系分成为四块区域, 从横轴的正向开始, 逆时针依次为每个象限的节点进行编号。 (2) 基于传统编号方案及编号思路 (1) , 对所有节点划分不同的所属带, 然后按照“S”形对所有有限元节点按照带的顺序依次进行编号。

2 象限划分式编号方法

根据特定网格图形的特点, 在传统方法的启发下进行改动, 其八层均匀网格图形的节点编号具体过程如下:

步骤1:首先将圆形网格图形按照平面坐标系划分为四个象限。其中横轴和纵轴的正半轴分属于第一与第二象限, π方向上的节点属于第三象限, 3π/2方向上的节点属于第四象限。

步骤2:借鉴传统的逆时针由内而外按层编号方法, 本方法亦采用逆时针顺序编号。但基于象限的编号顺序之上, 对每个象限内节点进行编号的规则, 亦采用逆时针由内而外的按层编号方法。

步骤3:以中心原点作为起始节点, 从第一象限开始标号, 则每个象限的节点数为72, 因此网格内所有节点数即为289, 记作total_number, 故最终将产生一个289阶的整体刚度矩阵, 而节点编号优化的最终目的, 正是使该矩阵的带宽有效降低, 以提高正问题计算效率。

步骤4:编写程序对每个节点进行编号。最终将所有节点信息放入total_number×3的矩阵中, 第三列为节点序号, 第一二列分别为该序号节点的横纵坐标值, 将该矩阵记为jnn。

3 启发式全新编号方法

基于原有网格的拓扑结构, 在总结本文已有编号方案的优缺点基础之上, 创新出了一种全新的编号方案。该方案不仅编程简单, 而且在降低整体刚度矩阵带宽值、提高电阻层析成像正问题计算效率方面均比传统编号方案理想。

步骤1:在原有拓扑结构的基础之上, 将图形网格所有有限元节点的信息放入total_pointnumber×3的矩阵str中。

步骤2:根据横坐标值对网格图形划分不同的节点所属条形带, 从左至右分别编号, 共16个带状区域。

步骤3:将所有节点按带分类。属同一横坐标值范围的条形区域内的节点, 划为同一类, 即同一带。

步骤4:根据带域编号对同一类节点进行排序。由于启发式全新方案采用“S”形顺序进行编号, 即若带号为奇数, 则其内节点按纵坐标从大到小依次排列;若带号为偶数, 则其内节点按纵坐标从小到大依次排列。

步骤5:将矩阵str中的节点信息按照带的序号以及步骤4中的规则进行重新编号, 并整理到矩阵pnn中, 则该矩阵中的数据即为在启发式全新编号方案下的节点编号信息。矩阵中, 第一列元素为所有节点的横坐标值, 第二列元素为所有节点的纵坐标值, 第三列元素为所有节点在传统编号方法下的序号值。

4 带宽优化结果比较

图1为传统有限元节点编号方案, 即“由内及外, 逆时针旋转”的编号方法及相应的整体刚度矩阵带宽效果图。图2为基于象限划分式和启发全新式的节点编号示意图, 图3为以上两种编号方法下的整体刚度矩阵带宽效果图。

由图1、图2、图3可知, 与传统编号方案比较, 按象限顺序进行有限元节点编号的方法对于降低整体刚度矩阵带宽及提高正问题计算效率无积极作用, 只是在理论上提供了一种新的网格节点编号思路;而基于该方法之上的全新式节点编号方法不仅在编号次序的选择上实现了创新, 而且在吸取了传统编号方法优点的基础上, 明显降低了整体刚度矩阵的带宽值, 实现了提高计算效率和缩减占用计算机内存空间两方面的创新。

为进一步验证启发式全新方案在有限元节点编号优化方面的效果, 现将同一拓扑结构下的节点数目和整体刚度矩阵带宽值列写如下, 见表1。

5 结论

在对电阻层析成像技术的研究中, 国内外相关研究人员较为常用的一个方法就是“以图形中心为初始编号节点, 逆时针方向按照层数由内到外”的编号规则。实验证明虽然该方法有很多优点, 但它在降低整体刚度矩阵带宽这一人们最关心的问题上并没有实现最优。本文在研究了传统方法并吸取其经典思想后, 首先提出了一种依象限划分网格进行有限元节点编号的方法。虽然该方法被证明在降低带宽方面效果不及传统编号方案, 但基于该方法的思想并吸取其优点, 最后创新出了一套全新的有限元节点编号方法, 该方案可有效降低矩阵带宽值, 从而可在一定程度上提高正问题计算效率。

摘要：有限元网格节点编号规则对整体刚度矩阵的带宽有着非常重要的影响, 而正问题计算效率在很大程度上取决于整体刚度矩阵的带宽值。为了减少电阻层析成像中有限元运算时间, 节省计算机的内存空间, 本文在仔细研究了目前相关学者常用的一些编号方法基础上, 提出了一种启发式全新“S”形编号规则, 在一定程度上降低了整体刚度矩阵带宽值, 减少了数据存储量, 大大提高了正问题的计算效率。

关键词：有限元,节点编号,带宽,电阻层析成像

参考文献

[1]肖理庆, 王化祥.改进遗传算法的ERT有限元网格节点编号优化[J].计算机工程与应用, 2011, 47 (7) :20-22.

[2]肖理庆, 邵晓根, 王琳琳等.电阻层析成像有限元仿真模型分析与设计[J].仪器仪表学报, 2008, 29 (2) :354-360.

[3]欧阳兴, 陈忠奎, 施法中.有限元网格节点编号[J].北京航空航天大学学报, 2002, 28 (3) :339-342.

[4]孟宪海, 李吉刚, 蒋丽等.PEBI网格节点编号优化方法研究[J].计算机工程与应用, 2010, 46 (24) :197-200.

[5]郭晓霞.四边形有限元网格重分技术研究及软件开发[D].太原:太原重型机械学院, 1999, 1-53.

层析成像 篇10

折射层析成像是基于对初至波进行射线追踪反演, 构建相应的速度层析成像图, 由此确定地质体中速度异常。与传统的折射波法相比, 折射层析成像最大的优点是在纵向梯度变化的速度层、强烈的水平速度变化地层、大倾角地层、隐伏层和任意地形起伏下, 都可以通过折射波初至获取复杂的地下结构, 给出较好的反演结果。

1 折射层析成像基本原理

折射层析成像法是通过对初至波进行射线追踪反演, 构建相应的速度层析成像图, 由此可以确定地质体中速度异常体的大小、位置、物性等参数。

(1) 初至拾取。在地震记录中实际观测到的初至包含直达波、回折波、折射波或三者组合的初至波, 主要在近地表层进行传播, 一般能量较强, 便于识别, 且可追踪性好, 其走时包含了近地表层介质的速度信息。初至的正确拾取是折射层析反演工作的关键, 起跳点的识别是初至拾取的难点所在, 关系着反演结果的可靠性。图1左侧窗口为初至拾取窗口, 右侧窗口为走时曲线显示窗口, 红色为根据地震记录拾取初至后的时距曲线, 绿色为根据观测系统某次正演计算得到的时距曲线, 两曲线通过不断正演, 反演计算后, 达到一定误差精度要求后便构建出相应的速度层析成像图。

(2) 速度模型的建立。层析反演过程中, 初始模型的选取至关重要, 初始模型与实际结果越接近, 层析反演计算速度越快, 效果越好。过于简化的模型可能使结构中有意义的信息被忽略, 复杂的模型可能使反演的不确定性增强, 同时可能引入虚假信息。在实际计算时, 初始模型应该根据测井资料或速度谱等资料来确定。实际反演过程中, 可能先验信息不足, 给定与实际结果相近的初始速度模型难度较大, 一般是通过给定常速初始模型进行反演计算。层析反演过程中, 需对速度加以约束限制, 速度约束值范围越接近实际速度的变化范围, 计算速度越快反演效果越好。在实际计算时, 应该根据测井资料以及速度谱等资料来确定速度约束值。

(3) 射线追踪正演。本次使用Geogiga公司的折射层析成像模块进行处理, 方法的关键是对模型进行射线追踪, 程序中使用的射线追踪方法基于惠更斯和费马原理, 联合使用最短路径法和弯曲射线法, 首先用最短路径法寻找全局最短路径, 然后再用弯曲射线法寻找局部最优, 该方法运算速度较快且准确性高。

(4) 迭代反演。折射层析成像的实质是反演, 即根据初至时间推断地下速度结构。从反演角度讲, 网格剖分较大时每个网格内通过的射线越多, 越利于反演, 但从正演角度讲, 网格剖分过大则会降低正演精度, 从而直接影响反演的精度。实际计算时, 网格剖分要同时兼具两者才能取得更好的反演效果。

2 工程实例分析

2.1 工程概况

某水电站左、右岸江面以上均分布Ⅱ级阶地, 河床总体上河谷相对开阔, 呈基本对称的“V”字型峡谷, 为斜向河谷, 下伏基岩为玄武岩、杏仁状玄武岩、花岗闪长岩等, 岩性复杂。

折射层析成像法主要用于探测覆盖层厚度和坝址岩体的波速分层情况, 文章以下坝址Z11剖面来进行举例说明。图2为Z11测线示意图。

2.2外业工作方法

此次外业数据采集使用美国SI公司的S-Land全数字化地震勘探数据采集系统, 单个采集站控制36个接收道, 28Hz数字化检波器。外业数据采集使用夯锤作为震源, 每个记录进行多次叠加, 确保每一个检波道的初至清晰。为提高解释精度, 采用纵测线观测系统, 并用13个以上不同炮检距激发, 保证目的层有足够射线。

2.3 折射层析资料处理

基于上章节所述的原理与处理流程, 在取得有效折射波地震资料后, 对折射波地震记录进行了初至拾取, 速度模型建立时底层速度是根据声波测井实测数值取值, 测区玄武岩的新鲜完整岩块的波速取值5.8km/s, 变质砂岩及泥岩的新鲜完整岩块波速取值4.5km/s。表层速度取1.0km/s左右, 速度模型最大深度设置为100m~120m。迭代反演网格设置为道间隔距一半 (即2.5m×1.5m) , 初至提取误差设置为0.5ms, 最大反演次数设置为12次, 有效的将反演单炮拟合误差控制在5ms以下。图3为处理流程的示意图。

2.4折射层析成果分析

Z11剖面表层主要由沙土、碎石、卵石、砾石、块石组成, 下伏基岩为砂泥岩。长度555米, 两个排列共激发42炮, 最远偏移距离为102.5米。由折射层析反演速度剖面推测覆盖层、全风化层速度0.75km/s~1.9km/s, 厚度24.5m~34.6m;强风化速度1.9km/s~2.7km/s, 厚度9.5m~34.6m;弱风化上段速度2.7km/s~3.55km/s, 厚度24.7m~45.6m;深度78.1m以下为弱风化下段, 速度大于3.55km/s, 解释结果与钻孔资料较吻合。

3结束语

影响折射层析反演结果的因素主要有观测系统、数据质量、正反演算法、目标区域的形态及性质等方面, 合理的观测系统可以获得研究区域较密集的射线覆盖, 为最终的成像打下良好的基础。正反演算法直接关系到地震层析成像的速度和精度, 网格划分、正反演所采用的计算方法、初始速度模型的选择等都是其主要的影响因素。目标区域的形态及性质决定射线的聚散程度和分分布情况, 在参数选择时, 应尽可能根据近地表调查资料得到的某些些参数的精确值, 或者从已知地质背景知识得到的某些参数的取值范范围和模型空间的参数分布特点, 然后建立相应的速度模型来进行初初至的射线追踪反演。

摘要：目前主要利用反射波法和折射波法来确定覆盖层厚度, 折射波法初至容易识别, 资料解释较容易, 但是传统折射解释方法却有着许多限制, 对于速度横向不均匀、下覆地层起伏变化较大或速度渐变、存在透射地形地质情况, 传统折射解释方法一般都不能给出较好的结果。地震折射层析成像有效地克服了一般折射解释方法的这些缺点。