穩(wěn)健收斂的時差頻差定位技術(shù)

房嘉奇馮大政 李 進(jìn)

(西安電子科技大學(xué)雷達(dá)信號處理國家重點實驗室 西安 710071)

穩(wěn)健收斂的時差頻差定位技術(shù)

房嘉奇*馮大政 李 進(jìn)

(西安電子科技大學(xué)雷達(dá)信號處理國家重點實驗室 西安 710071)

為實現(xiàn)對目標(biāo)位置和速度的精確定位，該文提出一種基于正則化理論的時差頻差定位技術(shù)。該算法首先利用最大似然方法確定目標(biāo)函數(shù)，然后通過傳統(tǒng)牛頓法對目標(biāo)位置和速度進(jìn)行迭代求解。眾所周知傳統(tǒng)牛頓法對初始值要求較高，較差初始值會導(dǎo)致Hess矩陣趨于病態(tài)，從而致使迭代發(fā)散，該文引入正則化理論修正Hess矩陣，使其更加合理，保證算法穩(wěn)健收斂。實驗結(jié)果表明：相對于傳統(tǒng)牛頓法，該文算法在初始值的選取上具有穩(wěn)健性，對誤差選取較大的初始值，仍能夠保證算法的收斂性；相對于現(xiàn)有閉合式定位方法，該文算法在噪聲較大時具有較好的定位精度，定位精度接近于Cramer-Rao界，具有廣泛的實用價值。

無源定位；到達(dá)時差；到達(dá)頻差；傳統(tǒng)牛頓法；Hess矩陣；正則化算法

1 引言

對未知目標(biāo)的無源定位問題，由于其在導(dǎo)航，跟蹤測軌，無線傳感器網(wǎng)絡(luò)，移動通信等方面的廣泛應(yīng)用，一直都是信號處理領(lǐng)域的研究熱點。目前常用的無源定位技術(shù)可以分為到達(dá)角度(Angle Of Arrival, AOA)、到達(dá)時間(Time Of Arrival, TOA)、到達(dá)時間差(Time-Difference-Of-Arrival, TDOA)、到達(dá)頻差(Frequency-Difference-Of-Arrival, FDOA)等方法。對于固定的目標(biāo)輻射源來說，時差(TDOA)定位技術(shù)利用目標(biāo)到兩個基站的到達(dá)時差信息確定一條以兩基站為焦點的雙曲線，通過多個基站的TDOA測量值確定多條雙曲線，其交點為目標(biāo)位置。時差定位技術(shù)由于不需要目標(biāo)與基站的合作關(guān)系且成本較低，因而在現(xiàn)實中得到廣泛的應(yīng)用。在現(xiàn)代定位問題中，目標(biāo)輻射源時常處于運動狀態(tài)，而接收機(jī)也多裝載于汽車，飛機(jī)，衛(wèi)星等移動平臺中。因此，需要在時差定位的基礎(chǔ)上引入到達(dá)頻差(FDOA)技術(shù)來聯(lián)合估計目標(biāo)的位置以及速度。

估計目標(biāo)的位置和速度并不簡單，因為其需要求解的時差，頻差等式是非線性的。最大似然估計能夠?qū)Ψ蔷€性問題能夠提供漸進(jìn)無偏解，但它需要大量的搜索，在實際中需要耗費大量的計算，難以實現(xiàn)。一個普遍的做法就是通過泰勒展開式(Taylor-Series, TS)將時差頻差非線性等式線性化[1]，這種方法要求初始值必須具備良好的精度才能夠保證迭代收斂?；谶@個缺點，一些閉合式解法被提出用于求解非線性的時差頻差方程[2?13]，其中廣泛使用的兩步WLS(Weighted Least Squares)算法[2?5]引入關(guān)于目標(biāo)位置的附加參數(shù)，采用二次最小二乘算法給出了定位方程組的非迭代解析解，特點是計算量小，在噪聲較小且服從高斯分布的環(huán)境下定位精度高，但當(dāng)噪聲測量誤差較大時，該算法的性能會顯著下降。近年來興起的SDP(Semi-Definite Programming)技術(shù)[14,15]利用凸優(yōu)化方法對目標(biāo)初始位置進(jìn)行求解，然后通過迭代確定目標(biāo)位置，該方法計算量較大，如果初始值不夠精確，目標(biāo)位置可能陷入局部最小解中。

