王亞菲，鐘世明，王海濤，歐吉坤

1. 中國科學院測量與地球物理研究所，湖北 武漢 430077； 2. 大地測量與地球動力學國家重點實驗室，湖北 武漢 430077； 3. 中國科學院大學，北京 100049

?

LEO衛(wèi)星軌道預(yù)報精度分析

王亞菲1,2,3，鐘世明1,2，王海濤1,2，歐吉坤1,2

1. 中國科學院測量與地球物理研究所，湖北 武漢 430077； 2. 大地測量與地球動力學國家重點實驗室，湖北 武漢 430077； 3. 中國科學院大學，北京 100049

利用動力學擬合法，以HY-2A衛(wèi)星為例分析不同擬合區(qū)間對預(yù)報不同弧長的軌道精度的影響?；贑NES提供的事后精密星歷和采用非差簡化動力學方法獲得的厘米級快速軌道兩種產(chǎn)品，進行衛(wèi)星軌道預(yù)報。結(jié)果表明：采用24 h和12 h擬合區(qū)間分別預(yù)報12 h軌道時，其3DRMS優(yōu)于6 dm；當預(yù)報24 h軌道時，3DRMS優(yōu)于1 m。

動力學擬合；Collocation積分；擬合區(qū)間；預(yù)報； 快速定軌

基于低軌衛(wèi)星的導(dǎo)航增強系統(tǒng)是目前國內(nèi)外研究的熱點[1-4]。從美國Iridium系統(tǒng)輔助GPS系統(tǒng)的研究成果來看，利用LEO衛(wèi)星對導(dǎo)航系統(tǒng)進行性能增強是一種有效的手段[3]。在此增強系統(tǒng)中，LEO衛(wèi)星具有“源”“端”“中繼”3種作用[4]。當LEO衛(wèi)星作為測距源時，與WAAS系統(tǒng)中GEO衛(wèi)星類似，需要發(fā)布自身的廣播星歷，而廣播星歷的生成則需要進行軌道預(yù)報。因此，LEO衛(wèi)星軌道的預(yù)報是該導(dǎo)航增強系統(tǒng)需要解決的重要問題之一。

目前，軌道預(yù)報方法主要有3種。第1種是分析法[5-6]，即在分析獲得衛(wèi)星t時刻平均根數(shù)的基礎(chǔ)上給出的一種軌道預(yù)報方法。該方法在空間環(huán)境監(jiān)測和實時跟蹤測量等對預(yù)報軌道精度要求不高的領(lǐng)域有著很好的應(yīng)用效果。第2種是基于多項式擬合的方法，在不考慮衛(wèi)星的受力情況下，利用多項式擬合衛(wèi)星的精密軌道得到衛(wèi)星的速度場，據(jù)此進行軌道預(yù)報。常用的有最小二乘曲線擬合法、切比雪夫多項式擬合法、拉格朗日多項式擬合法[7-8]。低軌衛(wèi)星由于軌道較低，受力情況相對于中高軌衛(wèi)星更為復(fù)雜，對于不涉及動力學模型的預(yù)報方法而言，只能進行短期(幾分鐘甚至更短)預(yù)報，預(yù)報精度為厘米甚至分米級[9-10]。第3種是基于動力學軌道擬合的方法[11-12]，用已知的衛(wèi)星位置坐標采用動力學方法擬合出初始軌道和動力學參數(shù)，在此基礎(chǔ)上通過積分進行軌道外推。該方法與前兩種方法相比，在達到分米級預(yù)報精度的前提下，能夠?qū)崿F(xiàn)更長弧長的軌道預(yù)報，滿足低軌衛(wèi)星增強導(dǎo)航衛(wèi)星系統(tǒng)性能的應(yīng)用[11]。

本文以HY-2A衛(wèi)星為例進行LEO衛(wèi)星軌道預(yù)報研究，結(jié)合動力學擬合及Collocation方法[12-15]，重點分析了擬合弧長及預(yù)報時長對軌道預(yù)報精度的影響。在此基礎(chǔ)上，基于超快速星歷產(chǎn)品實現(xiàn)HY-2A衛(wèi)星高精度快速定軌的同時開展軌道預(yù)報研究，對進一步研究基于LEO衛(wèi)星的導(dǎo)航增強系統(tǒng)具有很好的參考價值。

