王志超，曹起武，張 全

(1.沈陽工業(yè)大學信息學院，遼寧沈陽 110870；2.遼寧機電職業(yè)技術(shù)學院信息系，遼寧丹東 118009)

?

MCMC粒子濾波和復化Newton-cotes算法測算區(qū)域面積的方法

王志超1，曹起武2，張 全1

(1.沈陽工業(yè)大學信息學院，遼寧沈陽 110870；2.遼寧機電職業(yè)技術(shù)學院信息系，遼寧丹東 118009)

針對不規(guī)則區(qū)域面積測算中定位精度和面積計算精度兩方面不足，提出一種定位精度高、面積誤差小的面積測算新方法。其采用一種組合定位方法精確定位，即將差分GPS測量系統(tǒng)(DGPS)與馬爾可夫鏈蒙特卡羅(Markov chain Monte Carol，MCMC)粒子濾波相結(jié)合，再配合復化Newton-cotes算法，擬合邊界曲線并準確求得區(qū)域面積。將MCMC粒子濾波應用于DGPS定位數(shù)據(jù)處理，其既可處理非高斯分布噪聲，又解決粒子濾波(PF)的粒子退化問題，提高定位精度。將復化Newton-cotes算法應用于面積計算，其既避免高次插值的舍入誤差，又將面積區(qū)間進一步細分，提高面積計算精度。實驗結(jié)果表明，該新方法定位精度更高，面積誤差更小。

不規(guī)則區(qū)域；面積測算方法；差分GPS測量系統(tǒng)；馬爾可夫鏈蒙特卡羅；粒子濾波；復化Newton-cotes算法

0 引言

借助GPS 測算不規(guī)則區(qū)域面積的方法研究既是當今物聯(lián)網(wǎng)的一個具體應用，又是智能控制理論研究的課題之一[1-2]。GPS已經(jīng)在林業(yè)面積估算、農(nóng)業(yè)土地面積測量以及城市規(guī)劃等不規(guī)則區(qū)域面積測量方面得到了越來越廣泛的應用[3]。然而，一些測量不規(guī)則區(qū)域面積的傳統(tǒng)方法存在定位精度低、速度慢、面積計算誤差大等諸多問題[4]。以傳統(tǒng)方格法和數(shù)字電子地圖法為例，主要存在以下缺陷：

(1)環(huán)境干擾導致定位不準確[5]；

(2)直線代替曲線導致邊界擬合誤差大[6]；

(3)面積求算方法誤差大[7]。

綜上，提出一種不規(guī)則區(qū)域面積測算的新方法，其能較快速且準確地測算不規(guī)則區(qū)域的面積。

首先介紹了C/A偽碼差分GPS測量系統(tǒng)(DGPS)，并引用馬爾可夫鏈蒙特卡羅移動步驟(Markov Chain Monte Carlo，MCMC )對粒子濾波(Particle filter，PF)進行改進，提出DGPS與MCMC粒子濾波相結(jié)合的組合定位方法；隨后詳細介紹并推導復化牛頓-柯特斯算法(Compound Newton-cotes algorithm)，并給出其與改進方格法相結(jié)合求被測區(qū)域面積的詳細過程描述；最后分別對兩種算法進行實驗驗證，給出實驗結(jié)果。仿真結(jié)果表明：與傳統(tǒng)不規(guī)則區(qū)域面積測量方法相比，DGPS與MCMC粒子濾波組合定位能大幅提高GPS坐標定位精度，復化牛頓-柯特斯求積算法能提高求積精度。該方法是一種定位精度更高，計算誤差更小，抗干擾性更強的不規(guī)則區(qū)域面積測量新方法。

1 DGPS與MCMC粒子濾波組合定位

1.1 C/A偽碼DGPS定位

由于大氣條件等不定因素影響，大多單機GPS定位精度只能達到±10 m級的水平[8-9]，如此精度測量較小面積顯然無精度可言。故通常采用C/A偽碼差分GPS測量系統(tǒng)(DGPS)提高定位精度[10]。DGPS可根據(jù)基準站已知位置信息，計算偽距中誤差修正值，實時修正定位數(shù)據(jù)，有效消除常值誤差。但其差分精度隨基準站到用戶的距離增加而降低，而這種誤差是用任何差分法都不能消除的[11-12]。因此，采用DGPS與MCMC粒子濾波組合定位的方法來盡可能消除該誤差。