在噪聲較大的情況下，由于時差頻差等式的非線性特點，閉合式算法和SDP算法所求得的目標(biāo)位置，速度會存在較大誤差，而運用傳統(tǒng)牛頓法進(jìn)行迭代求解時，又會因為較差的初始值而導(dǎo)致迭代發(fā)散。本文在傳統(tǒng)牛頓法的基礎(chǔ)上，提出了一種基于正則化理論的時差頻差定位技術(shù)。由時差頻差測量信息可知，求解目標(biāo)位置不需要用到目標(biāo)的速度信息，而求解目標(biāo)的速度則需要用到目標(biāo)的位置信息，針對這一特點，本文利用最大似然估計建立目標(biāo)函數(shù)，首先運用傳統(tǒng)牛頓法對目標(biāo)位置進(jìn)行迭代求解，引入正則化理論修正迭代過程中因為初始值差導(dǎo)致的病態(tài)Hess矩陣，從而使算法在初始值的選取上具有穩(wěn)健性，保證迭代收斂。然后再次運用傳統(tǒng)牛頓法對目標(biāo)位置和速度進(jìn)行聯(lián)解，求得目標(biāo)的速度的同時進(jìn)一步精確目標(biāo)的位置。本文算法有以下優(yōu)點：(1)對目標(biāo)位置，速度初始值的選取具有強(qiáng)大的穩(wěn)健性，在初始值較差的情況下仍能準(zhǔn)確的搜索到目標(biāo)的位置，速度。(2)在測量噪聲較大時相對于閉合式算法具有更高的定位精度，定位誤差接近Cramer-Rao界。

2 系統(tǒng)模型

在多站無源時差頻差定位系統(tǒng)下，假設(shè)待測目標(biāo)輻射源位置x=[x,y,z]T，速度，第i個已知接收基站的位置si=[xi,yi,zi]T，速度=,i=1,2,…,M ,M為參與定位基站數(shù)目，在3維空間內(nèi)至少需要4個基站來產(chǎn)生3組TDOA和FDOA值來確定目標(biāo)的位置和速度。為了獲得目標(biāo)位置和速度的唯一解，接收基站不能設(shè)置在一個點

其中ir為目標(biāo)和第i個接收基站之間的距離?；蛘咭粭l直線上。不考慮非視距傳播的影響，以第1個接收基站為中心基站，則第i個基站和第1個基站之間的TDOA測量值為

c為電波傳播速度，ti1為第i個基站和第1個基站之間的TDOA測量值，di1為與其對應(yīng)的測量距離差(Range-Difference-Of-Arrival, RDOA),ni1為TDOA測量值中的噪聲，文中用RDOA來代替TDOA，因為其在本質(zhì)上只是相差了一個常數(shù)乘子c。對式(1)求導(dǎo)可得到第i個基站和第1個基站之間FDOA測量值。

測得的頻差需要在式(3)中乘以常數(shù)fc/c, fc表示目標(biāo)載波頻率，實際中忽略fc/c不影響定位精度。在通常情況下，檢測到得TDOA, FDOA測量值分別服從零均值，方差為,的高斯分布，并且二者之間相互獨立。將TDOA, FDOA測量等式方程組采用向量形式表達(dá)為

其中d=[d21,d31,…,dM1]T,f=[d.21,d.31,…,d.M1]T, r =[r1,r2,…,rM]T,r.=[r.1,r.2,…,r.M]T,G為第1列均為-1的1×(M?1)的列向量與維數(shù)為M-1的單位矩陣所合并的(M?1)×M的矩陣，nt=c·[n21,n31,…, nM1]T ,。利用最大似然估計準(zhǔn)則，令θ=[xT,x.T]T，目標(biāo)函數(shù)可寫為

其中Qt,Qf為TDOA, FDOA測量值的協(xié)方差矩陣。由式(1)和式(3)可知目標(biāo)位置信息僅由多個TDOA測量值獲得，與測量的頻差信息無關(guān)，而求解目標(biāo)的速度則需要利用目標(biāo)的位置信息，因此在求解式(7)之前，先根據(jù)式(8)：