1　動力學擬合方法

1.1動力學擬合基本原理

(1)

假設(shè)在LEO衛(wèi)星精密星歷中，對應(yīng)t時刻的位置為rt，則采用動力學擬合的觀測方程可以表示如下[18]

(2)

假定在一段弧長[t1,tm]內(nèi)，LEO衛(wèi)星精密星歷中的ti時刻衛(wèi)星位置及速度為ri，則可以構(gòu)建m維觀測方程

(3)

式中

(4)

(5)

1.2定軌和預(yù)報采用的攝動力模型

LEO衛(wèi)星由于軌道高度較低，其受力情況更為復(fù)雜。因此，在進行LEO衛(wèi)星軌道確定及預(yù)報過程中需要考慮地球非球形引力、N體攝動、固體潮、海潮、大氣阻力、太陽光壓、地球輻射壓及相對論效應(yīng)等引起的攝動[19-21]。其中，在采用非差簡化動力學法進行快速精密定軌時，大氣阻力通過設(shè)置偽隨機脈沖參數(shù)[22]進行補償(軌道預(yù)報過程，沒有設(shè)置偽隨機脈沖參數(shù))。本文在數(shù)據(jù)處理過程中所采用的動力學模型如表1所示。

表1　相關(guān)的動力學模型

2　基于事后精密星歷的軌道預(yù)報精度分析

為分析不同弧長軌道數(shù)據(jù)的擬合精度以及不同的擬合區(qū)間長度對不同弧長的預(yù)報影響，本文采用法國國家空間中心(CNES)提供的HY-2A衛(wèi)星上SLR、DORIS和星載GPS數(shù)據(jù)聯(lián)合定軌[23-24]生成的精密星歷進行方法測試(徑向誤差在1 cm左右，沿跡和法向上的誤差為2～3 cm，采樣間隔60 s)。利用Bernese5.2軟件對2014年沒有軌道機動情況下的DOY 60—128共計65 d(DOY 73和DOY 106出現(xiàn)軌道機動，暫不處理)的精密星歷進行動力學擬合及預(yù)報，并統(tǒng)計了軌道擬合及預(yù)報結(jié)果與參考軌道的位置偏差，預(yù)報方案如表2所示。

表2　LEO衛(wèi)星軌道擬合與預(yù)報方案

2.1衛(wèi)星軌道擬合精度分析

根據(jù)以上提出的擬合方案，采用動力學擬合法對HY-2A衛(wèi)星進行相應(yīng)的軌道擬合與預(yù)報數(shù)據(jù)處理，其不同擬合弧長的擬合精度(3DRMS值)如圖1所示。結(jié)果顯示：隨著擬合弧長的增加，擬合殘差逐漸增大。其中，24 h弧長擬合3DRMS值小于8 cm，2 h弧長擬合3DRMS在毫米級水平。

圖1　HY-2A衛(wèi)星不同弧長軌道擬合的3DRMS值Fig.1　3DRMS of HY-2A satellite different arc length orbit fitting

2.2衛(wèi)星軌道預(yù)報結(jié)果分析

根據(jù)2.1節(jié)擬合得到相關(guān)的動力學參數(shù)進行數(shù)據(jù)處理，分別基于不同擬合區(qū)間預(yù)報2 h、4 h、8 h、12 h和24 h弧長的HY-2A衛(wèi)星軌道。在預(yù)報過程中采用與軌道擬合時相同的動力學模型及參數(shù)，積分步長為3 min，式(4)中多項式階數(shù)q=10。圖2—圖6統(tǒng)計了不同擬合區(qū)間分別預(yù)報2 h、4 h、8 h、12 h和24 h弧長的軌道預(yù)報精度(3DRMS)，圖7為2014年DOY 64的不同擬合區(qū)間預(yù)報24 h弧長在R(徑向)、S(沿跡方向)、W(法向)3個方向上的差值時序圖。

(1) 圖2—圖6表明擬合區(qū)間對軌道的預(yù)報精度影響較大。當預(yù)報2 h弧長時，選擇2 h的擬合區(qū)間長度較為適宜，其3DRMS可達厘米級；當預(yù)報4 h、8 h、12 h和24 h弧長時，采用12 h或者24 h的擬合區(qū)間長度進行外推，其預(yù)報結(jié)果相對比較穩(wěn)定，精度較好，24 h弧長預(yù)報軌道的3DRMS達到分米級，滿足高精度應(yīng)用需求。

