Research on TV-L1 Optical Flow Model for Image Registration Based on Fractional-order Differentiation
-
摘要: 图像的非刚性配准在计算机视觉和医学图像分析中有着重要的作用.TV-L1（全变分L1范数、Total variation-L1）光流模型是解决非刚性配准问题的有效方法，但TV-L1光流模型的正则项是一阶导数，会导致纹理特征等具有弱导数性质的信息模糊.针对该问题，将G-L（Grünwald-Letnikov）分数阶引入TV-L1光流模型，提出基于G-L分数阶微分的TV-L1光流模型，并应用原始-对偶算法求解该模型.新的模型用G-L分数阶微分代替正则项中的一阶导数，由于分数阶微分比整数阶微分具有更好的细节描述能力，并能有效地、非线性地保留具有弱导数性质的纹理特征，从而提高图像的配准精度.另外，通过实验给出了配准精度与G-L分数阶模板参数之间的关系，从而为模板最佳参数的选取提供了依据.尽管不同类型的图像其最佳参数是不同的，但是其最佳配准阶次一般在1 ~2之间.理论分析和实验结果均表明，提出的新模型能够有效地提高图像配准的精度，适合于包含较多弱纹理和弱边缘信息的医学图像配准，该模型是TV-L1光流模型的重要延伸和推广.Abstract: In computer vision and medical image analysis, non-rigid image registration is a challenging task. TV-L1 (Total variation-L1) optical flow model has been proved to be an effective method in the field of non-rigid image registration. It can solve the problem of fuzzy edge caused by smooth displacement fields of Horn-Schunck, but its first-order derivative in regularization term leads to fuzzy texture information with weak derivative property. Aiming at the problem, this paper introduces G-L (Grünwald-Letnikov) fractional differentiation to TV-L1 optical flow model, and proposes a new TV-L1 optical flow model based on fractional differentiation, and then finds the solution of the model using primal-dual algorithm. In this paper we use Grünwald-Letnikov fractional order differential instead of the first-order derivative in the regularization term for its better ability of detail description than first-order's. Then we purposefully control to retain or inhibit the texture information with weak derivative nature, thus improving the registration accuracy. Experimental results show that the proposed method has a better registration accuracy in registration of texture information with weak derivative, and that it can be considered an important extension and generalization of TV-L1 optical flow modes.
-
瞬变源是一种偶发的短暂的非周期性的天文现象.从观测上, 其持续时标从数秒到数周甚至数年.目前已知的主要瞬变源为:超新星1、伽玛暴2、微引力透镜3、恒星被大质量黑洞所瓦解的潮汐瓦解事件4以及引力波的电磁对应体5等.瞬变源对于研究宇宙的起源和极端环境下的物理现象有着重要意义.宇宙加速膨胀现象的发现[1-3]正是通过对大样本瞬变源超新星的观测研究发现的.超新星、伽玛暴等瞬变源爆发时辐射的能量超过整个星系的光度, 这种极高能量的短时大爆发为研究极端物理环境下的物理现象提供了难得的观测条件.
1https://en.wikipedia.org/wiki/Supernova
2https://en.wikipedia.org/wiki/Gamma-ray\_burst
3https://en.wikipedia.org/wiki/Gravitational\_microlensing
4https://en.wikipedia.org/wiki/Tidal\_disruption\_event
5https://en.wikipedia.org/wiki/Gravitational\_wave
由于瞬变源是偶发天文事件, 要求瞬变源搜索的观测设备具有大视场(即单位时间内能观测到更大的天区)和高时间采样率(即对同一天区的高回访观测频率)的特点.我国建设中的瞬变源搜索设备地基广角相机阵(Ground wide angle camera, GWAC), 由36台直径为18 cm的广角相机组成, 每个广角相机配有4 k×4 k的CCD探测器, 整个相机阵的视场达到5 000平方度.每15秒产生一幅观测图像(10秒曝光+5秒读出数据), 即15秒会对原视场作一次回访.该设备对于未知短时标瞬变源的搜索具有重要的意义, 同时也对传统的瞬变源搜索技术提出了巨大挑战.相对于国际主流的瞬变源搜索设备, 地基广角相机阵设备在视场和采样频率上都将提高1~2个数量级.
经典的瞬变源搜索过程如图 1所示.主要原理是通过将观测图像(图 1(a))与模板图像(图 1(b))进行相减, 如果是一个瞬变源(即新出现的源), 那么在减完后的残差图像(图 1(c))中就是一个类似完整点源的像(如图 1(c)中的$o1$和$o2$), 而其他残缺的像斑则为相减过程中产生的噪声(如图 1(c)中的$n1$, $n2$, $n3$).如何将残差图像中的瞬变源从周围的噪声中自动快速地识别出来是本文要解决的关键问题.传统的识别方法是人眼识别, 由天文学家对所有观测图像进行逐幅识别.这种方法虽然正确率高, 但是效率非常低, 对于现代的大数据瞬变源巡天的处理是无法适用的.
随着数据处理技术的发展, 不同国家的天文学家曾尝试开展利用机器学习的自动分类方法进行瞬变源自动分类的研究.最早的开创性工作是2007年Bailey等[4]将监督式机器学习分类技术实验性地应用于超新星工厂(The nearby supernova factory, SNFactory)巡天项目; 2012年, Brink等[5]在Bailey等[4]的算法基础上, 开发了基于高精随机森林架构下的分类器来识别帕洛马(Palomar transient factory, PTF)巡天项目的瞬变源; 2013年, Brink等[5]利用递归方法对Bloom等[6]使用的特征参量进行了优化, 提高了分类性能. 2014年, Buisson等[7]基于由主成分分析法得到的特征参量对随机森林、K-近邻、支持向量机、神经网络、贝叶斯等多种算法进行综合比较, 结果表明随机森林算法在瞬变源识别中具有最好性能.最新工作是2015年Goldstein等[8]在文献[5-6]的基础上进行特征参量的添加和优化, 最后选取38个特征参量来进行训练分类器, 基于随机森林算法, 使得分类的正确率和处理效率较前面文献中的工作都有了较大的提高[8].
本文以文献[8]的算法为基础, 结合地基广角相机阵的数据特点, 提出基于等光度测量星像轮廓等新的特征参量, 使用实际星像轮廓仿真和构建较真实的训练样本算法; 加入基于实测数据分析的噪声过滤判据等方法, 实现一个优化的瞬变源快速自动识别系统.
本文的组织结构如下:第1节阐述瞬变源识别系统的特征参量及其提取, 着重描述了基于等光度测量星像轮廓的特征参量; 第2节描述训练样本的构建过程; 第3节和第4节阐述本识别系统的实现过程和具体的测试与验证; 第5节对本文提出的识别系统进行讨论与总结.
1. 瞬变源特征参量的研究
瞬变源的自动识别就是将瞬变源从残差图像提取的点源样本中识别出来.残差图像中的识别样本主要分为真实点源(瞬变源候选体)和噪声源两类.其实例的效果如图 1(c)中$o1$, $o2$和$n1$, $n2$, $n3$所示.从理论上分析, 残差图像中的真实点源应具有类似于观测图像中点源的能量分布轮廓, 即星像的能量从中心到边缘应具有仪器自身特点的点扩散函数分布(一般为类似于高斯函数分布); 而残差图像中的噪声则主要来自于个别像元的随机噪声、位置及轮廓匹配中的差异引起的残差像斑、饱和星留下的残差等.虽然噪声源的类型多种多样, 但都不具有真实点源所具有的能量分布和形状等主要特点, 因此点源的能量分布及形状相关的特征参量是进行自动分类的基本特征参量.如何提取特征参量来表达星像轮廓的特征直接会影响到分类模型的最后识别结果和数据处理速度.
分析文献[8]中的38个特征参量, 保留其中与星像轮廓相关的以及具有高权重值的18个特征参量.去除如星系相关星等、CCD编号等与GWAC项目及其科学目标无关的参量.同时, 为了提高处理速度, 取消部分耗时过长和权重较低的特征参量.保留的18个特征参量根据是否与星像轮廓相关分成两组参量, 即第Ⅱ和第Ⅲ组(详见表 1).第Ⅱ组特征参量的计算需要对图像进行预处理.第Ⅲ组的参量主要是利用点扩散函数(PSF)拟合方法来确定星像的轮廓.本文新引入13个特征参量记为第Ⅰ组参量, 主要是基于等光度方法测定星像轮廓参量, 取代原算法中的第Ⅲ组参量来获得更好的性能优化.
表 1 特征参量Table 1 Feature sets组号 序号 特征参量 参量描述 权重 排序 来源 Ⅰ 1 flux_radius2 20%能量处的像斑孔径大小(单位:像素) 0.1391 1 本文新参量 2 flux_radius1 10%能量处的像斑孔径大小(单位:像素) 0.0548 6 3 flux_aper 固定孔径($r$ = 2.5像元)的流量 0.0287 13 4 ISO 0 等光度区域0的面积(单位:像素平方) 0.0559 5 5 ISO 1 等光度区域1的面积(单位:像素平方) 0.0308 11 6 ISO 2 等光度区域2的面积(单位:像素平方) 0.0145 20 7 ISO 3 等光度区域3的面积(单位:像素平方) 0.0145 19 8 ISO 4 等光度区域4的面积(单位:像素平方) 0.0099 23 9 r_max_aper 最大像元光度流量与固定孔径流量之比 0.1056 3 10 r_aper_ISO 固定孔径流量与等光度流量之比 0.0295 12 11 r_aper_ISOCOR 固定孔径流量与修正等光度流量之比 0.0213 16 12 mag_err_aper 星等的均方根误差 0.0349 10 13 class_star 恒星与星系分类标识(取值: 0~1) 0.0072 25 Ⅱ 14 diffsum 在矩阵$R(d)$上, 以对象为中心构成的5×5矩阵中所有元素的和 0.0512 7 文献[8] 15 colmeds 在矩阵$B(d)$上, 每列元素中位数的最大值 0.0152 18 16 numneg 在矩阵$R(d)$上, 以对象为中心构成的7×7矩阵中负元素的个数 0.0086 24 17 a_image 长轴方向上的均方根, 来自SExtractor 0.0138 21 18 b_image 短轴方向上的均方根, 来自SExtractor 0.1152 2 19 ellipticity 1-b_image/a_image, 来自SExtractor 0.0908 4 20 flags SExtractor在矩阵$I(d)$上的提取标志, 来自SExtractor 0.0404 8 21 mag_aper 固定孔径的星等, 来自于SExtractor 0.0361 9 22 n2sig3 在矩阵$R(d)$上, 以对象为中心构成的5×5矩阵中元素值＜-2的个数 0.0273 14 23 n3sig3 在矩阵$R(d)$上, 以对象为中心构成的5×5矩阵中元素值＜-3的个数 0.0228 15 24 n3sig5 在矩阵$R(d)$上, 以对象为中心构成的7×7矩阵中元素值＜-3的个数 0.0188 17 25 n2sig5 在矩阵$R(d)$上, 以对象为中心构成的7×7矩阵中元素值＜-2的个数 0.0135 22 Ⅲ r_aper_psf (flux_aper+flux_psf)/flux_psf 0.148 文献[8] flux_ratio 矩阵$I(d)$上以对象为中心的5个像素上的流量值与矩阵$I(t)$上以对象为中心的5个像素上流量值的绝对值之比 0.037 n3sig3shift 矩阵$R(d)$上以对象为中心构成的5×5矩阵中元素≥ 3的个数与矩阵$R(t)$上以对象为中心构成的5×5矩阵中元素大于等于3的个数之差 0.019 n3sig5shift 矩阵$R(d)$上以对象为中心构成的7×7矩阵中元素≥ 3的个数与矩阵$R(t)$上以对象为中心构成的7×7矩阵中元素大于等于3的个数之差 0.018 n2sig3shift 矩阵$R(d)$上以对象为中心构成的5×5矩阵中元素≥ 2的个数与矩阵$R(t)$上以对象为中心构成的5×5矩阵中元素大于等于2的个数之差 0.014 n2sig5shift 矩阵$R(d)$上以对象为中心构成的7×7矩阵中元素≥ 2的个数与矩阵$R(t)$上以对象为中心构成的7×7矩阵中元素大于等于2的个数之差 0.012 1.1 残差图像的预处理
我们的实测数据6对比测试实验结果表明:对于地基广角相机阵做图像的预处理[5, 8], 同样能提高系统的识别正确率, 但预处理过程的基本参数选取与文献[8]不同, 具体处理过程如下:
6迷你地基广角相机阵(Mini-GWAC)是GWAC的先导项目, 由12个7厘米望远望组成的阵.其观测策略、科学目标和数据特点与GWAC相同.
1) 以残差图像中找到的星像目标为中心, 截取出$(2k+1)$像素× $(2k+1)$像素的窗口图像.我们的实测对比试验结果表明:对于GWAC图像采用$k = 15$ (即31像素×31像素)的窗口像图最为有效.等效于星像轮廓参量半高全宽(FWHM)的大约20倍为经验的合理参量.其中, 每个$(x, y)$处的像元响应量标记为$I_{x, y}$.
2）计算表 1中第Ⅱ组参量中涉及的矩阵R和B, 计算方法即文献[8]中的式(1)、式(4)和式(5), 具体表述如下:
$ C_{x, y} = \frac{1}{N_{u}} \sum\limits_{i=0}^1 \sum\limits_{j=0}^1I_{2x+i, 2y+j} $
(1) 矩阵$C_{x, y}$即为原图像每相邻4个像元的压缩.压缩处理提高了图像的对比度, 能更容易地区分出真实像斑与噪声的轮廓.对矩阵$C_{x, y}$做流量的归一化处理后, 即可得到矩阵R.其计算表达式为
$ R_{x, y}\approx\frac{1}{1.4826}\left[\frac{C_{x, y}-{\rm med}(C)}{{\rm med}(| C_{x, y}-{\rm med}(C)|)}\right] $
(2) 其中, ${\rm med}(\cdot)$为中值计算符.
矩阵B由未经过压缩的原窗口图像(GWAC采用: 31像素×31像素)的直接处理得到, 表示像元响应与中值的偏移量的归一化值.计算公式为
$ B_{x, y}=\frac{I_{x, y}-{\rm med}(I)}{\max(| I |)} $
(3) 3) 表 1第Ⅱ组中的相关参量可根据“参量描述”列中的说明, 由前面的相关矩阵公式计算得出.
1.2 特征参量的优化提取
基于GWAC数据的特点分析, 利用点扩散函数(Point spread function, PSF)拟合来描述星像轮廓主要有以下不足: 1)拟合所需的计算比较耗时; 2)对于拟合轮廓所用的数学模型有依赖性. GWAC的PSF轮廓主要由光学轮廓而非大气视宁度决定, 相比由大气视宁度占主导的望远镜, 其PSF比较复杂难以用简单的高斯模型拟合.因此, 通过引入新的基于等光度测量星像轮廓的特征参量取代原有的PSF拟合参量进行优化(表 1中第Ⅰ组参量).
图 2显示了三种不同的星像轮廓测量方法[9].等光度轮廓(ISO)的测量主要将光度水平相同处连结成一条线构成星像的轮廓区域.第Ⅰ组参量中的ISO 0~ISO 4表示5组处于不同光度水平的轮廓面积.不同光度$I_{i}$的计算公式为
$ I_{i}=S\ast\left( \frac{I_{p}}{S}\right)^{\frac{i}{8}} $
(4) 其中, S表示背景涨落标准方差($\sigma$)的5倍, $I_{p}$为像斑中最大的像元响应值. ISO 0~ISO 4的光度值计算取i=(0~4)由式(4)计算可得. ISOCOR表示将等光度计算的轮廓等效到高斯模型下的圆形轮廓.
除了等光度轮廓的测量, 第Ⅰ组新参量中还包含其他的轮廓辅助测量参量.例如参量1~3分别表示占流量20%和10%处的孔径大小, 以及孔径为2.5像元内的总流量; 参量9~11表示由不同测光方法得到的光度流量比值; 参量12表示星等的测量误差; 参量13表示恒星与星系的分类标识, 取值为0~1之间, 是一个与星像椭率相关的量.
2. 训练样本的构建
天文瞬变源自动识别系统使用监督式机器学习方法.即在大量训练样本类别已知的情况下, 通过机器学习训练分类器.数学表述为: $o[(v_{1}, v_{2}, v_{3}, \cdots), class]$, o表示一个对象, $v_{i}$表示对象第i个特征参量, 所有特征参量一起构成对象的特征参量空间, 而$class$则表示对象o的类别.
天文瞬变源是相对稀有事件, 实际观测瞬变源数据难以提供足够数量的训练样本集, 尤其对于无历史数据积累的刚建成的观测设备.利用仿真的方法构建训练样本是唯一可行的途径, 但如何使仿真的样本具有较高的真实性, 是要解决问题的关键点.像斑的轮廓是能正确分类的关键因素.因此在仿真重构训练样本时, 采用真实的星像作为星像轮廓模板, 而仿真调整的参数仅仅是像斑的位置和像斑的响应流量.
主要通过两种方法仿真瞬变源: 1)从去除背景的图像中选出一批(约400颗/幅图)从亮到暗不等的星, 作为星像轮廓模板.选取星像轮廓模板时需要保证这些星像不受周围星的干扰, 相对比较孤立.将这些星像轮廓模板按随机位置(或有规则排列)叠加到原始观测图像中构建出含有瞬变源的仿真观测图像. 2)从去除背景的图像中选出一颗较为孤立的星, 以10倍半高全宽7的窗口从图像中裁剪出来, 作为星像轮廓模板.然后, 对该模板的流量进行仿真重构.仿真重构的流量从饱和星等开始一直到极限星等附近(最暗星仿真到$2.5\, \sigma$).最后按照随机方式(或者有规则排列)将这些仿真的星撒回到实际观测的图像中构建出含有瞬变源的仿真图像.
7半高全宽(FWHM):二维高斯函数拟合计算.
以上两种方法仿真瞬变源的整个过程如图 3所示, 得到的仿真观测图像与实测的模板图像相减得到残差图像, 然后通过特征参量的提取得到前述的$v_{1}, v_{2}, v_{3}, \cdots, v_{i}$, 而$class$可由仿真瞬变源注入时的位置信息, 利用位置搜索得到相应的瞬变源分类信息.最后加入肉眼识别去除因饱和星带来的干扰.考虑饱星和变星等干扰因素后, 我们的样本污染程度＜5%, 而根据以前的研究[6], 样本的受污染程度＜10%都是可靠的.
以上两种方法各有优点.方法1能够仿真构建出多种真实轮廓, 因为一幅图像在不同位置星像轮廓会有一些细小的差异.方法2能根据需要仿真出任何不同亮度的目标, 便于对探测极限附近星像的仿真和探测能力评估.两种方法相辅相成, 共同完成训练样本的有效构建.
3. 系统的实现
天文瞬变源候选体搜索的数据处理流程如图 4所示.主要包括如下过程: 1)图像相减过程, 即在完成观测图像与模板图的轮廓与流量匹配后进行相减获取残差图像; 2)对残差图像完成点源提取; 3)点源提取的星表和残差图像输入到自动识别系统(虚线框部分).自动识别系统主要完成图像预处理与特征参量提取、数值过滤器、自动分类器. 4)输出自动识别出的瞬变源候选体.
本识别系统的工程化实现主要基于pyth-on(2.7)和机器学习处理包python-sklearn以及其他相关的天文数据处理包pyfits等.除前述算法外, 还加入一个数值过滤器处理模块, 主要实现对亮星相减后残差噪声的定向去除.判据的物理意义是相减后的残差图像像斑的相应孔径内存在若干个光度小于(接近)零或者大大小于背景噪声水平的像元, 则认为是噪声而非真实瞬变源.判断标准主要通过对实测数据的处理分析总结得出, 数学描述如下:
$ \begin{align} \begin{cases} {\rm len}(Flux15_{x, y}==1e-30)>10\\[1mm] {\rm len}(Flux8_{x, y}==1e-30)>3\\[1mm] {\rm len}(Flux15_{x, y}<{\rm med}(Flux15_{x, y})-6\sigma)>5\\[1mm] {\rm len}(Flux8_{x, y}<{\rm med}(Flux8_{x, y})-4\sigma)>3 \end{cases} \end{align} $
(5) 其中, $Flux15_{x, y}$和$Flux8_{x, y}$分别表示以像斑中心为中心分别截取的$15\times15$和$8\times8$像元大小的窗口图像. len$(\cdot)$表示统计满足条件像元数目算符, med$(\cdot)$为中值计算算符, $\sigma$为窗口图像背景的标准方差.式(5)中的4个判据条件只要满足其一即被证伪, 从瞬变源候选体中排除出去.
Buisson等[7]的系统性对比测试表明针对瞬变源识别的数据和处理特点, 随机森林(Random forest, RF)[10-12]算法具有较优秀的表现.随机森林是利用多棵决策树对样本进行训练和预测的一种分类器.在预测某一个测试样本类别时, 由随机森林中的所有树共同投票决定, 样本的类别取决于投票数多少.该算法具有训练速度快、容易实现并行化, 能够快速处理高维数据、可以处理离散型变量(分类)和连续型变量(回归)、分类器训练完成后能够给出特征参量的重要性信息、预测时能够给出测试集中每个实例属于不同类别的概率等优点.因此, 采用随机森林算法实现样本的训练与分类.主要的设置参数取值如表 2所示.
表 2 随机森林主要参数Table 2 The main parameters of random forest超参数名称 取值 描述 n_estimators 100 随机森林中树的个数 criterion entropy 决定树中节点是否进行分割的决策函数 n_jobs -1 随机森林中并行训练树的个数, -1表示并行训练树的个数等于计算机CPU的核数 max_features 5 训练节点时无放回随机抽取的最大特征维数 min_samples_split 3 训练分割节点时需要的最少样本数 max_depth unlimited 随机森林中树的最大深度 4. 自动识别系统的测试与验证
为了测试与验证本系统的识别正确率及在处理速度上的表现能力, 主要通过两种途径: 1)在Mini-GWAC的实测数据中加入仿真瞬变源的方法. 2) Mini-GWAC在实际观测中的性能表现.
测试途径1.采用1 200幅与训练样本不同的分别来自不同观测夜的数据.在每一幅实测的观测图像中注入大约400颗不同亮度的仿真瞬变源, 仿真瞬变源通过实际星像轮廓模板的仿真方法得出.然后按照图 4的流程执行处理.为了测试数值过滤器的性能, 对过滤器执行过滤的源不做真实剔除而是仅做标识便于前后对比.对比测试主要分成3组:测试A组为本系统, 采用特征参量为Ⅰ+Ⅱ组; 测试B组采用特征参量Ⅱ+Ⅲ组; 测试C组仅采用特征参量第Ⅱ组.测试平台的硬件设备CPU为: Core i7 2 600 K, 内存15 GB.软件系统为Scientific Linux 6.0版本.
测试的结果与分析: 1)数据处理速度: A、B、C三组处理每幅图像的平均时间为9.7 s, 98 s, 8.4 s.三组数据的比较表明, A组(本系统)较B组(文献[8]的主要参量)提速近10倍.从总的分析来看, 主要的数据处理耗时来自特征参量的提取.差异的主要原因分析:等光度轮廓测量较PSF做拟合的测量方法处理速度更快.另外, B组部分参量需要对模板图像进行操作, 从过程来说更为复杂, 因而增加了数据处理时间. C组仅用一组特征参量因此耗时最短. 2)筛选样本的正确检出率与错误检出率:图 5表示A(Ⅰ+Ⅱ)、B(Ⅱ+Ⅲ)、C(Ⅱ)三组方法在对不同信噪比(点源测光提取时的信噪比)瞬变源的正确检出率.结果表明, 信噪比>20的所有源都能被三组方法100%正确识别.信噪比等于=14时A和B能保持一致的10%识别正确率.随着信噪降低, A方法相对B方法在低信噪比降低过程中, 体现出更高的识别正确率, 在信噪比为3.4时8仍有85%以上的正确识别率.错误检出率(不是瞬变源当成瞬变源的数目占总检出数目的比率): A方法为8.6%, B方法为3.6%, C方法为6.4%. A的错检率最高, 当加入数值滤波器后, A的错检率会降到和B同等水平, 接近3%.结果分析: PSF拟合方法对星像轮廓的描述对于低信噪比不敏感, 而等光度轮廓测量不依赖拟合而是直接测量对于低信噪比部份更为敏感, 同时对于是否圆形形状不敏感, 导致对部分亮星留下的残差的错误识别.当加数值过滤处理后, 能很好地消除这方面噪声引起的错误检出.
8在此信噪比下即使肉眼辨星也有些困难.
测试途径2.本系统应用于Mini-GWAC实际观测的实时处理测试.经过大约半年的应用测试, 结果表明, 能实时地完成瞬变源的自动快速识别, 通过与星表交叉法找瞬变源在线处理结果交叉验证, 对比测试交叉的正确率在$99%$以上.交叉测试结果表明, 本文的识别方法在暗弱目标及有背景亮星干扰的情况下, 具有更好的筛选能力.图 6为实际观测中, 本文的快速自动识别方法发现的一个真实的瞬变源耀星.
5. 讨论与总结
针对现代大数据瞬变源巡天要求快速自动搜索瞬变源的技术需求, 结合我国在建中的地基广角相机阵的数据特点, 抓住星像轮廓的光度分布是自动识别的关键特征参量, 通过优化研究并开发了天文瞬变源自动识别系统.
本识别系统主要通过引入新的13个包括等光度测量星像轮廓的特征参量取代国际主流算法(Goldstein等[8])中的PSF拟合方法测量星像轮廓的参量.同时, 去除与模板相关的特征参量, 降低了数据处理的复杂度, 提高了数据处理速度.与Goldstein等[8]的算法相比, 处理速度提高近10倍.而正确检出率具有相同水平, 尤其在低信噪比处, 等光度测光星像轮廓参量比拟合法测轮廓参量更为敏感.从理论上分析也支持了这一测试结果.
基于天文瞬变源是相对稀有事件, 难以获取足够数量的训练样本.利用仿真方法构建训练样本集, 即采用真实的星像作为星像轮廓模板, 而仿真调整的参数仅是像斑的位置和像斑的响应流量, 从而实现较真实的数据仿真.研究表明, 对于由饱和星(较亮星)相减留下的部分残差噪声被误识别成瞬变源, 导致错误检出率较文献[8]的算法高出5%$, 表明等光度测量轮廓法相对于PSF拟合法对于星像轮廓是否是圆形形状不敏感.我们通过引入数值过滤器专门对此类噪声进行滤除, 最后的错误检出率能控制到与文献[8]的相同水平.
本识别系统已成功应用于我国已建成的迷你地基广角相机阵的实际数据在线处理.通过与星表交叉法找瞬变源在线处理结果的交叉验证, 表明正确率在99%$以上.在正式的地基广角相机阵(GWAC)建成以后, 只需要根据GWAC数据特点重新构建分类训练器, 便可快速实现系统的移植.因而, 本系统对于其他类似的天文大视场, 要求快速实时处理的天文瞬变源识别的项目也具有应用与参考价值.
致谢: 感谢GWAC项目组天文观测与工程维护人员在本系统开发测试过程中提供实测数据及在硬件设备支撑上给予的帮助.
-
图 4 验图像(其中第一行为参考图像, 第二行为与参考图像对应的浮动图像)
Fig. 4 Experimental images (The first line are reference images, the second line are floating images corresponding to reference images)
图 5 Lena图像实验(第一行为浮动图像和配准后图像, 第二行为差值图像. (a)为浮动图像, (b) $\sim$ (e)为配准后的图像; (b) TV-L$^{1}$方法; (c)本文方法($\alpha=1.2, k=1$); (d)本文方法($\alpha=1.2, k=2$); (e)本文方法($\alpha=1.2, k=3$); (f) $\sim$ (j)分别为第一行图像与参考图像(图 4(a))的差值图像)
Fig. 5 Lena image (The first line is floating image and registered image, the second line is difference image. (a) Floating image, (b) $\sim$ (e) are registered images, (b) TV-L$^{1}$, (c) Our method ($\alpha=1.2, k=1$), (d) Our method ($\alpha=1.2, k=2$), (e) Our method ($\alpha=1.2, k=3$), (f) $\sim$ (j) are difference images)
图 7 Brain1图像实验(第一行为浮动图像和配准后图像, 第二行为差值图像. (a)为浮动图像, (b) $\sim$ (e)为配准后的图像; (b) TV-L$^{1}$方法; (c)本文方法($\alpha=1.3, k=1$); (d)本文方法($\alpha=1.3, k=2$); (e)本文方法($\alpha=1.3, k=3$); (f) $\sim$ (j)分别为第一行图像与参考图像(图 4(b))的差值图像)
Fig. 7 Brain1 image (The first line is floating image and registered image, the second line is difference image. (a) Floating image, (b) $\sim$ (e) are registered images, (b) TV-L$^{1}$, (c) Our method ($\alpha=1.3, k=1$), (d) Our method ($\alpha=1.3, k=2$), (e) Our method ($\alpha=1.3, k=3$), (f) $\sim$ (j) are difference images)
图 8 Brain2图像实验(第一行为浮动图像和配准后图像, 第二行为差分图像. (a)为浮动图像, (b) $\sim$ (e)为配准后的图像; (b) TV-L$^{1}$方法; (c)本文方法($\alpha=1.3$, $k=1$); (d)本文方法($\alpha=1.3$, $k=2$); (e)本文方法($\alpha=1.3$, $k=3$); (f) $\sim$ (j)分别为第一行图像与参考图像(图 4(c))的差值图像)
Fig. 8 Brain2 image. The first line is floating image and registered image, the second line is difference image ((a) floating image, (b) $\sim$ (e) are registered images, (b) TV-L$^{1}$, (c) Our method ($\alpha=1.3$, $k=1$), (d) Our method ($\alpha=1.3$, $k=2$), (e) Our method ($\alpha=1.3$, $k=3$), (f) $\sim$ (j) are difference images)
表 1 参考图像和配准图像的均方误差(MSE)
Table 1 Mean square error (MSE) of reference image and registered
输入图片 配准前 TV-L $^{1}$ 本文算法($k$ = 1) 本文算法($k$ = 2) 本文算法($k$ = 3) Lena ($\alpha$ = 1.2) 669.33 14.86 10.36 9.17 11.56 Brain1 ($\alpha$ = 1.3) 295.85 27.51 18.20 15.93 20.22 Brain2 ($\alpha$ = 1.3) 813.77 31.02 11.47 9.95 17.89 
下载: 导出CSV
表 2 峰值信噪比(PSNR)
Table 2 Peak signal to noise ratio (PSNR)
输入图片 配准前 TV-L $^{1}$ 本文算法($k$ = 1) 本文算法($k$ = 2) 本文算法($k$ = 3) Lena ($\alpha$ = 1.2) 19.32 35.55 38.18 38.88 37.50 Brain1 ($\alpha$ = 1.3) 22.78 33.73 35.34 36.17 34.89 Brain2 ($\alpha$ = 1.3) 19.03 31.21 37.73 38.15 35.60
下载: 导出CSV
-
[1] 蒲亦非, 王卫星.数字图像的分数阶微分掩模及其数值运算规则.自动化学报, 2007, 33(11):1128-1135 http://www.aas.net.cn/CN/abstract/abstract13448.shtmlPu Yi-Fei, Wang Wei-Xing. Fractional differential masks of digital image and their numerical implementation algorithms. Acta Automatica Sinica, 2007, 33(11):1128-1135 http://www.aas.net.cn/CN/abstract/abstract13448.shtml [2] 陈青, 刘金平, 唐朝晖, 李建奇, 吴敏.基于分数阶微分的图像边缘细节检测与提取.电子学报, 2013, 41(10):1873-1880 doi: 10.3969/j.issn.0372-2112.2013.10.001Chen Qing, Liu Jin-Ping, Tang Zhao-Hui, Li Jian-Qi, Wu Min. Detection and extraction of image edge curves and detailed features using fractional differentiation. Acta Electronica Sinica, 2013, 41(10):1873-1880 doi: 10.3969/j.issn.0372-2112.2013.10.001 [3] Pu Y F, Siarry P, Zhou J L, Liu Y G, Zhang N, Huang G, Liu Y Z. Fractional partial differential equation denoising models for texture image. Science China Information Sciences, 2014, 57(7):1-19 doi: 10.1007/s11432-014-5112-x [4] Liu J, Chen S C, Tan X Y. Fractional order singular value decomposition representation for face recognition. Pattern Recognition, 2008, 41(1):378-395 doi: 10.1016/j.patcog.2007.03.027 [5] Zhang Y, Pu Y F, Hu J R, Zhou J L. A class of fractional-order variational image inpainting models. Applied Mathematics and Information Sciences, 2012, 6(2):299-306 http://www.researchgate.net/publication/260210952_A_Class_of_Fractional-Order_Variational_Image_Inpainting_Models [6] Ren Z M. Adaptive active contour model driven by fractional order fitting energy. Signal Processing, 2015, 117:138-150 doi: 10.1016/j.sigpro.2015.05.009 [7] 薛鹏, 杨佩, 曹祝楼, 贾大宇, 董恩清.基于平衡系数的Active Demons非刚性配准算法.自动化学报, 2016, 42(9):1389-1400 http://www.aas.net.cn/CN/abstract/abstract18927.shtmlXue Peng, Yang Pei, Cao Zhu-Lou, Jia Da-Yu, Dong En-Qing. Active demons non-rigid registration algorithm based on balance coefficient. Acta Automatica Sinica, 2016, 42(9):1389-1400 http://www.aas.net.cn/CN/abstract/abstract18927.shtml [8] 张桂梅, 曹红洋, 储珺, 曾接贤.基于Nyström低阶近似和谱特征的图像非刚性配准.自动化学报, 2015, 41(2):429-438 http://www.aas.net.cn/CN/abstract/abstract18621.shtmlZhang Gui-Mei, Cao Hong-Yang, Chu Jun, Zeng Jie-Xian. Non-rigid image registration based on low-rank Nyström approximation and spectral feature. Acta Automatica Sinica, 2015, 41(2):429-438 http://www.aas.net.cn/CN/abstract/abstract18621.shtml [9] 闫德勤, 刘彩凤, 刘胜蓝, 刘德山.大形变微分同胚图像配准快速算法.自动化学报, 2015, 41(8):1461-1470 http://www.aas.net.cn/CN/abstract/abstract18720.shtmlYan De-Qin, Liu Cai-Feng, Liu Sheng-Lan, Liu De-Shan. A fast image registration algorithm for diffeomorphic image with large deformation. Acta Automatica Sinica, 2015, 41(8):1461-1470 http://www.aas.net.cn/CN/abstract/abstract18720.shtml [10] Thirion J P. Image matching as a diffusion process:an analogy with Maxwell's demons. Medical Image Analysis, 1998, 2(3):243-260 doi: 10.1016/S1361-8415(98)80022-4 [11] Wang H, Dong L, O'Daniel J, Mohan R, Garden A A, Ang K K, Kuban D A, Bonnen M, Chang J Y, Cheung R. Validation of an accelerated 'demons' algorithm for deformable image registration in radiation therapy. Physics in Medicine and Biology, 2005, 50(12):2887-2905 doi: 10.1088/0031-9155/50/12/011 [12] Palos G, Betrouni N, Coulanges M, Vermandel M, Devlaminck V, Rousseau J. Multimodal matching by maximisation of mutual information and optical flow technique. In:Proceedings of the 26th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (IEMBS). San Francisco, CA, USA:IEEE, 2004. 1679-1682 http://europepmc.org/abstract/MED/17272026 [13] Pock T, Urschler M, Zach C, Beichel R, Bischof H. A duality based algorithm for TV-L1-optical-flow image registration. Medical Image Computing and Computer-Assisted Intervention-MICCAI 2007. Berlin Heidelberg:Springer-Verlag, 2007. 511-518 [14] Pérez J S, Meinhardt-Llopis E, Facciolo G. TV-L1 optical flow estimation. Image Processing on Line, 2013, 3:137-150 doi: 10.5201/ipol [15] Yip S S F, Coroller T P, Sanford N N, Huynh E, Mamon H, Aerts H J W L, Berbeco R I. Use of registration-based contour propagation in texture analysis for esophageal cancer pathologic response prediction. Physics in Medicine and Biology, 2016, 61(2):906-922 doi: 10.1088/0031-9155/61/2/906 [16] Horn B K, Schunck B G. Determining optical flow. In:Technical Symposium East. Washington, D.C., USA:International Society for Optics and Photonics, 1981. 319-331 [17] Rudin L I, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D:Nonlinear Phenomena, 1992, 60(1-4):259-268 doi: 10.1016/0167-2789(92)90242-F [18] Chambolle A, Pock T. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 2011, 40(1):120-145 doi: 10.1007/s10851-010-0251-1 [19] Cafagna D. Fractional calculus:a mathematical tool from the past for present engineers. IEEE Industrial Electronics Magazine, 2007, 1(2):35-40 doi: 10.1109/MIE.2007.901479 期刊类型引用(1)
1. 祝杰，周丹，郑立新，曹建军，药新雨，陈国平，于涌，葛健，唐正宏，潘翔，杨臣威，姜鹏. 基于漂移扫描CCD技术的南极时域天文观测阵原型机的设计与实现. 中国科学:物理学 力学 天文学. 2024(08): 115-127 . 
百度学术其他类型引用(1)
-
下载:
计量
- 文章访问数: 2686
- HTML全文浏览量: 316
- PDF下载量: 871
- 被引次数: 2
