平码五不中公式规律
  • / 18
  • 下载费用:30 金币  

一种极低照度条件下光子计数激光三维计算成像方法.pdf

关 键 ?#21097;?/dt>
一种 照度 条件下 光子 计数 激光 三维 计算 成像 方法
  专利查询网所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
摘要
申请专利号:

CN201611088456.3

申请日:

2016.11.30

公开号:

CN106683168A

公开日:

2017.05.17

当前法律状态:

实审

有效性:

审中

法?#19978;?#24773;: 实质审查的生效IPC(主分类):G06T 17/00申请日:20161130|||公开
IPC分类号: G06T17/00 主分类号: G06T17/00
申请人: ?#26412;?#33322;天自动控制研究所
发明人: 李少军; 靳松直; 张聪; 高仕博; 蔡伟
地址: 100854 ?#26412;?#24066;海淀区永定路50号
优?#28909;ǎ?/td>
专利代理机构: ?#26412;?#21531;恒知识产权代理事务所(普通合伙) 11466 代理人: 张璐;?#30772;?#34892;
PDF完整版下载: PDF下载
法律状态
申请(专利)号:

CN201611088456.3

授权公告号:

|||

法律状态公告日:

2017.06.09|||2017.05.17

法律状态类型:

实质审查的生效|||公开

摘要

本发明公开了一种极低照度条件下光子计数激光三维计算成像方法,包括:同时测量得到所有像素的光子计数数据;构建光子数量与反射率的测量模型,并构建第一负对数似然函数;根据总变分半范数和第一负对数似然函数,构建第一目标函数;逐点迭代求解第一目标函数的凸优化问题,得到反射率图像;根据空间相关性构造二元假设检验统计量,筛选背景光子噪声;建立光子到达时间与测量距离的测量模型,并构建第二负对数似然函数;根据总变分半范数和第二负对数似然函数,构建第二目标函数;逐点迭代求解第二目标函数的凸优化问题,得到深度图像。通过使用本发明所提供的方法,可以在光子计数模式下同时获取高质量的激光三维距离像和强度像。

权利要求书