1.2 MCMC粒子濾波

粒子濾波是基于貝葉斯估計的濾波方法，其基本思想是：由系統(tǒng)狀態(tài)向量經(jīng)驗條件分布，在其狀態(tài)空間中產(chǎn)生一組隨機樣本(即粒子)，隨后由觀測量情況不斷調(diào)整粒子的位置和權(quán)重，根據(jù)調(diào)整后粒子信息，修正最初經(jīng)驗條件分布，如果有足夠多的粒子，則經(jīng)修正后的經(jīng)驗條件分布將收斂于系統(tǒng)狀態(tài)向量的真實條件分布，此時，系統(tǒng)狀態(tài)向量的估計值就可由這些粒子的均值得到[13]。

通常GPS動態(tài)定位的觀測方程均為非線性的，為使定位數(shù)據(jù)平滑輸出，提高數(shù)據(jù)精度，常用卡爾曼濾波等濾波方法[14]。然而，這些濾波方法要求動力學模型噪聲和觀測噪聲為獨立不相關(guān)的高斯白噪聲，且要求觀測模型和動力學模型為線性模型，而線性化過程會引入舍入誤差[15]，因此這些濾波算法不具備最優(yōu)性。而粒子濾波(PF)算法可彌補該不足，其無須對觀測模型線性化處理，且不要求動力學模型噪聲和觀測噪聲為高斯白噪聲，而直接通過蒙特卡羅模擬方法遞推貝葉斯濾波，因此可用于由非線性觀測模型表示的動態(tài)濾波系統(tǒng)[16]。

但PF經(jīng)多次迭代后，易出現(xiàn)粒子退化現(xiàn)象[17]，它使粒子喪失多樣性，導致采樣枯竭，降低算法性能。為解決該問題，在PF重采樣后引入馬爾可夫鏈蒙特卡羅移動步驟，增加粒子的多樣性，提高算法的性能[18]。仿真結(jié)果表明：MCMC粒子濾波能有效提高GPS動態(tài)定位的精度，其濾波效果優(yōu)于一般PF和UKF等濾波算法[19]。

1.3 MCMC粒子濾波算法描述

設(shè)動態(tài)濾波系統(tǒng)的狀態(tài)方程和觀測方程可描述為

(1)

式中：Xk為狀態(tài)向量;fk為狀態(tài)轉(zhuǎn)移函數(shù);hk為觀測向量和狀態(tài)向量間的傳遞函數(shù);Zk為測量向量;nk為觀測噪聲;υk-1為系統(tǒng)噪聲[20]。

設(shè)第k-1時刻有一組后驗粒子集{xk-1(i),ωk-1(i);i=1,2,…,N}，其中粒子數(shù)為N，xk-1(i) 表示第k-1時刻第i個粒子，ωk-k-1(i)表示第k-1時刻第i個粒子的權(quán)重。

(1)初始化粒子集，令k=0:

(2)當k=1,2,…時，執(zhí)行以下步驟：

①狀態(tài)預測

根據(jù)系統(tǒng)的狀態(tài)方程計算p(Xk|Xk-1)，由計算得到的概率密度p(Xk|Xk-1)抽取第時刻的先驗粒子：{Xk|k-1(i);i=1,2…,N}～P(Xk|Xk-1)。

(2)

歸一化權(quán)值為

(3)

(4)

(c)估計

計算系統(tǒng)此時的狀態(tài)估計值：

(5)

(1)在區(qū)間上以均勻概率抽樣得門限值u,u～U[0，1]。

(2)由重要性概率密度函數(shù)p(xk|xk-1(i))抽樣得xk|k-1(i),即xk|k-1(i)～p(xk|xk-1)。

1.4 MCMC粒子濾波處理GPS定位數(shù)據(jù)算法描述

將上述MCMC粒子濾波應用于GPS定位數(shù)據(jù)的處理算法描述如下：

1.4.1 初始化粒子集，k=0

1.4.2 Fork=1,2,…循環(huán)執(zhí)行以下步驟

(1)根據(jù)系統(tǒng)狀態(tài)方程抽取第k時刻的先驗粒子{Xk|k-1(i):i=1,2,…,N}～p(Xk|Xk-1)。

(2)更新

