一種激光雷達全波形信號分析方法
【專利摘要】本發明提出一種激光雷達全波形信號分析方法。本發明方法在每次更新操作時在回波信號的數量、入射到激光雷達探測器的噪聲信號、每個回波信號的峰值位置和每個回波信號的峰值幅中隨機選擇任意一個參量繼續進行更新操作；參量采用Metropolis算法思想進行更新，參量采用Metropolis算法或者模擬回火算法思想進行更新，參量和分別采用模擬回火算法思想進行更新。本發明能夠快速處理特性未知目標的全波形信號，且不受全波形信號信噪比的影響。
【專利說明】一種激光雷達全波形信號分析方法
【技術領域】
[0001]本發明屬于激光雷達【技術領域】，具體涉及一種激光雷達全波形信號分析方法。
【背景技術】
[0002]傳統的激光雷達通常只能提取激光光斑內的一個距離特征，無法處理光斑內多個距離特征目標的回波信號。采用具有完整地記錄接收信號功能的全波形激光雷達可以得到這種包含多個回波信號的全波形信號，為全波形信號分析提供了基礎，所謂全波形信號由激光光斑印記內所有特征距離的回波信號共同產生。為了有效解析單個激光光斑印記內包含的距離特征，需要采用全波形信號分析技術提取出其中包含的回波信號的峰值位置、數量、峰值幅度、噪聲共4個關鍵的全波形參量，分別標記為β、Β。其中，回波數量和峰值位置是全波形參量是最為重要的兩個參數，用于標記距離特征的數量和相互位置關系并完成全波形信號分解，峰值位置即準確距離上的光子飛行時間。
[0003]全波形信號y可以看成是接收信號F的后驗抽樣結果，以能量在時間軸上的能量分布形式來表示。全波形信號分析的目的就是利用全波形信號I求解全波形參量tpk、β、B，其過程可以看作是利用后驗信息對先驗信息進行估計。其中，后驗信息是全波形信號y，而相對應的先驗信息是全波形參量tpk、β、Β。通常可以采用后驗分布概率密度函數來評價全波形參量估計值的準確程度。
[0004]典型的全波形分析策略主要包括采用LM算法的非線性最小二乘波形擬合方法和采用EM算法及其改進算法的全波形參量最大似然值估計方法。前者需要事先預知目標的有效數據模型，只適用于單一特性的場景，例如:植被、建筑等；后者的收斂性能強烈依賴于初始值的預設，無法實現完全的非人工干預全波形參量提取，有時存在失效的風險。
[0005]近來，一種基于貝葉斯估計框架的可逆跳躍馬爾可夫蒙特卡羅算法(RJMCMC)(詳見 IEEE Transactions on Pattern Analysis and Machine Intelligence2007, 29 (12), 2170-2180)被提出作為一種全波形分析策略。在此策略中，其關鍵之處是采用RJMCMC算法按照固定的順序以循環形式逐個進行待解全波形參量h、k、β和B，并由后驗分布概率密度函數評估所有由參量近似解構成的全局近似解(tpk、β、Β)，挑選最為合適的全局近似解(tpk、β、Β)作為全局最優解。該方法允許在更新過程中待解全波形參量可以在不同尺度的預設分布模型之間進行跳躍完成對分布空間的探測而不依賴先驗信息，從而實現對目標特性未知的全波形信號分析。但該方法對分布空間的快速探測主要依賴于有效的預設分布，而且各參量的轉移概率通常較低，從而導致較慢的運算速度。
【發明內容】
[0006]本發明提出一種激光雷達全波形信號分析方法，能夠快速處理特性未知目標的全波形信號，且不受全波形信號信噪比的影響。
[0007]為了解決上述技術問題，本發明提供一種激光雷達全波形信號分析方法，包括以下步驟:[0008]步驟一、輸入一個激光雷達全波形信號y，明確待求參量為:‘、k、β、B，其中:
[0009]k為回波信號的數量，其取值范圍為:{0，1，2Kkmax}，kmax是預設的k的最大值；
[0010]B為入射到激光雷達探測器的噪聲信號，其取值范圍為:(0，max(y))，max(y)是全波形信號I的最大值；
[0011]h為每個回波信號中的光子準確飛行時間構成的k維向量，其元素對應k個回波信號的峰值位置，to的所有元素的取值范圍為:(0，ifflax)，Ifflax是全波形信號y的長度；
[0012]β為每個回波信號的峰值幅度構成的k維向量，其元素對應k個回波信號的峰值幅度，β的所有元素的取值范圍為:(0，max (y))；
[0013]從各待求參量的取值范圍內選取任意值分別賦值給參量h、k、β和B，完成各待求參量的初始化；在完成初始化的待求參量中選取任意一個參量，以被選中的參量開始進行一次更新操作；
[0014]步驟二、將更新操作后的所有待求參量組成一個全局近似解(tpk、β、Β)，在每次更新操作后組成的全局近似解(h、k、β、B)中隨機選擇任意一個參量繼續進行更新操作；重復本步驟，直至更新操作次數的總和達到預先設定的更新操作次數的最大值后停止更新操作；
[0015]步驟三、將每一次更新操作后獲得的所有全局近似解分別代入后驗分布概率密度函數，求出每一個全局近似解的后驗分布概率值，將值最大的后驗分布概率值對應的全局近似解作為全局最優解；所述后驗分布概率密度函數如公式(I)所示:
【權利要求】
1.一種激光雷達全波形信號分析方法，其特征在于，包括以下步驟: 步驟一、輸入一個激光雷達全波形信號1，明確待求參量為Apk、β、B，其中: k為回波信號的數量，其取值范圍為:{O，1，2Kkmax}，kmax是預設的k的最大值； B為入射到激光雷達探測器的噪聲信號，其取值范圍為:(O，max (y))，max (y)是全波形信號I的最大值； h為每個回波信號中的光子準確飛行時間構成的k維向量，其元素對應k個回波信號的峰值位置，t0的所有元素的取值范圍為:(0，ifflax)，Ifflax是全波形信號y的長度； β為每個回波信號的峰 值幅度構成的k維向量，其元素對應k個回波信號的峰值幅度，β的所有元素的取值范圍為:(0，max(y))； 從各待求參量的取值范圍內選取任意值分別賦值給參量tpk、β和B，完成各待求參量的初始化；在完成初始化的待求參量中選取任意一個參量，以被選中的參量開始進行一次更新操作； 步驟二、將更新操作后的所有待求參量組成一個全局近似解(h、k、β、B)，在每次更新操作后組成的全局近似解(tpk、β、B)中隨機選擇任意一個參量繼續進行更新操作；重復本步驟，直至更新操作次數的總和達到預先設定的更新操作次數的最大值后停止更新操作； 步驟三、將每一次更新操作后獲得的所有全局近似解分別代入后驗分布概率密度函數，求出每一個全局近似解的后驗分布概率值，將值最大的后驗分布概率值對應的全局近似解作為全局最優解；所述后驗分布概率密度函數如公式(1)所示:
2.如權利要求1所述的激光雷達全波形信號分析方法，其特征在于，當一次更新操作選定的參量為k時，可以進一步在k'和k''這兩種狀態中隨機選擇k經過本次更新操作后可能達到的狀態，然后根據隨機選擇的結果，按照k'和k''分別對應的不同的更新方法進行本次更新操作；或者在k'、k' '、h、β和B中選擇任意一個參量,將該參量對應的更新方法作為本次更新操作的方法，其中V = k+1, , = k-1 ； 所述隨機選擇的，k經過本次更新操作后可能達到的狀態為k'時，本次更新操的計算過程為: 2.1.1、若k = kmax，則各參量Vk、β和B保持不變，跳出本次更新操作，重新選擇參量進行下一次更新操作；否則繼續步驟2.1.2 ; 2.1.2、設〖/表示h經過本次更新操作后可能到達的狀態，且t0'為k,維向量，t0'的前k個元素依次為h的前k個元素； 設β '表示β經過本次更新操作后可能到達的狀態，且β '為k'維向量，β，的前k個元素依次為h的前k個元素； 在區間(0，imax)中選擇任意一個服從均勻分布的隨機數賦值給t/的第k'個元素，在區間(0，max(y))中選擇任意一個服從均勻分布的隨機數賦值給β /的第k'個元素； 2.1.3、將y、h和β代入公式(2)求解依據y、、和β預估的回波信號的能量分布{f(i, β j, t0J)}, j是不大于k的任意非負整數，β j和、分別是向量β和h的第j個元素；將y、t/和β'代入公式(2)求解依據y't/和β'預估的回波信號的能量分布{f(i, β/ ,t0/ )}，m是不大于k'的任意非負整數，β/和分別是向量β'和t/的第m個元素； 然后將{f(i，和B代入公式(3)求解預估的入射到激光雷達探測器的信號
,to/' )}、k'和B代入公式(3)求解預估的入射到激光雷達探測器的信號F(i，β '，tQ'，k'，B); 再將F(i, β，t0, k, B)和y代入公式(4)求解全波形信號y在信號F(i, β，t0, k, B)的聯合分布概率Uyltci, k，β，Β)，將F(i，β ' ,t/，k'，B)和y代入公式(4)求解全波形信號y在信號F(i，i3'，tQ'，k'，B)的聯合分布概率L(y|t0'，k' , β ' ,B)； 最后將L(y I tQ, k, β，B)、L(y 11/，k' , β '，B)代入公式(5)求解k'的轉移概率a (k/，k)，公式(3)、(4)、(5)具體如下所示:
3.如權利要求1所述的激光雷達全波形信號分析方法，其特征在于，當一次更新操作選定的參量為B時，本次更新操作的過程為: `3.1、設B'是B可能更新到的狀態，在區間(0，max(y))選擇一個服從伽瑪分布的隨機數賦值給B'； `3.2、首先將y、、和β代入公式(2)求解依據y、、和β預估的回波信號的能量分布{f(i, β j, t0J)}, j是不大于k的任意非負整數，β j和、是β和h的第j個元素； 然后將{f(i，和B代入公式(3)求解依據{f(i，β和B預估的入射到激光雷達探測器的信號F(i，β,^Ι?,Β)，將{f(i，β和B'代入公式(3)求解依據{f(i，i3j，tQj)}、k和B'預估的入射到激光雷達探測器的信號F(i，i3，tQ，k，B')； 再將F(i，β , t0, k, B )和y代入公式(4)求解y在信號F (i，β , t0, k, B)的聯合分布概率L(y|t0, k, i3，B)JfF(i，i3，tQ，k，B')和 y 代入公式(4)求解 y 在信號 F(i，i3，tQ，k，B')的聯合分布概率L(y|tQ，k，β，Β'); 最后將L(y|tQ，k，β，Β)、L(y|tQ，k，β，Β ')代入公式(7)求解B'的轉移概率α (B' ,B)；
α (B' ，B) = L(y I tQ, k, β，B' ) -L(y 110, k, β , B) (7) `3.3、若α (B' ,B) >0，則將B'賦值給B，其它參量t^k、β保持不變，完成本次更新操作，在由Vk、β和B組成的新的全局近似解(tpk、β、B)中隨機選擇任意參量進行下一次更新操作； ` 3.4若α (B' ,B) <0，則進一步判斷是否到達最大嘗試更新次數Nbi，其中步驟3.1至步驟3.3是一次嘗試更新，若達到最大嘗試更新次數Nbi，則所有參量tpk、β和B保持不變，跳出本次更新操作，重新選擇參量進行下一次更新操作；若未到達最大更新次數Nbi，則跳到步驟3.1繼續本次更新操作。
4.如權利要求1所述的激光雷達全波形信號分析方法，其特征在于，當一次更新操作選定的參量為B時，本次更新操作的過程為: `4.1、根據模擬回火算法的要求，設定函數T = g(T0, u), I < T < Tmax，且T同時滿足公式(8)和(9)； T' > O (8) (Τ, )2_2Τ(Τ, )2+t, / ⑴2 ≤ οO) 公式(8)和(9)中，T'和T''分別為T關于u的一次導函數和二次導函數，初始化此函數為 T = TQ(TQ> I)； `4.2、設Tv = g(T0, V), V是u的范圍內的任意值,設置在Tv = g(T0, V)下的最大嘗試更新次數Nb2 ； `4.3、設B'是B可能更新到的狀態，在區間(0，max(y))選擇一個服從伽瑪分布的隨機數賦值給B'；.4.4、首先將y、、和β代入公式(2)求解依據y、、和β預估的回波信號的能量分布{f(i, β j, t0J)}, j是不大于k的任意非負整數，β j和、是β和h的第j個元素； 然后將{f(i，和B代入公式(3)求解預估的入射到激光雷達探測器的信號F(i,B)，將{f(i，和B'代入公式(3)求解預估的入射到激光雷達探測器的信號 F(i，i3，tQ，k，B'); 再將F(i，β , t0, k, B)和y代入公式(4)求解y在信號F (i，β , t0, k, B)的聯合分布概率L(y|t0, k, i3，B)JfF(i，i3，tQ，k，B')和 y 代入公式(4)求解 y 在信號 F(i，i3，tQ，k，B')的聯合分布概率L(y|tQ，k，β，Β'); 最后將L(y|tQ，k，β，Β)、L(y|tQ，k，β，Β ')代入公式(7)求解B'的轉移概率α (B' ,B)； .4.5、若α (B' ,B) >0，則將B'賦值給B，其它參量t^k和β保持不變，完成本次更新操作，在由Vk、β和B組成的新的全局近似解(tpk、β、B)中隨機選擇任意參量進行下一次更新操作； .4.6、若α (B' ,B) <0，則計算提高轉移概率η (B'，B)，計算方式如公式(10)所示， η (B' ，B)=exp(a(B' , B)/Tv) (10) .4.7、若η (B' , B) > ε , ε為在區間(0，I)的一個滿足均勻分布的隨機數，則將K賦值給B，其它參量h、k、β保持不變，完成本次更新操作，在由tpk、β和B組成的新的全局近似解(tpk、β、B)中隨機選擇任意參量進行下一次更新操作； .4.8、若η (B' , B) ^ ε ,則進一步判斷是否到達在Tv = g(T0, v)條件下最大嘗試更新次數Nb2,其中步驟4.3至步驟4.7是一次嘗試更新,若未到達在Tv = g(T0, V)條件下最大嘗試更新次數NB2，則跳到步驟4.3繼續本次更新操作；若到達在Tv = g(T0, v)條件下最大嘗試更新次數Nb2，則按照公式(11)所示方式調節T:
Tv+1 = g(T0, v+1) (11) 式(11)中，Tv+1為調節后的T ； .4.9、若Tv+1 < Tmax，轉到步驟4.2繼續本次更新操作；若Tv+1 ^ Tmax，則參量h、k、β和B均保持不變，跳出本次更新操作，重新選擇參量進行下一次更新操作。
5.如權利要求1所述的激光雷達全波形信號分析方法，其特征在于，當一次更新操作選定的參量為h時，本次更新操作的過程為: . 5.1、根據模擬回火算法的要求，設定函數T = g(T0, u)，且滿足公式(8) (9);其中，限制u的取值范圍使得I < T < Tmax，T'和Ti ,分別為T關于u的一次導函數和二次導函數，初始化此函數T = T0 (T0 > I)； .5.2、設Tv = g(T0, V), V是u范圍內的任意值,設置在Tv = g(T0, V)下的最大嘗試更新次數義:. 5.3、設k維向量t/是h可能更新到的狀態，在區間(0，ifflax)選擇k個服從均勻分布的隨機數賦值給k維向量t/的所有元素； .5.4、將7、‘和β代入公式(2)求解依據y、、和β預估的回波信號的能量分布{f(i,和β代入公式( 2)求解依據yU/和β預估的回波信號的能量分布{f(i, β j, t0J; )}, j是不大于k的任意非負整數,hj、β」和分別是向量?ρβ和t/的第j個元素；然后將{f(i，和B代入公式(3)求解預估的入射到激光雷達探測器的信號F(i, Itc^k, B),將{f(i, β j, t0J; )}、k和B代入公式(3)求解預估的入射到激光雷達探測器的信號F (i，β，V , k, B)； 再將F(i，β , t0, k, B)和y代入公式(4)求解y在信號F (i，β , t0, k, B)的聯合分布概率L(y|t0, k, i3，B)JfF(i，β，ν，k，B)和 y 代入公式(4)求解 y 在信號 F(i，β，t/，k，B)的聯合分布概率Uyltc/ , k, β , B)； 最后將L (y I tQ, k, β，B)、L (y 11。'，k, β，B)代入公式(12 )求解tQ'的轉移概率α (V，t0)；
6.如權利要求1所述的激光雷達全波形信號分析方法，其特征在于，當一次更新操作選定的參量為β時，本次更新操作的過程為: .6.1、根據模擬回火算法的要求，設定函數T = g(T0, u)，且滿足公式(8) (9);其中，限制u的取值范圍使得I < T < Tmax，T'和Ti ,分別為T關于u的一次導函數和二次導函數，初始化此函數T = T0 (T0 > I)； . 6.2、設Tv = g(T0, V), V是u范圍內任意值,設置在Tv = g(T0, V)下的最大嘗試更新次數\ ； .6.3、設k維向量t/是tQ可能更新到的狀態,在區間(0，max (y))選擇k個服從均勻分布的隨機數賦值給β ,的所有個元素； .6.4、首先將y、、和β代入公式(2)求解依據y、、和β預估的回波信號的能量分布{f(i,β，代入公式(2)求解依據y、、和β，預估的回波信號的能量分布{f(i, β / , t0J)}, j是不大于k的任意非負整數，t0J> β」和β/分別是向量h、β和β ,的第j個元素； 然后將{f(i，i^，tw)}、k和B代入公式(3)求解入射到激光雷達探測器的信號F(i, β，將{f(i，β/，、)}、k和B代入公式(3)求解入射到激光雷達探測器的信號 F(i，i3' , t0, k, B)； 再將F(i, β , t0, k, B)和y分別代入公式(4)求解y在信號F(i, β , t0, k, B)的聯合分布概率Uyltci, k，i3，B)JfF(i，β ' , t0, k, B)和y分別代入公式(4)求解y在信號F(i, β '，tQ, k, B)的聯合分布概率 L(y I tQ, k, β ' ,B); 最后將L (y I tQ, k, β，B)、L (y I tQ, k, β ' , B)代入公式(15 )求解β '的轉移概率α (β ; , β)；
α (β ' , β ) = L (y 110, k, & 1 , B) -L (y 110, k, β , B) (15) .6.5、若α(β'，3)>0，則將3'賦值給β，其它參量h、k和B保持不變，完成本次更新操作，在由h、k、β和B組成的新的全局近似解(h、k、β、B)中隨機選擇任意參量進行下一次更新操作； . 6.6、若0(0'，β≥0，則計算提高的轉移概率η (β '，β)，計算方式如公式(16)所示，
rI ( β ' , β ) = exp ( α ( β ' , β )/Τν) (16) . 6.7、若η (β '，β ) > ε，ε為區間(0，I)的一個滿足均勻分布隨機數，則將β '賦值給β，其它參量Vk和B保持不變，完成本次更新操作，在由νΚβ和B組成的新的全局近似解(tpk、β、B)中隨機選擇任意參量進行下一次更新操作； . 6.8、若η (β '，β)≤ε ,則進一步判斷是否到達在Tv = g (Tci, V)條件下最大嘗試更新.次數Ne,其中步驟6.3至步驟6.7是一次嘗試更新,若未到達在Tv = g(T0, V)條件下最大嘗試更新次數Ne，則轉到步驟6.3繼續本次更新操作；若到達在Tv = g (Ttl, V)條件下最大嘗試更新次數Ne，則按照公式(11)所示方式調節T ； . 6.9、若1；+1 < Tmax，其中Tv+1是調節后的T，轉到步驟6.2繼續本次更新操作；若Tv+1 ^ Tmax，參量tpk、β和B均保持不變，跳出本次更新操作，重新選擇參量進行下一次更新操作。
【文檔編號】G01S7/493GK103926575SQ201410077558
【公開日】2014年7月16日 申請日期:2014年3月5日 優先權日:2014年3月5日
【發明者】何偉基, 尹文也, 王瑋, 陳錢, 顧國華, 張聞文, 錢惟賢, 隋修寶, 任侃, 路東明, 于雪蓮 申請人:南京理工大學