1.一种极低照度条件下光子计数激光三维计算成像方法,其特征在于,该方法包括如
下步骤:
同时测量得到所有像素的光子计数数据;
构建光子数量与反射率的测量模型,并构建给定光子数量条件下反射率的第一负对数
似然函数;
根据总变分半范数和第一负对数似然函数,构建第一目标函数Gα(x,y);
逐点迭代求解第一目标函数Gα(x,y)的凸优化问题,得到并输出反射率图像;
根据空间相关性构造二元假设检验统计量,筛选背景光子噪声;
建立光子到达时间与测量距离的测量模型,并构建光子达到时间条件下测量距离的第
二负对数似然函数;
根据总变分半范数和第二负对数似然函数,构建第二目标函数JR(x,y);
逐点迭代求解第二目标函数JR(x,y)的凸优化问题,得到并输出深度图像。
2.根据权利要求1所述的方法,其特征在于,所述光子计数数据包括:
每个像素接收到的光子数量和每个光子的到达时间。
3.根据权利要求2所述的方法,其特征在于,所述同时测量得到所有像素的光子计数数
据包括:
设定激光探测器阵列的驻留时间Tw;在所述驻留时间内,使用泛光激光脉冲持续照射某
一预设目标场景;
激光探测器阵列中的每个像素同时检测接收到的光子数量,并测量每个光子到达时
间,从而得到所有像素的光子计数数据。
4.根据权利要求1所述的方法,其特征在于,所述第一负对数似然函数为:
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>L</mi> <mo>&lsqb;</mo> <mi>&alpha;</mi> <mo>|</mo> <mi>n</mi> <mo>&rsqb;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <mi>log</mi> <mi> </mi> <msubsup> <mi>C</mi> <mi>n</mi> <mi>N</mi> </msubsup> <msub> <mi>P</mi> <mn>0</mn> </msub> <msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mrow> <mi>N</mi> <mo>-</mo> <mi>n</mi> </mrow> </msup> <msup> <mrow> <mo>&lsqb;</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>P</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&rsqb;</mo> </mrow> <mi>n</mi> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>=</mo> <mrow> <mo>(</mo> <mrow> <mi>N</mi> <mo>-</mo> <mi>n</mi> </mrow> <mo>)</mo> </mrow> <mi>a</mi> <mrow> <mo>(</mo> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mi>A</mi> <mo>-</mo> <mi>n</mi> <mi> </mi> <mi>log</mi> <mrow> <mo>{</mo> <mrow> <mn>1</mn> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mrow> <mo>&lsqb;</mo> <mrow> <mi>a</mi> <mrow> <mo>(</mo> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mi>A</mi> <mo>+</mo> <mi>B</mi> </mrow> <mo>&rsqb;</mo> </mrow> </mrow> </msup> </mrow> <mo>}</mo> </mrow> </mrow> </mtd> </mtr> </mtable> <mo>;</mo> </mrow>
其中,log(·)为自然对数,P0(x,y)为没有检测到光子的概率,N为激光脉冲照射的次
数,n为检测接收到的光子数量,α(x,y)为反射率,A表示平均光子数量,B表示背景光子平均
到达速率。
5.根据权利要求4所述的方法,其特征在于,所述逐点迭代求解第一目标函数Gα(x,y)的
凸优化问题,得到并输出反射率图像包括:
采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(α),迭代求解第一目标函数Gα
(x,y)的最优解的最小值,即迭代求解如下公式的最小值:
<mrow> <msub> <mi>&alpha;</mi> <mrow> <mi>M</mi> <mi>L</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>argmin</mi> <mrow> <mi>a</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&GreaterEqual;</mo> <mn>0</mn> </mrow> </munder> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>&beta;</mi> <mi>&alpha;</mi> </msub> <mo>)</mo> </mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>W</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>H</mi> </munderover> <mi>L</mi> <mo>&lsqb;</mo> <mi>&alpha;</mi> <mo>|</mo> <mi>n</mi> <mo>&rsqb;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&beta;</mi> <mi>&alpha;</mi> </msub> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>&alpha;</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中,αML(x,y)是第一目标函数Gα(x,y)的最优解,βa是惩罚系数,惩罚函数pl(α)是凸函
数;H、W分别表示图像的行和列;
惩罚函数pl(α)定义为反射率{α(x,y)}的总变分半范数:
<mrow> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>&alpha;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>W</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>H</mi> </munderover> <mo>|</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>x</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>|</mo> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>W</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>H</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mo>|</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>;</mo> </mrow>
在k次迭代,H(α)的二阶Taylor级数?#24179;麳k(z):
<mrow> <msup> <mi>H</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>;</mo> </mrow>
其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z与迭
代值zk的二范数的平方;
k+1次迭代公式如下:
<mrow> <msup> <mi>z</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <munder> <mrow> <mi>arg</mi> <mi>min</mi> </mrow> <mrow> <mi>z</mi> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <mi>W</mi> <mo>&times;</mo> <mi>H</mi> </mrow> </msup> </mrow> </munder> <mo>{</mo> <msup> <mi>H</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&beta;</mi> <mo>&CenterDot;</mo> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>}</mo> <mo>;</mo> </mrow>
其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;
凸目标函数Gα(x,y)的最优解简化为:
<mrow> <msup> <mi>z</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <munder> <mi>argmin</mi> <mrow> <mi>z</mi> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <mi>W</mi> <mi>H</mi> </mrow> </msup> </mrow> </munder> <mo>{</mo> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>s</mi> <mi>k</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <mi>&beta;</mi> <mo>&CenterDot;</mo> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mo>}</mo> <mo>;</mo> </mrow>
<mrow> <msup> <mi>s</mi> <mi>k</mi> </msup> <mo>=</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>-</mo> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
然后,通过如下的步骤得到反射率图像:
步骤A1,进行初始化,包括:
步骤A1a,设每个像素反射率初始值为
步骤A1b,对第一负对数似然函数求导,计算梯度
<mrow> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msub> <mi>&alpha;</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mrow> <mo>(</mo> <mi>N</mi> <mo>-</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>A</mi> <mo>+</mo> <mi>n</mi> <mi>A</mi> <mfrac> <mrow> <mi>exp</mi> <mo>&lsqb;</mo> <mo>-</mo> <msub> <mi>a</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mi>A</mi> <mo>-</mo> <mi>B</mi> <mo>&rsqb;</mo> </mrow> <mrow> <mn>1</mn> <mo>-</mo> <mi>exp</mi> <mo>&lsqb;</mo> <mo>-</mo> <msub> <mi>a</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&CenterDot;</mo> <mi>A</mi> <mo>-</mo> <mi>B</mi> <mo>&rsqb;</mo> </mrow> </mfrac> <mo>;</mo> </mrow>
步骤A1c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;
步骤A2,开始进行迭代,k=0,……,M;
步骤A3,迭代求解简化后的凸目标函数Gα(x,y)的最优解,得到αk+1;
步骤A4,判断αk+1是否满足迭代下降准则:
<mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>&le;</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mn>1</mn> <mo>&le;</mo> <mi>i</mi> <mo>&le;</mo> <mi>k</mi> </mrow> </munder> <mo>{</mo> <mi>&Phi;</mi> <mrow> <mo>(</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&epsiv;&sigma;</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>&alpha;</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>}</mo> <mo>;</mo> </mrow>
如果否,则αk+1←ηαk+1,然后再返回执行步骤A3;
如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:
<mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>&alpha;</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>/</mo> <mo>|</mo> <mo>|</mo> <msub> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>&le;</mo> <mi>&eta;</mi> <mo>;</mo> </mrow>
其中,η>0是很小的常数。
6.根据权利要求5所述的方法,其特征在于,所述根据空间相关性构造二元假设检验统
计量包括如下步骤:
对于每一个像素位置(x,y),根据八邻域像素检测到的光子到达时间的中值的平均值
构造二元假设检验统计量H(x,y);
根据二元假设检验统计量H(x,y)设置筛选光子准则。
7.根据权利要求6所述的方法,其特征在于,所述二元假设检验统计量H(x,y)通过如下
的公式计算得到:
<mrow> <mi>H</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>tm</mi> <mn>1</mn> </msub> <mo>+</mo> <msub> <mi>tm</mi> <mn>2</mn> </msub> <mo>+</mo> <msub> <mi>tm</mi> <mn>3</mn> </msub> <mo>+</mo> <msub> <mi>tm</mi> <mn>4</mn> </msub> <mo>+</mo> <mo>+</mo> <msub> <mi>tm</mi> <mn>5</mn> </msub> <mo>+</mo> <mo>+</mo> <msub> <mi>tm</mi> <mn>6</mn> </msub> <mo>+</mo> <mo>+</mo> <msub> <mi>tm</mi> <mn>7</mn> </msub> <mo>+</mo> <mo>+</mo> <msub> <mi>tm</mi> <mn>8</mn> </msub> </mrow> <mn>8</mn> </mfrac> <mo>;</mo> </mrow>
其中,八邻域像素检测到的光子到达时间分别为:
{tj(x-1,y-1),j=1,…,k1}、{tj(x-1,y),j=1,…,k2}、{tj(x-1,y+1),j=1,…,k3}、{tj
(x,y-1),j=1,…,k4}、{tj(x,y+1),j=1,…,k5}、{tj(x+1,y-1),j=1,…,k6}、{tj(x+1,y),
j=1,…,k7}、{tj(x+1,y+1),j=1,…,k8};
对应到达时间中值分别为tm1、tm2、tm3、tm4、tm5、tm6、tm7、tm8;
k1,k2,k3,k4,k5,k6,k7,k8分别表?#26223;?#37051;域像素接收到的光子数量。
8.根据权利要求7所述的方法,其特征在于,所述筛选光子准则为:
<mrow> <mo>|</mo> <msub> <mi>t</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>H</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>|</mo> <mo>&lt;</mo> <mn>2</mn> <msub> <mi>T</mi> <mi>p</mi> </msub> <mfrac> <mi>B</mi> <mrow> <msub> <mi>&alpha;</mi> <mrow> <mi>M</mi> <mi>L</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mi>S</mi> <mo>+</mo> <mi>B</mi> </mrow> </mfrac> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中,tj(x,y)为光子的到达时间,H(x,y)为二元假设检验统计量,αML(x,y)是反射率的
最大似然估计,Tp是激光脉冲宽度的均方根。
9.根据权利要求8所述的方法,其特征在于,所述第二负对数似然函数为:
<mrow> <mi>L</mi> <mo>&lsqb;</mo> <mi>R</mi> <mo>|</mo> <msub> <mi>t</mi> <mi>j</mi> </msub> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> <mo>&rsqb;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </munderover> <mi>l</mi> <mi>o</mi> <mi>g</mi> <mo>&lsqb;</mo> <mi>g</mi> <mo>&lsqb;</mo> <msub> <mi>t</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mn>2</mn> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>c</mi> <mo>&rsqb;</mo> <mo>&rsqb;</mo> <mo>;</mo> </mrow>
其中,L[R(x,y)|tj(x,y),j=1,…,n(x,y)]为第二负对数似然函数,当n(x,y)=0时,L
[R(x,y)|tj(x,y),j=1,…,n(x,y)]=0;g(·)为激光脉冲的波形。
10.根据权利要求9所述的方法,其特征在于,所述逐点迭代求解第二目标函数JR(x,y)
的凸优化问题,得到并输出深度图像包括如下步骤:
采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(R),迭代求解第二目标函数JR
(x,y)的最优解的最小值,即迭代求解如下公式的最小值:
<mrow> <msub> <mi>R</mi> <mrow> <mi>M</mi> <mi>L</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mi>argmin</mi> <mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>&Element;</mo> <mo>&lsqb;</mo> <mn>0</mn> <mo>,</mo> <mi>c</mi> <mi>T</mi> <mi>r</mi> <mo>/</mo> <mn>2</mn> <mo>)</mo> </mrow> </munder> <mrow> <mo>(</mo> <mn>1</mn> <mo>-</mo> <msub> <mi>&beta;</mi> <mi>R</mi> </msub> <mo>)</mo> </mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>W</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>H</mi> </munderover> <mi>L</mi> <mo>&lsqb;</mo> <mi>R</mi> <mo>|</mo> <msub> <mi>t</mi> <mi>j</mi> </msub> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> <mo>&rsqb;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&beta;</mi> <mi>R</mi> </msub> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>R</mi> <mo>)</mo> </mrow> <mo>;</mo> </mrow>
其中,RML(x,y)是第二目标函数JR(x,y)的最优解,βR是惩罚系数,惩罚函数pl(R)是凸函
数,H、W分别表示图像的行和列;
惩罚函数pl(R)定义为深度图像{R(x,y)}的总变分半范数:
<mrow> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>R</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>W</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>y</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>H</mi> </munderover> <mo>|</mo> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>|</mo> <mo>+</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>x</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>W</mi> </munderover> <munderover> <mo>&Sigma;</mo> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>H</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <mo>|</mo> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>R</mi> <mrow> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <mo>;</mo> </mrow>
在k次迭代,H(R)的二阶Taylor级数?#24179;麳k(z):
<mrow> <msup> <mi>H</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>;</mo> </mrow>
其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z与迭
代值zk的二范数的平方;
k+1次迭代公式如下:
<mrow> <msup> <mi>z</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <munder> <mrow> <mi>arg</mi> <mi>min</mi> </mrow> <mrow> <mi>z</mi> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <mi>W</mi> <mo>&times;</mo> <mi>H</mi> </mrow> </msup> </mrow> </munder> <mo>{</mo> <msup> <mi>H</mi> <mi>k</mi> </msup> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>&beta;</mi> <mo>&CenterDot;</mo> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>}</mo> <mo>;</mo> </mrow>
其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;
凸目标函数JR(x,y)的最优解简化为:
<mrow> <msup> <mi>z</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <munder> <mi>argmin</mi> <mrow> <mi>z</mi> <mo>&Element;</mo> <msup> <mi>R</mi> <mrow> <mi>W</mi> <mi>H</mi> </mrow> </msup> </mrow> </munder> <mo>{</mo> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>z</mi> <mo>-</mo> <msup> <mi>s</mi> <mi>k</mi> </msup> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>/</mo> <mn>2</mn> <mo>+</mo> <mi>&beta;</mi> <mo>&CenterDot;</mo> <mi>p</mi> <mi>l</mi> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mo>}</mo> <mo>;</mo> </mrow>
<mrow> <msup> <mi>s</mi> <mi>k</mi> </msup> <mo>=</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>-</mo> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msup> <mi>z</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>/</mo> <msub> <mi>&sigma;</mi> <mi>k</mi> </msub> <mo>;</mo> </mrow>
然后,通过如下所述的步骤得到深度图像:
步骤B1,进行初始化,包括:
步骤B1a,经过筛选后的像素,初始深度值为其中,d为距离分辨
率,δ为测量偏差;未经过筛选后的像素
步骤B1b,对第二负对数似然函数求导,计算梯度
<mrow> <mo>&dtri;</mo> <mi>H</mi> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>-</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>n</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> </mrow> </munderover> <mi>d</mi> <mi>v</mi> <mo>&lsqb;</mo> <msub> <mi>t</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>-</mo> <mn>2</mn> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>/</mo> <mi>c</mi> <mo>&rsqb;</mo> <mo>/</mo> <msub> <mi>dt</mi> <mi>j</mi> </msub> <mo>;</mo> </mrow>
步骤B1c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;
步骤B2,开始进行迭代,k=0,……,M;
步骤B3,迭代求解简化后的凸目标函数JR(x,y)的最优解,得到Rk+1;
步骤B4,判断Rk+1是否满足迭代下降准则:
<mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>&le;</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mn>1</mn> <mo>&le;</mo> <mi>i</mi> <mo>&le;</mo> <mi>k</mi> </mrow> </munder> <mo>{</mo> <mi>&Phi;</mi> <mrow> <mo>(</mo> <msub> <mi>R</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&epsiv;&sigma;</mi> <mi>k</mi> </msub> </mrow> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <msub> <mi>R</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>R</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>}</mo> <mo>;</mo> </mrow>
如果否,则Rk+1←ηRk+1,然后再返回执行步骤B3;
如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:
<mrow> <mo>|</mo> <mo>|</mo> <msub> <mi>R</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>-</mo> <msub> <mi>R</mi> <mi>k</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>/</mo> <mo>|</mo> <mo>|</mo> <msub> <mi>R</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>&le;</mo> <mi>&eta;</mi> <mo>;</mo> </mrow>
其中,η>0是很小的常数。

说明书

一种极低照度条件下光子计数激光三维计算成像方法

技术领域

本发明涉及激光三维计算成像技术,特别涉及一种极低照度条件下光子计数激光
三维计算成像方法。

背景技术

现有技术中的激光三维成像探测技术既可实时获取目标的深度图像和激光反射
率图像,从而同时生成目标场景的三维图像,又可实时自主测量相对姿态与运动参数,为视
觉操控提供实时的三维测量与感知信息,是一种新型的自主三维视觉测量技术。现有技术
中的激光三维成像遥感技术,可以直接快速生成目标场景的高分辨率三维影像或地形图,
不需要地面控制点的配合,是一种新型定量遥感技术。

然而,现有技术中的激光三维成像探测和激光三维成像遥感技术等存在的共性难
点问题在于:

1)远距离探测导致激光回波能量很低,激光探测器阵列接收到的光辐射只有光子
量级照度;

2)需要在极低照度条件下同时生成深度图像和激光反射率图像,或三维立体测绘
和光谱图像;