對目標(biāo)位置進(jìn)行求解，將結(jié)果作為目標(biāo)位置的初始值，再通過式(7)對目標(biāo)速度進(jìn)行求解，同時再次完善目標(biāo)位置。本文算法這么做的目的在于：通過分步計算解決了式(7)中的位置，速度雙重變量的初始值選取難題(位置或速度任何一個變量的初始值較差都可能導(dǎo)致迭代發(fā)散)，使算法對位置，速度初始值的選取都具有強(qiáng)大的穩(wěn)定性，保證迭代收斂。在每一次迭代中通過求解目標(biāo)函數(shù)的局部最小二乘解(LS)來改進(jìn)估計值，直到搜索到正確解上。每次迭代的改變量為

每次迭代后新的迭代點xk+1可表示為

3 正則化算法

式(7)和式(8)中的梯度向量和Hess矩陣?f1(θ)，?2f1(θ)和?f2(x),?2f2(x)可表示為

本文利用傳統(tǒng)牛頓法對目標(biāo)位置，速度信息進(jìn)行求解，牛頓法是一種需要初始估計值的迭代算法，

由式(2)，式(4)可以得到

本文算法先利用式(8)求解目標(biāo)位置，然后將求得的目標(biāo)位置代入式(7)求解目標(biāo)位置和速度。因為目標(biāo)的速度信息與目標(biāo)的位置信息有很大關(guān)系，而一般情況下FDOA噪聲功率又比TDOA噪聲功率小的多，因此如果式(8)求得的目標(biāo)位置信息較為精確，給出的速度初始值即使誤差較大，式(7)仍然能夠通過牛頓法快速準(zhǔn)確地搜索出目標(biāo)的位置及速度。 因此，本文將關(guān)注的重點集中在對目標(biāo)位置的求解上，只要能夠?qū)δ繕?biāo)位置進(jìn)行精確估計，目標(biāo)的速度也能夠隨之被精確算出。

在求解式(8)時，傳統(tǒng)的牛頓算法必須要求目標(biāo)位置具有良好的初始值估計，當(dāng)初始值和目標(biāo)真實位置相距較近時，Hess矩陣為正定矩陣，其迭代方向為下降方向，因此傳統(tǒng)牛頓法能夠準(zhǔn)確地搜索到得到目標(biāo)的真實位置。然而，當(dāng)目標(biāo)初始值和目標(biāo)真值相距較遠(yuǎn)時，Hess矩陣趨于病態(tài)，從而導(dǎo)致其逆矩陣不準(zhǔn)確，致使算法迭代發(fā)散。導(dǎo)致Hess矩陣病態(tài)的原因有兩點：(1)Hess矩陣的特征值逐漸降低且趨近于零。(2)矩陣的條件數(shù)，即最大特征值和最小特征值之間的比率較大。本文通過公式量化分析Hess矩陣病態(tài)原因并且提出改進(jìn)方法，根據(jù)式(13)和式(14)令A(yù)=?2f2(x), b =??f2(x)， 因為Hess矩陣為對稱矩陣，對矩陣A做SVD分解可得

其中U為單位正交矩陣UTU=In, Σ=diag(σ1, σ2,…,σn),σi和向量ui為特征值和與其相對應(yīng)的特征向量，由式(9)本文得到迭代的改變量：

因為式(22)中小特征值σi所對應(yīng)的傅里葉系數(shù)|b|的下降速度并不如這些特征值的下降速度快，因此每次迭代的改變量Δx受小特征值所對應(yīng)項的影響較大。作為結(jié)果，Δx會表現(xiàn)出許多符號改變以及隨機(jī)特征，為了保證迭代中Δx的正確性，就需要減低或者濾出Δx中小特征值所對應(yīng)的那部分輸出的貢獻(xiàn)。本文引用正則化理論知識[16]，對病態(tài)的Hess矩陣進(jìn)行修正，將式(22)的問題轉(zhuǎn)化為

令式(23)的代價函數(shù)最小來求得修正解Δx，正則化的方法有很多種，本文選擇常用的對角加載DSVD方法來修正Hess矩陣，則式(22)可修正為

修正后的矩陣變的更加合理，矩陣條件數(shù)降低，其逆矩陣不再出現(xiàn)符號變化以及隨機(jī)特征， 式(23)和式(24)中λ為重要的正則化參數(shù)，用來控制式(23)中前項(數(shù)據(jù)能量)與后項(穩(wěn)定能量)之間的比例。所有的正則化方法都需要確定合理的參數(shù)λ保證算法快速準(zhǔn)確地收斂，本文引入著名的L曲線(L-curve)理論[17]來確定λ，該方法將作為變量取對數(shù)畫圖，在圖中尋找最顯著的拐點，根據(jù)該拐點，求出其對應(yīng)的正則化參數(shù)λ。實驗結(jié)果驗證了L-curve理論在本文算法中的有效性。本文算法步驟為：

