[發明專利]一種光子計數目標高效率高分辨力深度重建方法有效

申請號：	202110484526.1	申請日：	2021-04-30
公開（公告）號：	CN113203475B	公開（公告）日：	2022-08-19
發明（設計）人：	沈姍姍;孫肖琳;蘇適;何瑞清;冒添翼	申請（專利權）人：	南京工業職業技術大學
主分類號：	G01J1/44	分類號：	G01J1/44;H04N5/335
代理公司：	南京燦爛知識產權代理有限公司 32356	代理人：	朱妃
地址：	210023 江蘇省***	國省代碼：	江蘇;32
權利要求書：	查看更多	說明書：	查看更多
摘要：
搜索關鍵詞：	一種光子計數目標高效率分辨力深度重建方法
鉆瓜網技術展會專利詞庫專利權人專利榜在售專利公布日期熱門專利

【權利要求書】：

1.一種光子計數目標高效率高分辨力深度重建方法，基于一體化光纖式偽隨機碼幅度調制深度矯正裝置的成像系統，其特征在于，包括以下步驟：

S1、將成像系統的準直器固定在二維導軌平臺上，去掉A單光子探測器、固定衰減器和分路器，偽隨機序列發生器產生一路同步高電平脈沖，該脈沖產生時間和開始發送偽隨機序列的時間一致，將同步高電平脈沖連接光子時間到達記錄儀，作為開始計時的觸發信號，偽隨機序列發生器驅動激光器發送激光脈沖，激光脈沖信號通過成像系統后得到單點擴頻光子計數值分布，表達為系統沖激響應、噪聲光子和攜帶目標特征的目標光子計數值分布的卷積積分；

其中，激光脈沖信號通過成像系統后得到單點擴頻光子計數值分布，表達為系統沖激響應、噪聲光子和攜帶目標特征的目標光子計數值分布的卷積積分，具體為，

激光脈沖發射后，先后通過大氣介質、目標和接收系統，在短距離下在到達接收系統前，脈沖回波信號s(n)由目標沖激響應決定；

通過接收系統后，目標沖激響應h_o(n)和探測器的沖激響應h_d(n)是造成接收的目標信號產生瞬時失真的主要因素，則偽隨機光脈沖碰到目標后，往返一次，回波到達接收系統，探測到的單點擴頻光子計數值分布R(n)為，

式中，R(n)為攜帶目標信息并且被系統沖激響應模糊的單點擴頻光子計數值分布，C(n)為互相關函數，為系統沖激響應，表征系統自身產生的瞬時干擾，f(n)＝h_o(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、采用高斯擬合法擬合平均的系統沖激響應，

設平均的系統沖激響應數據(o_i,H_i)i＝1,2,3,…,n，用高斯函數描述為，

式中，H_max、o_max和S為待估參數，H_max代表高斯曲線的峰高，o_max代表高斯曲線峰位置，S代表半高全寬；

將(3)式兩邊取自然對數，化為，

令InH_i＝z_i,

則(4)式化為二次多項式擬合函數，

考慮全部數據和測量誤差ε_i，并以矩陣形式表示如下，

簡記為Z_n×1＝O_n×3B_3×1+E_n×1，

根據最小二乘原理，求得擬合常數b₀﹑b₁﹑b₂構成的矩陣B的廣義最小二乘解為，

B＝(O^TO)^-1O^TZ (8)

進而根據(5)式，求出待估參數H_max﹑o_max和S，代入(3)式，從而得到擬合的高斯函數；

S3、恢復成像系統結構到初始狀態，基于最佳檢測的運算機制，通過構建合適的稀疏基，高效反演目標的擴頻光子計數值分布，具體為，

恢復成像系統結構，對目標成像，獲得接收路光子到達點時間，

S3.1、根據步驟S2.2的序列重構方法，重構接收路光子到達時間點，一個周期的序列長度為n_a，接收的經過延時后的光子到達時間信號記為

S3.2、生成數據壓縮比α，令M＝αn_a，以此為參量，隨機生成測量矩陣

S3.3、基于最佳檢測的運算機制，生成稀疏基Ψ，

基于最佳檢測的運算機制，采用一路偽隨機序列p(n)作為發送的參考碼型，另一路為接收的經過延時后的光子到達時間點，重構后得到序列x(n)，兩個序列最佳檢測的運算表示為，

其中，n_a為一個周期的序列長度，m為時間單元，

擴頻光子計數值分布R(m)的峰值對應橫坐標即為時間延時，其代表目標所處的距離信息；基于最佳檢測的運算方法，以參考序列為模板序列，擬采用包含時間位移的Hankel矩陣結構，構造稀疏基Ψ，該稀疏基的列元素包含所有可能的參考序列的時間位移，分以下兩種情況建立，

1)當積分的時間t_a和序列的周期t_p相等，即模板序列和接收序列的長度相等，則構建的稀疏基表示為，