3)非扫描激光凝视成像,需要满足动态条件下可靠性、?#39318;?#24615;要求。

现有技术中的激光显微断层成像仪存在的主要难点问题是:需要在极低照度条件
下同时获取高质量的立体图像和激光反射率图像(反射率信息),满足对活体组织结构动
态、定量、三维的显微观测需求,也要求非扫描激光凝视成像。

目前,现有技术中已有的激光成像探测体制主要可以分为以下几类:

1)基于盖革模式雪崩光电二极管(Avalanche Photodiode in Geiger Mode,GM-
APD)激光探测器阵列的激光凝视成像探测体制。该体制是解决此类激光成像探测的有效方
法,但是,通过时间相关方法检测激光回波的飞行时间,只能产生三维距离像,而不能产生
强度像;而?#36965;?#24403;背景光不连续时,即使有光子计数过程中内在的泊松噪声,每个像素检测
到几十个光子,也足以精确的距离成像,当由背景光产生的光子检测导致?#26412;?#31163;时,距离成
像需要与反射成像相当的光子数量,需要几百个光子才能产生精确的强度像。

2)基于像增强CCD的激光强度成像体制仅接收激光回波的强度,不测量激光脉冲
的飞行时间,只能产生强度像;在没有背景光照射的条件下,精确的反射成像需要几百个光
子才能产生精确的强度像,当存在背景光照射?#26412;?#30830;的反射成像需要的激光回波强度更
大。