①更新權(quán)值：測量結(jié)束后，由系統(tǒng)的觀測方程及式(6)計算粒子權(quán)值ωk(i)：

ωk(i)=ωk-1(i)p(Zk|Xk(i))i=1,…,N

(6)

將權(quán)值歸一化為

(7)

(3)估算此刻系統(tǒng)狀態(tài)估計值:

1.5 組合定位

綜上所述，MCMC粒子濾波在消除系統(tǒng)隨機誤差，保證動態(tài)性能的前提下，減小GPS定位數(shù)據(jù)的誤差；而DGPS技術(shù)可有效消除電離層折射等常值誤差[24]。若將這2種技術(shù)組合起來，既可減少由非正確預測引起的隨機誤差，又可減少在GPS定位時所產(chǎn)生的常值誤差，從而大大提高GPS定位精度。實驗證明：使用單機定位的精度誤差為2.2 m (2σ)，而2種技術(shù)組合使用后其定位誤差減少到0.8 m(2σ)以內(nèi)[25]。

2 復化Newton-cotes算法

2.1 算法推導

本方法采用復化Newton-cotes算法計算區(qū)域面積。其基本思想是通過復化公式等分目標區(qū)間，并在每一個小區(qū)間內(nèi)均使用一次Newton-cotes算法，即通過插值法用一個簡單且易于積分的函數(shù)逼近目標函數(shù)fi(x),以計算積分I(fi),最終將每份面積I(fi)加和得到該區(qū)域面積I(f)。該算法通過復化公式降低插值階數(shù)，從而減小舍入誤差，同時又將目標區(qū)間進一步細化，精化計算面積。

(1)插值法求擬合曲線方程：由n次插值多項式：

(8)

擬合函數(shù)f(x)。

由f(x)≈Ln(x)，得

(9)

(2)Newton-cotes公式：采用等距節(jié)點情形的Newton-cotes算法。構(gòu)造插值型求積公式：

(10)

對比式(9)，得

(11)

(12)

則當n=4時，式(10)可化為

(13)

式(13)即為柯特斯公式，它具有5次代數(shù)精度。

由式(12)可推出：當n增加到8以上時出現(xiàn)負系數(shù)，這會增大舍入誤差，從而無法保證收斂性和穩(wěn)定性。故采用高次插值不僅計算復雜，且易出現(xiàn)龍格現(xiàn)象，因此在實際計算中不能用高階的Newton-cotes公式[26]。為在算法上進一步提高面積計算精度，則需要加密求積節(jié)點，但由上述問題不能無限制提高求積公式的階數(shù)。故最終在Newton-cotes公式中引用復化公式，其既避免了由高次插值帶來的舍入誤差，又達到了對面積區(qū)間進一步細分的目的，從而大幅提高了面積的計算精度。

(14)

2.2 算法描述

復化Newton-cotes算法描述如下：

(1)擬合未知曲線：輸入節(jié)點初值x0,x1,…,xn及其對應的函數(shù)值f(xk)(k=0,1,…,n)，由插值公式(8)求得未知曲線擬合方程F(x)。

(3)計算面積：由式(14)計算未知區(qū)域面積。

3 不規(guī)則區(qū)域面積測算新方法

綜上，結(jié)合上述兩方法各自優(yōu)點，將傳統(tǒng)方格法進行改進，提出不規(guī)則區(qū)域面積測算新方法：

(1)先由C/A偽碼差分系統(tǒng)(DGPS)獲得被測區(qū)域的邊界坐標，隨后配合MCMC粒子濾波算法平滑濾波輸出，進一步提高定位坐標精度，在滿足定位精度的條件下可將被測區(qū)域邊界以間隔為1 m的點描標記出來，并分別找出上、下、左、右四個方向的極值點，以確定該不規(guī)則區(qū)域所在的最小外接矩形，如圖1下半部分不規(guī)則區(qū)域外接矩形所示。

圖1 分割、擬合方法示意圖

(2)隨后，在確定劃分方格的大小時，為了便于使用復化Newton-cotes算法，這里采用大小為40 m×40 m的小方格將該外接矩形坐標平面分割為若干等大的小方格，如圖1下半部分劃分的小方格所示。