其中，Ψ(:,i)代表稀疏基的第i列向量，zeros代表零向量，[]^T代表矩陣的轉置，p(i:j)為從第i個元素到第j個元素的參考向量；

2)當積分時間t_a大于序列的周期t_p，向量p右側補零，直至p和x長度相等；

S3.4、計算稀疏基的相干性δ(Ψ)

采用δ(Ψ)表示稀疏基的相干性，最劣的相干性即為兩個不同稀疏基元素之間內積的絕對值的最大值，δ(Ψ)的取值范圍為[0，1]，計算式(11)以防止出現相干的稀疏基元素對，

式中，Ψ_i為第i列元素構成的向量，Ψ_j為第j列元素構成的向量；

S3.5、重構擴頻光子計數值分布

設測量信號為為接收的經過延時后的光子到達時間信號，為變換系數，即待恢復的擴頻光子計數值分布，其描述光子在時間軸上的分布特征，其具有尖銳的峰值，在時間分布上具有稀疏特性，即||R||₀≤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計算R_k＝(A^TA+ρI)^-1[A^Ty+ρ(z_k-1-u_k-1)]，

其中R_k為第k次迭代計算得到的擴頻光子計數值分布，即式(1)中的R(n)，ρ為懲罰系數，I為單位矩陣，()^-1為矩陣的逆，()^T為矩陣的轉置，z_k-1和u_k-1均為第k-1次迭代的長度為2n_a-1的向量；

S3.5.3計算

其中prox為關于變量Z的L1范數近端算子，具有顯式解；

S3.5.4計算u_k＝u_k-1+R_k-z_k；

S3.6、判斷δ(Ψ)＜ε是否成立，如果小于設定數值ε，則k＝k+1，繼續采用S3.5的步驟迭代更新R，如果大于設定數值ε，則停止迭代，輸出R；

S4、采用基于貝葉斯理論和最大似然估計的迭代重建算法，引入矢量外推加速算法，估計最終攜帶目標特征的目標光子計數值分布，

具體為，

S4.1、初始化擬合的沖激響應﹑最大迭代次數k_max和加速參數α；

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)是迭代計算公式的一般形式，其中，ψ(f_k)表示線性搜索的迭代算法，對于一維擴頻光子計數值分布，f_k表示第k次迭代計算時產生的目標光子計數值的估計值，f_k+1表示第k+1次迭代計算時產生的目標光子計數值的估計值；

S4.3、基于差分表示的估計值變化引入

根據式(19)從現有的估計值f_k-1產生估計值f_k的迭代出發，給出一階泰勒矢量，

f_k＝ψ(f_k-1)＝f_k-1+g_k (20)

設f_k逐漸收斂到第s次迭代計算時產生的估計值f_s，在第k次迭代中引入的變化g_k朝著解的梯度方向▽m(f)生成一條通向多維空間的解，目標光子計數值的求解過程中，需要最大化的函數可局部近似為多維二次函數m(f)，則有，

m(f)＝-(f-f_s)^TA(f-f_s) (21)

▽m(f)＝-2A(f-f_s) (22)

式(22)中，()^T表示轉置，A表示平方、對稱、正定矩陣，