3)基于线性雪崩光电二极管(APD)激光探测器阵列激光凝视成像探测体制,既可
以测量激光回波的飞行时间,还可测量激光回波的强度。但是,该体制强度成像需要连续的
光子流,不能工作在光子计数模式,灵敏度很低,激光回波利用效率不高,不能满足极低照
度条件下成像。

4)激光扫描成像体制,通过光机扫描或微扫描机构实现激光束逐点成像,从而获
得距离像、强度像。但是,该技术要求激光重复频率很高,存在运动光机结构,平台适应性
差,应用范围受限,因此不适用激光三维成像遥感,激光三维成像探测。

由此可知,现有技术中所要解决的技术问题是:在极低照度条件下、在很短驻留时
间内,每个GM-APD像素只能检测到少量、不连续的光子,这些光子远不足以清晰成像;在没
有背景光照射的条件下,清晰反射率成像需要几百个连续光子流,清晰距离成像则需要几
十个光子流,因此难以同时生成激光三维距离像和反射率图像。

发明内容

有鉴于此,本发明提供一种极低照度条件下光子计数激光三维计算成像方法,从
而可以在光子计数模式下同时获取高质量的激光三维距离像和强度像。

本发明的技术方案具体是这样实现的:

一种极低照度条件下光子计数激光三维计算成像方法,该方法包括如下步骤:

同时测量得到所有像素的光子计数数据;

构建光子数量与反射率的测量模型,并构建给定光子数量条件下反射率的第一负
对数似然函数;

根据总变分半范数和第一负对数似然函数,构建第一目标函数Gα(x,y);

逐点迭代求解第一目标函数Gα(x,y)的凸优化问题,得到并输出反射率图像;

根据空间相关性构造二元假设检验统计量,筛选背景光子噪声;

建立光子到达时间与测量距离的测量模型,并构建光子达到时间条件下测量距离
的第二负对数似然函数;

根据总变分半范数和第二负对数似然函数,构建第二目标函数JR(x,y);

逐点迭代求解第二目标函数JR(x,y)的凸优化问题,得到并输出深度图像。

较佳的,所述光子计数数据包括:

每个像素接收到的光子数量和每个光子的到达时间。

较佳的,所述同时测量得到所有像素的光子计数数据包括:

设定激光探测器阵列的驻留时间Tw;在所述驻留时间内,使用泛光激光脉冲持续
照射某一预设目标场景;

激光探测器阵列中的每个像素同时检测接收到的光子数量,并测量每个光子到达
时间,从而得到所有像素的光子计数数据。

较佳的,所述第一负对数似然函数为:


其中,log(·)为自然对数,P0(x,y)为没有检测到光子的概率,N为激光脉冲照射
的次数,n为检测接收到的光子数量,α(x,y)为反射率,A表示平均光子数量,B表示背景光子
平均到达速率。

较佳的,所述逐点迭代求解第一目标函数Gα(x,y)的凸优化问题,得到并输出反射
率图像包括:

采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(α),迭代求解第一目标函
数Gα(x,y)的最优解的最小值,即迭代求解如下公式的最小值:


其中,αML(x,y)是第一目标函数Gα(x,y)的最优解,βα是惩罚系数,惩罚函数pl(α)
是凸函数;H、W分别表示图像的行和列;

惩罚函数pl(α)定义为反射率{α(x,y)}的总变分半范数:


在k次迭代,H(α)的二阶Taylor级数?#24179;麳k(z):


其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z
与迭代值zk的二范数的平方;

k+1次迭代公式如下:


其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;

凸目标函数Gα(x,y)的最优解简化为:



然后,通过如下的步骤得到反射率图像:

步骤A1,进行初始化,包括:

步骤A1a,设每个像素反射率初始值为

步骤A1b,对第一负对数似然函数求导,计算梯度


步骤A1c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;

步骤A2,开始进行迭代,k=0,······,M;

步骤A3,迭代求解简化后的凸目标函数Gα(x,y)的最优解,得到αk+1;

步骤A4,判断αk+1是否满足迭代下降准则:


如果否,则αk+1←ηαk+1,然后再返回执行步骤A3;

如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:


其中,η>0是很小的常数。

较佳的,所述根据空间相关性构造二元假设检验统计量包括如下步骤:

对于每一个像素位置(x,y),根据八邻域像素检测到的光子到达时间的中值的平
均值构造二元假设检验统计量H(x,y);

根据二元假设检验统计量H(x,y)设置筛选光子准则。

较佳的,所述二元假设检验统计量H(x,y)通过如下的公式计算得到:


其中,八邻域像素检测到的光子到达时间分别为:

{tj(x-1,y-1),j=1,…,k1}、{tj(x-1,y),j=1,…,k2}、{tj(x-1,y+1),j=1,…,
k3}、

{tj(x,y-1),j=1,…,k4}、{tj(x,y+1),j=1,…,k5}、{tj(x+1,y-1),j=1,…,k6}、

{tj(x+1,y),j=1,…,k7}、{tj(x+1,y+1),j=1,…,k8};

对应到达时间中值分别为tm1、tm2、tm3、tm4、tm5、tm6、tm7、tm8;

k1,k2,k3,k4,k5,k6,k7,k8分别表?#26223;?#37051;域像素接收到的光子数量。

较佳的,所述筛选光子准则为:


其中,tj(x,y)为光子的到达时间,H(x,y)为二元假设检验统计量,αML(x,y)是反射
率的最大似然估计,Tp是激光脉冲宽度的均方根。

较佳的,所述第二负对数似然函数为:


其中,L[R(x,y)|tj(x,y),j=1,…,n(x,y)]为第二负对数似然函数,当n(x,y)=0
时,L[R(x,y)|tj(x,y),j=1,…,n(x,y)]=0;g(·)为激光脉冲的波形。

较佳的,所述逐点迭代求解第二目标函数JR(x,y)的凸优化问题,得到并输出深度
图像包括如下步骤:

采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(R),迭代求解第二目标函
数JR(x,y)的最优解的最小值,即迭代求解如下公式的最小值:


其中,RML(x,y)是第二目标函数JR(x,y)的最优解,βR是惩罚系数,惩罚函数pl(R)
是凸函数,H、W分别表示图像的行和列;

惩罚函数pl(R)定义为深度图像{R(x,y)}的总变分半范数:


在k次迭代,H(R)的二阶Taylor级数?#24179;麳k(z):


其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z
与迭代值zk的二范数的平方;

k+1次迭代公式如下:


其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;

凸目标函数JR(x,y)的最优解简化为:



然后,通过如下所述的步骤得到深度图像:

步骤B1,进行初始化,包括:

步骤B1a,经过筛选后的像素,初始深度值为其中,d为距离
分辨率,δ为测量偏差;未经过筛选后的像素

步骤B1b,对第二负对数似然函数求导,计算梯度


步骤B1c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;

步骤B2,开始进行迭代,k=0,······,M;

步骤B3,迭代求解简化后的凸目标函数JR(x,y)的最优解,得到Rk+1;

步骤B4,判断Rk+1是否满足迭代下降准则:


如果否,则Rk+1←ηRk+1,然后再返回执行步骤B3;

如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:


其中,η>0是很小的常数。

如上可见,本发明所提供的极低照度条件下光子计数激光三维计算成像方法,可
以精确测量在很短的驻留时间内每个GM-APD像素检测到的光子数量以及每个光子到达时
间,对测量光子计数数据进行筛选,剔除噪声光子数据,根据经过筛选的光子数量和到达时
间数据分别重构目标场景的反射率图像和三维距离像,实现了GM-APD探测器阵列同时生成
三维距离像和反射率图像,从而可以在光子计数模式下同时获取高质量激光三维距离像和
强度像。

附图说明

图1为本发明实施例中的极低照度条件下光子计数激光三维计算成像方法的流程
示意图。

具体实施方式

为使本发明的目的、技术方案及优点更加清楚明白,以下参照附图并举实施例,对
本发明进一步详细说明。

本实施例提供了一种极低照度条件下光子计数激光三维计算成像方法。

图1为本发明实施例中的极低照度条件下光子计数激光三维计算成像方法的流程
示意图。如图1所示,本发明实施例中的极低照度条件下光子计数激光三维计算成像方法主
要包括如下所述的步骤:

步骤11,同时测量得到所有像素的光子计数数据。

在本步骤中,将同时测量所有像素的光子计数数据,从而获取驻留时间内的光子
数据。

具体来说,较佳的,在本发明的具体实施例中,所述步骤11可以包括如下所述的步
骤:

步骤111,设定激光探测器阵列的驻留时间Tw;在所述驻留时间内,使用泛光激光
脉冲持续照射某一预设目标场景;

另外,较佳的,在本发明的具体实施例中,所述激光探测器阵列可以是:GM-APD激
光探测器阵列。

步骤112,激光探测器阵列中的每个像素同时检测接收到的光子数量,并测量每个
光子到达时间,从而得到所有像素的光子计数数据。

