姚宜斌,熊朝暉,張 豹,張 良,孔 建
1.武漢大學測繪學院,湖北 武漢 430079;2.武漢大學地球空間環(huán)境與大地測量教育部重點實驗室,湖北 武漢 430079;3.地球空間信息技術協(xié)同創(chuàng)新中心,湖北 武漢 430079;4.武漢大學中國南極測繪研究中心,湖北 武漢 430079
顧及設計矩陣誤差的AR模型新解法
姚宜斌1,2,3,熊朝暉1,張 豹1,張 良1,孔 建4
1.武漢大學測繪學院,湖北 武漢 430079;2.武漢大學地球空間環(huán)境與大地測量教育部重點實驗室,湖北 武漢 430079;3.地球空間信息技術協(xié)同創(chuàng)新中心,湖北 武漢 430079;4.武漢大學中國南極測繪研究中心,湖北 武漢 430079
在自回歸模型求解中,設計矩陣和觀測值均存在誤差,傳統(tǒng)的最小二乘法不能很好地解決這一問題。本文提出了一種顧及設計矩陣誤差的AR模型新解法,通過引入虛擬觀測值,使觀測向量與設計矩陣不僅同源而且?guī)д`差的元素個數(shù)相同,然后通過對觀測方程進行等價變換巧妙實現(xiàn)了在最小二乘框架下求解自回歸問題。利用模擬數(shù)據(jù)及實測數(shù)據(jù)分別對新算法進行了內(nèi)符合精度檢驗,并利用實測數(shù)據(jù)對新算法進行外符合精度檢驗,結果表明新算法得到的結果顯著優(yōu)于奇異值分解(singular value decomposition,SVD)解法及傳統(tǒng)最小二乘解法,驗證了算法的精度和有效性。
AR模型;設計矩陣誤差;整體最小二乘;虛擬觀測值;奇異值分解
平差的實際問題中常存在設計矩陣由含有誤差的觀測值構成,傳統(tǒng)最小二乘方法在求解此類模型的過程中,因忽略設計矩陣中實際存在的誤差而具有理論的不完備性。考慮設計矩陣誤差,將所有觀測值改正數(shù)平方和最小作為平差準則的求解方法最初由Adcock于1877年提出[1]。1980年,Golub和Van Loan提出了奇異值分解(SVD)算法,并將這類設計矩陣帶有誤差的問題命名為整體最小二乘方法(total least square,TLS)[2]。使用SVD解法雖能獲得穩(wěn)定數(shù)值解,但計算復雜度高[3],如果設計矩陣中含有誤差的元素重復出現(xiàn),使用SVD解法,同一觀測值在不同位置改正數(shù)將不一樣。另外,使用SVD解法估計參數(shù),難以給出以觀測值表達的未知參數(shù)平差值表達式,從而無法利用協(xié)方差傳播率進行精度評定[4]。近些年,文獻[4]提出一種顧及設計矩陣隨機誤差的最小二乘組合新解法,其將設計矩陣中所有含誤差的元素作為虛擬觀測值引入虛擬觀測方程求解模型,解決了TLS方法的精度評定問題。加權整體最小二乘亦有較大進展[5-10]。此外,文獻[11]研究了三維坐標轉換的通用整體最小二乘解法;文獻[12]對部分變量誤差模型提出了整體抗差最小二乘估計方法;文獻[13]對Partial EIV模型提出了一種新解法;文獻[14]提出了一種基于中位數(shù)法的抗差整體最小二乘估計方法。
AR模型是一類基本的時間序列模型,在工程實踐與導航定位等方面有著廣泛應用[15-22]。對于t階AR模型,任意時刻觀測值zm為自身最近t階滯后項的線性組合[23],如下
(1)
式(1)用矩陣形式表達可寫為
L=BX
(2)
由此可見,含誤差的觀測值規(guī)律地出現(xiàn)在設計矩陣不同位置,且觀測向量中元素亦重復出現(xiàn)在設計矩陣中。使用傳統(tǒng)最小二乘方法估計AR模型參數(shù),因忽略設計矩陣中存在含誤差的元素而理論缺乏嚴密性。AR模型中,觀測向量誤差與設計矩陣誤差同源,若使用SVD方法估計AR模型參數(shù),則多次出現(xiàn)的同一觀測值的改正數(shù)不再唯一。文獻[24]提出一種非線性的AR模型整體最小二乘迭代算法,但難以應用協(xié)方差傳播定律給出精度評定公式。
前述方法對求解AR模型存在一定缺陷,本文提出一種顧及設計矩陣誤差,適用于AR模型的整體最小二乘新解法,該方法通過引入未在觀測向量中出現(xiàn)且含誤差的觀測值作為虛擬觀測值,將設計矩陣對應的改正數(shù)與未知參數(shù)初值乘積改寫為由新設參數(shù)改正數(shù)及實際觀測值改正數(shù)表達的線性組合,使得觀測向量與設計矩陣中帶誤差的元素個數(shù)相同。由此將整體最小二乘問題轉化為經(jīng)典最小二乘問題,最終使用間接平差方法進行參數(shù)求解。本文方法有效地克服了AR模型中同一參數(shù)在不同位置改正數(shù)不一致的情況,并能依據(jù)協(xié)方差傳播定律進行精度評定。
1.1 算法推導
同時考慮設計矩陣B與觀測向量L誤差,式(2)可寫為
L+v=(B0+ΔB)(X0+x)=B0X0+
B0x+ΔBX0+ΔBx
(3)
忽略二階小項ΔBx,式(3)可改寫為
v=B0X0+B0x+ΔBX0-L
(4)
1.1.1 矩陣等價變換
矩陣ΔB中改正數(shù)規(guī)律地重復出現(xiàn),為將不同的改正數(shù)單獨取出,同時考慮未在設計矩陣中出現(xiàn)的改正數(shù)vt+n,將ΔBX0作如式(5)所示等價變換
(5)
1.1.2 引入虛擬觀測值
(6)
v=B0X0+B0x+B10xB+B20v-L
(7)
1.1.3 函數(shù)模型構建
將式(7)中右邊B20v項移至等式左邊,提取公因式v得(E-B20)v。由于矩陣B20中非零元素集中出現(xiàn)在主對角線下方,主對角線及其上方元素全為0,故行列式|E-B20|恒為1,即矩陣E-B20可逆,從而可得式(8)
v=(E-B20)-1B0x+(E-B20)-1B10xB-(E-B20)-1(L-B0X0)
(8)
虛擬觀測值誤差方程為
vB=xB
(9)
如果令
則有
vg=Bgxg-lg
(10)
1.2 解算流程
本文方法解算流程為:
(3) 根據(jù)式(10)計算vg,更新實際觀測值L估值及虛擬觀測值LB估值。
(4) 重復步驟(2)、(3)直到未知參數(shù)改正數(shù)x小于一定閾值停止迭代,得到最終的未知參數(shù)結果。
1.3 精度評定
在經(jīng)典最小二乘中,權由觀測值精度確定,而觀測值出現(xiàn)的次數(shù),則會在平差過程中,通過對系數(shù)矩陣的作用,進一步影響法方程進而影響平差結果。但是在整體最小二乘中,系數(shù)矩陣中的觀測值具有雙重身份,除了影響系數(shù)矩陣本身外,還作為觀測值參與平差,本文認為觀測值這部分影響與觀測值在系數(shù)矩陣中出現(xiàn)的次數(shù)也有關系。本文所定義的權陣P由兩部分構成,即由觀測值精度信息構建的權陣P′,以及一個與次數(shù)有關的因子陣K,其中P=P′K。為簡化問題,本文先行假定AR模型各期觀測值等精度,即P′為單位陣,P=K,若觀測值不等精度,則在按照精度信息得到的權陣P′的基礎上再乘以一個次數(shù)有關的因子陣K。
由于在AR模型中,設計矩陣中元素在不同位置中規(guī)律地重復出現(xiàn),觀測值亦重復出現(xiàn)于設計矩陣中,在使用整體最小二乘方法進行求解未知參數(shù)時,實際上將同一觀測值當作多次相同觀測建模,為此,在各期觀測值等精度的條件下,模型各觀測值權可由觀測值出現(xiàn)的次數(shù)確定,考慮引入虛擬觀測值,該算法權陣P定義如下
t+1,t,t-1,…,2,1)
(11)
使用上述方法經(jīng)過若干次迭代可以獲得未知參數(shù)X平差值,單位權中誤差σ計算如式(12)所示
(12)
式中,模型自由度f=(n+u)-(t+u)=n-t。
在本文算法中,參數(shù)X0=[φ10φ20…φt0]T的初值為最小二乘所得到的AR系數(shù)。需要注意的是,若初值精度不夠,則會導致迭代次數(shù)多,收斂速度慢,初值精度過差會導致迭代不能收斂。
2.1 模擬數(shù)據(jù)
為驗證算法的可行性與精度,依據(jù)式(13)模擬一個平穩(wěn)的AR(3)模型,其AR系數(shù)分別取為φ1=0.8、φ2=-0.5、φ3=-0.3,模擬數(shù)據(jù)如表1所示。
表1 AR(3)模型模擬數(shù)據(jù)Tab.1 Simulation date of AR(3) model
zt=φ1zt-1+φ2zt-2+φ3zt-3+et
(13)
式中,et為高斯白噪聲,表1模擬數(shù)據(jù)所附加噪聲et~N(0,1)。
對于表1模擬數(shù)據(jù),分別采用本文方法(迭代閾值為0.000 000 1)、SVD方法以及直接最小二乘3種方法求解參數(shù)及單位權中誤差,列于表2。
表23種不同方法所求參數(shù)值及單位權中誤差
Tab.2Theparametervalueandtheerrorofunitweightofthreemethods
參數(shù)SVDLS本文方法真值φ10.83630.77090.80430.8φ2-0.5461-0.4667-0.5084-0.5φ3-0.2725-0.3246-0.3068-0.3σ01.44551.41600.9042
從表2可知,本文方法所得到的參數(shù)估值較SVD方法和傳統(tǒng)最小二乘方法更接近于真值,由于考慮了設計矩陣誤差及同一觀測值出現(xiàn)在不同位置改正數(shù)應一致等事實,本文方法的單位權中誤差明顯小于最小二乘方法和SVD方法,并且更加接近虛擬觀測噪聲的中誤差。當對模擬數(shù)據(jù)不附加高斯白噪聲時,本文方法所得參數(shù)估值將與真值一致。
2.2 實測沉降數(shù)據(jù)
本文算法、普通最小二乘(LS)以及SVD 3種方法所得參數(shù)估值如表3所示。實測觀測數(shù)據(jù)及3種方法所得平差值如表4所示。圖1所示為本文算法、普通最小二乘(LS)以及SVD 3種方法平差結果與原始觀測值殘差絕對值對比。
表3中本文方法所求結果與SVD法較為接近,但本文方法所得單位權中誤差小于SVD法,直接最小二乘結果與本文結果相差較遠,查閱文獻[23]知,本文所引用實測沉降數(shù)據(jù)并非平穩(wěn)的自回歸時間序列,這可能是和直接最小二乘結果估值結果相差較大的原因。由表4及圖1可知,本文方法所得觀測值的平差值更契合原始觀測序列變化趨勢。新方法看起來誤差相對平穩(wěn),波動小一些。第14期的結果,其他方法較差,新方法較好,說明新方法由于理論嚴密,所以對于抑制較大誤差效果更好。
表33種方法所求沉降數(shù)據(jù)參數(shù)估值
Tab.3Theparametervalueofobservationdataofsettlementofthreemethods
參數(shù)SVDLS本文方法φ1-0.37480.0411-0.4863φ20.28980.32780.4264φ31.09010.63511.0670σ00.89540.76730.7320
表4 沉降數(shù)據(jù)實測值及3種方法平差值Tab.4 Adjustment of observation of three method and the original data mm
圖1 3種方法平差結果與原始觀測值比較Fig.1 Comparative results between original data and adjustment of observation of three value
2.3 預報精度分析
AR模型能利用前若干期觀測值預報后若干期觀測值,為檢驗本文算法的外符合精度,以表4所列前30期實測沉降觀測值為基礎,分別利用本文算法、普通最小二乘(LS)以及SVD 3種方法求解模型系數(shù)并預報第31期至36期沉降數(shù)據(jù)。依據(jù)前30期實測沉降數(shù)據(jù)所求參數(shù)結果及相應方法預報值如表5所示。圖2所示為本文算法、普通最小二乘(LS)以及SVD 3種方法預報結果與實測沉降觀測值殘差絕對值對比。
表5 3種方法預報結果Tab.5 Forecast result of three method
圖2 3種方法預測結果與實測沉降觀測值比較Fig.2 Comparative results between original data and forecast result of observation of three value
由表5及圖2可知,在第31期至36期數(shù)據(jù)中,由SVD方法所求參數(shù)的預測效果最差,本文算法整體上預測效果最優(yōu),在第34期預報中3種方法預測效果均差,可能是觀測噪聲較大所致。向后預報期數(shù)越小,本文算法所求系數(shù)的預報結果與普通最小二乘方法相比,更接近與實際觀測值,而在較遠期預報中,二者效果相當??傮w而言,在本文算法、普通最小二乘(LS)以及SVD 3種方法中,本文算法的外符合精度最優(yōu)。
在AR模型中,觀測向量與設計矩陣的誤差同源,觀測值規(guī)律地重復出現(xiàn)在設計矩陣中,且設計矩陣中自身元素亦重復出現(xiàn)。傳統(tǒng)最小二乘方法難以解決系數(shù)矩陣及觀測值向量皆帶誤差的數(shù)據(jù)處理問題。本文提出了一種考慮設計矩陣誤差的AR模型整體最小二乘新解法,引入僅在設計矩陣中出現(xiàn)且含誤差的元素作為虛擬觀測值,使觀測向量與設計矩陣中帶誤差的元素個數(shù)相同。然后巧妙地對觀測方程進行等價變換,實現(xiàn)了最小二乘框架下求解整體最小二乘問題,有效地克服了傳統(tǒng)SVD方法的理論缺陷且能應用協(xié)方差傳播定律給出未知參數(shù)的精度評定公式。最后通過對模擬數(shù)據(jù)及實測數(shù)據(jù)的驗證,發(fā)現(xiàn)本文方法具有比SVD方法及經(jīng)典最小二乘方法更高的精度及更優(yōu)的外符合精度。
[1] ADCOCK R J.Note on the Method of Least Squares[J].Analyst,1877,4(6):183-184.
[2] GOLUB G H,VAN LOAN C F.An Analysis of the Total Least Squares Problem[J].SIAM Journal on Numerical Analysis,1980,17(6):883-893.
[3] VAN HUFFEL S,ZHA Hongyuan.An Efficient Total Least Squares Algorithm Based on a Rank-revealing Two-sided Orthogonal Decomposition[J].Numerical Algorithms,1993,4(1):101-133.
[4] 姚宜斌,孔建.顧及設計矩陣隨機誤差的最小二乘組合新解法[J].武漢大學學報(信息科學版),2014,39(9):1028-1032.
YAO Yibin,KONG Jian.A New Combined LS Method Considering Random Errors of Design Matrix[J].Geomatics and Information Science of Wuhan University,2014,39(9):1028-1032.
[5] 曾文憲,方興,劉經(jīng)南,等.附有不等式約束的加權整體最小二乘算法[J].測繪學報,2014,43(10):1013-1018.DOI:10.13485/j.cnki.11-2089.2014.0173.
ZENG Wenxian,F(xiàn)ANG Xing,LIU Jingnan,et al.Weighted Total Least Squares Algorithm with Inequality Constraints[J].Acta Geodaetica et Cartographica Sinica,2014,43(10):1013-1018.DOI:10.13485/j.cnki.11-2089.2014.0173.
[6] SCHAFFRIN B,WIESER A.On Weighted Total Least-squares Adjustment for Linear Regression[J].Journal of Geodesy,2008,82(7):415-421.
[7] SCHUERMANS M,MARKOVSKY I,VAN HUFFEL S.An Adapted Version of the Element-wise Weighted Total Least Squares Method for Applications in Chemometrics[J].Chemometrics and Intelligent Laboratory Systems,2007,85(1):40-46.
[8] VAN HUFFEL S,VANDEWALLE J.Analysis and Properties of the Generalized Total Least Squares Problem AX≈B When Some or All Columns in A are Subject to Error[J].SIAM Journal on Matrix Analysis and Applications,1989,10(3):294-315.
[9] XU Peiliang,LIU Jingnan,SHI Chuang.Total Least Squares Adjustment in Partial Errors-in-variables Models:Algorithm and Statistical Analysis[J].Journal of Geodesy,2012,86(8):661-675.
[10] FANG X.Weighted Total Least Squares Solutions for Applications in Geodesy[D].Hannover,Germany:Leibniz University,2011.
[11] 方興,曾文憲,劉經(jīng)南,等.三維坐標轉換的通用整體最小二乘算法[J].測繪學報,2014,43(11):1139-1143.DOI:10.13485/j.cnki.11-2089.2014.0193.
FANG Xing,ZENG Wenxian,LIU Jingnan,et al.A General Total Least Squares Algorithm for Three-dimensional Coordinate Transformations[J].Acta Geodaetica et Cartographica Sinica,2014,43(11):1139-1143.DOI:10.13485/j.cnki.11-2089.2014.0193.
[12] 趙俊,歸慶明.部分變量誤差模型的整體抗差最小二乘估計[J].測繪學報,2016,45(5):552-559.DOI:10.11947/j.AGCS.2016.20150374.
ZHAO Jun,GUI Qingming.Total Robustified Least Squares Estimation in Partial Errors-in-variables Model[J].Acta Geodaetica et Cartographica Sinica,2016,45(5):552-559.DOI:10.11947/j.AGCS.2016.20150374.
[13] 王樂洋,余航,陳曉勇.Partial EIV模型的解法[J].測繪學報,2016,45(1):22-29.DOI:10.11947/j.AGCS.2016.20140560.
WANG Leyang,YU Hang,CHEN Xiaoyong.An Algorithm for Partial EIV Model[J].Acta Geodaetica et Cartographica Sinica,2016,45(1):22-29.DOI:10.11947/j.AGCS.2016.20140560.
[14] 陶葉青,高井祥,姚一飛.基于中位數(shù)法的抗差總體最小二乘估計[J].測繪學報,2016,45(3):297-301.DOI:10.11947/j.AGCS.2016.20150234.
TAO Yeqing,GAO Jingxiang,YAO Yifei.Solution for Robust Total Least Squares Estimation Based on Median Method[J].Acta Geodaetica et Cartographica Sinica,2016,45(3):297-301.DOI:10.11947/j.AGCS.2016.20150234.
[15] 吳富梅,楊元喜.基于高階AR模型的陀螺隨機漂移模型[J].測繪學報,2007,36(4):389-394.
WU Fumei,YANG Yuanxi.Gyroscope Random Drift Model Based on the Higher-order AR Model[J].Acta Geodaetica et Cartographica Sinica,2007,36(4):389-394.
[16] 潘國榮,劉大杰.顧及鄰近點變形因素項的動態(tài)模型辨識及預測[J].測繪學報,2001,30(1):32-35.
PAN Guorong,LIU Dajie.Dynamic Modeling Identification and Predication in Consideration of the Adjacent Point Deformation[J].Acta Geodaetica et Cartographica Sinica,2001,30(1):32-35.
[17] 楊元喜,崔先強.動態(tài)定位有色噪聲影響函數(shù)——以一階AR模型為例[J].測繪學報,2003,32(1):6-10.
YANG Yuanxi,CUI Xianqiang.Influence Functions of Colored Noises on Kinematic Positioning:Taking the AR Model of First Class as an Example[J].Acta Geodaetica et Cartographica Sinica,2003,32(1):6-10.
[18] 葉志偉,尹暉,張守建.AR模型譜在超導重力數(shù)據(jù)信號檢測中的分析研究[J].武漢大學學報(信息科學版),2007,32(6):536-539.
YE Zhiwei,YIN Hui,ZHANG Shoujian.Using AR Model Spectrum Algorithms to Detect Superconducting Gravimetric Signals[J].Geomatics and Information Science of Wuhan University,2007,32(6):536-539.
[19] 張昊,王琪潔,朱建軍,等.對錢德勒參數(shù)進行時變修正的CLS+AR模型在極移預測中的應用[J].武漢大學學報(信息科學版),2012,37(3):286-289.
ZHANG Hao,WANG Qijie,ZHU Jianjun,et al.Application of CLS+AR Model Polar Motion to Prediction Based on Time-varying Parameters Correction of Chandler Wobble[J].Geomatics and Information Science of Wuhan University,2012,37(3):286-289.
[20] 王樂洋,許才軍,魯鐵定.邊長變化反演應變參數(shù)的總體最小二乘方法[J].武漢大學學報(信息科學版),2010,35(2):181-184.
WANG Leyang,XU Caijun,LU Tieding.Inversion of Strain Parameter Using Distance Changes Based on Total Least Squares[J].Geomatics and Information Science of Wuhan University,2010,35(2):181-184.
[21] 魏二虎,殷志祥,李廣文,等.虛擬觀測值法在三維坐標轉換中的應用研究[J].武漢大學學報(信息科學版),2014,39(2):152-156.
WEI Erhu,YIN Zhixiang,LI Guangwen,et al.On 3D Coordinate Transformations with Virtual Observation Method[J].Geomatics and Information Science of Wuhan University,2014,39(2):152-156.
[22] 姚宜斌,黃書華,孔建,等.空間直線擬合的整體最小二乘算法[J].武漢大學學報(信息科學版),2014,39(5):571-574.
YAO Yibin,HUANG Shuhua,KONG Jian,et al.Total Least Squares Algorithm for Fitting Spatial Straight Lines[J].Geomatics and Information Science of Wuhan University,2014,39(5):571-574.
[23] CRYER J D,CHAN K S.時間序列分析及應用[M].潘紅宇,譯.北京:機械工業(yè)出版社,2011.
CRYER J D,CHAN K S.Time Series Analysis with Applications in R[M].PAN Hongyu,tran.Beijing:China Machine Press,2011.
[24] 姚宜斌,黃書華,陳家君.求解自回歸模型參數(shù)的整體最小二乘新方法[J].武漢大學學報(信息科學版),2014,39(12):1463-1466.
YAO Yibin,HUANG Shuhua,CHEN Jiajun.A New Method of TLS to Solving the Autoregressive Model Parameter[J].Geomatics and Information Science of Wuhan University,2014,39(12):1463-1466.
[25] 王新洲,陶本藻,邱衛(wèi)寧,等.高等測量平差[M].北京:測繪出版社,2013.
WANG Xinzhou,TAO Benzao,QIU Weining,et al.Advanced Surveying Adjustment[M].Beijing:Mapping Publishing Company,2013.
A New Method to Solving AR Model Parameters Considering Random Errors of Design Matrix
YAO Yibin1,2,3,XIONG Zhaohui1,ZHANG Bao1,ZHANG Liang1,KONG Jian4
1.School of Geodesy and Geomatics,Wuhan University,Wuhan 430079,China; 2.Key Laboratory of Geospace Environment and Geodesy,Ministry of Education,Wuhan University,Wuhan 430079,China; 3.Collaborative Innovation Center for Geospatial Technology,Wuhan 430079,China; 4.Chinese Antarctic Center of Surveying and Mapping,Wuhan 430079,China
The ordinary least square method could not solve the problem that the error exist both in design matrix and observation vector while compute parameter values of AR model.In this article, a new method is proposed which consider the random errors of design matrix.The source of design matrix and observation vector is same and the amount of parameters contain error can be equal by introducing virtual observation.Then, this problem could be solved under the framework of normal least square by equivalence transformation of observation equation.The result of this new method is superior to SVD method and normal least square method by simulation date and observation data which verify the feasibility and effectiveness of this method.
AR model;design matrix error;TLS;virtual observations;SVD method
The General Program of National Natural Science Foundation of China(Nos.41274022;41574028);Natural Science Foundation for Distinguished Young Scholars of Hubei Province of China(No.2015CFA036)
YAO Yibin(1976—),male,professor,majors in geodetic data processing,GNSS space environment science,etc.
XIONG Zhaohui
姚宜斌,熊朝暉,張豹,等.顧及設計矩陣誤差的AR模型新解法[J].測繪學報,2017,46(11):1795-1801.
10.11947/j.AGCS.2017.20170004.
YAO Yibin,XIONG Zhaohui,ZHANG Bao,et al.A New Method to Solving AR Model Parameters Considering Random Errors of Design Matrix[J].Acta Geodaetica et Cartographica Sinica,2017,46(11):1795-1801.DOI:10.11947/j.AGCS.2017.20170004.
P228
A
1001-1595(2017)11-1795-07
國家自然科學基金(41274022;41574028);湖北省杰出青年科學基金(2015CFA036)
(責任編輯:宋啟凡)
2017-01-03
修回日期:2017-08-18
姚宜斌(1976—),男,教授,主要從事測量數(shù)據(jù)處理理論與方法、GNSS空間環(huán)境學研究。
E-mail:ybyao@sgg.whu.edu.cn
熊朝暉
E-mail:cehui_xiong@whu.edu.cn