根據式(21)和(22)，采用差分形式給出第k次迭代引入的變化g_k，如式(23)所示，

g_k＝ε▽m(f_k)＝-2εA(f_k-f_s) (23)

式(19)中，根據矩陣對角化原理，通過矩陣A的特征值構成的對角矩陣Λ_A和特征向量S給出矩陣A，即，

A＝SΛ_AS^-1 (24)

設I為單位矩陣，Λ_B是矩陣B的特征值構成的對角矩陣，B＝I-2εA,Λ_B＝I-2εΛ_A，f₀表示估計值的初值，Λ_B^k為Λ_B的k次方，ε表示一個用來防止Λ_B出現負值的向量，根據式(20)、式(23)和式(24)，f_k為，

根據式(25)得到第k-1次迭代引入的變化g_k-1，k-δ-1次迭代引入的變化g_k-1-δ，k-δ次迭代的估計值f_k-δ，

g_k-1＝SΛ_B^k(I-Λ_B^-1)S^-1(f₀-f_s)+f_s (26)

f_k-δ＝SΛ_B^k-δS^-1(f₀-f_s)+f_s (27)

g_k-δ-1＝SΛ_B^k-δ(I-Λ_B^-1)S^-1(f₀-f_s) (28)

式(28)中，Λ_B^k-δ為Λ_B的k-δ次方，且有I-Λ_B^-δ≈δ(I-Λ_B^-1)；

S4.4、計算估計值之差

估計值f_k和估計值f_k-δ之差h_k,-δ表示為，

h_k,-δ＝f_k-f_k-δ＝SΛ_B^k(I-Λ_B^-δ)S^-1(f₀-f_s) (29)

對于第k+δ次迭代的估計值f_k+δ的預測為，

f_k+δ＝SΛ_B^k+δS^-1(f₀-f_s)+f_s (30)

第k+δ-1次迭代引入的變化g_k+δ-1為，

g_k+δ-1＝SΛ_B^k+δ(I-Λ_B^-1)S^-1(f₀-f_s) (31)

估計值f_k+δ和估計值f_k之差h_k,δ為，

h_k,δ＝f_k+δ-f_k＝SΛ_B^k+δ(I-Λ_B^-δ)S^-1(f₀-f_s) (32)

式(32)中，Λ_B^k+δ為Λ_B的k+δ次方，Λ_B^-δ為Λ_B的-δ次方，Λ_B^-1為Λ_B的-1次方，S^-1為S的-1次方；

S4.5、估計加速參數

基于泰勒級數的一階矢量外推法，通過最小二乘投影α_kh_k,-δ預測h_k,δ，得，

h_k,δ＝α_kh_k,-δ (33)

(h_k,δ-α_kh_k,-δ)^Th_k,-δ＝0 (34)

式(34)中，α_k為加速參數，

求解式(34)，得到，

又因為I-Λ_B^-δ≈δ(I-Λ_B^-1),將I-Λ_B^-δ≈δ(I-Λ_B^-1)代入式(32)可得，

h_k,δ≈δSΛ_B^k(I-Λ_B^-1)S^-1(f₀-f_s)＝δg_k-1 (36)

以此類推得到，

h_k,-δ≈δg_k-δ-1 (37)

將式(35)和式(37)代入式(35)，可得加速參數α_k為，

S4.6、估計攜帶目標特征的目標光子計數值

設第k次待估計的目標光子計數值為l_k，根據第k次估計中已經得到的f_k、f_k-1、α_k，代入式l_k＝f_k+α_k(f_k-f_k-1)，得到的l_k，k+1后，將l_k代入得到f_k+1，同理通過式(20)得到g_k+1，通過式(38)得到α_k+1，再代入l_k＝f_k+α_k(f_k-f_k-1)的變換式l_k+1＝f_k+1+α_k+1(f_k+1-f_k)，得到l_k+1，以此類推，對l_k進行循環更新，直到達到最大迭代次數k_max，獲得最終攜帶目標特征的目標光子計數值l_kmax；