例如,较佳的,在本发明的具体实施例中,所述像素可以是:GM-APD激光探测器阵
列中的GM-APD像素。

另外,在本发明的较佳实施例中,所述光子计数数据包括?#22909;?#20010;像素接收到的光子
数量和每个光子的到达时间,因此,在本发明的具体实施例中,可以设每个像素同时检测接
收到的光子数量为:{n(x,y)},每个光子的到达时间为:{tj(x,y),j=1,…,n(x,y)};其中,
tj(x,y)是相对于最近一个激光脉冲的时间。

在上述的步骤111、112中,由于为激光探测器阵列设定了驻留时间Tw,相当于为激
光探测器阵列中的每个像素设定了相同的驻留时间,从而可以保证面阵探测器中的所有像
素可以同步工作,适用于面阵探测器的并行成像。

当然,并不是每个像素都能产生有效的光子计数数据,有的像素可能没有检测到
任何光子,对于这类数据缺失,可以根据场景目标的空间相关性和反射率相关性来推断。

步骤12,构建光子数量与反射率的测量模型,并构建给定光子数量条件下反射率
的第一负对数似然函数。

在本步骤中,将先构建光子数量n(x,y)与反射率α(x,y)的测量模型,从而建立了
光子数量n(x,y)与反射率α(x,y)的数量关系;然后,即可根据测量模型构建给定光子数量
{n(x,y)}条件下反射率α(x,y)的负对数似然函数L[α(x,y)|n(x,y)],可称之为第一负对数
似然函数。

在相同的驻留时间内,相同数量的激光脉冲照射目标时,反射率越高的目标返回
的激光回波能量越大,单光子探测器检测到的光子数量越多;反射率越低的目标返回的激
光回波能量越小,单光子探测器检测到的光子数量越小。因此,目标像素点的激光反射率α
(x,y)可由检测到的光子数量n(x,y)进行编码,即每个像素的激光反射率α(x,y)是该像素
检测到的光子数量n(x,y)的函数。

根据单光子激光探测器的统计特性可知,一次激光脉冲照射,光子计数探测器最
多产生一次计数脉冲,其发生概率与背景光、信号光和反射率大小有关,且分布服从泊松分
布模型:

没有检测到光子的概率:P0(x,y)=exp{-[a(x,y)·A+B]} (1)

检测到光子的概率:P1(x,y)=1-P0(x,y) (2)

其中,A表示平均光子数量,B表示背景光子平均到达速率,α(x,y)为反射率。

对于N次激光脉冲照射,每个GM-APD像素检测到的光子数量n(x,y)服从二项分布,
该分布是反射率α(x,y)的函数:


因此,在本发明的技术方案中,可以构建光子数量n(x,y)与反射率α(x,y)的测量
模型。然后,构建给定光子数量{n(x,y)}条件下反射率α(x,y)的第一负对数似然函数L[α
(x,y)|n(x,y)],将反射率图像{α(x,y)}估计转换为逐点求解负对数似然函数与惩罚函数
构成的目标函数Gα(x,y)的最优解。

较佳的,在本发明的具体实施例中,所述第一负对数似然函数可以表示为:


其中,log(·)为自然对数。去掉与α(x,y)无关的常数,L[α(x,y)|n(x,y)]是反射
率α(x,y)的凸函数。

步骤13,根据总变分半范数和第一负对数似然函数,构建第一目标函数Gα(x,y)。

在本步骤中,可以由总变分半范数和所有像素位置的激光反射率{α(x,y)}的第一
负对数似然函数L[α(x,y)|n(x,y)]构建第一目标函数Gα(x,y)。

单光子探测器阵列获取的光子测量数据均为非负数,第一负对数似然函数L[α(x,
y)|n(x,y)]是光子数量n(x,y)的凸函数。因此,可以选择总变分半范数作为惩罚函数。第一
目标函数Gα(x,y)也是光子数量n(x,y)的凸函数,则估计反射率图像即是逐点求解凸函数Gα
(x,y)的最优解。

但是,测量数据的非负性在目标函数最优化中引入了一组不等式约束,导致最优
化求解变得困?#36873;?#22240;此,较佳的,在本发明的具体实施例中,针对非负约束的目标函数的凸
优化问题,可以使用一种二阶Taylor级数?#24179;?#20108;次序列目标函数的迭代数值求解算法。

较佳的,在本发明的具体实施例中,可以考虑反射率估计?#20132;?#32422;束,在第一目标函
数中选择总变分半范数作为惩罚函数,从所有像素检测到的光子数量{n(x,y)}估计反射率
图像{α(x,y)},即是逐点迭代求解负对数似然函数与惩罚函数构成的目标函数Gα(x,y)的
凸优化问题的最优解。

例如,较佳的,在本发明的具体实施例中,可以通过如下所述的公式得到第一目标
函数Gα(x,y)的最优解:


其中,αML(x,y)是最优解,βα是惩罚系数,惩罚函数pl(α)是凸函数。约束反射率估
计的非?#20132;?#24615;,H、W分别表示图像的行和列。

另外,惩罚函数pl(α)可以定义为反射率{α(x,y)}的总变分半范数:


步骤14,逐点迭代求解第一目标函数Gα(x,y)的凸优化问题,得到并输出反射率图
像。

测量数据的非负性在目标函数最优化中引入了一组不等式约束,导致最优化求解
变得困?#36873;?#22240;此,较佳的,在本发明的具体实施例中,针对非负约束的目标函数的凸优化问
题,可以使用一种二阶Taylor级数?#24179;?#20108;次序列目标函数的迭代数值求解算法。

例如,较佳的,在本发明的具体实施例中,所述步骤14可以通过如下所述的方法实
现:

采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(α),迭代求解第一目标函
数Gα(x,y)的最优解的最小值,即迭代求解公式(5)的最小值;

在k次迭代,H(α)的二阶Taylor级数?#24179;麳k(z):


其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z
与迭代值zk的二范数的平方;

k+1次迭代公式如下:


其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;

凸目标函数Gα(x,y)的最优解简化为:



然后,通过如下所述的步骤得到反射率图像:

步骤141,进行初始化,包括:

步骤141a,设每个像素反射率初始值为

步骤141b,对第一负对数似然函数,即公式(4)求导,计算梯度



步骤141c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;

步骤142,开始进行迭代,k=0,······,M;

步骤143,迭代求解简化后的凸目标函数Gα(x,y)的最优解,即迭代求解公式(17)、
(18),得到αk+1;

步骤144,判断αk+1是否满足迭代下降准则:


如果否(即αk+1不满足迭代下降准则),则αk+1←ηαk+1,然后再返回执行步骤143;

如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:


其中,η>0是很小的常数。

通过上述的方法,即可逐点迭代求解第一目标函数Gα(x,y)的凸优化问题,得到反
射率图像。

另外,较佳的,在本发明的具体实施例中,在得到反射率图像之后,即可输出该反
射率图像{α(x,y)},x=1,…,W,j=1,…,H。其中,H、W分别表示图像的行和列。

步骤15,根据空间相关性构造二元假设检验统计量,筛选背景光子噪声。

由于单光子探测器无法区分检测到的光子是来自目标返回的激光回波、背景光还
是暗计数,而背景光产生的光子以及探测器本身的暗计数导致光子到达时刻的负对数似然
函数呈现多模态,表现为非凸函数,是导致逐点最大似然估计失败的原因,因此检测到的光
子到达时间tj(x,y)必须在逐点似然估计距离之前进行背景噪声筛选,剔除噪声光子,即剔
除来自背景光的光子和探测器本身的暗计数,才能用于后续的三维结构估计。

在物理场景中三维空间结构的高度空间相关性隐含信号光产生的光子到达时刻
有非常小的条件方差,它们在时间上集中,在空间上相关;噪声光子的到达时刻是在时间区
间[0,Tr]上均匀分布,在空间位?#27809;?#30456;独立分布,有较高的方差。

因此,在本发明的技术方案中,可以根据空间相关性原理,构造二元假设检验统计
量H(x,y)来筛选背景光子噪声。

例如,较佳的,在本发明的具体实施例中,所述根据空间相关性构造二元假设检验
统计量可以包括如下所述的步骤:

步骤151,对于每一个像素位置(x,y),根据八邻域像素检测到的光子到达时间的
中值的平均值构造二元假设检验统计量H(x,y);

例如,较佳的,在本发明的具体实施例中,所述二元假设检验统计量H(x,y)可通过
如下所述的公式计算得到:



其中,八邻域像素检测到的光子到达时间分别为:

{tj(x-1,y-1),j=1,…,k1}、{tj(x-1,y),j=1,…,k2}、{tj(x-1,y+1),j=1,…,
k3}、

{tj(x,y-1),j=1,…,k4}、{tj(x,y+1),j=1,…,k5}、{tj(x+1,y-1),j=1,…,k6}、

{tj(x+1,y),j=1,…,k7}、{tj(x+1,y+1),j=1,…,k8}。

对应到达时间中值分别为tm1、tm2、tm3、tm4、tm5、tm6、tm7、tm8;

k1,k2,k3,k4,k5,k6,k7,k8分别表?#26223;?#37051;域像素接收到的光子数量。

因此可知,上述的二元假设检验统计量H(x,y)为八个邻域像素光子到达时间的中
值的平均值。

对于缺失光子数据的像素,H(x,y)=0。

步骤152,根据二元假设检验统计量H(x,y)设置筛选光子准则;

例如,较佳的,在本发明的具体实施例中,所述筛选光子准则可以是:


其中,αML(x,y)是反射率的最大似然估计,Tp是激光脉冲宽度的均方根。

因此,满足上述公式(7)的光子属于信号光,否则,该光子属于背景光。

在低照度条件下,单光子激光探测器测量大光子计数数据包含大量的背景光子或
暗计数,二元统计量H(x,y)是真实光子达到时间的较好近似值。单光子激光成像来说,α(x,
y)A+B<<1,激光脉冲宽度在纳秒级,公式(8)?#20918;?#39033;是很小的稳定值,不会导致二元统计量H
(x,y)筛选光子准则性能降低。

步骤16,建立光子到达时间与测量距离的测量模型,并构建光子达到时间条件下
测量距离的第二负对数似然函数。

在本步骤中,将建立光子到达时间tj(x,y)与测量距离R(x,y)的测量模型,即建立
光子到达时间tj(x,y)与深度图像R(x,y)的数量关系;然后,即可根据测量模型构建光子到
达时间tj(x,y)条件下测量距离R(x,y)的负对数似然函数L[R(x,y)|tj(x,y)],可称之为第
二负对数似然函数。

在真实世界中三维空间结构的高度空间相关性隐含信号光产生的光子到达时刻
具有非常小的条件方差,因此基于空间相关性构造二元假设检验统计量,审查光子到达时
刻tj(x,y)是否满足时间条件约束。

单光子探测器在一次激光脉冲持续期间只发生一次光子计数事件,只提供第一个
检测到的光子的计时信息,此光子计数过程具有独立增量特性。因此,光子到达时间tj(x,
y)可用返回光信号的时间平移波形来表征,而平移时间由测量距离R(x,y)?#33539;ǎ?#20174;而可以
建立光子到达时间tj(x,y)与测量距离R(x,y)的测量模型。

例如,较佳的,在本发明的具体实施例中,对于经过筛选后的光子数据,在激光脉
冲持续期间光子到达时间t(x,y)∈[0,Tr)的概率为:



其中,第一项中g(·)表示光信号的时间平移波形,系数为常数,第二项为常数。因
此,光子到达时间tj(x,y)的测量模型是测量距离R(x,y)的函数。

然后,可根据光子到达时间tj(x,y)与测量距离R(x,y)的测量模型,构造光子达到
时间tj(x,y)与距离R(x,y)的第二负对数似然函数L[R(x,y)|tj(x,y)],将三维距离像{R(x,
y))}估计转换为逐点求解由负对数似然函数与惩罚函数构成的目标函数JR(x,y)的最优
解。

