基于Bayesian估計(jì)的X射線脈沖星微弱信號檢測*

張 華 許錄平 謝 強(qiáng) 羅 楠

(西安電子科技大學(xué)電子工程學(xué)院，西安 710071)(2010年4月8日收到;2010年7月10日收到修改稿)

基于Bayesian估計(jì)的X射線脈沖星微弱信號檢測*

張 華 許錄平 謝 強(qiáng) 羅 楠

(西安電子科技大學(xué)電子工程學(xué)院，西安 710071)(2010年4月8日收到;2010年7月10日收到修改稿)

累積輪廓、流量和周期是X射線脈沖星輻射信號的三個(gè)重要特征，將其應(yīng)用于X射線脈沖星信號檢測中，提出了一種基于Bayesian估計(jì)的X射線脈沖星周期輻射信號時(shí)域檢測方法.該方法以非脈沖區(qū)噪聲觀測為先驗(yàn)知識，利用X射線脈沖星輻射信號的泊松分布模型推導(dǎo)了信號概率密度分布函數(shù)，以該函數(shù)的累積分布函數(shù)為判據(jù)，對X射線脈沖星微弱信號進(jìn)行檢測，并提取位相偏移量.利用仿真數(shù)據(jù)和RXTE衛(wèi)星的實(shí)測數(shù)據(jù)進(jìn)行實(shí)驗(yàn)驗(yàn)證，結(jié)果表明:本文方法性能優(yōu)于同類的基于高斯分布模型的檢測方法，在檢測信號的同時(shí)能在一定精度下給出信號位相偏移值.

脈沖星，Bayesian估計(jì)，位相測量，時(shí)域檢測

PACS:97.60.Gb，98.70.Qy，96.60.tk

1.引 言

脈沖星是一類高速穩(wěn)定自旋的中子星，其輻射信號覆蓋了射電、紅外、可見光、紫外、X射線及γ射線波段［1］.其中X射線波段集中了輻射信號的大部分能量，易于小型化設(shè)備探測和處理，使基于X射線脈沖星的航天器定時(shí)、定姿和定位成為可能［2］.脈沖星距太陽系通常都有幾千光年，空間X射線探測器探測到的X射線脈沖星輻射信號非常微弱且湮沒在背景噪聲中，因此X射線脈沖星信號檢測就成為應(yīng)用中的一大難題.常用的檢測方法以基于FFT的頻域周期檢測法為主［3］，有時(shí)還需要將脈沖星信號進(jìn)行周期累加得到脈沖星累積脈沖輪廓模型，再提取平移不變特征量進(jìn)行檢測［4，5］.但基于FFT的頻域檢測方法僅在高斯白噪聲前提下適用，而且當(dāng)脈沖星的輻射脈沖較窄時(shí)，信號能量分散于諧波中，若噪聲包含頻率成分與諧波接近，則諧波能量得到增強(qiáng)，這將導(dǎo)致信號被誤檢測，比如，利用該方法對ROSAT望遠(yuǎn)鏡觀測的光變曲線數(shù)據(jù)進(jìn)行周期檢測就存在失效情況［6］;此外，這種頻域方法計(jì)算量大，這就限制了其在計(jì)算能力有限的星載計(jì)算機(jī)中的實(shí)現(xiàn).相對于這種頻域方法而言，時(shí)域方法可以避免頻域變換時(shí)能量分散于諧波中，因此對信息利用更為充分，可以得到更優(yōu)的檢測性能［6］.文獻(xiàn)［7］最早提出了利用統(tǒng)計(jì)推理方法對感興趣區(qū)域脈沖信號進(jìn)行檢測的時(shí)域方法，但其前提是信號和噪聲呈高斯分布.本文對這一時(shí)域方法進(jìn)行拓展，根據(jù)X射線脈沖星輻射信號的泊松分布模型，提出一種新的基于 Bayesian估計(jì)［8］的 X射線脈沖星輻射信號時(shí)域檢測方法，討論了其在位相測量中的作用，并通過實(shí)驗(yàn)說明了該方法的有效性.

2.基于Bayesian估計(jì)的X射線脈沖星輻射信號檢測

2.1.脈沖星輻射信號劃分及其建模