(3)接下來要使用改進的方格法對區(qū)域面積進行計算。通過上述劃分方法，可以得到與被測區(qū)域無重疊、部分重疊和完全重疊這三類小方格。其中，與被測區(qū)域無重疊的方格在計算面積時不予考慮；與被測區(qū)域完全重疊的M個方格其重疊面積已知，記為S=1 600 m2；而對于那些與被測區(qū)域部分重疊的N個方格，要分別使用復化Newton-cotes算法進行面積計算：邊界曲線f(x)如圖1上半部分所示，先將該方格某邊按等間隔細分10份，再分別于每一份小區(qū)間中使用一次Newton-cotes公式(即再次細分4份，此時該方格相當于被細分為40份，每個分割節(jié)點間距為1 m，符合上述坐標定位精度)，最終可由復化Newton-cotes公式求得與區(qū)域重疊面積大小，記為Si，如圖1上半部分中，邊界曲線f(x)與坐標軸所圍面積。

4 實驗與仿真

4.1 MCMC粒子濾波算法MATLAB仿真

4.1.1 仿真參數(shù)設(shè)定

為驗證MCMC粒子濾波優(yōu)于UKF、PF等傳統(tǒng)濾波算法，對同一狀態(tài)模型分別采用UKF、PF和MCMC粒子濾波算法進行比較。MATLAB仿真參數(shù)取值：時間步長T取5 s，采樣間隔取τ= 1 s，濾波時間取50 s，粒子濾波的粒子數(shù)取500。 模型的狀態(tài)初值x0=0.1,初始狀態(tài)估計值：Xu=Xpf=Xmc=1,初始估計方差值：Pu=Ppf=Pmc=5 。設(shè)過程系統(tǒng)噪聲和觀測噪聲均為白噪聲，則過程系統(tǒng)噪聲均方差：Vu=Vpf=Vmc=1,觀測噪聲均方差：Nu=Npf=Nmc=0.1。

4.1.2 仿真結(jié)果

通過MCMC粒子濾波算法得到的狀態(tài)估計值與真實值關(guān)系如圖2所示。由圖2可以看出，經(jīng)MCMC粒子濾波算法濾波得到的估計值均分布在真實值的95%置信區(qū)間之中，且其很接近真實值。

圖2 MCMC粒子濾波狀態(tài)向量估計值

由圖3可以看出，由UKF、PF和MCMC粒子濾波算法得到的濾波結(jié)果雖然都能逼近真值，但MCMC粒子濾波算法逼近真值的程度會更好一些。

圖3 各算法的狀態(tài)向量估計值與真實值

由圖4可以看出，分別采用UKF、PF和MCMC粒子濾波時，通過MCMC粒子濾波算法得到的狀態(tài)估計值誤差始終很小且相對穩(wěn)定，沒有出現(xiàn)類似于UKF和PF的狀態(tài)估計值誤差那樣發(fā)散的趨勢，因此，也可說明MCMC粒子濾波算法優(yōu)于UKF和PF等算法。

圖4 各算法的狀態(tài)向量估計值誤差

仿真得到的各算法估計誤差均分值如表1所示。

表1 各算法估計誤差均方值

由仿真結(jié)果可知，采用的MCMC粒子濾波作為一種非線性的濾波算法，避免了對狀態(tài)方程的線性化，減少大量線性化誤差，并可有效抑制傳統(tǒng)粒子濾波(PF)算法的粒子退化問題，增加粒子多樣性，提高濾波精度。

4.2 復化Newton-cotes算法實驗

4.2.1 實驗步驟

為了驗證復化Newton-cotes算法在面積計算精度上優(yōu)于傳統(tǒng)求積算法，本文引用張汶賢文中的一個實例區(qū)域[27],分別應用梯形法、Simpson 1/3、Simpson 3/8和復化Newton-cotes算法進行面積測算。具體步驟如下：

(1)首先，對于參考面積已知的圖例區(qū)域以1 m為間隔分別得到區(qū)域邊界各坐標點，圖例區(qū)域如圖5所示。

圖5 圖例區(qū)域

(2)對于該圖例區(qū)域，分別采用梯形法、Simpson 1/3法、Simpson 3/8法以及復化Newton-cotes算法計算其面積，對比各面積結(jié)果及相對誤差。

4.2.2 實驗結(jié)果

由上述4種求積算法計算所得的面積值及相對誤差如表2所示。其中參考面積為430.78 m2。

表2 計算結(jié)果對比