例如,较佳的,在本发明的具体实施例中,可以根据驻留时间内光子到达时间的测
量模型,给定光子到达时间数据{tj(x,y),j=1,…,n(x,y)}条件下测量距离R(x,y)的第二
负对数似然函数L[R(x,y)|tj(x,y),j=1,…,n(x,y)]:


当n(x,y)=0时,L[R|tj,j=1,…,n](x,y)=0,对距离估计没有贡献。

激光脉冲的波形g(·)通常近似为:

g(t)∝exp[-v(t)] (11)

其中,t是激光脉冲持续时间,v(t)是激光脉冲波形的近似函数。


第二负对数似然函数L[R(x,y)|tj(x,y)]可以记为H(R)。

另外,在本发明的技术方案中,上述的步骤16可以与步骤12同时执行,也可以按照
预设执行顺序执行(例如,先执行步骤12再执行步骤16,或者相反),本发明对此并不做限
定。

步骤17,根据总变分半范数和第二负对数似然函数,构建第二目标函数JR(x,y)。

单光子探测器阵列获取的光子测量数据均为非负数,第二负对数似然函数L[R(x,
y)|tj(x,y)]是光子到达时间tj(x,y)的凸函数。因此,可以选择总变分半范数作为惩罚函
数。第二目标函数JR(x,y)也是光子到达时间tj(x,y)的凸函数,则估计三维距离像即是逐点
求解凸函数JR(x,y)的最优解。

经过二元假设检验筛选后的光子数据,才可用于深度图像{R(x,y)}估计。由筛选
后的像素位置的负对数似然函数,构建第二目标函数JR(x,y),考虑空间相关性约束,在目
标函数中选择总变分半范数作为惩罚函数。从所有像素检测到的光子到达时间{tj(x,y)}
估计深度图像{R(x,y)},即是逐点迭代求解负对数似然函数与惩罚函数构成的第二目标函
数JR(x,y)的凸优化问题的最优解。

例如,较佳的,在本发明的具体实施例中,可以通过如下所述的公式得到第二目标
函数JR(x,y)的最优解:



其中,RML(x,y)是最优解,βR是惩罚系数,惩罚函数pl(R)是凸函数。约束空间距离
估计的非?#20132;?#24615;,H、W分别表示图像的行和列。

另外,惩罚函数pl(R)可以定义为深度图像{R(x,y)}的总变分半范数:


步骤18,逐点迭代求解第二目标函数JR(x,y)的凸优化问题,得到并输出三维距离
像(即深度图像)。

测量数据的非负性在目标函数最优化中引入了一组不等式约束,导致最优化求解
变得困?#36873;?#22240;此,较佳的,在本发明的具体实施例中,针对非负约束的目标函数的凸优化问
题,可以使用一种二阶Taylor级数?#24179;?#20108;次序列目标函数的迭代数值求解算法。具体过程
与步骤14中的相关处理过程相类似。

例如,较佳的,在本发明的具体实施例中,所述步骤18可以通过如下所述的方法实
现:

采用序列二阶Taylor级数?#24179;?#21487;分解的对数似然函数H(R),迭代求解第二目标函
数JR(x,y)的最优解的最小值,即迭代求解公式(13)的最小值;

在k次迭代,H(R)的二阶Taylor级数?#24179;麳k(z):


其中,为负对数似然函数H(·)在zk处的一阶导数,参数σk>0,为向量z
与迭代值zk的二范数的平方;

k+1次迭代公式如下:


其中,z≥0,惩罚函数pl(α)为总变分半范数,?#29575;?#24809;罚系数;

凸目标函数JR(x,y)的最优解简化为:



然后,通过如下所述的步骤得到深度图像:

步骤181,进行初始化,包括:

步骤181a,经过筛选后的像素,初始深度值为其中,d为距
离分辨率,δ为测量偏差;未经过筛选后的像素

步骤181b,对第二负对数似然函数,即公式(12)求导,计算梯度


步骤181c,选择参数σk∈[σmin,σmax]、梯度下降步长η>0、最大迭代次数M∈Z+;

步骤182,开始进行迭代,k=0,······,M;

步骤183,迭代求解简化后的凸目标函数JR(x,y)的最优解,即迭代求解公式(17)、
(18),得到Rk+1;

步骤184,判断Rk+1是否满足迭代下降准则:


如果否(即Rk+1不满足迭代下降准则),则Rk+1←ηRk+1,然后再返回执行步骤183;

如果是,则k←k+1,然后进行下一次算法迭代,直到算法迭代满足收敛准则:


其中,η>0是很小的常数。

综上可知,本发明所提供的极低照度条件下光子计数激光三维计算成像方法,可
以精确测量在很短的驻留时间内每个GM-APD像素检测到的光子数量以及每个光子到达时
间,对测量光子计数数据进行筛选,剔除噪声光子数据,根据经过筛选的光子数量和到达时
间数据分别重构目标场景的反射率图像和三维距离像,实现了GM-APD探测器阵列同时生成
三维距离像和反射率图像,从而可以在光子计数模式下(极低照度条件、很短驻留时间)同
时获取高质量激光三维距离像和强度像。

该方法在激光成像探测、激光遥感和激光显微断层成像仪等领域有相当大的应用
价值。

以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精
神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明保护的范围之内。

关于本文
本文标题:一种极低照度条件下光子计数激光三维计算成像方法.pdf
链接地址:http://www.pqiex.tw/p-6079630.html
关于我们 - 网站声明 - 网?#38236;?#22270; - 资源地图 - 友情链接 - 网站客服 - 联系我们

[email protected] 2017-2018 zhuanlichaxun.net网?#26223;?#26435;所有
经营许可证编号:粤ICP备17046363号-1 
 


收起
展开
平码五不中公式规律 捕鱼大师稳赢版最新131 罗永浩锤子手机赚钱吗 安徽快三投注技巧 安徽时时计划软件qq 创业如何赚钱 手机版的奔驰宝马游戏 棋牌赢钱游戏 艇1000赢3万公式 蟹蟹优选赚钱 河北快三计划大神