脈沖星的脈沖周期和輪廓特征在X射線探測器探測的原始信號上表現(xiàn)為X射線光子流的平均強(qiáng)度變化，為了同時(shí)利用觀測數(shù)據(jù)所包含的脈沖星輻射信號的輪廓信息、強(qiáng)度信息和周期信息，根據(jù)累積輪廓的結(jié)構(gòu)，將脈沖星信號分為信號段和噪聲段.以噪聲段的噪聲觀測作為先驗(yàn)知識，利用Bayesian估計(jì)推導(dǎo)出信號段信號強(qiáng)度的后驗(yàn)概率分布.不失一般性，以脈沖星B1937+21為例，如圖1，該脈沖星一段觀測信號累積出的脈沖輪廓包含一個(gè)主脈沖和一個(gè)次脈沖，組成信號段，其在一個(gè)脈沖周期內(nèi)持續(xù)的時(shí)間分別設(shè)為 WP1，WP2，令 WP=WP1+WP2;其他部分組成噪聲段，在一個(gè)脈沖周期內(nèi)持續(xù)的時(shí)間分別設(shè)為 WN1，WN2和 WN3，令 WN=WN1+WN2+WN3;WP+WN=T.

圖1 脈沖星信號劃分

空間X射線光子到達(dá)探測器具有隨機(jī)性，若將單位時(shí)間內(nèi)到達(dá)探測器的光子數(shù)視為隨機(jī)事件，可用泊松分布對該事件進(jìn)行建模［9，10］.由于 X射線脈沖星輻射信號強(qiáng)度周期性變化，光子到達(dá)探測器的密度隨時(shí)間分布是不均勻的，定義落在區(qū)間(t1，t2)內(nèi)的光子數(shù)等于k的概率為［9］

式中，λ(t)表示 t時(shí)刻光子流量密度，λ(t)=λs(t)+λn，其中λs(t)為信號輻射強(qiáng)度，當(dāng)脈沖周期為T時(shí)，有λs(t)=λs(t+T)，λn為各種噪聲強(qiáng)度的總和，包括背景輻射噪聲、系統(tǒng)熱噪聲等，簡單起見視噪聲平均強(qiáng)度為常數(shù);t2－t1=ΔT為觀測時(shí)間.

2.2.X射線脈沖星信號后驗(yàn)概率表示

為了簡化分析，假設(shè):1)背景輻射噪聲平穩(wěn)，信號段的背景噪聲分布與噪聲段一致，且可用噪聲段觀測數(shù)據(jù)建模;2)一次分析時(shí)間為一個(gè)脈沖周期，在整個(gè)觀測期間假設(shè)脈沖周期不變;3)在背景輻射噪聲強(qiáng)度未知時(shí)假設(shè)其服從均勻分布.定義:Cg，Cb和Cp分別為信號段光子數(shù)、噪聲段光子數(shù)和一個(gè)脈沖周期時(shí)間觀測的光子總數(shù)，則Cp=Cb+Cg，Cg=定義:μp，μb和μg分別為信號光子平均計(jì)數(shù)率、背景光子平均計(jì)數(shù)率和信號段總的光子計(jì)數(shù)率，有μg=μp+μb.對X射線脈沖星信號分析的目的即在已知Cb和Cg的前提下，估計(jì)后驗(yàn)概率 P(μp/Cb，Cg).由查普曼-科爾莫格洛夫方程有

在平均背景輻射強(qiáng)度已知時(shí)，即μb已知時(shí)，對于泊松分布模型而言，可認(rèn)為 P(μp/μb)=P(μg= μp+μb)，得到

其中 μg僅與 Cg有關(guān)，有

由X射線脈沖星輻射的泊松分布模型可得

根 據(jù) Bayesian 原 理 有 P(μg/Cg) =，假設(shè) P(Cg)服從均勻分布，有E［P(μg)］WP=E［P(Cg)］，近似得到

在假設(shè)1)前提下，用噪聲段的觀測估計(jì)μb，

與(6)式同理得到

將(3)，(4)，(6)，(7)和(8)式代入(2)式得到

其中μg=μp+μb，根據(jù)二項(xiàng)式定理有(μgWP)Cg=，代入(9)式得到

再根據(jù) Gamma積分［11］

對(10)式積分運(yùn)算得到

(12)式中 P(μp/Cb，Cg)表示在已知觀測量Cb，Cg時(shí)該觀測中存在信號的可能性，當(dāng)觀測到的信號段總計(jì)數(shù)Cg趨近無窮大時(shí)，P(μp/Cb，Cg)趨近于1.