(2) 圖7表明當擬合區(qū)間長度一定時，隨著預(yù)報弧長的增加，相比較于R方向及 W方向，S方向的預(yù)報精度最差，下降最快[25-26]。

3　基于超快速精密定軌結(jié)果的軌道預(yù)報精度分析

由于CNES發(fā)布的HY-2A衛(wèi)星星歷是一種事后精密軌道，存在較大時間延遲，在此基礎(chǔ)上的軌道預(yù)報在實際應(yīng)用中存在一定限制。因此，需要對HY-2A衛(wèi)星實現(xiàn)快速精密定軌，以便獲得較小延遲下的HY-2A精密軌道。本文采用中科院測量與地球物理研究所(WHIGG)iGMAS數(shù)據(jù)分析中心提供的GPS超快速產(chǎn)品(時延3～9個小時，星歷采樣率為15 min，鐘差采樣率為30 s)對HY-2A衛(wèi)星進行快速精密定軌及預(yù)報。

圖2　HY-2A不同弧長預(yù)報2 h軌道3DRMSFig.2　3DRMS of 2 h prediction orbit with HY-2A satellite different arc length

圖3　HY-2A不同弧長預(yù)報4 h軌道3DRMSFig.3　3DRMS of 4 h prediction orbit with HY-2A satellite different arc length

圖4　HY-2A不同弧長預(yù)報8 h軌道3DRMSFig.4　3DRMS of 8 h prediction orbit with HY-2A satellite different arc length

圖5　HY-2A不同弧長預(yù)報12 h軌道3DRMSFig.5　3DRMS of 12 h prediction orbit with HY-2A satellite different arc length

圖6　HY-2A不同弧長預(yù)報24 h軌道3DRMSFig.6　3DRMS of 24 h prediction orbit with HY-2A satellite different arc length

圖7　HY-2A不同弧長預(yù)報24 h的差值時序圖Fig.7　Error of 24 h prediction orbit with HY-2A satellite different arc length

圖8統(tǒng)計的是基于2015年DOY 209—342的快速星歷獲得HY-2A衛(wèi)星非差簡化動力學軌道與CNES提供的參考軌道在R、S、W和三維位置比較的RMS。從中可以得到，快速精密軌道與參考軌道的差值在徑向、沿跡和法向上的平均值分別為1.41 cm、3.30 cm、2.29 cm，3DRMS的平均值為4.27 cm，與文獻[27]中結(jié)論相一致。

圖8　采用WHIGG數(shù)據(jù)分析中心超快速星歷定軌結(jié)果Fig.8　RMS of HY-2A satellite orbit by using WHIGG analysis center GPS ephemeris

根據(jù)2.2節(jié)得到的初步結(jié)論，本節(jié)分別采用HY-2A衛(wèi)星12 h和24 h弧長的快速定軌結(jié)果進行動力學擬合并預(yù)報軌道。圖9表示的是兩種不同擬合弧長(12 h和24 h)分別進行軌道預(yù)報得到2015年DOY 330—343軌道結(jié)果與參考軌道的三維位置誤差。從圖中可以得到，24 h和12 h弧長分別預(yù)報12 h軌道時，其3DRMS優(yōu)于6 dm；當預(yù)報24 h軌道時，3DRMS優(yōu)于1 m。

圖9　利用iGMAS超快速產(chǎn)品計算的不同弧長精密軌道預(yù)報12 h和24 h軌道比較結(jié)果Fig.9　3DRMS of prediction orbit with 24 h or 12 h rapid precise orbit

4　結(jié)　論

本文基于動力學擬合法以HY-2A衛(wèi)星為例開展了軌道預(yù)報研究，分析了不同擬合區(qū)間長度對不同弧長預(yù)報精度的影響。在此基礎(chǔ)上，提出利用GPS超快速產(chǎn)品在實現(xiàn)快速精密定軌的同時，通過采用動力學擬合法可以獲得較高精度的預(yù)報軌道。研究結(jié)果表明：

(1) 對于要求預(yù)報軌道超過2 h弧長的應(yīng)用而言，建議采用12 h或者24 h擬合區(qū)間進行軌道預(yù)報，其預(yù)報精度相對較為穩(wěn)定。

