CN113203475B - 一种光子计数目标高效率高分辨力深度重建方法 - Google Patents
一种光子计数目标高效率高分辨力深度重建方法 Download PDFInfo
- Publication number
- CN113203475B CN113203475B CN202110484526.1A CN202110484526A CN113203475B CN 113203475 B CN113203475 B CN 113203475B CN 202110484526 A CN202110484526 A CN 202110484526A CN 113203475 B CN113203475 B CN 113203475B
- Authority
- CN
- China
- Prior art keywords
- value
- target
- photon
- formula
- photon counting
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 238000009826 distribution Methods 0.000 claims abstract description 100
- 230000004044 response Effects 0.000 claims abstract description 55
- 238000001228 spectrum Methods 0.000 claims abstract description 49
- 238000003384 imaging method Methods 0.000 claims abstract description 48
- 230000001360 synchronised effect Effects 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 53
- 239000013598 vector Substances 0.000 claims description 36
- 238000001514 detection method Methods 0.000 claims description 29
- 238000005259 measurement Methods 0.000 claims description 28
- 238000004364 calculation method Methods 0.000 claims description 22
- 230000001133 acceleration Effects 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 12
- 230000003287 optical effect Effects 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 12
- 230000007246 mechanism Effects 0.000 claims description 11
- 230000003111 delayed effect Effects 0.000 claims description 10
- 238000007476 Maximum Likelihood Methods 0.000 claims description 9
- 108700041286 delta Proteins 0.000 claims description 9
- 230000010354 integration Effects 0.000 claims description 8
- 150000001875 compounds Chemical class 0.000 claims description 6
- 238000013213 extrapolation Methods 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 5
- 238000006073 displacement reaction Methods 0.000 claims description 5
- 239000013307 optical fiber Substances 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 238000012935 Averaging Methods 0.000 claims description 3
- 230000015556 catabolic process Effects 0.000 claims description 3
- 230000001427 coherent effect Effects 0.000 claims description 3
- 239000002131 composite material Substances 0.000 claims description 3
- 238000005314 correlation function Methods 0.000 claims description 3
- 238000013144 data compression Methods 0.000 claims description 3
- 238000006731 degradation reaction Methods 0.000 claims description 3
- 238000012886 linear function Methods 0.000 claims description 3
- 238000010606 normalization Methods 0.000 claims description 3
- 238000012887 quadratic function Methods 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 229910001316 Ag alloy Inorganic materials 0.000 claims description 2
- BQCADISMDOOEFD-UHFFFAOYSA-N Silver Chemical compound [Ag] BQCADISMDOOEFD-UHFFFAOYSA-N 0.000 claims description 2
- 125000004122 cyclic group Chemical group 0.000 claims description 2
- 230000001747 exhibiting effect Effects 0.000 claims 1
- 238000011084 recovery Methods 0.000 description 9
- 238000000926 separation method Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000007906 compression Methods 0.000 description 3
- 230000006835 compression Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 238000005265 energy consumption Methods 0.000 description 2
- 230000036314 physical performance Effects 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 238000013459 approach Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000004304 visual acuity Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J1/00—Photometry, e.g. photographic exposure meter
- G01J1/42—Photometry, e.g. photographic exposure meter using electric radiation detectors
- G01J1/44—Electric circuits
-
- H—ELECTRICITY
- H04—ELECTRIC COMMUNICATION TECHNIQUE
- H04N—PICTORIAL COMMUNICATION, e.g. TELEVISION
- H04N25/00—Circuitry of solid-state image sensors [SSIS]; Control thereof
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J1/00—Photometry, e.g. photographic exposure meter
- G01J1/42—Photometry, e.g. photographic exposure meter using electric radiation detectors
- G01J1/44—Electric circuits
- G01J2001/4413—Type
- G01J2001/442—Single-photon detection or photon counting
Landscapes
- Engineering & Computer Science (AREA)
- Multimedia (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明公开了一种光子计数目标高效率高分辨力深度重建方法，包括以下步骤，S1、将成像系统的准直器固定在二维导轨平台上，去掉A单光子探测器、固定衰减器和分路器，伪随机序列发生器产生一路同步高电平脉冲，该脉冲产生时间和开始发送伪随机序列的时间一致，将同步高电平脉冲连接光子时间到达记录仪，作为开始计时的触发信号，伪随机序列发生器驱动激光器发送激光脉冲，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统的冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分。省去单独的噪声去除步骤，以最快的速度方向，加速估计下一次迭代的数值，运算速度快于线性搜索的迭代算法。
Description
技术领域
本发明涉及时间相关单光子计数深度成像技术领域，具体为一种光子计数目标高效率高分辨力深度重建方法。
背景技术
近年来，光子计数型伪随机码幅度调制二维成像方法的研究已经广泛展开，成像过程中，存在以下两个问题：
1)通过两路光子到达时间点的最佳检测运算，得到扩频光子计数值分布，这就是传统的基于最佳检测的相关接收技术。为了得到高精度的扩频光子计数值分布，往往通过延长积分时间来提高信噪比，以至于增大了采集的数据量、加重数据处理负担、降低成像速度和信号恢复效率。在特殊环境下，存在后向散射噪声光的干扰，而信号光衰减加剧，光子计数率很低，通常为每秒几个或更少，由于信息量的不足，难以恢复具有尖锐峰值的扩频光子计数值分布。
2)构建系统的单光子探测器、时间到达记录仪等设备引入系统的瞬时干扰，该瞬时干扰可采用系统的冲激响应来描述。以卷积积分为理论依据，即输入信号通过线性系统冲激响应的移位、倍乘、叠加运算。信号通过系统后计算获得扩频光子计数值分布，该波形包含系统的瞬时干扰项(Tamer F.Refaat，Syed Ismail，M.Nurul Abedin，ScottM.Spuler，Shane D.Mayor，and Upendra N.Singh.Lidar backscatter signal recoveryfrom phototransistor systematic effect by deconvolution[J].Appl.Opt，2008，47:5281-5288.Joel F.Campbell，Bing Lin，Amin R.Nehrir，F.Wallace Harrison，andMichael D.Obland.Super-resolution technique for CW lidar using Fouriertransform reordering and Richardson–Lucy deconvolution[J].Opt.Lett，2014，39:6981-6984.)，导致失真的扩频光子计数值分布，使得距离分辨率往往无法达到伪随机码发送速率决定的距离分辨率，例如，码发送速率为1GHz，则最小距离单元为一个码元宽度1ns，对应距离为15厘米。实际探测中，可分辨的最小距离达不到15厘米。针对复杂的目标场景，单个光斑内往往包含多个距离的反射光子，当目标在横向和纵向距离上较为接近时，无法分辨多个距离，造成目标探测的遗漏，目标成像的细节性不高等问题，如图1所示，图1-(a)表示在激光光斑内，纵向方向上存在多个距离分布，探测器接收到两种距离的返回回波光子，纵向间隔距离分布比较密集时，目标的光子统计直方图会产生重叠如图1-(c)上方，采用峰值搜索法提取扩频光子计数值分布的峰值时，往往会遗漏部分被模糊的距离；光斑大小不变，在纵向分辨能力的限制下，图1-(b)中，当对横向方向上两个距离分布比较密集的目标成像时，无法分辨两个目标，如图1-(c)右下方所示。
发明内容
本发明的目的在于提供一种光子计数目标高效率高分辨力深度重建方法，基于一体化光纤式伪随机码幅度调制偏移校正装置(申请号为201710325775.X)的成像系统，后期数据处理中，采用压缩感知的方法高校恢复扩频光子计数值分布，省去单独的噪声去除步骤，以最快的速度方向，加速估计下一次迭代的数值，运算速度快于线性搜索的迭代算法。
为达到上述目的，根据本发明的一个方面，本发明提供如下技术方案：
一种光子计数目标高效率高分辨力深度重建方法，基于一体化光纤式伪随机码幅度调制深度矫正装置的成像系统，其特征在于，包括以下步骤：
S1、将成像系统的准直器固定在二维导轨平台上，去掉A单光子探测器、固定衰减器和分路器，伪随机序列发生器产生一路同步高电平脉冲，该脉冲产生时间和开始发送伪随机序列的时间一致，将同步高电平脉冲连接光子时间到达记录仪，作为开始计时的触发信号，伪随机序列发生器驱动激光器发送激光脉冲，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分；
其中，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分，具体为，
激光脉冲发射后，先后通过大气介质、目标和接收系统，在短距离下在到达接收系统前，脉冲回波信号s(n)由目标冲激响应决定；
通过接收系统后，目标冲激响应ho(n)和探测器的冲激响应hd(n)是造成接收的目标信号产生瞬时失真的主要因素，则伪随机光脉冲碰到目标后，往返一次，回波到达接收系统，探测到的单点扩频光子计数值分布R(n)为，
式中，R(n)为携带目标信息并且被系统冲激响应模糊的单点扩频光子计数值分布，C(n)为互相关函数，
为系统冲激响应，表征系统自身产生的瞬时干扰，f(n)＝ho(n)为携带目标特征的目标光子计数值分布，B为噪声光子；
S2、改变S1中的成像系统结构，两个准直器面对面发射并接收光子，获取平均的系统冲激响应并拟合，具体为，
S2.1、在成像系统上，断开可调衰减器和光环行器的连接，将准直器B和可调衰减器连接，作为收准直器端，将连接光环行器的准直器A作为发准直器端，将准直器A和准直器B面对面放置，通过微调和粗调准直器的位置，使得从准直器A中发出的光信号导入接收的准直器B中，再由B单光子探测器接收后，得到接收路光子时间到达点；
S2.2、将matlab的序列重构方法嵌入Labview，基于Labview，获取接收路的光子时间到达点，并重构0﹑1序列，表示为y(n)，伪随机序列发生器中存储的模板序列p(n)与y(n)进行快速傅里叶变换，得到扩频光子计数值分布，即系统冲激响应H，其代表系统的瞬时干扰，
式中，
表示傅里叶变换，
为傅里叶逆变换；
S2.3、测量若干次的系统冲激响应并作平均，得到平均的系统冲激响应；
S2.4、采用高斯拟合法拟合平均的系统冲激响应，
设平均的系统冲激响应数据(oi,Hi)i＝1,2,3,…,n，用高斯函数描述为，
式中，Hmax、omax和S为待估参数，Hmax代表高斯曲线的峰高，omax代表高斯曲线峰位置，S代表半高全宽；
将(3)式两边取自然对数，化为，
令
则(4)式化为二次多项式拟合函数，
考虑全部数据和测量误差εi，并以矩阵形式表示如下，
简记为Zn×1＝On×3B3×1+En×1，
根据最小二乘原理，求得拟合常数b0﹑b1﹑b2构成的矩阵B的广义最小二乘解为，
B＝(OTO)-1OTZ (8)
进而根据(5)式，求出待估参数Hmax﹑omax和S，代入(3)式，从而得到拟合的高斯函数；
S3、恢复成像系统结构到初始状态，基于最佳检测的运算机制，通过构建合适的稀疏基，高效反演目标的扩频光子计数值分布，具体为，
恢复成像系统结构，对目标成像，获得接收路光子到达点时间，
S3.1、根据步骤S2.2的序列重构方法，重构接收路光子到达时间点，一个周期的序列长度为na，接收的经过延时后的光子到达时间信号记为
S3.2、生成数据压缩比α，令M＝αna，以此为参量，随机生成测量矩阵
S3.3、基于最佳检测的运算机制，生成稀疏基Ψ，
基于最佳检测的运算机制，采用一路伪随机序列p(n)作为发送的参考码型，另一路为接收的经过延时后的光子到达时间点，重构后得到序列x(n)，两个序列最佳检测的运算表示为，
其中，na为一个周期的序列长度，m为时间单元，
扩频光子计数值分布R(m)的峰值对应横坐标即为时间延时，其代表目标所处的距离信息；基于最佳检测的运算方法，以参考序列为模板序列，拟采用包含时间位移的Hankel矩阵结构，构造稀疏基Ψ，
该稀疏基的列元素包含所有可能的参考序列的时间位移，分以下两种情况建立，
1)当积分的时间ta和序列的周期tp相等，即模板序列和接收序列的长度相等，则构建的稀疏基表示为，
其中，Ψ(:,i)代表稀疏基的第i列向量，zeros代表零向量，[]T代表矩阵的转置，p(i:j)为从第i个元素到第j个元素的参考向量；
3)当积分时间ta大于序列的周期tp，向量p右侧补零，直至p和x长度相等；
S3.4、计算稀疏基的相干性δ(Ψ)
采用δ(Ψ)表示稀疏基的相干性，最劣的相干性即为两个不同稀疏基元素之间内积的绝对值的最大值，δ(Ψ)的取值范围为[0，1]，计算式(11)以防止出现相干的稀疏基元素对，
式中，Ψi为第i列元素构成的向量，Ψj为第j列元素构成的向量；
S3.5、重构扩频光子计数值分布
设测量信号为
为接收的经过延时后的光子到达时间信号，
为变换系数，即待恢复的扩频光子计数值分布，其描述光子在时间轴上的分布特征，其具有尖锐的峰值，在时间分布上具有稀疏特性，即||R||0≤K，是K稀疏的，则有
在稀疏基
的表示下是K稀疏的，基于压缩感知的测量过程理解为对经过延时后的光子到达时间信号的M次独立测量，每次以线性函数
的形式测量接收数据的线性投影
M个测量向量构成的测量矩阵为
其将高维信号x映射到低维，则测量过程表述为
y＝Φx＝Φ(ΨR) (12)
设v表示噪声，带噪观测情况下，式(12)表示为，
y＝AR+v (13)
其中A＝ΦΨ；
S3.5.1式(13)转化为以下lasso标准问题的求解，λ为代价参数，
S3.5.2计算Rk＝(ATA+ρI)-1[ATy+ρ(zk-1-uk-1)]，
其中Rk为第k次迭代计算得到的扩频光子计数值分布，即式(1)中的R(n)，ρ为惩罚系数，I为单位矩阵，()-1为矩阵的逆，()T为矩阵的转置，zk-1和uk-1均为第k-1次迭代的长度为2na-1的向量；
S3.5.3计算
其中prox为关于变量Z的L1范数近端算子，具有显式解；
S3.5.4计算uk＝uk-1+Rk-zk；
S3.6、判断δ(Ψ)＜ε是否成立，如果小于设定数值ε，则k＝k+1，继续采用S3.5的步骤迭代更新R，如果大于设定数值ε，则停止迭代，输出R；
S4、采用基于贝叶斯理论和最大似然估计的迭代重建算法，引入矢量外推加速算法，估计最终携带目标特征的目标光子计数值分布，
具体为，
S4.1、初始化拟合的冲激响应﹑最大迭代次数kmax和加速参数α；
S4.2、设携带目标特征的目标光子计数值分布为f，简称为目标光子计数值，由贝叶斯模型可得，
式(14)中，p(f/R)为已知扩频光子计数值分布情况下，目标光子计数值的概率分布情况，p(R/f)为已知目标光子计数值的情况下，扩频光子计数值概率分布情况，p(f)为目标光子计数概率分布；当p(f/R)取得最大数值，则恢复的目标信号接近理想的波形，
式(14)的分母部分为常量，通过最大似然估计法求解其最大值，
ML(f)＝max[p(R/f)] (15)
式(15)中，ML(f)表示采用最大似然估计的方法估计f，max为最大值，
设待恢复的波形，即目标光子计数值f，其与系统的冲激响应H的卷积表示为
T表示每个像素的积分时间，则联合概率密度分布L(R/f)为，
式(16)等号两边求对数得到，
对等式(17)两边求导，并令其为0，根据退化函数归一化的性质，可得，
引入迭代运算，则有，
式(19)是迭代计算公式的一般形式，其中，ψ(fk)表示线性搜索的迭代算法，对于一维扩频光子计数值分布，fk表示第k次迭代计算时产生的目标光子计数值的估计值，fk+1表示第k+1次迭代计算时产生的目标光子计数值的估计值；
S4.3、基于差分表示的估计值变化引入
根据式(19)从现有的估计值fk-1产生估计值fk的迭代出发，给出一阶泰勒矢量，
fk＝ψ(fk-1)＝fk-1+gk (20)
设fk逐渐收敛到第s次迭代计算时产生的估计值fs，在第k次迭代中引入的变化gk朝着解的梯度方向▽m(f)生成一条通向多维空间的解，目标光子计数值的求解过程中，需要最大化的函数可局部近似为多维二次函数m(f)，则有，
m(f)＝-(f-fs)TA(f-fs) (21)
式(22)中，()T表示转置，A表示平方、对称、正定矩阵，
根据式(21)和(22)，采用差分形式给出第k次迭代引入的变化gk，如式(23)所示，
式(19)中，根据矩阵对角化原理，通过矩阵A的特征值构成的对角矩阵ΛA和特征向量S给出矩阵A，即，
A＝SΛAS-1 (24)
设I为单位矩阵，ΛB是矩阵B的特征值构成的对角矩阵，B＝I-2εA,ΛB＝I-2εΛA，f0表示估计值的初值，ΛB k为ΛB的k次方，ε表示一个用来防止ΛB出现负值的向量，根据式(20)、式(23)和式(24)，fk为，
根据式(25)得到第k-1次迭代引入的变化gk-1，k-δ-1次迭代引入的变化gk-1-δ，k-δ次迭代的估计值fk-δ，
gk-1＝SΛB k(I-ΛB -1)S-1(f0-fs)+fs (26)
fk-δ＝SΛB k-δS-1(f0-fs)+fs (27)
gk-δ-1＝SΛB k-δ(I-ΛB -1)S-1(f0-fs) (28)
式(28)中，ΛB k-δ为ΛB的k-δ次方，且有I-ΛB -δ≈δ(I-ΛB -1)；
S4.4、计算估计值之差
估计值fk和估计值fk-δ之差hk,-δ表示为，
hk,-δ＝fk-fk-δ＝SΛB k(I-ΛB -δ)S-1(f0-fs) (29)
对于第k+δ次迭代的估计值fk+δ的预测为，
fk+δ＝SΛB k+δS-1(f0-fs)+fs (30)
第k+δ-1次迭代引入的变化gk+δ-1为，
gk+δ-1＝SΛB k+δ(I-ΛB -1)S-1(f0-fs) (31)
估计值fk+δ和估计值fk之差hk,δ为，
hk,δ＝fk+δ-fk＝SΛB k+δ(I-ΛB -δ)S-1(f0-fs) (32)
式(32)中，ΛB k+δ为ΛB的k+δ次方，ΛB -δ为ΛB的-δ次方，ΛB -1为ΛB的-1次方，S-1为S的-1次方；
S4.5、估计加速参数
基于泰勒级数的一阶矢量外推法，通过最小二乘投影αkhk,-δ预测hk,δ，得，
hk,δ＝αkhk,-δ (33)
(hk,δ-αkhk,-δ)Thk,-δ＝0 (34)
式(34)中，αk为加速参数，
求解式(34)，得到，
又因为I-ΛB -δ≈δ(I-ΛB -1),将I-ΛB -δ≈δ(I-ΛB -1)代入式(32)可得，
hk,δ≈δSΛB k(I-ΛB -1)S-1(f0-fs)＝δgk-1 (36)
以此类推得到，
hk,-δ≈δgk-δ-1 (37)
将式(35)和式(37)代入式(35)，可得加速参数αk为，
S4.6、估计携带目标特征的目标光子计数值
设第k次待估计的目标光子计数值为lk，根据第k次估计中已经得到的fk、fk-1、αk，代入式lk＝fk+αk(fk-fk-1)，得到的lk，k+1后，将lk代入
得到fk+1，同理通过式(20)得到gk+1，通过式(38)得到αk+1，再代入lk＝fk+αk(fk-fk-1)的变换式lk+1＝fk+1+αk+1(fk+1-fk)，得到lk+1，以此类推，对lk进行循环更新，直到达到最大迭代次数kmax，获得最终携带目标特征的目标光子计数值lkmax；
S5、根据最终携带目标特征的目标光子计数值分布，通过残差去除﹑峰值搜索和质心拟合方法提取距离数值，具体为，
S5.1、通过最终携带目标特征的目标光子计数值lkmax，最高幅度值选取阈值β，提取出满足lkmax＜β对应的所有k值，对应k的lkmax幅度值置为0，得到的数据记为lkmax_opt，即第k个深度单元对应的目标光子计数值；
S5.2、搜索lkmax_opt中的峰值，并且提取峰值对应的k值；
S5.3、在搜索到的峰值附近提取L个非零数值点，采用质心拟合算法计算距离数值d，
式(39)中，τ为最小的距离分辨单元；
S6、逐点重构距离数值，得到重构的二维深度图像，具体为，
将准直器固定在二维导轨平台上，扫描获得每个像素点的光子到达时间点，采用步骤S3～步骤S5所述的方法，逐点重构距离数值，得到重构的二维深度图像。
与现有技术相比，本发明具有的优势是：提出将压缩感知信号恢复方法引入基于最佳检测的单光子计数成像方法中，确立稀疏表示的信息高效恢复技术。通过稀疏基稀疏表达隐藏的延时信息，利用最佳检测的运算核心结构来构建稀疏基，并提出双相干性约束，进一步保证信息恢复的效率。该方法突破成像系统中前端接收系统的物理性能极限，省去冗余数据的存取，降低计算负担。丰富了压缩感知理论和最佳检测理论的科学内涵，为信息的高效恢复提供新的应用价值。并且本发明省去单独的噪声去除步骤，以最快的速度方向，加速估计下一次迭代的数值，运算速度快于线性搜索的迭代算法，具备节约硬件资源、降低能源损耗、提高成像效率的优势；提高了深度成像的横向和纵向分辨能力。
附图说明
图1为本发明成像目标结构示意图；
图2为实时数据处理流程图，图2-(a)Labview接收序列数据存取，图2-(b)Matlab实时接收序列重构算法；
图3为高效率高分辨力成像方法流程；
图4为系统冲激响应图，其中实线为测量的平均的系统冲激响应，虚线为拟合的冲激响应；
图5为最小均方误差对应最大迭代次数的确定图。
具体实施方式
下面结合说明书附图，对本发明作进一步的说明。
与现有技术相比，本发明具有的优势是：提出将压缩感知信号恢复方法引入基于最佳检测的单光子计数成像方法中，确立结构化稀疏表示的信息高效恢复技术。利用最佳检测的算法核心结构来构建稀疏基，通过稀疏基稀疏表达隐藏的延时信息，并提出双相干性约束，进一步保证信息恢复的效率。该方法突破成像系统中前端接收系统的物理性能极限，省去冗余数据的存取，降低计算负担。丰富了压缩感知理论和最佳检测的理论的科学内涵，为信息的高效恢复提供新的应用价值。
并且本发明省去单独的噪声去除步骤，以最快的速度方向，加速估计下一次迭代的数值，运算速度快于线性搜索的迭代算法，具备节约硬件资源、降低能源损耗、提高成像效率的优势；提高了深度成像的横向和纵向分辨能力。
一种光子计数目标高效率高分辨力深度重建方法，基于一体化光纤式伪随机码幅度调制偏移校正装置(申请号为201710325775.X)的成像系统，包括以下步骤：
S1、将成像系统的准直器固定在二维导轨平台上，去掉A单光子探测器、固定衰减器和分路器，伪随机序列发生器产生一路同步高电平脉冲，该脉冲产生时间和开始发送伪随机序列的时间一致，将同步高电平脉冲连接光子时间到达记录仪，作为开始计时的触发信号，伪随机序列发生器驱动激光器发送激光脉冲，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统的冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分，
其中，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统的冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分，具体为，
激光脉冲发射后，先后通过大气介质、目标和接收系统，在短距离下在到达接收系统前，脉冲回波信号s(n)由目标冲激响应决定；
通过接收系统后，目标冲激响应ho(n)和探测器的冲激响应hd(n)是造成接收的目标信号产生瞬时失真的主要因素，则伪随机光脉冲碰到目标后，往返一次，回波到达接收系统，探测到的单点扩频光子计数值分布R(n)为，
式中，R(n)为携带目标信息并且被系统的冲激响应模糊的单点扩频光子计数值分布，C(n)为互相关函数，
为系统的冲激响应，表征系统自身产生的瞬时干扰，f(n)＝ho(n)为携带目标特征的目标光子计数值分布，B为噪声光子。
S2、改变S1中的成像系统结构，将两个准直器面对面发射并接收光子，获取平均的系统冲激响应并拟合，具体为，
S2.1、在一体化光纤式伪随机码幅度调制深度矫正装置的成像系统上，断开可调衰减器和光环行器的连接，将准直器B和可调衰减器连接，作为收准直器端，将连接光环行器的准直器A作为发准直器端，将准直器A和准直器B面对面放置，通过微调和粗调准直器的位置，使得从准直器A中发出的光信号导入接收的准直器B中，再由B单光子探测器接收后，得到接收路光子时间到达点；
S2.2、基于Labview，获取接收路的光子时间到达点如图2(a)，并采用matlab重构0﹑1序列如图2(b)，从获取的时间到达点中重构发序列码型，根据码发送速率和系统最小深度分辨率，比如码发送速率为2.5GHz，则码元宽度为400ps，时间相关记录仪分辨率为4ps，则时间分辨因子为400ps/4ps＝100，根据该通道的第一个时间到达点，通过时间分辨因子确定第一个‘1’码元的位置，再依次确定第二个、第三个的‘1’码元。接着，判断两个‘1’码元，即两个时间到达点之间序列零的个数，最后通过序列的对比和补零，得到长度一致的收发序列，基于Matlab实现重构算法，获得一定积分时间下的序列码型。通过循环发送序列，一路是已经生成好的模板序列，一路是首尾相接的接收序列码型。伪随机码是循环发送的，通过接收路的光子时间到达点重构0﹑1序列，表示为y(n)，模板序列p(n)与y(n)进行快速傅里叶变换，得到扩频光子计数值分布，即系统冲激响应H，其代表系统的瞬时干扰，
式中，
表示傅里叶变换，
为傅里叶逆变换；
S2.3、测量50次的系统冲激响应并作平均，得到平均的系统冲激响应；
S2.4、采用高斯拟合法拟合平均的系统冲激响应，
设平均的系统冲激响应数据(oi,Hi)i＝1,2,3,…,n，用高斯函数描述为，
式中，Hmax、omax和S为待估参数，Hmax代表高斯曲线的峰高，omax代表高斯曲线峰位置，S代表半高全宽；
将(3)式两边取自然对数，化为，
令
则(4)式化为二次多项式拟合函数，
考虑全部数据和测量误差εi，并以矩阵形式表示如下，
简记为Zn×1＝On×3B3×1+En×1，
根据最小二乘原理，求得拟合常数b0﹑b1﹑b2构成的矩阵B的广义最小二乘解为，
B＝(OTO)-1OTZ (8)
进而根据(5)式，求出待估参数Hmax﹑omax和S，代入(3)式，从而得到拟合的高斯函数。图4实线为测量得到的平均的系统冲激响应，采用高斯拟合法拟合平均的系统冲激响应如图4虚线，拟合曲线的半高全宽(FWHM)约为1.4ns，意味着仅仅采用峰值搜索法可分辨的目标的最小纵向间隔约为21厘米。再加上背景噪声、复杂的目标属性等外界条件，使得两个距离峰值不断融合，距离波形被模糊，有必要引入信号处理算法，高效去除系统的瞬时干扰，提高距离分辨率。
获得了拟合的冲激响应后，二维图像的重构流程如图3所示。
S3、恢复成像系统结构并成像，基于最佳检测的运算机制，通过稀疏基表示方法，高效反演目标的单点扩频光子计数值分布，具体为，
恢复成像系统结构，对目标成像，获得接收路光子到达点时间，
S3.1、根据步骤s2.2的序列重构方法，重构接收路光子到达时间点，一个周期的序列长度为na，记为
S3.2、生成数据压缩比α，令M＝αna，以此为参量，随机生成测量矩阵
S3.3、基于最佳检测的运算机制核心机制，生成稀疏基Ψ，
基于最佳检测的运算机制，采用一路伪随机序列p(n)作为发送的参考码型，另一路为接收的经过延时后的光子到达时间点，重构后得到序列x(n)，两个序列最佳检测的运算如下，
其中，na为一个周期的序列长度，m为时间单元，
扩频光子计数值分布R(m)的峰值对应横坐标即为时间延时，其代表目标所处的距离信息，其核心机制是接收序列和参考序列的各次位移相乘再求和的过程。基于最佳检测的运算方法，以参考序列为模板序列，拟采用包含时间位移的Hankel矩阵结构，构造参考序列结构化的稀疏基，参考序列结构化的稀疏基的列元素包含所有可能的参考序列的时间位移，其分以下两种情况建立：
1)如果积分时间ta和序列的周期tp相等，即模板序列和接收序列的长度相等，则构建的稀疏基表示为，
其中，Ψ(:,i)代表第i列向量，zeros代表零向量，[]T代表矩阵的转置，p(i:j)为从第i个元素到第j个元素的参考向量；
4)当积分时间ta大于序列的周期tp，向量p右侧补零，直至p和x长度相等；
S3.4、计算稀疏基的相干性δ(Ψ)
采用δ(Ψ)表示稀疏基的相干性，最劣的相干性即为两个不同稀疏基元素之间内积的绝对值的最大值，δ(Ψ)的取值范围为[0，1]，计算式(11)以防止出现相干的稀疏基元素对，
式中，Ψi为第i列元素构成的向量，Ψj为第j列元素构成的向量；
S3.5、重构扩频光子计数值分布
设测量信号为
为接收的经过延时后的光子到达时间信号，
为变换系数，即待恢复的扩频光子计数值分布，其描述光子在时间轴上的分布特征，其具有尖锐的峰值，在时间分布上具有稀疏特性，即||R||0≤K，是K稀疏的，则有
在稀疏基
的表示下是K稀疏的，基于压缩感知的测量过程理解为对经过延时后的光子到达时间信号的M次独立测量，每次以线性函数
的形式测量接收数据的线性投影
M个测量向量构成的测量矩阵为
其将高维信号x映射到低维，则测量过程表述为
y＝Φx＝Φ(ΨR) (12)
设v表示噪声，带噪观测情况下，式(12)表示为，
y＝AR+v (13)
其中A＝ΦΨ；
S3.5.1式(13)转化为以下lasso标准问题的求解，λ为代价参数，
S3.5.2计算Rk＝(ATA+ρI)-1[ATy+ρ(zk-1-uk-1)]
其中Rk为第k次迭代计算得到的扩频光子计数值分布，即式(1)中的R(n)，ρ为惩罚系数，I为单位矩阵，()-1为矩阵的逆，()T为矩阵的转置，zk-1和uk-1均为第k-1次迭代的长度为2na-1的向量；
S3.5.3计算
其中prox为关于变量Z的L1范数近端算子，具有显式解；
S3.5.4计算uk＝uk-1+Rk-zk；
S3.6、判断δ(Ψ)＜ε是否小于ε，如果小于ε，则k＝k+1，继续采用S3.5的步骤迭代更新R，如果大于设定数值，则停止迭代，输出R。
S4、采用基于贝叶斯理论和最大似然估计的迭代重建算法，引入矢量外推加速算法，估计最终携带目标特征的目标光子计数值分布，具体为，
S4.1、初始化拟合的冲激响应﹑最大迭代次数kmax和加速参数α；
S4.2、设携带目标特征的目标光子计数值分布为f，简称为目标光子计数值，由贝叶斯模型可得，
式(14)中，p(f/R)为已知扩频光子计数值分布情况下，目标光子计数值的概率分布情况，p(R/f)为已知目标光子计数值的情况下，扩频光子计数值概率分布情况，p(f)为目标光子计数概率分布；当p(f/R)取得最大数值，则恢复的目标信号接近理想的波形，
式(14)的分母部分为常量，通过最大似然估计法求解其最大值，
ML(f)＝max[p(R/f)] (15)
式(15)中，ML(f)表示采用最大似然估计的方法估计f，max为最大值，
设待恢复的波形，即目标光子计数值f，其与系统的冲激响应H的卷积表示为
T表示每个像素的积分时间，则联合概率密度分布L(R/f)为，
式(16)等号两边求对数得到，
对等式(17)两边求导，并令其为0，根据退化函数归一化的性质，可得，
引入迭代运算，则有，
式(19)是迭代计算公式的一般形式，其中，ψ(fk)表示线性搜索的迭代算法，对于一维扩频光子计数值分布，fk表示第k次迭代计算时产生的目标光子计数值的估计值，fk+1表示第k+1次迭代计算时产生的目标光子计数值的估计值；
S4.3、基于差分表示的估计值变化引入
根据式(19)从现有的估计值fk-1产生估计值fk的迭代出发，给出一阶泰勒矢量，
fk＝ψ(fk-1)＝fk-1+gk (20)
设fk逐渐收敛到第s次迭代计算时产生的估计值fs，在第k次迭代中引入的变化gk朝着解的梯度方向▽m(f)生成一条通向多维空间的解，目标光子计数值的求解过程中，需要最大化的函数可局部近似为多维二次函数m(f)，则有，
m(f)＝-(f-fs)TA(f-fs) (21)
式(22)中，()T表示转置，A表示平方、对称、正定矩阵，
根据式(21)和(22)，采用差分形式给出第k次迭代引入的变化gk，如式(23)所示，
式(19)中，根据矩阵对角化原理，通过矩阵A的特征值构成的对角矩阵ΛA和特征向量S给出矩阵A，即，
A＝SΛAS-1 (24)
设I为单位矩阵，ΛB是矩阵B的特征值构成的对角矩阵，B＝I-2εA,ΛB＝I-2εΛA，f0表示估计值的初值，ΛB k为ΛB的k次方，ε表示一个用来防止ΛB出现负值的向量，根据式(20)、式(23)和式(24)，fk为，
根据式(25)得到第k-1次迭代引入的变化gk-1，k-δ-1次迭代引入的变化gk-1-δ，k-δ次迭代的估计值fk-δ，
gk-1＝SΛB k(I-ΛB -1)S-1(f0-fs)+fs (26)
fk-δ＝SΛB k-δS-1(f0-fs)+fs (27)
gk-δ-1＝SΛB k-δ(I-ΛB -1)S-1(f0-fs) (28)
式(28)中，ΛB k-δ为ΛB的k-δ次方，且有I-ΛB -δ≈δ(I-ΛB -1)；
S4.4、计算估计值之差
估计值fk和估计值fk-δ之差hk,-δ表示为，
hk,-δ＝fk-fk-δ＝SΛB k(I-ΛB -δ)S-1(f0-fs) (29)
对于第k+δ次迭代的估计值fk+δ的预测为，
fk+δ＝SΛB k+δS-1(f0-fs)+fs (30)
第k+δ-1次迭代引入的变化gk+δ-1为，
gk+δ-1＝SΛB k+δ(I-ΛB -1)S-1(f0-fs) (31)
估计值fk+δ和估计值fk之差hk,δ为，
hk,δ＝fk+δ-fk＝SΛB k+δ(I-ΛB -δ)S-1(f0-fs) (32)
式(32)中，ΛB k+δ为ΛB的k+δ次方，ΛB -δ为ΛB的-δ次方，ΛB -1为ΛB的-1次方，S-1为S的-1次方；
S4.5、估计加速参数
基于泰勒级数的一阶矢量外推法，通过最小二乘投影αkhk,-δ预测hk,δ，得，
hk,δ＝αkhk,-δ (33)
(hk,δ-αkhk,-δ)Thk,-δ＝0 (34)
式(34)中，αk为加速参数，
求解式(34)，得到，
又因为I-ΛB -δ≈δ(I-ΛB -1),将I-ΛB -δ≈δ(I-ΛB -1)代入式(32)可得，
hk,δ≈δSΛB k(I-ΛB -1)S-1(f0-fs)＝δgk-1 (36)
以此类推得到，
hk,-δ≈δgk-δ-1 (37)
将式(35)和式(37)代入式(35)，可得加速参数αk为，
S4.6、估计携带目标特征的目标光子计数值
设第k次待估计的目标光子计数值为lk，根据第k次估计中已经得到的fk、fk-1、αk，代入式lk＝fk+αk(fk-fk-1)，得到的lk，k+1后，将lk代入
得到fk+1，同理通过式(20)得到gk+1，通过式(38)得到αk+1，再代入lk＝fk+αk(fk-fk-1)的变换式lk+1＝fk+1+αk+1(fk+1-fk)，得到lk+1，以此类推，对lk进行循环更新，直到达到最大迭代次数kmax，获得最终携带目标特征的目标光子计数值lkmax。
S5、根据最终携带目标特征的目标光子计数值分布，通过残差去除﹑峰值搜索和质心拟合方法提取距离数值，具体为，
S5.1、通过最终携带目标特征的目标光子计数值lkmax，最高幅度值选取阈值β，提取出满足lkmax＜β对应的所有k值，对应k的lkmax幅度值置为0，得到的数据记为lkmax_opt，即第k个深度单元对应的目标光子计数值；
S5.2、搜索lkmax_opt中的峰值，并且提取峰值对应的k值；
S5.3、在搜索到的峰值附近提取L个非零数值点，采用质心拟合算法计算距离数值d，
式(39)中，τ为最小的距离分辨单元。
S6、逐点重构距离数值，得到重构的二维深度图像，具体为，
将准直器固定在二维导轨平台上，扫描获得每个像素点的光子到达时间点，采用步骤S3～步骤S6所述的方法，如图3所示，逐点重构距离数值，得到重构的二维深度图像。
为了描述单点测距的性能，定义m次测量的单点距离均方误差为：
其中，
为第j次估计的距离值，dactual为实际距离数值，Tc为码元宽度，c为光速。
定义m次测量的单点距离标准差为：
其中，
为估计值的均值。
为了初始化最大迭代次数，需要定义纵向距离间隔均方误差(RMSE)为：
实验中，首先，对7米处的白墙成像，每个像素点获取的信号光子计数值为5400个，噪声信号光子计数值为600个，总共为6000个。分别采用50％、30％、10％三种压缩比，则每个像素点获取的光子计数值分别为3000个、1800个和600个，采用图3中带噪声高效重构方法得到扩频光子计数分布。测量50次，得到该分布的峰值对应距离数值的RMSEsingle和STDsingle如表1所示。同等条件下，采用最佳检测的方法重构的扩频光子计数分布的峰值对应距离数值的RMSEsingle和STDsingle分别为3.1cm和2.7cm。也即采用50％的光子计数值即可达到最佳检测的方法的距离精确度，如表1所示。
其次，对放置在白墙前面的积木成一维距离像，积木和白墙间隔8厘米，成像目标如图1(a)所示。将激光打在积木边缘，根据图3高效率高分辨力成像方法流程，采用步骤3-6，先通过带噪声高效重构方法得到扩频光子计数值分布，以最小的RMSE得到最大迭代次数kmax为5次，如图5所示，根据步骤S2计算得到拟合的冲激响应，再采用快速迭代重构目标光子计数值方法得到高分辨力距离像。表2给出8厘米间隔，50次测量的间隔均方误差。表3为在同等条件下，采用传统的最佳检测获得扩频光子计数值分布，再采用Ancombe变换和快速反卷积算法提高分辨能力。比较表2和3，本发明提出的方法只需要获取60％的光子计数值，即可以达到原方法估计的距离间隔均方误差，省去了单独的去噪声算法，降低了数据存取复杂度，提高运算速度和运算效率。
表1带噪高效重建方法距离估计值统计量
| 压缩比 | RMSE<sub>single</sub>[cm] | STD<sub>single</sub>[cm] |
| 50％ | 3 | 2.1 |
| 30％ | 9 | 8.6 |
| 10％ | 14.4 | 12.7 |
表2高效率高分辨力间隔估计值统计量
| 压缩比 | RMSE[cm] | STD[cm] |
| 60％ | 3.9 | 2.9 |
| 50％ | 8.5 | 6.3 |
| 30％ | 11.5 | 10.9 |
| 10％ | 20.8 | 18.9 |
表3传统最佳检测高分辨力间隔估计值统计量
| RMSE[cm] | STD[cm] |
| 3.8 | 2.7 |
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (1)
1.一种光子计数目标高效率高分辨力深度重建方法，基于一体化光纤式伪随机码幅度调制深度矫正装置的成像系统，其特征在于，包括以下步骤：
S1、将成像系统的准直器固定在二维导轨平台上，去掉A单光子探测器、固定衰减器和分路器，伪随机序列发生器产生一路同步高电平脉冲，该脉冲产生时间和开始发送伪随机序列的时间一致，将同步高电平脉冲连接光子时间到达记录仪，作为开始计时的触发信号，伪随机序列发生器驱动激光器发送激光脉冲，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分；
其中，激光脉冲信号通过成像系统后得到单点扩频光子计数值分布，表达为系统冲激响应、噪声光子和携带目标特征的目标光子计数值分布的卷积积分，具体为，
激光脉冲发射后，先后通过大气介质、目标和接收系统，在短距离下在到达接收系统前，脉冲回波信号s(n)由目标冲激响应决定；
通过接收系统后，目标冲激响应ho(n)和探测器的冲激响应hd(n)是造成接收的目标信号产生瞬时失真的主要因素，则伪随机光脉冲碰到目标后，往返一次，回波到达接收系统，探测到的单点扩频光子计数值分布R(n)为，
式中，R(n)为携带目标信息并且被系统冲激响应模糊的单点扩频光子计数值分布，C(n)为互相关函数，
为系统冲激响应，表征系统自身产生的瞬时干扰，f(n)＝ho(n)为携带目标特征的目标光子计数值分布，B为噪声光子；
S2、改变S1中的成像系统结构，两个准直器面对面发射并接收光子，获取平均的系统冲激响应并拟合，具体为，
S2.1、在成像系统上，断开可调衰减器和光环行器的连接，将准直器B和可调衰减器连接，作为收准直器端，将连接光环行器的准直器A作为发准直器端，将准直器A和准直器B面对面放置，通过微调和粗调准直器的位置，使得从准直器A中发出的光信号导入接收的准直器B中，再由B单光子探测器接收后，得到接收路光子时间到达点；
S2.2、将matlab的序列重构方法嵌入Labview，基于Labview，获取接收路的光子时间到达点，并重构0﹑1序列，表示为y(n)，伪随机序列发生器中存储的模板序列p(n)与y(n)进行快速傅里叶变换，得到扩频光子计数值分布，即系统冲激响应H，其代表系统的瞬时干扰，
式中，
表示傅里叶变换，
为傅里叶逆变换；
S2.3、测量若干次的系统冲激响应并作平均，得到平均的系统冲激响应；
S2.4、采用高斯拟合法拟合平均的系统冲激响应，
设平均的系统冲激响应数据(oi,Hi)i＝1,2,3,…,n，用高斯函数描述为，
式中，Hmax、omax和S为待估参数，Hmax代表高斯曲线的峰高，omax代表高斯曲线峰位置，S代表半高全宽；
将(3)式两边取自然对数，化为，
令InHi＝zi,
则(4)式化为二次多项式拟合函数，
考虑全部数据和测量误差εi，并以矩阵形式表示如下，
简记为Zn×1＝On×3B3×1+En×1，
根据最小二乘原理，求得拟合常数b0﹑b1﹑b2构成的矩阵B的广义最小二乘解为，
B＝(OTO)-1OTZ (8)
进而根据(5)式，求出待估参数Hmax﹑omax和S，代入(3)式，从而得到拟合的高斯函数；
S3、恢复成像系统结构到初始状态，基于最佳检测的运算机制，通过构建合适的稀疏基，高效反演目标的扩频光子计数值分布，具体为，
恢复成像系统结构，对目标成像，获得接收路光子到达点时间，
S3.1、根据步骤S2.2的序列重构方法，重构接收路光子到达时间点，一个周期的序列长度为na，接收的经过延时后的光子到达时间信号记为
S3.2、生成数据压缩比α，令M＝αna，以此为参量，随机生成测量矩阵
S3.3、基于最佳检测的运算机制，生成稀疏基Ψ，
基于最佳检测的运算机制，采用一路伪随机序列p(n)作为发送的参考码型，另一路为接收的经过延时后的光子到达时间点，重构后得到序列x(n)，两个序列最佳检测的运算表示为，
其中，na为一个周期的序列长度，m为时间单元，
扩频光子计数值分布R(m)的峰值对应横坐标即为时间延时，其代表目标所处的距离信息；基于最佳检测的运算方法，以参考序列为模板序列，拟采用包含时间位移的Hankel矩阵结构，构造稀疏基Ψ，
该稀疏基的列元素包含所有可能的参考序列的时间位移，分以下两种情况建立，
1)当积分的时间ta和序列的周期tp相等，即模板序列和接收序列的长度相等，则构建的稀疏基表示为，
其中，Ψ(:,i)代表稀疏基的第i列向量，zeros代表零向量，[]T代表矩阵的转置，p(i:j)为从第i个元素到第j个元素的参考向量；
2)当积分时间ta大于序列的周期tp，向量p右侧补零，直至p和x长度相等；
S3.4、计算稀疏基的相干性δ(Ψ)
采用δ(Ψ)表示稀疏基的相干性，最劣的相干性即为两个不同稀疏基元素之间内积的绝对值的最大值，δ(Ψ)的取值范围为[0，1]，计算式(11)以防止出现相干的稀疏基元素对，
式中，Ψi为第i列元素构成的向量，Ψj为第j列元素构成的向量；
S3.5、重构扩频光子计数值分布
设测量信号为
为接收的经过延时后的光子到达时间信号，
为变换系数，即待恢复的扩频光子计数值分布，其描述光子在时间轴上的分布特征，其具有尖锐的峰值，在时间分布上具有稀疏特性，即||R||0≤K，是K稀疏的，则有
在稀疏基
的表示下是K稀疏的，基于压缩感知的测量过程理解为对经过延时后的光子到达时间信号的M次独立测量，每次以线性函数
的形式测量接收数据的线性投影
M个测量向量构成的测量矩阵为
其将高维信号x映射到低维，则测量过程表述为
y＝Φx＝Φ(ΨR) (12)
设v表示噪声，带噪观测情况下，式(12)表示为，
y＝AR+v (13)
其中A＝ΦΨ；
S3.5.1式(13)转化为以下lasso标准问题的求解，λ为代价参数，
S3.5.2计算Rk＝(ATA+ρI)-1[ATy+ρ(zk-1-uk-1)]，
其中Rk为第k次迭代计算得到的扩频光子计数值分布，即式(1)中的R(n)，ρ为惩罚系数，I为单位矩阵，()-1为矩阵的逆，()T为矩阵的转置，zk-1和uk-1均为第k-1次迭代的长度为2na-1的向量；
S3.5.3计算
其中prox为关于变量Z的L1范数近端算子，具有显式解；
S3.5.4计算uk＝uk-1+Rk-zk；
S3.6、判断δ(Ψ)＜ε是否成立，如果小于设定数值ε，则k＝k+1，继续采用S3.5的步骤迭代更新R，如果大于设定数值ε，则停止迭代，输出R；
S4、采用基于贝叶斯理论和最大似然估计的迭代重建算法，引入矢量外推加速算法，估计最终携带目标特征的目标光子计数值分布，
具体为，
S4.1、初始化拟合的冲激响应﹑最大迭代次数kmax和加速参数α；
S4.2、设携带目标特征的目标光子计数值分布为f，简称为目标光子计数值，由贝叶斯模型可得，
式(14)中，p(f/R)为已知扩频光子计数值分布情况下，目标光子计数值的概率分布情况，p(R/f)为已知目标光子计数值的情况下，扩频光子计数值概率分布情况，p(f)为目标光子计数概率分布；当p(f/R)取得最大数值，则恢复的目标信号接近理想的波形，
式(14)的分母部分为常量，通过最大似然估计法求解其最大值，
ML(f)＝max[p(R/f)] (15)
式(15)中，ML(f)表示采用最大似然估计的方法估计f，max为最大值，
设待恢复的波形，即目标光子计数值f，其与系统的冲激响应H的卷积表示为
T表示每个像素的积分时间，则联合概率密度分布L(R/f)为，
式(16)等号两边求对数得到，
对等式(17)两边求导，并令其为0，根据退化函数归一化的性质，可得，
引入迭代运算，则有，
式(19)是迭代计算公式的一般形式，其中，ψ(fk)表示线性搜索的迭代算法，对于一维扩频光子计数值分布，fk表示第k次迭代计算时产生的目标光子计数值的估计值，fk+1表示第k+1次迭代计算时产生的目标光子计数值的估计值；
S4.3、基于差分表示的估计值变化引入
根据式(19)从现有的估计值fk-1产生估计值fk的迭代出发，给出一阶泰勒矢量，
fk＝ψ(fk-1)＝fk-1+gk (20)
设fk逐渐收敛到第s次迭代计算时产生的估计值fs，在第k次迭代中引入的变化gk朝着解的梯度方向▽m(f)生成一条通向多维空间的解，目标光子计数值的求解过程中，需要最大化的函数可局部近似为多维二次函数m(f)，则有，
m(f)＝-(f-fs)TA(f-fs) (21)
▽m(f)＝-2A(f-fs) (22)
式(22)中，()T表示转置，A表示平方、对称、正定矩阵，
根据式(21)和(22)，采用差分形式给出第k次迭代引入的变化gk，如式(23)所示，
gk＝ε▽m(fk)＝-2εA(fk-fs) (23)
式(19)中，根据矩阵对角化原理，通过矩阵A的特征值构成的对角矩阵ΛA和特征向量S给出矩阵A，即，
A＝SΛAS-1 (24)
设I为单位矩阵，ΛB是矩阵B的特征值构成的对角矩阵，B＝I-2εA,ΛB＝I-2εΛA，f0表示估计值的初值，ΛB k为ΛB的k次方，ε表示一个用来防止ΛB出现负值的向量，根据式(20)、式(23)和式(24)，fk为，
根据式(25)得到第k-1次迭代引入的变化gk-1，k-δ-1次迭代引入的变化gk-1-δ，k-δ次迭代的估计值fk-δ，
gk-1＝SΛB k(I-ΛB -1)S-1(f0-fs)+fs (26)
fk-δ＝SΛB k-δS-1(f0-fs)+fs (27)
gk-δ-1＝SΛB k-δ(I-ΛB -1)S-1(f0-fs) (28)
式(28)中，ΛB k-δ为ΛB的k-δ次方，且有I-ΛB -δ≈δ(I-ΛB -1)；
S4.4、计算估计值之差
估计值fk和估计值fk-δ之差hk,-δ表示为，
hk,-δ＝fk-fk-δ＝SΛB k(I-ΛB -δ)S-1(f0-fs) (29)
对于第k+δ次迭代的估计值fk+δ的预测为，
fk+δ＝SΛB k+δS-1(f0-fs)+fs (30)
第k+δ-1次迭代引入的变化gk+δ-1为，
gk+δ-1＝SΛB k+δ(I-ΛB -1)S-1(f0-fs) (31)
估计值fk+δ和估计值fk之差hk,δ为，
hk,δ＝fk+δ-fk＝SΛB k+δ(I-ΛB -δ)S-1(f0-fs) (32)
式(32)中，ΛB k+δ为ΛB的k+δ次方，ΛB -δ为ΛB的-δ次方，ΛB -1为ΛB的-1次方，S-1为S的-1次方；
S4.5、估计加速参数
基于泰勒级数的一阶矢量外推法，通过最小二乘投影αkhk,-δ预测hk,δ，得，
hk,δ＝αkhk,-δ (33)
(hk,δ-αkhk,-δ)Thk,-δ＝0 (34)
式(34)中，αk为加速参数，
求解式(34)，得到，
又因为I-ΛB -δ≈δ(I-ΛB -1),将I-ΛB -δ≈δ(I-ΛB -1)代入式(32)可得，
hk,δ≈δSΛB k(I-ΛB -1)S-1(f0-fs)＝δgk-1 (36)
以此类推得到，
hk,-δ≈δgk-δ-1 (37)
将式(35)和式(37)代入式(35)，可得加速参数αk为，
S4.6、估计携带目标特征的目标光子计数值
设第k次待估计的目标光子计数值为lk，根据第k次估计中已经得到的fk、fk-1、αk，代入式lk＝fk+αk(fk-fk-1)，得到的lk，k+1后，将lk代入
得到fk+1，同理通过式(20)得到gk+1，通过式(38)得到αk+1，再代入lk＝fk+αk(fk-fk-1)的变换式lk+1＝fk+1+αk+1(fk+1-fk)，得到lk+1，以此类推，对lk进行循环更新，直到达到最大迭代次数kmax，获得最终携带目标特征的目标光子计数值lkmax；
S5、根据最终携带目标特征的目标光子计数值分布，通过残差去除﹑峰值搜索和质心拟合方法提取距离数值，具体为，
S5.1、通过最终携带目标特征的目标光子计数值lkmax，最高幅度值选取阈值β，提取出满足lkmax＜β对应的所有k值，对应k的lkmax幅度值置为0，得到的数据记为lkmax_opt，即第k个深度单元对应的目标光子计数值；
S5.2、搜索lkmax_opt中的峰值，并且提取峰值对应的k值；
S5.3、在搜索到的峰值附近提取L个非零数值点，采用质心拟合算法计算距离数值d，
式(39)中，τ为最小的距离分辨单元；
S6、逐点重构距离数值，得到重构的二维深度图像，具体为，
将准直器固定在二维导轨平台上，扫描获得每个像素点的光子到达时间点，采用步骤S3～步骤S5所述的方法，逐点重构距离数值，得到重构的二维深度图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110484526.1A CN113203475B (zh) | 2021-04-30 | 2021-04-30 | 一种光子计数目标高效率高分辨力深度重建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110484526.1A CN113203475B (zh) | 2021-04-30 | 2021-04-30 | 一种光子计数目标高效率高分辨力深度重建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113203475A CN113203475A (zh) | 2021-08-03 |
| CN113203475B true CN113203475B (zh) | 2022-08-19 |
Family
ID=77029693
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110484526.1A Active CN113203475B (zh) | 2021-04-30 | 2021-04-30 | 一种光子计数目标高效率高分辨力深度重建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113203475B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115902835B (zh) * | 2021-09-30 | 2024-02-27 | 深圳市速腾聚创科技有限公司 | 一种雷达数据收发装置、测距方法及激光雷达 |
| CN115001587B (zh) * | 2022-06-03 | 2024-05-28 | 上海交通大学 | 单光子傅立叶变换信号分析方法 |
| CN114884575A (zh) * | 2022-06-20 | 2022-08-09 | 济南量子技术研究院 | 一种可靠性增强的双向qkd系统及其光纤链路监测方法 |
| CN116776901B (zh) * | 2023-08-25 | 2024-04-30 | 深圳市爱德泰科技有限公司 | 一种应用于电力通信机房的光纤配线架标签管理系统 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107015233A (zh) * | 2017-05-10 | 2017-08-04 | 南京理工大学紫金学院 | 一体化光纤式伪随机码幅度调制偏移校正装置 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20110260036A1 (en) * | 2010-02-22 | 2011-10-27 | Baraniuk Richard G | Temporally- And Spatially-Resolved Single Photon Counting Using Compressive Sensing For Debug Of Integrated Circuits, Lidar And Other Applications |
| KR101209908B1 (ko) * | 2011-08-04 | 2012-12-11 | 광주과학기술원 | 희소 신호 전송 방법 및 장치, 그리고 희소 신호 복구 방법 및 장치 |
-
2021
- 2021-04-30 CN CN202110484526.1A patent/CN113203475B/zh active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107015233A (zh) * | 2017-05-10 | 2017-08-04 | 南京理工大学紫金学院 | 一体化光纤式伪随机码幅度调制偏移校正装置 |
Non-Patent Citations (1)
| Title |
|---|
| 沈姗姗 等.光子计数深度获取系统误差机理研究和矫正.《光子学报》.2017, * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113203475A (zh) | 2021-08-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113203475B (zh) | 一种光子计数目标高效率高分辨力深度重建方法 | |
| JP5994024B2 (ja) | Ｐｅｔイメージングにおける単一イベントリスト型のデータ同期方法及びシステム | |
| Berdyugina | Surface imaging by the Occamian approach. Basic principles, simulations, and tests | |
| CN101509972B (zh) | 基于高分辨目标距离像修正相关矩阵的宽带雷达检测方法 | |
| CN110688763A (zh) | 一种基于脉冲型ToF相机深度和光强图像的多径效应补偿方法 | |
| CN110095766B (zh) | 基于非均匀重采样技术的机动目标相干积累检测方法 | |
| CN105375931B (zh) | 一种基于卡尔曼滤波的复杂环境下信号重构方法 | |
| US11415697B2 (en) | Real-time image formation from Geiger-mode LADAR | |
| Castorena et al. | Compressive sampling of lidar: Full-waveforms as signals of finite rate of innovation | |
| Castorena et al. | Modeling lidar scene sparsity using compressive sensing | |
| CN115508788A (zh) | 基于稳健噪底估计的信道化检测方法 | |
| CN109270344B (zh) | 脉冲丢失下的相参脉冲信号频率估计方法 | |
| CN103985093B (zh) | 基于多随机测量迭代像素判决的压缩感知稳健重构方法 | |
| Chen et al. | An Efficient Phase Unwrapping Method Based on Unscented Kalman Filter | |
| CN113109777B (zh) | 基于stokes矢量分解的宽带极化雷达目标检测方法 | |
| CN113049098B (zh) | 基于Richardson–Lucy反卷积的穿透散射介质高分辨率成像方法 | |
| Hu et al. | A Joint TDOA, FDOA and Doppler Rate Parameters Estimation Method and Its Performance Analysis | |
| CN109344783A (zh) | 一种合成信号的确定方法及系统 | |
| Yang et al. | A fast Wigner Hough transform algorithm for parameter estimation of low probability of intercept radar polyphase coded signals | |
| Li et al. | An Unfolded Deep Network for SAR Imaging Based on General Regularization and S-TLS Model | |
| US9762217B2 (en) | Random sampler adapted to one-dimension slow-varying signal | |
| CN113687325B (zh) | 基于lp和hrrp模型的遮蔽小目标检测方法 | |
| Zhang et al. | Seismic data recovery with curvelet bivariate shrinkage function based on compressed sensing | |
| CN113219429B (zh) | 基于多测量压缩感知下的捷变频雷达高速目标重构方法 | |
| CN113702970A (zh) | 一种基于2d-fomp的二维联合稀疏成像算法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |