姚海洋1)2) 王海燕1)2)? 張之琛1)2) 申曉紅1)2)
1)(海洋聲學(xué)信息感知工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室(西北工業(yè)大學(xué)),西安 710072)
2)(西北工業(yè)大學(xué)航海學(xué)院,西安 710072)
一種基于廣義Duffing振子的水中弱目標(biāo)檢測(cè)方法?
姚海洋1)2) 王海燕1)2)? 張之琛1)2) 申曉紅1)2)
1)(海洋聲學(xué)信息感知工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室(西北工業(yè)大學(xué)),西安 710072)
2)(西北工業(yè)大學(xué)航海學(xué)院,西安 710072)
(2016年11月24日收到;2017年3月26日收到修改稿)
海洋環(huán)境中,在水下目標(biāo)的線譜頻率未知或者目標(biāo)輻射噪聲的連續(xù)譜很弱時(shí),很難實(shí)現(xiàn)水中弱目標(biāo)的準(zhǔn)確檢測(cè),本文提出基于廣義Duffing振子檢測(cè)系統(tǒng)的水下目標(biāo)輻射噪聲檢測(cè)方法.通過(guò)對(duì)傳統(tǒng)周期擾動(dòng)的Duffing振子信號(hào)檢測(cè)系統(tǒng)的分析和推廣,提出了一種可輸入非周期、非平穩(wěn)信號(hào)的廣義Duffing振子檢測(cè)系統(tǒng),可檢測(cè)輸入的無(wú)先驗(yàn)信息目標(biāo)信號(hào).為實(shí)現(xiàn)廣義Duffing振子系統(tǒng)運(yùn)動(dòng)狀態(tài)的精確、有效判斷,提出了一種相空間圖形的離散分布列計(jì)算方法,通過(guò)類(lèi)網(wǎng)格函數(shù)實(shí)現(xiàn)了利用統(tǒng)計(jì)復(fù)雜度對(duì)系統(tǒng)輸出的嵌入式表征,從而實(shí)現(xiàn)了無(wú)先驗(yàn)信息時(shí)的水中弱目標(biāo)的嵌入式檢測(cè).相同條件下與傳統(tǒng)檢測(cè)方法仿真對(duì)比可知,本文提出的方法可以檢測(cè)到更低信噪比下的目標(biāo),并能滿足水中檢測(cè)實(shí)時(shí)性要求.
水中目標(biāo)檢測(cè),廣義Duffing振子系統(tǒng),統(tǒng)計(jì)復(fù)雜度,相空間
對(duì)水中遠(yuǎn)距離弱目標(biāo)的被動(dòng)檢測(cè)是通過(guò)提取目標(biāo)(指水中機(jī)動(dòng)目標(biāo),如船舶、潛艇、水下航行器等)輻射噪聲中的特征來(lái)實(shí)現(xiàn),目標(biāo)輻射噪聲頻譜具有特殊的線譜和連續(xù)譜,一般通過(guò)提取線譜和連續(xù)譜特征實(shí)現(xiàn)目標(biāo)檢測(cè).線譜是目標(biāo)輻射噪聲譜中的一個(gè)重要特征,李啟虎等[1,2]從理論方法和數(shù)值分析角度分別探討了自相關(guān)檢測(cè)、快速傅里葉變換方法、自適應(yīng)線譜增強(qiáng)等幾種方法的優(yōu)劣,相比之下分段快速傅里葉變換線譜檢測(cè)具有較好的效果,且對(duì)頻率漂移現(xiàn)象有較好的寬容性.Antoni和Hanson[3]提出利用循環(huán)頻率分析方法提取線譜頻率,并進(jìn)行了理論推導(dǎo)分析.周關(guān)林等[4]研究了在海洋混響背景下利用隨機(jī)共振原理檢測(cè)微弱線譜的方法,仿真顯示,經(jīng)隨機(jī)共振處理后,從信號(hào)頻譜圖可以清晰地看到原本非常微弱的線譜信息.基于傳統(tǒng)Duffing振子的一系列目標(biāo)輻射噪聲檢測(cè)方法研究結(jié)果表明,在待測(cè)線譜已知的條件下,混沌振子可以實(shí)現(xiàn)高斯白噪聲背景下信噪比低至?25 dB的微弱線譜檢測(cè)[5].以上方法均要求線譜頻率等參數(shù)已知,無(wú)此先驗(yàn)信息時(shí),結(jié)果會(huì)受到很大影響.叢超等[6]利用間歇Duffing混沌振子實(shí)現(xiàn)了線譜頻率未知的目標(biāo)輻射噪聲微弱線譜的檢測(cè),然而其多系統(tǒng)搜索方式導(dǎo)致資源消耗巨大,難以工程實(shí)現(xiàn),同時(shí)其最終基于圖像變化的判別方式無(wú)法實(shí)現(xiàn)嵌入式表征.
連續(xù)譜是目標(biāo)輻射噪聲頻譜的另一重要特征和組成部分,具有單獨(dú)的譜峰,一般通過(guò)檢測(cè)此譜峰實(shí)現(xiàn)檢測(cè).張曉勇和羅來(lái)源[7]從頻率與能量分布的角度出發(fā),利用瞬時(shí)頻率分析方法,實(shí)現(xiàn)了利用艦船輻射噪聲連續(xù)譜進(jìn)行目標(biāo)檢測(cè).然而,這種基于能量的檢測(cè)方法在遠(yuǎn)距離下結(jié)果會(huì)受到較大影響.同時(shí),海洋環(huán)境噪聲也并非一般情況所假設(shè)的復(fù)高斯隨機(jī)過(guò)程[8,9].
隨著人們對(duì)水下目標(biāo)輻射噪聲的深入研究,發(fā)現(xiàn)目標(biāo)輻射噪聲中含有混沌行為.Hinich等[10]利用雙譜對(duì)海洋聲信號(hào)進(jìn)行特征分析,證實(shí)以目標(biāo)輻射噪聲為主的信號(hào)中含有環(huán)境噪聲所不具備的混沌分量.基于目標(biāo)輻射噪聲信號(hào)中的特殊混沌成分,可進(jìn)行目標(biāo)檢測(cè).
在水下目標(biāo)的線譜頻率未知或者目標(biāo)輻射噪聲的連續(xù)譜很弱時(shí),上述方法很難實(shí)現(xiàn)對(duì)水中弱目標(biāo)的檢測(cè).針對(duì)這種問(wèn)題,本文通過(guò)對(duì)傳統(tǒng)周期擾動(dòng)的Duffing振子信號(hào)檢測(cè)系統(tǒng)進(jìn)行分析和推廣,提出了一種可直接輸入非周期、非平穩(wěn)信號(hào)的廣義Duffing振子檢測(cè)系統(tǒng),目標(biāo)輻射噪聲信號(hào)與其他信號(hào)輸入廣義Duffing振子系統(tǒng)后輸出相空間存在明顯差異.通過(guò)一種針對(duì)圖形的離散分布列計(jì)算方法,實(shí)現(xiàn)了利用統(tǒng)計(jì)復(fù)雜度對(duì)相空間差異的嵌入式表征,進(jìn)而實(shí)現(xiàn)了無(wú)目標(biāo)輻射噪聲線譜頻率等先驗(yàn)信息條件下水中目標(biāo)的低信噪比檢測(cè),分析了算法的實(shí)時(shí)性,得出了系統(tǒng)的最低檢測(cè)信噪比.
2.1 Duffing振子檢測(cè)系統(tǒng)
常規(guī)Duffing振子檢測(cè)系統(tǒng)描述為
式中k為阻尼系數(shù),?αx(t)+βx3(t)為非線性恢復(fù)力,rcos(ωt)為內(nèi)策動(dòng)力,ω為內(nèi)策動(dòng)力頻率,r為內(nèi)策動(dòng)力幅值,x為振動(dòng)幅值,y為振動(dòng)速率.
在穩(wěn)定內(nèi)策動(dòng)力下,取α,β,ω為常數(shù),固定k時(shí),隨著r從0逐漸增加,系統(tǒng)輸出相空間將依次出現(xiàn):同宿軌道周期狀態(tài)、分叉狀態(tài)、混沌狀態(tài)和大尺度周期狀態(tài).
基于Duffing振子檢測(cè)系統(tǒng)的信號(hào)檢測(cè)方法的基本思想是:將待測(cè)信號(hào)作為特定參數(shù)輸入臨界狀態(tài)下的檢測(cè)系統(tǒng),利用系統(tǒng)所擁有的豐富的非線性動(dòng)力學(xué)特點(diǎn),辨識(shí)其運(yùn)動(dòng)狀態(tài),根據(jù)運(yùn)動(dòng)狀態(tài)的變化,最終達(dá)到微弱信號(hào)檢測(cè)的目的[11].
由以上表述可知,Duffing振子檢測(cè)系統(tǒng)會(huì)在穩(wěn)定擾動(dòng)下發(fā)生規(guī)律變化.那么,將待測(cè)信號(hào)以加項(xiàng)形式輸入(1)式表征的系統(tǒng),待測(cè)信號(hào)中的成分構(gòu)成決定了系統(tǒng)輸出相空間的形態(tài),一般的周期信號(hào)(如單頻信號(hào)等)和隨機(jī)信號(hào)(如高斯白噪聲信號(hào)、色噪聲信號(hào))輸入系統(tǒng)總會(huì)圍繞上述四種狀態(tài)運(yùn)動(dòng).前期研究表明,特定的復(fù)雜待測(cè)信號(hào)會(huì)使得輸出相空間脫離以上四種標(biāo)準(zhǔn)形態(tài),本文據(jù)此進(jìn)行目標(biāo)信號(hào)的檢測(cè).
2.2 統(tǒng)計(jì)復(fù)雜度
以輸出相空間形態(tài)變化為依據(jù)進(jìn)行系統(tǒng)運(yùn)動(dòng)狀態(tài)判斷,為實(shí)現(xiàn)運(yùn)動(dòng)狀態(tài)判斷的嵌入式表述,引入以下統(tǒng)計(jì)復(fù)雜度的概念.統(tǒng)計(jì)復(fù)雜度是刻畫(huà)與系統(tǒng)時(shí)間信號(hào)相關(guān)聯(lián)的概率分布的函數(shù),通過(guò)標(biāo)準(zhǔn)Shannon熵和失衡度表示[12].
Shannon熵是經(jīng)典信息論中的測(cè)度,表示概率分布為P={pi,i=1,···,N}的物理過(guò)程的不確定程度:
其中,pi為概率值,N為隨機(jī)變量個(gè)數(shù).
系統(tǒng)處于完全隨機(jī)狀態(tài),即概率分布為均勻分布Pe={1/N,···,1/N}時(shí),Shannon熵的值達(dá)到最大S[Pe]=lnN.由此,標(biāo)準(zhǔn)Shannon熵為
其中Smax=S[Pe]=lnN,顯然,0 6 HS[P]6 1.
失衡度K用于度量系統(tǒng)任一狀態(tài)T時(shí)的概率分布P與均勻分布Pe之間的距離DS,表示為
其中K0是歸一化常數(shù),則0 6 K 6 1.DS通常選用Jensen-Shannon散度JS進(jìn)行刻畫(huà),對(duì)于概率空間中任意兩個(gè)分布P1和P2,其表示為
那么,失衡度為
其中,歸一化常數(shù)K0為JS[P,Pe]取最大值時(shí)的倒數(shù).
由(3)式所示的標(biāo)準(zhǔn)Shannon熵和(6)式所示的失衡度,可得下述統(tǒng)計(jì)復(fù)雜度:
其中P為樣本概率分布.
由于標(biāo)準(zhǔn)Shannon熵反映系統(tǒng)的無(wú)序性,失衡度描述系統(tǒng)時(shí)間信號(hào)在概率空間中的結(jié)構(gòu)特性,故(7)式表示的統(tǒng)計(jì)復(fù)雜度同時(shí)反映了系統(tǒng)的這兩種性質(zhì),其計(jì)算關(guān)鍵是找出信號(hào)的概率分布P.在計(jì)算輸出相空間的統(tǒng)計(jì)復(fù)雜度后,根據(jù)相空間的有界性和差異性,提出一種針對(duì)圖形的離散分布列計(jì)算方法.
由于目標(biāo)輻射噪聲為非周期、非平穩(wěn)信號(hào)[13],而常規(guī)Duffing振子檢測(cè)系統(tǒng)輸入為周期或平穩(wěn)隨機(jī)信號(hào).本文將常規(guī)Duffing振子檢測(cè)系統(tǒng)進(jìn)行推廣,稱非周期、非平穩(wěn)信號(hào)輸入下的Duffing振子檢測(cè)系統(tǒng)為廣義系統(tǒng).
將待測(cè)信號(hào)以加項(xiàng)形式輸入到(1)式表示的檢測(cè)系統(tǒng),描述為
其中g(shù)(t)=Acos(ω′t+ φ)+n′(t),n′(t)為高斯白噪聲信號(hào),A為目標(biāo)信號(hào)幅值,ω′為目標(biāo)信號(hào)頻率,φ為目標(biāo)信號(hào)相位.(8)式表示傳統(tǒng)周期擾動(dòng)下的常規(guī)Duffing振子檢測(cè)系統(tǒng).
設(shè)目標(biāo)輻射噪聲信號(hào)為s(t),海洋環(huán)境噪聲為n(t),檢測(cè)系統(tǒng)的輸入信號(hào)(即待測(cè)信號(hào))可表示為
顯然,(9)式表示非周期、非平穩(wěn)信號(hào).
根據(jù)(8)式,對(duì)此周期信號(hào)擾動(dòng)的檢測(cè)系統(tǒng)進(jìn)行推廣,將(9)式代入以代替g(t),得到
由于(8)式所示系統(tǒng)的檢測(cè)原理為對(duì)周期信號(hào)敏感和對(duì)平穩(wěn)噪聲信號(hào)免疫,當(dāng)輸入非周期、非平穩(wěn)信號(hào)時(shí),其原本具有的一些特殊的性質(zhì)和限制將發(fā)生改變.下面以系統(tǒng)(8)為出發(fā)點(diǎn),分情況討論推廣系統(tǒng)(10)的性質(zhì)和可行性.
情況1當(dāng)f(t)=g(t)時(shí),即(10)式表示的檢測(cè)系統(tǒng)等價(jià)于(8)式表示的系統(tǒng),其特性分析如下.
為了能檢測(cè)到更低幅值的單頻信號(hào),參數(shù)r的取值需盡可能精確.采用定步長(zhǎng)四階龍格-庫(kù)塔方法進(jìn)行數(shù)值計(jì)算,步長(zhǎng)設(shè)定為0.001,取k=0.5,α =1,β =1,ω =1,初值為(1,1). 利用Melnikov方法[14]計(jì)算得知,r=0.863時(shí),(1)式所示的Duffing振子系統(tǒng)處于混沌臨界狀態(tài),分別取不同精度r值,r=0.8,0.86,0.863,得到臨界狀態(tài)系統(tǒng)相空間形態(tài),如圖1所示.
由圖1所示的系統(tǒng)臨界狀態(tài)相空間形態(tài)可以看出,內(nèi)策動(dòng)力幅值r取值越精確,相空間顯示出的混沌形態(tài)越不明顯,越趨于大尺度周期狀態(tài).在r=0.8和0.863兩種參數(shù)時(shí)分別輸入單頻信號(hào)、高斯白噪聲信號(hào)及單頻+高斯白噪聲信號(hào)三種情況下的系統(tǒng)相空間如圖2和圖3所示.由圖3(b)和圖3(c)可以看出,單頻信號(hào)加背景噪聲的輸入系統(tǒng)狀態(tài)判斷非常困難.且由仿真及文獻(xiàn)[14]可得結(jié)論1:Duffing振子檢測(cè)系統(tǒng)對(duì)高斯白噪聲和色噪聲信號(hào)免疫,僅使得相空間曲線變得不光滑.
仿真結(jié)果表明,周期信號(hào)加平穩(wěn)噪聲輸入下的Duffing振子檢測(cè)系統(tǒng)存在以下問(wèn)題:
圖1 臨界狀態(tài)下不同內(nèi)策動(dòng)力幅值相空間 (a)r=0.8;(b)r=0.86;(c)r=0.863Fig.1.Output phase spaces under critical state with di ff erent amplitude of policy dynamic:(a)r=0.8;(b)r=0.86;(c)r=0.863.
圖2 內(nèi)策動(dòng)力幅值r=0.8時(shí)檢測(cè)過(guò)程相空間 (a)輸入單頻信號(hào);(b)輸入高斯白噪聲信號(hào);(c)輸入單頻信號(hào)+高斯白噪聲信號(hào)Fig.2.Output phase spaces in detection when r=0.8:(a)Single-frequency signal;(b)white Gaussian noise signal;(c)single-frequency signal+white Gaussian noise signal.
圖3 參數(shù)r=0.863時(shí)檢測(cè)過(guò)程相空間 (a)輸入單頻信號(hào);(b)輸入高斯白噪聲信號(hào);(c)輸入單頻信號(hào)+高斯白噪聲信號(hào)Fig.3.Output phase spaces in detection when r=0.863:(a)Single-frequency signal;(b)white Gaussian noise signal;(c)single-frequency signal+white Gaussian noise signal.
1)低信噪比下,對(duì)參數(shù)r取值越精確系統(tǒng)輸出相空間差異越不明顯,判斷系統(tǒng)運(yùn)動(dòng)狀態(tài)越困難;
2)系統(tǒng)(8)進(jìn)行單頻信號(hào)檢測(cè)時(shí),需滿足內(nèi)策動(dòng)力rcos(ωt)的頻率ω與待測(cè)信號(hào)g(t)=Acos(ω′t+ φ)+n′(t)的頻率ω′很接近,即
其中ε為常數(shù),|ε|<0.03ω;
3)待測(cè)信號(hào)g(t)需要與內(nèi)策動(dòng)力進(jìn)行相位同步,形成正激勵(lì);
4)低信噪比下,隨著隨機(jī)噪聲的方差變大,通過(guò)輸出相空間清晰地進(jìn)行運(yùn)動(dòng)狀態(tài)判斷,定步長(zhǎng)四階龍格-庫(kù)塔方法的步長(zhǎng)需變小,計(jì)算量大大增加,檢測(cè)系統(tǒng)的實(shí)時(shí)性變差.
以上四個(gè)問(wèn)題為(8)式表示的Duffing振子檢測(cè)系統(tǒng)所固有的困難,下面將會(huì)討論廣義Duffing振子系統(tǒng)是否存在這些問(wèn)題.
情況2當(dāng)f(t)=s(t)+n(t)時(shí),(10)式表示廣義Duffing振子檢測(cè)系統(tǒng),其特性分析如下.
由結(jié)論1可知,高斯白噪聲和色噪聲對(duì)相空間形態(tài)無(wú)影響,于是,令輸入為f(t)=s(t).為方便觀察,取r=0.8,使得系統(tǒng)進(jìn)入混沌臨界狀態(tài),如圖1(a)所示.
混沌信號(hào)為一類(lèi)具有代表性的非周期、非平穩(wěn)信號(hào),對(duì)輸入混沌信號(hào)時(shí)的廣義Duffing振子檢測(cè)系統(tǒng)的性質(zhì)進(jìn)行分析.令s(t)為混沌特性信號(hào),利用奇異吸引子產(chǎn)生混沌特性信號(hào),Lorenz吸引子表示為[15]
取初值(x0,y0,z0)=(0,0,10?10);σ,b,c均為常系數(shù),σ=10,b=?8/3,c=28.利用定步長(zhǎng)四階龍格-庫(kù)塔方法進(jìn)行計(jì)算,取步長(zhǎng)為0.001,得到Lorenz吸引子信號(hào)W=(x,y,z),其中x,y,z均為一維向量.
由于Duffing振子輸入初值不能為零,取sL(t)=z,令f(t)=sL(t).同時(shí),為觀察改變輸入信號(hào)相位時(shí)系統(tǒng)輸出相空間的變化,對(duì)輸入信號(hào)進(jìn)行平移.分別輸入
其中m為遠(yuǎn)小于n的整數(shù).廣義Duffing振子系統(tǒng)輸出相空間形態(tài)如圖4所示.
圖4 輸入Lorenz信號(hào)后系統(tǒng)輸出相空間 (a)相空間1;(b)相空間2Fig.4.Output phase spaces when the input is Lorenz signal:(a)Phase space one;(b)phase space two.
輸入不同相位的Lorenz信號(hào),廣義Duffing振子系統(tǒng)輸出相空間具有兩種形態(tài):一種是相空間內(nèi)部向焦點(diǎn)運(yùn)動(dòng)軌跡發(fā)生變化,如圖4(a)所示;另一種是相空間左右鞍點(diǎn)軌道數(shù)比重發(fā)生變化,如圖4(b)所示.兩種相空間形態(tài)無(wú)固定規(guī)律,系統(tǒng)運(yùn)動(dòng)狀態(tài)沒(méi)有發(fā)生變化,仍為混沌狀態(tài).同時(shí),將2倍、1.5倍、1/2、1/10的原信號(hào)能量的信號(hào)輸入,得到的相空間依然如圖4所示.于是得到結(jié)論2:Lorenz吸引子產(chǎn)生的混沌特性信號(hào)輸入廣義Du ffing振子系統(tǒng)時(shí),相空間形態(tài)變化,系統(tǒng)運(yùn)動(dòng)狀態(tài)不發(fā)生改變,且不受能量變化影響.
推廣結(jié)論2到任意混沌信號(hào),分別考察與Lorenz系統(tǒng)相空間運(yùn)動(dòng)方式不同的Henon吸引子[16]和與Lorenz系統(tǒng)不拓?fù)涞葍r(jià)的Chen吸引子[17]產(chǎn)生的信號(hào)輸入廣義Duffing振子系統(tǒng)時(shí)輸出相空間的變化.仿真顯示,Henon吸引子信號(hào)和Chen吸引子信號(hào)輸入廣義Duffing系統(tǒng)后,輸出相空間與如圖4所示的Lorenz吸引子信號(hào)的輸出相空間形態(tài)相似,僅左右鞍點(diǎn)軌道占比不同.
故而,非平穩(wěn)的混沌特性信號(hào)輸入廣義Du ffing振子系統(tǒng)可以得到有效的相空間輸出,且可能會(huì)使得左右鞍點(diǎn)軌道運(yùn)動(dòng)發(fā)生變化,系統(tǒng)運(yùn)動(dòng)狀態(tài)不變.
情況3當(dāng)f(t)=s(t)+n(t)時(shí),(10)式表示廣義Duffing振子檢測(cè)系統(tǒng),為了解更加復(fù)雜的非周期、非平穩(wěn)信號(hào)輸入系統(tǒng)后的運(yùn)動(dòng)狀態(tài)變化,設(shè)計(jì)以下四種輸入形式:
其中sL(t),sH(t),sC(t)分別表示Lorenz吸引子信號(hào)、Henon吸引子信號(hào)和Chen吸引子信號(hào);λ,a,b,c為常數(shù),0 6 λ 6 1,a+b+c=1.此種情況下輸出相空間形態(tài)依然如圖4所示.
情況4當(dāng)f(t)=s(t)+n(t)時(shí),(10)式表示廣義Duffing振子檢測(cè)系統(tǒng),綜合分析如下.
若輸入信號(hào)同時(shí)含有單頻信號(hào)、混沌特性信號(hào)和隨機(jī)噪聲信號(hào),即
其 中s1(t)= Acos(ωt);s2(t)為 形 式 如(13)—(16)式的混沌特性信號(hào);s3(t)為隨機(jī)噪聲信號(hào);,,是值為0或1的常數(shù), 且,不同時(shí)為零.此種情況下輸出相空間形態(tài)依然如圖4所示,若存在s3(t)則會(huì)使得相空間變得粗糙.
因此,非周期、非平穩(wěn)信號(hào)輸入廣義Duffing振子檢測(cè)系統(tǒng)可以得到有效相空間輸出,并且對(duì)臨界狀態(tài)下的系統(tǒng)運(yùn)動(dòng)狀態(tài)沒(méi)有產(chǎn)生如單頻信號(hào)的擾動(dòng).
綜合以上分析可知,改變輸入信號(hào)組成模式有可能使得相空間形態(tài)發(fā)生一些改變,但一般的加性非周期、非平穩(wěn)信號(hào)不會(huì)對(duì)運(yùn)動(dòng)狀態(tài)產(chǎn)生影響.試驗(yàn)證明,目標(biāo)輻射噪聲輸入系統(tǒng)時(shí),會(huì)產(chǎn)生如圖5所示的特殊輸出相空間,依據(jù)如圖1(a)、圖4所示相空間與如圖5所示相空間差異可實(shí)現(xiàn)目標(biāo)檢測(cè).
圖5 特殊輸出相空間 (a)相空間1;(b)相空間2Fig.5.Special output phase spaces:(a)Phase space one;(b)phase space two.
非周期、非平穩(wěn)信號(hào)輸入廣義Duffing振子系統(tǒng)時(shí),情況1中所述常規(guī)系統(tǒng)(8)存在的四種問(wèn)題均可進(jìn)行推廣.
1)相空間形態(tài)發(fā)生巨大變化,不再是密閉的橢圓形態(tài),更易判斷系統(tǒng)運(yùn)動(dòng)狀態(tài),且由于不涉及單頻擾動(dòng),檢測(cè)精度與參數(shù)r設(shè)置無(wú)關(guān),而取決于系統(tǒng)對(duì)目標(biāo)信號(hào)復(fù)雜度的敏感度,可設(shè)臨界參數(shù)為r=0.8.
2)輸入信號(hào)改變,廣義Duffing振子系統(tǒng)相空間變化主要依賴于輸入信號(hào)中的復(fù)雜的非線性成分,故而不存在頻率接近的問(wèn)題.
3)輸入信號(hào)發(fā)生變化,無(wú)需使得輸入信號(hào)對(duì)策動(dòng)力產(chǎn)生正激勵(lì),故而無(wú)需進(jìn)行相位同步.
4)檢測(cè)系統(tǒng)的實(shí)時(shí)性將在下面一節(jié)給出.
綜上所述,廣義Duffing振子系統(tǒng)是可行的,并且可以實(shí)現(xiàn)對(duì)水下目標(biāo)輻射噪聲的有效檢測(cè).
為實(shí)現(xiàn)嵌入式檢測(cè),需借助2.2節(jié)中所述的統(tǒng)計(jì)復(fù)雜度,計(jì)算系統(tǒng)輸出的統(tǒng)計(jì)復(fù)雜度需要得到輸出時(shí)間序列的概率分布.從已知時(shí)間信號(hào)中提取概率分布的一種有效方法為Bandt-Pompe算法[18],然而B(niǎo)andt-Pompe算法針對(duì)時(shí)間序列進(jìn)行重排列時(shí),由于嵌入維數(shù)DM是確定的,無(wú)法靈活調(diào)整計(jì)算精度,一些情況下過(guò)于精確的計(jì)算會(huì)使得實(shí)時(shí)性降低,影響檢測(cè)性能.
本文直接針對(duì)輸出相空間提出一種基于類(lèi)網(wǎng)格函數(shù)的相圖離散分布列計(jì)算方法,可達(dá)到相空間離散分布列提取的目的,并可以通過(guò)控制類(lèi)網(wǎng)格函數(shù)進(jìn)行精度控制.
通過(guò)搜索檢測(cè)系統(tǒng)輸出(x,y),分別得到x的最大值xmax與最小值xmin,y的最大值ymax與最小值ymin,則
其中xscale為相空間x軸方向上的相圖范圍,yscale為相空間y軸方向上的相圖范圍.
根據(jù)(18)和(19)式可計(jì)算輸出相空間的網(wǎng)格步長(zhǎng):
其中hx為x軸上的步長(zhǎng),hy為y軸上的步長(zhǎng),n為精度控制參數(shù)即x軸劃分網(wǎng)格數(shù),DM為重構(gòu)維數(shù)(Duffing振子重構(gòu)維數(shù)為2).
由(20)和(21)式所示步長(zhǎng),分別以x=xmin,y=ymin為初值,可得到以下分段區(qū)間:
其中i=0,···,n ? 1;j=0,···,DM? 1.
本文對(duì)類(lèi)網(wǎng)格函數(shù)定義為
根據(jù)本文研究的非周期、非平穩(wěn)信號(hào)的相空間特性,建立基于類(lèi)網(wǎng)格函數(shù)的輸出相空間圖形的分布列計(jì)算公式:
其中pl為時(shí)間序列的分布列,#表示{}中對(duì)應(yīng)數(shù)對(duì)的數(shù)量,N為數(shù)對(duì)(x,y)的總個(gè)數(shù),l=1,···,n × DM.
依據(jù)概率統(tǒng)計(jì)理論,顯然,該算法需滿足n×DM?N,利用Grassberger-Procaccia(G-P)算法[19]得到嵌入維數(shù)DM.由于引入了嵌入維數(shù)DM,本方法也適用于其他類(lèi)型的混沌系統(tǒng)輸出分布列計(jì)算.
本方法利用類(lèi)網(wǎng)格函數(shù),推導(dǎo)出了輸出相空間圖形的分布列計(jì)算公式,可計(jì)算出系統(tǒng)輸出相空間圖形的離散分布列,其中n的取值決定了分布的準(zhǔn)確性,取值越大,劃分的網(wǎng)格越多,越接近真實(shí)分布.進(jìn)而可以計(jì)算出Shannon熵和統(tǒng)計(jì)復(fù)雜度,對(duì)系統(tǒng)輸出相空間進(jìn)行量化描述,實(shí)現(xiàn)了目標(biāo)的嵌入式檢測(cè).
5.1 檢測(cè)方法海上實(shí)驗(yàn)分析
為了驗(yàn)證以上檢測(cè)方法的可行性,以四種艦船作為目標(biāo)(記為1—4型目標(biāo)),1型和2型為以不同航速航行的300噸級(jí)艦船,3型和4型為同航速3000噸級(jí)和萬(wàn)噸級(jí)艦船,聲吶收集信號(hào)時(shí)的采樣頻率fs為48000 Hz,降采樣獲取100—8000 Hz的目標(biāo)輻射噪聲,樣本長(zhǎng)度N=5000,1型艦船輻射噪聲時(shí)域波形如圖6所示.將信號(hào)歸一化,記為sJ(t),分別取
即將信號(hào)進(jìn)行連續(xù)三點(diǎn)平移(取n=5000,即截取t=n/fs=0.104 s長(zhǎng)度的信號(hào)),得到的平移后的信號(hào)分別輸入(10)式所示的廣義Duffing振子檢測(cè)系統(tǒng),得到目標(biāo)存在時(shí)廣義Duffing振子檢測(cè)系統(tǒng)輸出相空間形態(tài)如圖7所示(僅顯示1型輻射噪聲輸出形態(tài)).
圖6 1型艦船輻射噪聲時(shí)域形態(tài)Fig.6.Time domain waveform of ship-radiated signal type one.
圖7 1型艦船輻射噪聲輸出相空間形態(tài) (a)相空間1;(b)相空間2;(c)相空間3Fig.7. Output phase spaces when the input is ship-radiated signal type one:(a)Phase space one;(b)phase space two;(c)phase space three.
其他三種艦船的輸出相空間也均為如圖7所示的特殊形態(tài),只是兩種特殊形態(tài)出現(xiàn)的次數(shù)和排列順序不同.
由以上廣義Duffing振子系統(tǒng)輸出相空間可知,連續(xù)相移的艦船輻射噪聲信號(hào)輸入后,所有輸出相空間均呈特殊形態(tài).為增加檢測(cè)準(zhǔn)確性,本文采用連續(xù)三點(diǎn)相移輸入的策略進(jìn)行目標(biāo)檢測(cè).
一定范圍(方圓5 km)內(nèi)無(wú)目標(biāo)經(jīng)過(guò)時(shí),采集到如圖8所示的未經(jīng)歸一化的海洋環(huán)境噪聲時(shí)域波形.
圖8 海洋環(huán)境噪聲時(shí)域波形Fig.8.Time domain waveform of ocean ambient noise.
海洋環(huán)境噪聲輸入后的系統(tǒng)輸出相空間如圖9所示.
由圖9可知,此段海洋環(huán)境噪聲信號(hào)輸入廣義Duffing振子檢測(cè)系統(tǒng)之后,不改變系統(tǒng)運(yùn)動(dòng)狀態(tài),僅對(duì)左右鞍點(diǎn)運(yùn)動(dòng)偏重有影響.同時(shí),在后續(xù)仿真中,設(shè)計(jì)如下試驗(yàn):
試驗(yàn)1輸入2倍、1.5倍、1/2、1/10的原海洋環(huán)境噪聲能量的信號(hào);
試驗(yàn)2取原海洋環(huán)境噪聲信號(hào)中不同時(shí)段部分分別輸入;
試驗(yàn)3對(duì)其他海域3—5級(jí)海況下收集到的環(huán)境噪聲數(shù)據(jù)分別做如試驗(yàn)1和2的處理.
上述三種試驗(yàn)均得到類(lèi)似圖9所示的相空間形態(tài),只是左右鞍點(diǎn)運(yùn)動(dòng)高度略有差別.
從以上結(jié)果可知,利用這種相空間形態(tài)的區(qū)別進(jìn)行海洋環(huán)境噪聲背景下艦船輻射噪聲的檢測(cè)具有一定的普適性.
根據(jù)2.2節(jié)所述過(guò)程及第4節(jié)的算法,可得到如圖1(a)所示的混沌臨界狀態(tài)和如圖2(a)所示的大尺度周期狀態(tài)相空間形態(tài)的統(tǒng)計(jì)復(fù)雜度.由G-P算法[19]得到Duffing振子重構(gòu)維數(shù)為2,即DM=2,取精度參數(shù)n=5,即將相空間分割為2×5的小塊(有限個(gè)數(shù)分割的前提下,均為離散型隨機(jī)變量,分割塊數(shù)越多,對(duì)應(yīng)的分布越精確,分割塊數(shù)為無(wú)窮時(shí),得到概率密度曲線).進(jìn)一步得到圖1(a)和圖2(a)對(duì)應(yīng)的分布列.根據(jù)表1和表2所列分布列,由2.2節(jié)可得混沌臨界狀態(tài)統(tǒng)計(jì)復(fù)雜度為0.0169,大尺度周期狀態(tài)統(tǒng)計(jì)復(fù)雜度為0.0349.
圖9 海洋環(huán)境噪聲輸入時(shí)系統(tǒng)輸出相空間 (a)相空間1;(b)相空間2Fig.9.Output phase spaces when the input is ocean ambient noise:(a)Phase space one;(b)phase space two.
表1 混沌臨界狀態(tài)相空間分布列Table 1.Phase space distribution sequence in the chaotic critical state.
表2 大尺度周期狀態(tài)相空間分布列Table 2.Phase space distribution sequence in the large-scale periodic state.
表3 特殊相空間1分布列Table 3.Phase space distribution sequence of special phase space type one.
表4 特殊相空間2分布列Table 4.Phase space distribution sequence of special phase space type two.
同理可得出圖7所示的兩種特殊形態(tài)相空間的分布列.根據(jù)表3和表4所列分布列,可得兩種分布的統(tǒng)計(jì)復(fù)雜度分別為0.1907和0.2800.
由以上計(jì)算可知,特殊相空間統(tǒng)計(jì)復(fù)雜度值與常規(guī)相空間值差別較大,將閾值定為0.05,得到輸入不同信號(hào)時(shí)的檢測(cè)效果示意,如圖10所示.
圖10 目標(biāo)檢測(cè)效果Fig.10.Results of target detection.
圖10中,虛線為閾值大小,圈為檢測(cè)值即統(tǒng)計(jì)復(fù)雜度值.圖中1,2檢測(cè)點(diǎn)系統(tǒng)處于臨界狀態(tài)為檢測(cè)準(zhǔn)備時(shí)刻;3—5檢測(cè)點(diǎn)為第一個(gè)檢測(cè)時(shí)段,第3檢測(cè)點(diǎn)存在單頻干擾信號(hào)進(jìn)入系統(tǒng),使得系統(tǒng)進(jìn)入了大尺度周期狀態(tài)且于第5檢測(cè)點(diǎn)干擾消失(此處假設(shè)單頻信號(hào)內(nèi)策動(dòng)力同步等理想情況,否則依然為混沌狀態(tài));6,7檢測(cè)點(diǎn)為檢測(cè)間隔時(shí)段(同11,12與16,17),可調(diào)節(jié);8—10檢測(cè)點(diǎn)為第二個(gè)檢測(cè)時(shí)段,檢測(cè)結(jié)果為存在目標(biāo);13—15時(shí)刻為第三個(gè)檢測(cè)時(shí)段,檢測(cè)結(jié)果為存在目標(biāo);18—20時(shí)刻為第四個(gè)檢測(cè)時(shí)段,檢測(cè)結(jié)果為無(wú)目標(biāo).
由1—4型目標(biāo)輻射噪聲特殊輸出相空間可以看出,不同目標(biāo)輻射噪聲兩種特殊輸出相空間排列次序不同、出現(xiàn)次數(shù)不同,由于兩種特殊相空間統(tǒng)計(jì)復(fù)雜度值不一樣,導(dǎo)致檢測(cè)點(diǎn)起伏變化不同,通過(guò)增加輸入信號(hào)平移次數(shù),突出這種區(qū)別,可實(shí)現(xiàn)目標(biāo)識(shí)別.
5.2 檢測(cè)性能分析
根據(jù)以下信噪比的計(jì)算公式,可得利用廣義Duffing振子檢測(cè)目標(biāo)輻射噪聲的最低檢測(cè)信噪比:
其中Psignal為信號(hào)的平均功率,Pnoise為噪聲的平均功率.當(dāng)艦船輻射噪聲信號(hào)能量過(guò)低時(shí),將無(wú)法對(duì)廣義Duffing振子系統(tǒng)形成擾動(dòng),產(chǎn)生如圖7所示的特殊相空間,按(25)式可得到最低檢測(cè)信噪比為?9.133 dB.考慮到Duffing振子檢測(cè)系統(tǒng)對(duì)高斯白噪聲的免疫性,若將海洋環(huán)境噪聲假定為高斯白噪聲,則基于廣義Duffing振子系統(tǒng)的檢測(cè)方法可得到最低檢測(cè)信噪比為?27.55 dB.由于本文提出的方法并非以統(tǒng)計(jì)理論為基礎(chǔ)得出的,故而無(wú)法得出在不同虛警概率下的檢測(cè)概率曲線.
利用能量方法,依據(jù)奈曼-皮爾遜準(zhǔn)則,分別取虛警概率Pf為0.1,0.01,0.001,0.0001時(shí),得到其檢測(cè)概率與信噪比關(guān)系如圖11所示.由計(jì)算可知,能量檢測(cè)在虛警概率0.1、檢測(cè)概率PD>98%時(shí),檢測(cè)信噪比為?11 dB左右.
在高斯白噪聲背景下基于廣義Duffing振子系統(tǒng)的最低檢測(cè)信噪比(?27.55 dB)遠(yuǎn)低于能量檢測(cè)的信噪比(?11 dB).然而,由于四階龍格-庫(kù)塔方法的計(jì)算方式等的影響,廣義Duffing振子檢測(cè)方法并非很好地適用于背景噪聲為高斯白噪聲的情況,而更適用于上述海洋環(huán)境噪聲下的目標(biāo)檢測(cè).
圖11 (網(wǎng)刊彩色)能量檢測(cè)方法檢測(cè)概率曲線Fig.11.(color online)Detection probability curves of energy detection method.
需要說(shuō)明的是,聲納接收到的目標(biāo)輻射噪聲是寬帶信號(hào),且含有海洋環(huán)境噪聲成分,故而實(shí)際最低檢測(cè)信噪比通過(guò)實(shí)驗(yàn)計(jì)算得到的還要低一些.
下面分析廣義Duffing振子檢測(cè)方法的實(shí)時(shí)性.上述方法主要涉及的計(jì)算包括歸一化、四階龍格-庫(kù)塔方法解廣義Duffing振子方程、相圖離散分布列方法和統(tǒng)計(jì)復(fù)雜度的計(jì)算.其中,歸一化和統(tǒng)計(jì)復(fù)雜度涉及簡(jiǎn)單的乘法和除法,相圖離散分布列方法主要為分段、搜索和除法,而龍格-庫(kù)塔方法比較復(fù)雜,占據(jù)了主要的運(yùn)算時(shí)間.將實(shí)時(shí)性以運(yùn)算時(shí)間衡量,分別輸入樣本點(diǎn)數(shù)為3000—5000時(shí),運(yùn)算量如表5所列.
計(jì)算機(jī)配置不同、數(shù)據(jù)不同等會(huì)得到不同的運(yùn)算時(shí)間,因而僅利用表5觀察實(shí)時(shí)性趨勢(shì).首先,由表5可以看出總體耗時(shí)均在10 s以下,對(duì)于水中目標(biāo)航行速度不可能高的實(shí)際情況,算法耗時(shí)基本達(dá)到水中檢測(cè)的實(shí)時(shí)性要求.其次,輸入信號(hào)點(diǎn)數(shù)對(duì)耗時(shí)沒(méi)有很大影響,主要取決于龍格-庫(kù)塔方法步長(zhǎng)的選擇,考慮到相空間的完整性,選取步長(zhǎng)為0.001.最后,表5也驗(yàn)證了廣義Duffing振子方法的計(jì)算量主要集中于龍格-庫(kù)塔方法.
表5 廣義Duffing振子檢測(cè)方法實(shí)時(shí)性Table 5.Real-time analysis of generalized Duffing oscillation detection method.
本文提出了一種基于廣義Duffing振子檢測(cè)系統(tǒng)的真實(shí)海洋環(huán)境噪聲背景下目標(biāo)輻射噪聲檢測(cè)方法,利用相圖離散分布列方法計(jì)算統(tǒng)計(jì)復(fù)雜度對(duì)相空間進(jìn)行表征,實(shí)現(xiàn)了線譜頻率未知或者目標(biāo)輻射噪聲的連續(xù)譜很弱時(shí)的目標(biāo)檢測(cè),給出了檢測(cè)系統(tǒng)的最低檢測(cè)信噪比,算法的實(shí)時(shí)性可滿足實(shí)際要求.與常規(guī)能量方法對(duì)比可知,廣義Duffing振子檢測(cè)方法可以檢測(cè)到更低信噪比下的目標(biāo)輻射噪聲信號(hào).本文方法的特點(diǎn)如下:
1)通過(guò)對(duì)傳統(tǒng)周期擾動(dòng)的Duffing振子信號(hào)檢測(cè)系統(tǒng)進(jìn)行分析和推廣,提出一種針對(duì)非周期、非平穩(wěn)輸入的廣義Duffing振子檢測(cè)系統(tǒng),可將無(wú)先驗(yàn)信息的復(fù)雜接收信號(hào)直接輸入系統(tǒng);
2)發(fā)現(xiàn)了一種目標(biāo)輻射噪聲輸入廣義Duffing振子系統(tǒng)后的特殊形態(tài)相空間,并分析了不同信號(hào)輸入后的相空間差異,通過(guò)表征特殊相空間和一般相空間差異可實(shí)現(xiàn)目標(biāo)檢測(cè)和識(shí)別;
3)提出了一種相空間圖形的離散分布列計(jì)算方法,利用類(lèi)網(wǎng)格函數(shù)推導(dǎo)出了相空間分布列計(jì)算公式,得到了系統(tǒng)輸出相空間的分布列,實(shí)現(xiàn)了利用統(tǒng)計(jì)復(fù)雜度對(duì)相空間的嵌入式表征,從而實(shí)現(xiàn)了水下真實(shí)環(huán)境下線譜頻率未知或者接收的目標(biāo)輻射信號(hào)的連續(xù)譜很弱時(shí)目標(biāo)的嵌入式檢測(cè).
[1]Li Q H,Li M,Yang X T 2008 Acta Acoust.33 193(in Chinese)[李啟虎,李敏,楊秀庭 2008聲學(xué)學(xué)報(bào) 33 193]
[2]Li Q H,Li M,Yang X T 2008 Acta Acoust.33 289(in Chinese)[李啟虎,李敏,楊秀庭 2008聲學(xué)學(xué)報(bào) 33 289]
[3]Antoni J,Hanson D 2012 IEEE J.Ocean.Eng.37 478
[4]Zhou G L,Li G H,Cheng J 2009 Comput.Simulat.7 337(in Chinese)[周關(guān)林,李鋼虎,成靜 2009計(jì)算機(jī)仿真7 337]
[5]Liu H B,Wu D W,Jin W,Wang Y Q 2013 Acta Phys.Sin.62 050501(in Chinese)[劉海波,吳德偉,金偉,王永慶2013物理學(xué)報(bào)62 050501]
[6]Cong C,Li X K,Song Y 2014 Acta Phys.Sin.63 064301(in Chinese)[叢超,李秀坤,宋揚(yáng) 2014物理學(xué)報(bào)63 064301]
[7]Zhang X Y,Luo L Y 2015 Acta Acoust.40 511(in Chinese)[張曉勇,羅來(lái)源 2015聲學(xué)學(xué)報(bào) 40 511]
[8]Kedar S 2011 Comptes Rendus Geosci.343 548
[9]Cato D H 2012 AIP Conf.Proc.1495 242
[10]Hinich M J,Marandino D,Sullivan E J 1989 J.Acoust.Soc.Am.85 1512
[11]Lai Z H,Leng Y G,Sun J Q,Fan S B 2012 Acta Phys.Sin.61 050503(in Chinese)[賴志慧,冷永剛,孫建橋,范勝波2012物理學(xué)報(bào)61 050503]
[12]He M J,Xu W,Sun Z K 2014 Sci.Sin.:Phys.Mech.Astron.44 981(in Chinese)[何美娟,徐偉,孫中奎 2014中國(guó)科學(xué):物理學(xué)力學(xué)天文學(xué)44 981]
[13]Bao F,Li C,Wang X,Wang Q,Du S 2010 J.Acoust.Soc.Am.128 206
[14]Li Y,Shi Y W,Ma H T,Yang B J 2004 Acta Electron.Sin.32 87(in Chinese)[李月,石要武,馬海濤,楊寶俊2004電子學(xué)報(bào)32 87]
[15]Lorenz E N 1995 The Essence of Chaos(Seattle:University of Washington Press)p186
[16]Chen G,Lai D 1998 Int.J.Bifurcat.Chaos 8 1585
[17]Xue Y J,Yang S Y 2003 Chaos Soliton.Fract.17 717
[18]Bandt C,Pompe B 2002 Phys.Rev.Lett.88 174102
[19]Grassberger P,Procaccia I 2004 The Theory of Chaotic Attractors(New York:Springer)pp189–208
PACS:43.60.–c,05.45.–aDOI:10.7498/aps.66.124302
A method of detecting underwater weak target based on generalized Duffing oscillator?
Yao Hai-Yang1)2)Wang Hai-Yan1)2)?Zhang Zhi-Chen1)2)Shen Xiao-Hong1)2)
1)(Key Laboratory of Ocean Acoustics and Sensing(Northwestern Polytechnical University),Ministry of Industry and Information Technology,Xi’an 710072,China)
2)(School of Marine Science and Technology,Northwestern Polytechnical University,Xi’an 710072,China)
24 November 2016;revised manuscript
26 March 2017)
In the marine environment,when the line spectra of underwater target radiated signal are unknown or the continuous spectra are weak,it is extremely hard to accurately detect the underwater weak target.The line spectrum based method commonly requires spectrum information for detection,and the continuous spectrum based energy method is hard to achieve accurate detection in long distance.In this paper,an underwater radiated noise detection method based on generalized Duffing oscillator detection system is proposed.Firstly,a generalized Duffing oscillator detection system for non-periodic and non-stationary input signal is proposed through deducing the traditional Duffing oscillator detection system that is perturbed by periodic signal.And the proposed generalized Duffing oscillator detection system is able to detect the signals of targets without needing prior information.Secondly,when the target radiated signal(nonperiodic and non-stationary signal)is input into the generalized Duffing oscillator,a special form of output phase space(a special state of motion)is discovered and the di ff erences in output phase space among di ff erent input signals(periodic stationary signals,nonperiodic non-stationary signals and the target radiated signals)are analyzed.It is found that the special phase space has di ff erent form from the output phase spaces of other kinds of signals;accordingly the underwater targets can be detected through the representation of the di ff erence between special phase space and ordinary phase space.Thirdly,a discrete distribution sequence calculation method based on phase space is proposed for the precise and efficient judgment of system motion.The proposed calculation method de fi nes a similar-grid function,based on which,the distribution sequence calculation method of output phase space is deduced,and the distribution sequences of di ff erent kinds of output phase spaces are calculated.The method realizes an embedded expression of system output by using the statistical complexity,therefore achieving the embedded underwater target detection when the line spectra of underwater target radiated signal are unknown or the continuous spectra are weak.The analysis result indicates that the method is of low-computation.Finally,the experimental results in the sea are described and the lowest signalto-noise ratio(SNR)of the method is calculated to be?9.133 dB.Simulation and experimental results have shown that the proposed method can detect target successfully in a lower SNR than traditional detection method,and the real-time performance can meet the demand for underwater detection.The method in this paper provides new ideas and ways of thinking for underwater target detection,and has very important reference value for low SNR long-distance target detection under real condition.
underwater target detection,generalized Duffing oscillator system,statistical complexity,phase space
10.7498/aps.66.124302
?國(guó)家自然科學(xué)基金(批準(zhǔn)號(hào):61571365)資助的課題.
?通信作者.E-mail:hyang@mail.nwpu.edu.cn
?2017中國(guó)物理學(xué)會(huì)Chinese Physical Society
http://wulixb.iphy.ac.cn
*Project supported by the National Natural Science Foundation of China(Grant No.61571365).
?Corresponding author.E-mail:hyang@mail.nwpu.edu.cn