對比以上4種求積算法對同一目標區(qū)域面積的計算結(jié)果可知，在定位坐標精度一定的條件下，采用復化的Newton-cotes算法的計算誤差最小，精度最高。雖然高次插值能增大舍入誤差，且容易出現(xiàn)龍格現(xiàn)象，但通過改進方格法和復化的低次Newton-cotes算法相結(jié)合，對積分區(qū)間進行多次細分，既避免了由高次插值帶來的舍入誤差，又進一步達到了細分積分區(qū)間的目的，大幅提高了面積計算的精度。

5 結(jié)束語

綜上所述，測量不規(guī)則區(qū)域面積的新方法在坐標定位方面，采用集DGPS與MCMC粒子濾波于一體的組合定位，由Matlab仿真結(jié)果知，在改善濾波效果的同時，其修正定位數(shù)據(jù)偏差，使得定位精度大幅提高；在面積計算方面，采用復化Newton-cotes算法計算數(shù)值積分，由實際測算結(jié)果知，其既避免由高次插值帶來舍入誤差，又達到對面積區(qū)間進一步細分的目的，從而大幅提高面積的計算精度。經(jīng)實驗測算，該方法與其他傳統(tǒng)方法相比具有精度高、速度快、誤差小、抗干擾性強等優(yōu)越性，可以廣泛地應用于林業(yè)、農(nóng)業(yè)、城市規(guī)劃以及水利建設(shè)等諸多領(lǐng)域中。同時，本研究將進一步探索基于GPS的不規(guī)則自由曲面面積的測量方法[28-29]。

[1] 張全法，韓要軒，楊海彬，等.用面陣CCD測量不規(guī)則平面物體的面積. 儀表技術(shù)與傳感器，2000(2):31-34.

[2] 關(guān)勝曉.不規(guī)則形體面積的CCD測量研究.儀表技術(shù)與傳感器，1998(11):36-39.

[3] ANTHONY C, JINLING W, CHRIS R. Bridging GPS outages in the agricultural environment using virtualite measurements. Location and Navigation Symposium, 2008, 7(6):497- 504.

[4] 張海濤，楊敏華.不規(guī)則區(qū)域面積的矩形四等分割計算法.測繪與空間地理信息，2011，34(6):272-274.

[5] 馮仲科，劉永霞，王小昆，等.林地面積量算方法的比較研究.北京林業(yè)大學學報，2004，5(2):97-101.

[6] YONG H, HAIHONG Y, HUI F. Study on improving GPS measurement accuracy. Instrumentation and Measurement Technology Conference, Canada, 2005.

[7] 宋鄂平，張前勇.CASS在林權(quán)調(diào)查面積量算中的應用.湖北民族學院學報(自然科學版)，2009，27(1):105-108.

[8] MICHIEL D. GPS based bistatic radar - beyond specular reflection. 24th Digital Avionics Systems Conference, Williamsburg, VA, USA, 2005.

[9] GUPTA S, LOW M, TAN C. The new approach of simulate GPS measurements. Position Location and Navigation Symposium, 1998, 7(1):236-242.

[10] 馮仲科，梁長秀，南永天.超小圖斑面積的RTD GPS量測.林業(yè)資源管理，2001，5(3):70-72.

[11] SHUPING G, ZHANGHAO Z, MATTHEW T. GPS Spoofing Based Time Stamp Attack on Real Time Wide Area Monitoring in Smart Grid . IEEE Smart Grid Comm 2012 Symposium - Cyber Security and Privacy, Nepal, 2012.

[12] ABDEL-SALAM M. Precise point positioning using un-differenced code and carrier phase observations. IEEE International Conference, Canada, 2005.

[13] 田世君，陳俊，皮亦鳴.粒子濾波在高動態(tài)GPS定位中的應用.測繪學報，2007，36(3)：274-278.

[14] JULIER S, UHLMANN J. A new extension of the Kalman filter to nonlinear systems. Proceedings of SPIE, 1997, 2(3)：182-193.

[15] JULIER S, UHLMANN J. A new method for the nonlinear transformation of means and covariance in filters and estimators .IEEE Transactions on Automatic Control, 2000 , 45(3):477-482.

[16] 王爾申，蔡明，龐濤.MCMC粒子濾波的GPS定位數(shù)據(jù)處理算法.數(shù)據(jù)采集與處理，2013，28(2):213-218.