(2) 利用WHIGG數(shù)據(jù)分析中心的超快速產(chǎn)品在實現(xiàn)LEO衛(wèi)星快速精密定軌后，采用動力學擬合方法進行軌道預(yù)報，可以獲得與采用事后精密星歷進行軌道預(yù)報接近的精度，預(yù)報12 h弧長軌道時，3DRMS優(yōu)于6 dm；預(yù)報24 h弧長軌道時，3DRMS優(yōu)于1 m，可以滿足較高精度的應(yīng)用需求。

(3) LEO衛(wèi)星軌道預(yù)報精度很大程度上取決于S方向，且隨著預(yù)報弧長增加，其預(yù)報精度下降很快。因此提高S方向的預(yù)報精度，有利于提高低軌衛(wèi)星的預(yù)報精度。

致謝：感謝中國科學院測量與地球物理研究所iGMAS分析中心提供的GPS超快速產(chǎn)品。

[1]鄧琳, 李廣俠, 田世偉, 等. 基于LEO增強的COMPASS導(dǎo)航系統(tǒng)抗干擾能力研究[J]. 軍事通信技術(shù), 2012, 33(2): 65-69, 80. DENG Lin, LI Guangxia, TIAN Shiwei, et al. Capabilities of LEO Enhanced COMPASS Navigation System for Anti-jamming[J]. Journal of Military Communication Technology, 2012, 33(2): 65-69, 80.

[2]RABINOWITZ M, PARKSINSON B W, SPILKER J J. Some Capabilities of a Joint GPS-LEO Navigation System[C]∥Proceedings of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation. Salt Lake City: ION-GPS, 2000: 255-256.

[3]趙金峰, 于笑, 馮少棟, 等. 基于LEO增強的COMPASS系統(tǒng)設(shè)計及性能分析[J]. 電訊技術(shù), 2013, 53(2): 131-135.

ZHAO Jinfeng, YU Xiao, FENG Shaodong, et al. Design and Performance Analysis of LEO Satellites Enhanced COMPASS System[J]. Telecommunication Engineering, 2013, 53(2): 131-135.

[4]田世偉. 從“傳統(tǒng)”到“協(xié)同: 衛(wèi)星導(dǎo)航增強理論與方法研究[D]. 南京: 中國人民解放軍理工大學, 2015.

TIAN Shiwei. From “Conventional” to “Cooperative”: Research on Theories and Methods on Satellite Navigation Augmentation[D]. Nanjing: PLA University of Science and Technology, 2015.

[5]劉林, 王彥榮. 衛(wèi)星軌道預(yù)報的一種分析方法[J]. 天文學報, 2005, 46(3): 307-313. LIU Lin, WANG Yanrong. An Analytical Method for Satellite Orbit Forecasting[J]. Acta Astronomica Sinica, 2005, 46(3): 307-313.

[6]柳青松, 劉峰, 巨濤. 基于GPS的低軌衛(wèi)星簡化解析法軌道預(yù)報[J]. 全球定位系統(tǒng), 2010, 35(2): 1-4. LIU Qingsong, LIU Feng, JU Tao. A Simplified Analytical Method of LEO Satellite Orbit Forecasting Based on GPS[J]. GNSS World of China, 2010, 35(2): 1-4.

[7]魏二虎, 柴華. GPS精密星歷插值方法的比較研究[J]. 全球定位系統(tǒng), 2006, 31(5): 13-15, 20.

WEI Erhu, CHAI Hua. The Comparison Research on GPS Precise Ephemeris Interpolation Methods[J]. GNSS World of China, 2006, 31(5): 13-15, 20.

[8]張守建, 李建成, 邢樂林, 等. 兩種IGS精密星歷插值方法的比較分析[J]. 大地測量與地球動力學, 2007, 27(2): 80-83.

ZHANG Shoujian, LI Jiancheng, XING Lelin, et al. Comparative Analysis on Two Methods for IGS Precise Ephemeris Interpolation[J]. Journal of Geodesy and Geodynamics, 2007, 27(2): 80-83.