步驟 1 給定初始值x(1),允許誤差Δ＞0,置k=1, l=1。

步驟 2 由式(13)，式(14)計算?f2(x(k))和?2f2(x(k))。

步驟 4 置k:=k+1，轉(zhuǎn)步驟2。

步驟 5 利用步驟3的輸出x作為本次迭代的位置初始值，令θ=[xT,x.T]T，由式(11)，式(12)計算?f1(θ(l)) 和?2f1(θ(l))。

步驟 7 置l:=l+1，轉(zhuǎn)步驟5。

4 仿真實驗

為了驗證本文方法的有效性，實驗中采用5個基站來估計目標(biāo)位置和速度，設(shè)置基站1為中心基站，接收基站位置及速度見表1，近距離目標(biāo)和遠(yuǎn)距離目標(biāo)位置分別為(280, 325, 375) m或(2800, 3250, 2750) m，目標(biāo)速度均設(shè)為(-20, 15, 40) m/s。目標(biāo)位置和速度的均方誤差定義為RMSE(x)=。實驗中假設(shè)TDOA, FDOA測量噪聲是不相關(guān)的，其協(xié)方差矩陣可表示為σ2R和0.1σ2R,σ2為時差測量中噪聲的方差，R為對角元素為1其余元素為0.5的矩陣。實驗結(jié)果為5000次蒙特卡洛實驗的平均數(shù)據(jù)。

表1 接收基站位置及速度

圖1，圖2給出了本文算法在好初始值和差初始值兩種情況下同傳統(tǒng)牛頓法的比較結(jié)果，本組實驗中固定σ2=1。圖1中目標(biāo)位置，速度初始值分別給定為(300, 300, 400) m和(50, 50, 50) m/s，圖2中目標(biāo)位置，速度初始值分別給定為(1000, 1000, 800) m和(2000, 2000, 2000) m/s，圖2實驗特別選擇一個很差的初始值來驗證本文算法的穩(wěn)健性。從圖1中可以看出，當(dāng)目標(biāo)位置，速度初始值選取較好時，本文算法和傳統(tǒng)牛頓法都能準(zhǔn)確的得到目標(biāo)位置，速度的正確解，但當(dāng)目標(biāo)位置，速度初始值較差時(圖2)，傳統(tǒng)牛頓法迭代發(fā)散，無法得到目標(biāo)的真實值，但本文算法憑借其強(qiáng)大的穩(wěn)定性，依然能夠通過迭代獲得目標(biāo)位置，速度的正確解，同時在圖2中可以看到，只要求得目標(biāo)的位置信息足夠精確，目標(biāo)的速度初始值即使偏差很大，依然能夠快速準(zhǔn)確的收斂。本文算法和傳統(tǒng)牛頓法主要的計算量在于矩陣的求逆運算，因此其計算復(fù)雜度分別約為O(M1×33+M2×63)和O(N×63)，其中M1, M2表示本文算法步驟1和步驟2的迭代次數(shù)，N表示傳統(tǒng)牛頓法的迭代次數(shù)。由圖1可以看到本文算法計算量約為O(13×33+4×63)，傳統(tǒng)牛頓法計算量約為O(5×63)，本文算法在初始值良好的情況下計算量要略多于傳統(tǒng)牛頓算法，其原因在于當(dāng)初始值足夠好時仍然對Hess矩陣進(jìn)行修正，阻礙了算法的收斂速度，但本文算法在兼顧了速度，準(zhǔn)度的情況下，又具有良好的穩(wěn)定性，從而適用范圍更廣，尤其是噪聲較大的環(huán)境。