2.3.基于后驗(yàn)概率累積分布函數(shù)的X射線脈沖星信號檢測算法

通過對(12)式的后驗(yàn)概率密度函數(shù)積分得到其累積分布函數(shù)轉(zhuǎn)4);

4)若 j≤k，s(t)=s(t+jΔt)轉(zhuǎn)2)，否則轉(zhuǎn)5);

5)求 R(n)(0≤n≤k)的方差δ，搜索是否存在R(n)>3δ，若存在說明檢測到信號.

2.4.基于信號檢出率的位相偏移測量算法

在圖1所示的脈沖星信號劃分中，噪聲段 WN

從2.1節(jié)和2.2節(jié)的分析知道，Cb，Cg分別反映了噪聲和含噪聲的信號的先驗(yàn)知識，(13)式中 PCb，Cg則反映了從先驗(yàn)知識Cb，Cg中獲取的信號的信息.先驗(yàn)知識越準(zhǔn)確，Cg中包含的信號信息越多，PCb，Cg越大，當(dāng)PCb，Cg大于或等于某一置信水平時(shí)，就認(rèn)為探測到了信號.對于X射線探測器接收到的信號而言，認(rèn)為其是以采樣時(shí)間為單位的光子計(jì)數(shù)率的時(shí)間序列，用 s(nΔt)表示該時(shí)間序列，s(n1Δt，n2Δt)表示序列中的某一段，T表示脈沖星周期，Δt表示采樣間隔，T=kΔt，k 為每周期采樣數(shù)，以 PCb，Cg為判據(jù)，檢測序列中是否包含特定脈沖星信號的步驟為

1)對某顆脈沖星，按圖1所示的方法將其輪廓劃分為信號段WP和噪聲段WN，設(shè)閾值水平為L，令i=0，j=0，HS=0，HB=0，待檢測序列長度為 mT(m為整周期數(shù));

2)對 s(ikΔt，(i+1)kΔt)計(jì)算 PCb，Cg，若 PCb，Cg≥ L，HS=HS+1，否則 HB=HB+1;

3)若 i

式中s∈{0，k－1}表示起始位相，P表示觀測整周期數(shù).具體實(shí)現(xiàn)方法為在2.3節(jié)的步驟5)中提取r=max

n{R(n)}.由于時(shí)域采樣的離散性，時(shí)域位相提取的精度受最小采樣間隔的限制.Taylor FFT等頻域方法［5］可以擺脫這種限制，使位相測量精度僅與累積輪廓信噪比有關(guān)，但頻域方法中諧波分量會導(dǎo)致多解情況，造成真實(shí)解判定的混亂.因此，位相時(shí)域解最少可以作為頻域解搜索過程中的初值，為更精確的位相提取提供有力的參考.

3.脈沖星特征參數(shù)對檢測性能影響的理論分析

為了分析各參數(shù)對新方法檢測性能影響，以(13)式所表示的概率分布函數(shù) PCb，Cg為分析對象，分析結(jié)果如圖2所示:Cg，Cb和 β分別表示每脈沖周期信號段光子計(jì)數(shù)、噪聲段光子計(jì)數(shù)和累積輪廓占空比;圖 2(a)以 Cg=0，1，…，5 ph(光子)和 Cb=3 ph為例，說明 PCb，Cg隨 β 變化情況;圖 2(b)以Cb=0，1，…，5 ph 和 Cg=3ph 為例，說明 PCb，Cg隨 β變化情況;圖 2(c)以 β =0.5，Cg=0，1，…，5 ph 為例，說明 PCb，Cg隨 Cb變化情況;圖 2(d)以 β =0.5，Cb=0，1，…，5 ph 為例，說明 PCb，Cg隨 Cg變化情況.

圖2 概率分布函數(shù)與各參數(shù)的關(guān)系 (a)Cb=3 ph，PCb，Cg與 β 關(guān)系;(b)Cg=3ph，PCb，Cg與 β 關(guān)系;(c)PCb，Cg與 Cb關(guān)系;(d)PCb，Cg與 Cg關(guān)系