[17] 秦臻.改進的粒子濾波及其在GPS動態(tài)定位中的應用.全球定位系統(tǒng)，2010，5(1):25-28.

[18] 高靜，李善姬，邵奎軍.MCMC粒子濾波算法研究及應用.電子測試，2009，12(1):19-22.

[19] LIU J, CHEN R .Sequential Monte Carlo methods for dynamic systems. J Amer Statist Assor, 1998, 9 (3):1032-1044.

[20] SANJEEV M, MASKELL S. A tutorial on particle filters for online nonlinear /non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 2002, 50(2):174-188.

[21] LIU J, COLAJANNI. Dynamic load balancing on Web-server systems. IEEE Internet Computing, 1999, 3(3):28-39.

[22] IKOMA N, ICHIMURA N, HIGUCHI T, et al. Particle filter based method for maneuvering target tracking. In Proceedings of the IEEE International Workshop of Intelligent Signal Processing, Budapest,2010 .

[23] 趙琳，聶琦，高偉.基于MCMC方法的正則粒子濾波算法及其應用.儀器儀表學報，2008，29(10)：2156-2162.

[24] HU H, JING Z, LI A, et al. An MCMC-based particle filter for tracking target in glint noise environment. Journal of Systems Engineering and Electronics, 2005, 12(2)：305-309.

[25] SPALL J C .Estimation via Markov chain Monte Carlo. IEEE Control Systems Magazine, 2003, 23(2):34-45.

[26] 王金銘，劉艷秋，陳欣.數(shù)值分析.大連：大連理工大學出版社，2007.

[27] 張汶賢.計算不規(guī)則圖形面積的一個新公式.江漢石油學院學報，1989，11(4):83-88.

[28] 胡丹，李敏，常江.復雜輪廓曲面非接觸三維數(shù)據(jù)測量的研究. 儀表技術(shù)與傳感器，2002(7):23-25.

[29] 王菽芳，邱自學，袁江，等.集成激光位移傳感器與線性編碼器的自由曲面測量方法及系統(tǒng). 儀表技術(shù)與傳感器，2010(7)：4-9.

作者簡介：王志超(1990—)，碩士，主要研究領(lǐng)域為高精度不規(guī)則圖形區(qū)域面積測量、中央空調(diào)控制器設(shè)計等。 E-mail：694428961@qq.com 曹起武(1979—)，副教授，碩士，主要研究領(lǐng)域為智能控制、軟件與網(wǎng)絡(luò)等。E-mail：caoqiwu@163.com

Region Area Measurement Method Based on MCMCParticle Filter and Compound Newton-cotes Algorithm

WANG Zhi-chao1, CAO Qi-wu2, ZHANG Quan1

(1.School of Information Engineering, Shenyang University of Technology, Shenyang 110870, China;2. Department of Information, Liaoning Mechanical and Electrical Polytechnic, Dandong 118009, China)

To improve disadvantages of irregular-area measurement in positioning accuracy and area calculation accuracy, a new area measuring method which has characteristics of high precision and tiny area error was proposed in this paper. It combined the differential GPS measurement system and Markov chain Monte Carol particle filter to locate. And it fitted the boundary curve and calculated areas through the compound Newton-cotes algorithm. The new method used MCMC particle filter to process GPS data. The filter algorithm processed the non-Gaussian distribution noise and improved particle degradation. And the new method used compound Newton-cotes algorithm to calculate the area. The quadrature algorithm avoided the rounding error brought by the high-order interpolation, and further subdivided the area. Therefore, the new method improved positioning accuracy and area calculation accuracy. Experimental results show that the new method has characteristics of high accuracy and tiny errors.

irregular regions; area measurement methods; DGPS; Markov chain Monte Carlo (MCMC); particle filter (PF); compound Newton-cotes algorithm

馬俊(1990—)，碩士研究生，主要研究方向為電力電子和嵌入式系統(tǒng)。E-mail:mark41@163.com 李春茂(1963—)，教授，博士學位，主要研究領(lǐng)域為控制系統(tǒng)與電工理論。E-mail:chunmao@126.com

遼寧省科技廳基金項目(20130125)

2014-02-24 收修改稿日期：2015-03-20

P218

A

1002-1841(2015)06-0121-06