[9]高鵬, 喬學軍, 范城城. HY-2衛(wèi)星精密軌道擬合與外推的兩種方法比較[J]. 海洋測繪, 2013, 33(4): 58-61. GAO Peng, QIAO Xuejun, FAN Chengcheng. Comparison of Two Methods of HY-2 Satellite Precision Orbit Fitting and Extrapolation[J]. Hydrographic Surveying and Charting, 2013, 33(4): 58-61.

[10]張如偉, 劉根友. 低軌衛(wèi)星軌道擬合及預(yù)報方法研究[J]. 大地測量與地球動力學, 2008, 28(4): 115-120.

ZHANG Ruwei, LIU Genyou. Discussion on Orbit Fitting and Orbit Forecasting of Low Earth Orbit Satellites[J]. Journal of Geodesy and Geodynamics, 2008, 28(4): 115-120.

[11]JOERGER M, GRATTON L, PERVAN B, et al. Analysis of Iridium-Augmented GPS for Floating Carrier Phase Positioning[J]. Navigation, 2010, 57(2): 137-160.

[12]BEUTLER G. Methods of Celestial Mechanics: Physical, Mathematical, and Numerical Principles[M]. Berlin Heidelberg: Springer, 2005.

[13]李征航, 龔曉穎, 劉萬科. 配置法在衛(wèi)星軌道積分中的應(yīng)用[J]. 武漢大學學報(信息科學版), 2010, 35(3): 253-256.

LI Zhenghang, GONG Xiaoying, LIU Wanke. Application of Collocation Method to Numerical Integration of Satellite Orbit[J]. Geomatics and Information Science of Wuhan University, 2010, 35(3): 253-256.

[14]BEUTLER G, BOCK H, DACH R, et al. Bernese GPS Software Version 5.0[M]. Bern: Astronomical Institute, University of Bern, 2004.

[15]龔曉穎, 蔡詩響, 李征航. 小步長低階Collocation方法軌道積分性能研究[J]. 宇航學報, 2013, 34(5): 642-648. GONG Xiaoying, CAI Shixiang, LI Zhenghang. Performance Research of Small-Steplength and Low-Order Collocation Method in Orbit Determination[J]. Journal of Astronautics, 2013, 34(5): 642-648.

[16]宋小勇, 毛悅, 賈小林. BERNESE光壓模型參數(shù)的統(tǒng)計分析[J]. 測繪科學, 2009, 34(3): 25-27.

SONG Xiaoyong, MAO Yue, JIA Xiaolin. Analysis of the Parameter of BERNESE Radiation Force Model by Statistic Method[J]. Science of Surveying and Mapping, 2009, 34(3): 25-27.

[17]陳俊平, 王解先. GPS定軌中的太陽輻射壓模型[J]. 天文學報, 2006, 47(3): 310-319.

CHEN Junping, WANG Jiexian. Solar Radiation Pressure Models for the GPS Satellites[J]. Acta Astronomica Sinica, 2006, 47(3): 310-319.

[18]何仕強, 吳斌, 陳俊平. 利用IGS星歷預(yù)報GPS衛(wèi)星軌道[J]. 中國科學院上海天文臺年刊, 2011(32): 25-34.

HE Shiqiang, WU Bin, CHEN Junping. Orbit Predictions for GPS Satellites Using IGS Ephemeris[J]. Annals of Shanghai Observation Academia Sinica, 2011(32): 25-34.

[19]李濟生. 人造衛(wèi)星精密軌道確定[M]. 北京: 解放軍出版社, 1995.

LI Jisheng. Satellite Precise Orbit Determination[M]. Beijing: The PLA Press, 1995.

[20]趙齊樂. GPS導(dǎo)航星座及低軌衛(wèi)星的精密定軌理論和軟件研究[D]. 武漢: 武漢大學, 2004.

ZHAO Qile. Research on Precise Orbit Determination Theory and Software of GPS Navigation Constellation and LEO Satellite[D]. Wuhan: Wuhan University, 2004.

[21]周建, 呂志偉. 地球重力場模型對低軌衛(wèi)星軌道積分的影響[J]. 測繪工程, 2012, 21(2): 14-16.

ZHOU Jian, Lü Zhiwei. The Effects of Earth’s Gravity Field Models to LEO Orbit Integration[J]. Engineering of Surveying and Mapping, 2012, 21(2): 14-16.