利用累積分布 PCb，Cg作為判據(jù)，對于圖2(a)而言，β 一定時(shí)，PCb，Cg在 Cg=0，1，…，5 ph 時(shí)的各曲線之間的差值越大，說明檢測靈敏度越高，檢測性能越好;通過觀察該差值隨β分布可以發(fā)現(xiàn)，其僅在β∈［0.2，0.8］時(shí)較大，說明本文方法不適用于占空比過大或過小的脈沖星信號檢測，而且該差值在β∈［0.2，0.5］時(shí)的值大于 β∈［0.5，0.8］的值，說明本文方法對信號段較窄的脈沖星檢測性能更好些;對于圖 2(b)而言，β 一定時(shí) PCb，Cg在 Cb=0，1，…，5 ph的各曲線之間的差值越大，說明檢測靈敏度受噪聲影響越大，與圖2(a)的分析方法相同，可以發(fā)現(xiàn)本文方法檢測信號段較窄的脈沖星信號時(shí)對噪聲的抑制能力強(qiáng).從圖2(a)和圖2(b)中還發(fā)現(xiàn)，當(dāng) β一定時(shí)，Cg越大檢測性能越好;相反，Cb越大檢測性能越差，這與一般物理規(guī)律是符合的.圖2(c)和圖2(d)進(jìn)一步說明了這一結(jié)論.如圖2(c)所示，隨著 Cb增加，PCb，Cg減小，Cg=0，1，…，5 ph 的 PCb，Cg各曲線間的差也在減小，說明檢測靈敏度降低;反之，如圖2(d)所示，對Cb一定的信號而言，隨著Cg增強(qiáng)，PCb，Cg增大，則檢測性能得到增強(qiáng)，當(dāng) Cb越小時(shí)，檢測性能隨Cg增強(qiáng)越快.

4.實(shí)驗(yàn)與分析

4.1.仿真實(shí)驗(yàn)及分析

以輻射強(qiáng)度較弱的 B1937+21和 B1821－24兩顆脈沖星為例，主要參數(shù)如表1.設(shè)X射線探測器有效探測面積為1 m2，平均背景輻射強(qiáng)度為5×102ph/cm2.s－1，根 據(jù) 脈 沖 星 周 期、標(biāo) 準(zhǔn) 累 積 輪廓［12—14］、平均流量和每周期采樣次數(shù)計(jì)算每兩次采樣之間信號流量δ(t)和背景輻射流量α，單周期仿真信號表示為:nΔt))+ α}，其中 S為采樣點(diǎn)數(shù)，poissrnd(·)為泊松隨機(jī)數(shù)生成函數(shù).設(shè)探測置信水平為95%(虛警率水平5%)，做1000次Monte Carlo實(shí)驗(yàn)，分別統(tǒng)計(jì)兩顆脈沖星信號檢出率和虛警率與累積周期數(shù)關(guān)系，結(jié)果如圖3(a).可以看到，隨著累積周期數(shù)的增加，檢測性能逐漸提高，且對 B1821-24檢測性能提高更快.需要說明的是，由于周期累加過程中，噪聲也增強(qiáng)了，根據(jù)圖2(c)的結(jié)果，為增加靈敏度，累加后的信號應(yīng)消除部分累加噪聲，這里減去最小值并乘以系數(shù)0.02.為了說明輪廓、流量和周期參數(shù)差異對檢出率的影響，將B1821-24的各個(gè)參數(shù)分別修改為與B1937+21一致，再進(jìn)行仿真，結(jié)果如圖3(b).從圖3(b)分析，當(dāng)脈沖星形狀、信號流量和周期改變時(shí)，本文方法在檢測性能上存在差異，對于脈沖星B1821－24而言，輪廓占空比影響比較大，但這僅是個(gè)例，不同的脈沖星受各參數(shù)的影響在不同條件下是不同的，圖2的分析可以說明這一點(diǎn).

為了評價(jià)本文方法檢測性能，與同類的基于高斯分布的檢測方法對比.基于高斯分布的檢測方法假設(shè)信號和背景噪聲服從高斯分布，其噪聲抽取方法為［7］

圖3 檢出率與累積周期的關(guān)系 (a)檢出率比較;(b)參數(shù)差異影響

此時(shí)判決閾值定義為