圖3，圖4比較了本文算法同傳統(tǒng)牛頓法的迭代收斂概率，在沒有先驗信息的情況下，目標(biāo)的初始值不能隨機(jī)選取，實驗中通過目前廣泛使用的WLS算法獲得目標(biāo)位置和速度的解作為初始值。從圖中可以看出，當(dāng)噪聲較小時由WLS算法得到的目標(biāo)初始值比較精確，本文算法和傳統(tǒng)牛頓法收斂概率較高；然而當(dāng)噪聲較大時，由于定位問題的非線性特征，WLS算法獲得的目標(biāo)位置，速度初始值與真實值的誤差逐漸變大，遇到較差的初始值時，傳統(tǒng)牛頓法很容易迭代發(fā)散，因此收斂概率降低。而本文算法憑借強(qiáng)大的穩(wěn)定性，仍然能夠保證在噪聲功率較大情況下的較高的收斂概率。

圖5，圖6為本文算法在近距離和遠(yuǎn)距離目標(biāo)兩種情況下同WLS算法[2]，多維標(biāo)度(Multi-Dimensional Scaling, MDS)算法[6]以及克拉美羅界(CRLB)的比較結(jié)果，圖中可以看出隨著測量噪聲誤差的增加，由于定位問題的非線性特點，閉合式算法WLS, MDS受到門檻效應(yīng)的影響，所求得的目標(biāo)位置，速度解逐漸變差，其誤差逐漸變大，不斷偏離CRLB界限，而本文算法在噪聲較大的情況下，相比于WLS, MDS算法具有更高的定位精度，位置和速度的定位精度達(dá)到了CRLB界限。

5 結(jié)束語

對于目標(biāo)的無源時差頻差定位問題，本文在傳統(tǒng)牛頓法的基礎(chǔ)上引入了正則化理論來修正由較差初始值而帶來的病態(tài)Hess矩陣，通過兩次牛頓算法精確的搜索目標(biāo)的位置，速度。分析和仿真結(jié)果表明，本文算法與傳統(tǒng)牛頓法相比，解決了長期存在的初始值問題，在初始值較差的情況下仍能夠?qū)δ繕?biāo)進(jìn)行精確定位，算法具有強(qiáng)大的收斂性，穩(wěn)定性。與WLS, MDS等閉合式算法相比，盡管其計算量相對較大(迭代算法特性)，但在噪聲較大的情況下定位精度更高，定位誤差更接近于Cramer-Rao界。本文算法核心思想在于對傳統(tǒng)牛頓法中的Hess矩陣進(jìn)行修正，對于其他的觀測模型TOA, AOA，需要用到傳統(tǒng)牛頓法求解目標(biāo)位置時，該算法均可使用，因此本文算法具有良好的實用價值和廣泛的適用范圍。

圖1 好初始值下本文算法與傳統(tǒng)牛頓法性能比較

圖2 差初始值下本文算法 與傳統(tǒng)牛頓法性能比較

圖3 目標(biāo)位置為(280, 325, 375) m時本文算法與傳統(tǒng)牛頓法收斂概率比較

圖4 目標(biāo)位置為(2800, 3250, 2750) m時本文算法與傳統(tǒng)牛頓法收斂概率比較

圖5 目標(biāo)位置為(280, 325, 375) m時本 文算法與WLS, MDS算法性能比較

圖6 目標(biāo)位置為(2800, 3250, 2750) m時 本文算法與WLS, MDS算法性能比較

[1] Kovavisaruch L and Ho K C. Modified Taylor-series method for source and receiver localization using TDOA measurements with erroneous receiver positions[C]. IEEE International Symposium on Circuits and Systems, Kobe, Japan, 2005: 2295-2298.

[2] Ho K C and Xu W. An accurate algebraic solution for moving source location using TDOA and FDOA measurements[J]. IEEE Transactions on Signal Processing, 2004, 52(9): 2453-2463.

[3] Sun M, Yang L, and Ho K C. Accurate sequential selflocalization of sensor nodes in closed-form[J]. Signal Processing, 2012, 92(12): 2940-2951.

[4] Ho K C. Bias reduction for an explicit solution of source localization using TDOA[J]. IEEE Transactions on Signal Processing, 2012, 60(5): 2101-2114.

[5] Wang Y and Ho K C. TDOA Source Localization in the presence of synchronization clock bias and sensor poition erors[J]. IEEE Transactions on Signal Processing, 2013, 61(18): 4532-4544.

[6] Wei H W, Peng R, Wan Q, et al.. Multidimensional scaling analysis for passive moving target localization with TDOA and FDOA measurements[J]. IEEE Transactions on Signal Processing, 2010, 58(3): 1677-1688.