[22]韓保民, 朱秀英, 柳林濤, 等. 偽隨機脈沖估計及其在簡化動力學定軌中的應(yīng)用[J]. 武漢大學學報(信息科學版), 2007, 32(5): 466-469. HAN Baomin, ZHU Xiuying, LIU Lintao, et al. Estimation of Pseudo-Stochastic Pulses and Their Applications in Reduced-Dynamic Orbit Determination[J]. Geomatics and Information Science of Wuhan University, 2007, 32(5): 466-469.

[23]ZHAO Gang, ZHOU Xuhua, WU Bin. Precise orbit Determination of Haiyang-2 Using Satellite Laser Ranging[J]. Chinese Science Bulletin, 2013, 58(6): 589-597.

[24]郭靖, 趙齊樂, 李敏, 等. 利用星載GPS觀測數(shù)據(jù)確定海洋2A衛(wèi)星cm級精密軌道[J]. 武漢大學學報(信息科學版), 2013, 38(1): 52-55.

GUO Jing, ZHAO Qile, LI Min, et al. Centimeter Level Orbit Determination for HY2A Using GPS Data[J]. Geomatics and Information Science of Wuhan University, 2013, 38(1): 52-55.

[25]周建華, 楊龍, 徐波, 等. 一種導(dǎo)航衛(wèi)星中長期軌道預(yù)報方法[J]. 測繪學報, 2011, 40(S1): 39-45. ZHOU Jianhua, YANG Long, XU Bo, et al. A Long Duration Orbit Prediction Method for Navigation Satellite[J]. Acta Geodaetica et Catographica Sinica, 2011, 40(S1): 39-45.

[27]郭南男, 周旭華, 吳斌. 利用星載GPS數(shù)據(jù)進行海洋2A衛(wèi)星快速精密定軌[J]. 宇航學報, 2015, 36(7): 797-803.

GUO Nannan, ZHOU Xuhua, WU Bin. Rapid Precise Orbit Determination for Haiyang-2A Using On-board GPS Data[J]. Journal of Astronautics, 2015, 36(7): 797-803.

(責任編輯：陳品馨)

Precision Analysis of LEO Satellite Orbit Prediction

WANG Yafei1,2,3,ZHONG Shiming1,2,WANG Haitao1,2,OU Jikun1,2

1. Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China; 2. State Key Laboratory of Geodesy and Earth’s Dynamics, Wuhan 430077, China; 3. University of Chinese Academy of Science, Beijing 100049, China

The accuracy of different fitting intervals on the predicted orbit of different arc lengths was analyzed by using the dynamic fitting method and taken HY-2A satellite as an example. According to two orbit products, one was the precise orbit obtained by CNES and the other was rapid orbit computed by adopting zero-different reduced dynamic method, obtaining the precise predicted orbit. By the results, we can get that the 3DRMS is near 6 dm by using 12 h or 24 h fitting orbit to forecast 12 h arc; the 3DRMS is near 1 m when forecasting 24 h arc.

dynamic fitting; Collocation algorithm; fitting interval; prediction; rapid orbit determination

The National Natural Science Foundation of China(Nos. 41474029；41574015)；State Key Laboratory of Geodesy and Earth’s Dynamic Project(No. SKLGED2013-4-2-Z)

WANG Yafei(1989—), female, postgraduate，majors in GNSS high-precise data processing.

ZHONG Shiming

王亞菲，鐘世明，王海濤，等.LEO衛(wèi)星軌道預(yù)報精度分析[J].測繪學報，2016,45(9)：1035-1041.

10.11947/j.AGCS.2016.20160045.

WANG Yafei, ZHONG Shiming, WANG Haitao, et al.Precision Analysis of LEO Satellite Orbit Prediction[J]. Acta Geodaetica et Cartographica Sinica,2016,45(9):1035-1041. DOI:10.11947/j.AGCS.2016.20160045.

P228

A

1001-1595(2016)09-1035-07

國家自然科學基金(41474029；41574015)；大地測量與地球動力學國家重點實驗室自主項目(SKLGED2013-4-2-Z)

2016-01-28

王亞菲(1989—)，女，碩士，研究方向為GNSS高精度數(shù)據(jù)處理。

E-mail： wangyafei14@mails.ucas.ac.cn

鐘世明

E-mail： smzhong@whigg.ac.cn

修回日期： 2016-06-17