其中，δ0為信號方差，K為系數(shù)，這里令 K=2，當(dāng)μgauss≤Lc，表示檢測到信號，否則表示出現(xiàn)虛警.設(shè)累積周期為10000，以B1937+21為例，信號檢出率與虛警率的差和背景輻射強(qiáng)度的關(guān)系如圖4所示.當(dāng)背景噪聲較小時(shí)，基于高斯分布檢測方法的性能與本文方法的性能相當(dāng)，但隨著背景計(jì)數(shù)率的增加，基于高斯分布檢測方法性能迅速下降，本文方法性能優(yōu)于基于高斯分布的檢測方法.

圖4 與高斯檢測方法性能的比較

對于觀測時(shí)間內(nèi)可能出現(xiàn)多個(gè)信號的情況，由于新方法引入了周期和輪廓特征信息，對多個(gè)信號間的干擾具有一定的免疫能力.如圖6(a)所示，設(shè)信號和背景噪聲輻射強(qiáng)度均為 5×10－2ph/cm2·s－1，在僅包含 B1821-24信號的一段觀測中滑動檢測B1937-21信號，然后疊加 B1937-21信號采用同樣的方法檢測，包含B1937-21信號可以觀察到明顯的尖峰，說明檢測到信號.同理，如圖6(b)所示，在B1937-21信號中疊加 B1821-24信號，得到類似的檢測結(jié)果.

圖5 位相測量的均方差

表1 X射線脈沖星參數(shù):周期、X射線流量密度、每周期采樣次數(shù)、占空比

本文方法利用了X射線脈沖星的輪廓信息，可以提取出觀測信號的位相，不過所提取位相的精度受到時(shí)域采樣間隔的限制，精度上不如頻域方法［5］，但頻域方法受諧波分量的影響，不可避免出現(xiàn)多解情況，因此時(shí)域位相解至少可以作為頻域方法初值，提高搜索真實(shí)位相解的速度.對于新方法在位相測量的性能，本文利用B1821-24和B1937+21的仿真數(shù)據(jù)做了仿真分析.設(shè)背景噪聲輻射強(qiáng)度均為 5 ×10－2ph/cm2·s－1，仿真中先隨機(jī)移動仿真數(shù)據(jù)的初始位相，然后利用本文方法滑動檢測初始位相，記錄下真實(shí)位相與檢測位相的差，多次測量求均方差，該均方差與累積周期數(shù)的關(guān)系如圖5所示.可見隨著疊加次數(shù)的增加，整體上位相測量方差在減小，在信噪比相同時(shí)，B1937+21位相測量均方差較小，這與B1937+21采樣率高的原因有直接關(guān)系.

圖6 混合多脈沖星信號情況下的檢測 (a)檢測B1937-21;(b)檢測B1821-24

4.2.基于RXTE實(shí)測數(shù)據(jù)的實(shí)驗(yàn)及分析

從HEASARC數(shù)據(jù)庫的Public Data Access連接進(jìn)入，以B1855+09為關(guān)鍵字搜索出以該脈沖星為觀測對象的觀測號，然后從其FTP［15］上下載與該觀測號對應(yīng)的數(shù)據(jù)，利用FTOOLS工具處理得到該脈沖星X射線光變曲線.經(jīng)過修正，剔除掉不良數(shù)據(jù)，得到256個(gè)周期長的信號，對該信號經(jīng)累加去均值歸一化后的輪廓如圖7所示.由于該段數(shù)據(jù)背景計(jì)數(shù)率較大，為了應(yīng)用本文方法，先逐周期疊加，再減去最小值，乘以系數(shù) 0.02，求 PCb，Cg，然后將信號循環(huán)右移1rad位相，用同樣的方法再計(jì)算PCb，Cg，結(jié)果如圖8所示，可見:位相無誤時(shí)，隨著累加周期數(shù)的增加，累積分布PCb，Cg逐漸提高，而當(dāng)位相偏移1rad時(shí)，PCb，Cg隨著累加周期數(shù)的變化并不明顯，說明了本方法在檢測信號和辨別位相時(shí)的有效性.為了進(jìn)一步說明本文方法在位相測量上作用，將原信號位相右移πrad，然后從零位相開始滑動檢測，PCb，Cg隨滑動位相變化情況如圖9所示，可見位相對 PCb，Cg有直接影響，準(zhǔn)確的位相能最大化 PCb，Cg，從而通過PCb，Cg可以提取出位相.

圖7 歸一化累積輪廓

圖8 分布與累積周期數(shù)的關(guān)系