[7] Yu H, Huang G, Gao J, et al.. An efficient constrained weighted least squares algorithm for moving source location using TDOA and FDOA measurements[J]. IEEE Transactions on Wireless Communications, 2012, 11(1): 44-47.

[8] Deng B, Xiong J, and Xia C. The feasibility analysis of aerial moving target location based on TDOA/FDOA/DSF Measurements[J]. Procedia Engineering, 2012, 29: 785-790.

[9] Lin L, So H C, Chan F K W, et al.. A new constrained weighted least squares algorithm for TDOA-based localization[J]. Signal Processing, 2013, 93(11): 2872-2878.

[10] Xu B, Sun G, Yu R, et al.. High-accuracy TDOA-based localization without time synchronization[J]. IEEE Transactions on Parallel and Distributed Systems, 2013, 24(8): 1567-1576.

[11] Meng W, Xie L, and Xiao W. Decentralized TDOA sensor pairing in multihop wireless sensor networks[J]. IEEE Signal Processing Letters, 2013, 20(2): 181-184.

[12] Huang J and Wan Q. Analysis of TDOA and TDOA/SS based geolocation techniques in a non-line-of-sight environment[J]. Journal of Communications and Networks, 2012, 14(5): 533-539.

[13] Hara S, Anzai D, Yabu T, et al.. A perturbation analysis on the performance of TOA and TDOA localization in mixed LOS/NLOS environments[J]. IEEE Transactions on Communications, 2013, 61(2): 679-689.

[14] Yang K, Wang G, and Luo Z Q. Efficient convex relaxation methods for robust target localization by a sensor network using time differences of arrivals[J]. IEEE Transactions on Signal Processing, 2009, 57(7): 2775-2784.

[15] Wang G and Chen H. An importance sampling method for TDOA-based source localization[J]. IEEE Transactions on Wireless Communications, 2011, 10(5): 1560-1568.

[16] Hansen P C. Regularization tools: a matlab package for analysis and solution of discrete ill-posed problems[J]. Numerical Algorithms, 1994, 6(1): 1-35.

[17] Hansen P C, Jensen T K, and Rodriguez G. An adaptive pruning algorithm for the discrete L-curve criterion[J]. Journal of Computational and Applied Mathematics, 2007, 198(2): 483-492.

房嘉奇： 男，1984年生，博士生，研究方向為無源定位與跟蹤.

馮大政： 男，1959年生，教授，博士生導(dǎo)師，研究方向為無源定位、雷達(dá)成像、陣列信號處理、盲信號處理、神經(jīng)網(wǎng)絡(luò)等.

李 進(jìn)： 男，1985年生，博士生，研究方向為盲信號處理.

A Robustly Convergent Algorithm for Source Localization Using Time Difference of Arrival and Frequency Difference of Arrival

Fang Jia-qi Feng Da-zheng Li Jin
(National Laboratory of Radar Signal Processing, Xidian University, Xi'an 710071, China)

To pursue accurate source location and velocity, this paper proposes a method based on the Regularization theory to solve the source localization problem utilizing Time-Difference-Of-Arrival (TDOA) and Frequency-Difference-Of-Arrival (FDOA). The proposed algorithm determines the objective function using the maximum likelihood estimator, and then uses classical Newton method to estimate the source position and velocity in an iterative way. It is known that the Newton method requires a good initial value, and a bad initial value can cause an ill-posed Hess matrix which leads to the iteration divergence. This paper introduces the Regularization theory to modify the Hess matrix to make it more proper, which ensures the iteration convergence. The experiment results show that compared with the classical Newton method, the proposed algorithm is robust to the initial value, and is still able to ensure its convergence even with an inaccurate initial value of large error. Compared with some other closed-form source location methods, the proposed algorithm has better location accuracy in large noise levels which can achieve the Cramer-Rao bound. The proposed algorithm can be widely applied in practice.

Passive location; Time-Difference-Of-Arrival (TDOA); Frequency-Difference-Of-Arrival (FDOA); Newton method; Hess matrix; Regularization algorithm

TN911.7

： A

：1009-5896(2015)04-0798-06

10.11999/JEIT140560

2014-04-30收到，2014-11-24改回

國家自然科學(xué)基金(61271293)資助課題

*通信作者：房嘉奇 fangjiaqi123@hotmail.com