圖9 PCb，Cg隨位相變化情況

5.結(jié) 論

基于Bayesian估計(jì)的X射線脈沖星信號檢測方法利用X射線脈沖星信號輪廓特征、周期特征和強(qiáng)度等信息作為信號檢測的基本參數(shù)，在時(shí)域?qū)射線脈沖星信號檢測，避免了頻域方法中諧波分量引起的能量分散導(dǎo)致的誤檢測，在檢測性能上，優(yōu)于同類的基于高斯分布模型的檢測方法.利用RXTE的實(shí)測數(shù)據(jù)也說明了該方法的有效性.在位相測量方面，該方法能在一定精度要求下給出位相偏移值，雖然受時(shí)域采樣精度的限制，性能上不如Taylor FFT等頻域方法，但至少能作為頻域方法的初值，減少頻域方法多解情況引起的誤判，提高頻域方法位相搜索的速度.

［1］ Pei Y J，Jun Z 2008 Chin.Phys.B 17 1

［2］ Sheikh S I 2005 Ph.D.Dissertation(Maryland:University of Maryland)

［3］ Seward F D，Harnden F R，Helfand D J 1984 Astrophysical Journal 287 L19

［4］ Xie Z H，Xu L P，Ni G R 2007 J.Infrared Millim.Waves 26 187(in Chinese)［謝振華、許錄平、倪廣仁 2008紅外與毫米波學(xué)報(bào) 26 187］

［5］ Xie Z H，Xu L P，Ni G R 2007 Acta Phys.Sin.26 187(in Chinese)［謝振華、許錄平、倪廣仁 2007物理學(xué)報(bào) 26 187］

［6］ Gregory P C 1996 Astrophysical Journal 10 1

［7］ Currie L A 1968 Analytical Chemistry 40 586

［8］ Gregory P 2003 Bayesian Logical Data Analysis for the Physical Sciences(U.K.:Cambridge Univ)p5

［9］ Sala J，Urruela A，Villares X 2004 European Space Agency Advanced Concepts Team ARIADNA Study June 4202 23

［10］ Li J X，Ke X Z 2010 Acta Phys.Sin.59 8304(in Chinese)［李建勛、柯熙政 2010物理學(xué)報(bào) 59 8304］

［11］ Young J M K A 2009 IEEE Trans.Nucl.Sci.56 1278

［12］ Epn.http://www.jb.man.ac.uk/research/pulsar/resources/epn/browser.html.

［13］ Kuiper L， Hermsen W 2003 proceedings of AGILE Science Workshop Rome，Italy，Dec，2003 p11

［14］ Nicastro L，Cusumano G，Lohmer O，Kramer M，Kuiper L，Hermsen W，Mineo T，Becker W 2004 Astron.Astrophys 413 1065

［15］ Rxte.FTP:legacy.gsfc.nasa.gov.

PACS:97.60.Gb，98.70.Qy，96.60.tk

X-ray pulsar weak signal detection based on Bayesian estimation*

Zhang Hua Xu Lu-Ping Xie Qiang Luo Nan
(School of Electronic Engineering，Xidian University Xi’an 710071，China)(Received 8 April 2010;revised manuscript received 10 July 2010)

Integrated profile，flux and period are three major characteristics of pulsar.Taking advantage of these characteristics，a time domain detecting method is presented based on Bayesian estimation to detect pulsar periodic radiation signals.The signal probability distribution is deduced based on the Poisson distribution model of X-ray pulsar，with noise observation of non-pulse region used as the priori knowledge.The cumulative distribution function of the signal probability is used as a criterion to detect weak signals of X-ray pulsars and extract the phase offset.The RXTE data and the simulation data are used for experimental verification，and the results show that the method outperforms a similar method which is based on the Gauss distribution model，in addition，it can give the phase offset in a certain accuracy.

pulsars，Bayesian estimation，phase measurement，time domain detection

*國家高技術(shù)研究發(fā)展計(jì)劃(批準(zhǔn)號:2007AA12Z323)，國家自然科學(xué)基金(批準(zhǔn)號:60772139)資助的課題.

E-mail:zhanghua@mail.xidian.edu.cn

*Project supported by the National High Technology Research and Development Program of China(Grant No.2007AA12Z323)，the National Natural Science Foundation of China(Grant No.60772139).

E-mail:zhanghua@mail.xidian.edu.cn