文超斌，王躍鋼田 琦郭志斌左朝陽(yáng)楊家勝

（1.第二炮兵工程大學(xué) 304教研室，西安 710025；2.中國(guó)人民解放軍96124部隊(duì)，通化 134000）

局部連續(xù)場(chǎng)下重力輔助導(dǎo)航模型構(gòu)建算法

文超斌1，2，王躍鋼1，田 琦1，郭志斌1，左朝陽(yáng)1，楊家勝1

（1.第二炮兵工程大學(xué) 304教研室，西安 710025；2.中國(guó)人民解放軍96124部隊(duì)，通化 134000）

為了研究高精度的重力輔助導(dǎo)航模型以克服傳統(tǒng)模型的局限，必須建立精度高且具有良好解析性質(zhì)的局部重力異常場(chǎng)解析模型，同時(shí)考慮模型誤差方程中的重力精準(zhǔn)補(bǔ)償問(wèn)題。針對(duì)二維高斯樣條函數(shù)逼近局部重力異常場(chǎng)中的局部支撐參數(shù)選擇問(wèn)題，通過(guò)對(duì)所涉及的系數(shù)矩陣、解誤差、插值模型精度評(píng)估等問(wèn)題進(jìn)行分析，提出了一種新的最優(yōu)局部支撐參數(shù)計(jì)算方法；基于此提出了一種高精度的基于高斯插值的重力輔助導(dǎo)航模型構(gòu)建算法，該精準(zhǔn)模型補(bǔ)償了重力擾動(dòng)矢量、標(biāo)準(zhǔn)重力值誤差、厄特弗斯修正計(jì)算值對(duì)導(dǎo)航模型的影響。實(shí)驗(yàn)結(jié)果表明利用新型重力輔助導(dǎo)航模型構(gòu)建算法，可使輔助導(dǎo)航系統(tǒng)位置精度提高1倍左右，姿態(tài)、速度精度提高1～2倍，定位誤差保持在100～200 m。

重力輔助導(dǎo)航；高斯樣條插值；重力誤差補(bǔ)償；導(dǎo)航模型

重力匹配輔助慣性導(dǎo)航作為一種利用地球重力場(chǎng)特征信息結(jié)合匹配算法對(duì)對(duì)慣性導(dǎo)航系統(tǒng)位置信息進(jìn)行初次估計(jì)，而后采用導(dǎo)航數(shù)學(xué)模型方程結(jié)合一定的濾波算法對(duì)慣性導(dǎo)航系統(tǒng)位置信息進(jìn)行進(jìn)一步精準(zhǔn)修正的新技術(shù)，近年來(lái)得到了快速的發(fā)展。目前關(guān)于重力匹配算法的研究文獻(xiàn)較多，按照匹配算法設(shè)計(jì)原理可分為序列相關(guān)匹配方法和遞推濾波方法兩種，序列相關(guān)匹配方法主要包括最近等值線迭代算法（ICCP），相關(guān)極值分析算法兩大類(lèi)[1-5]；獲取高精度、高分辨率重力異常圖是進(jìn)行準(zhǔn)確匹配的關(guān)鍵所在。高斯樣條函數(shù)具有局部支撐的良好特性，用其作為樣條基函數(shù)對(duì)計(jì)算區(qū)域重力異常進(jìn)行二維整體逼近，能在保證重力異常局部特性不失真的前提下獲取計(jì)算區(qū)域重力異常的統(tǒng)一解析式等優(yōu)點(diǎn)，近年來(lái)有關(guān)學(xué)者進(jìn)行了相關(guān)應(yīng)用研究。但針對(duì)該問(wèn)題中涉及到的局部支撐參數(shù)計(jì)算問(wèn)題分析并不多，局部支撐參數(shù)過(guò)大或過(guò)小都將導(dǎo)致很大的逼近誤差，太大或太小均失去插值意義，有必要對(duì)其進(jìn)行尋優(yōu)，研究相應(yīng)的尋優(yōu)算法[6-7]。

目前重力輔助導(dǎo)航系統(tǒng)狀態(tài)方程建模中通過(guò)地球正常重力模型得到的正常重力矢量進(jìn)行重力補(bǔ)償，未考慮重力異常擾動(dòng)矢量、標(biāo)準(zhǔn)重力值誤差（標(biāo)準(zhǔn)重力值誤差由測(cè)量位置維度誤差解算引起）和厄特弗斯修正計(jì)算值三類(lèi)誤差影響，這些誤差都將引起帶來(lái)導(dǎo)航誤差，隨著慣性元件精度的不斷提高，慣性元件帶來(lái)的系統(tǒng)誤差與上述三類(lèi)誤差帶來(lái)的系統(tǒng)誤差處在同一個(gè)量級(jí)，系統(tǒng)誤差不能被忽略，有必要采用精度更高的重力輔助導(dǎo)航系統(tǒng)狀態(tài)方程模型對(duì)正常重力矢量進(jìn)行補(bǔ)償[8-12]。

所以，本文首先針對(duì)二維高斯樣條函數(shù)逼近局部重力異常場(chǎng)中的局部支撐參數(shù)選擇問(wèn)題，通過(guò)對(duì)所涉及的插值系數(shù)矩陣可逆性進(jìn)行研究給出了局部支撐參數(shù)的上限值，而后通過(guò)對(duì)插值系數(shù)矩陣解誤差研究，給出了矩陣條件數(shù)和局部支撐參數(shù)定性關(guān)系，最后結(jié)合插值模型精度評(píng)估理論，提出了一種新的高斯基插值函數(shù)最優(yōu)局部支撐參數(shù)計(jì)算方法，接著局部連續(xù)重力異常場(chǎng)提出了一種考慮重力擾動(dòng)矢量、標(biāo)準(zhǔn)重力值誤差、厄特弗斯修正計(jì)算值的精準(zhǔn)重力輔助導(dǎo)航模型構(gòu)建算法，最終為重力輔助慣性導(dǎo)航理論的高精度實(shí)踐化應(yīng)用打基礎(chǔ)。

1 高斯基插值函數(shù)最優(yōu)局部支撐參數(shù)計(jì)算

1.1 二維高斯樣條優(yōu)化插值數(shù)學(xué)模型

式中，Δx、Δy為格網(wǎng)點(diǎn)分辨率[7]。這樣符合上述條件的xy平面二維高斯基函數(shù)可寫(xiě)成矩陣形式：

插值條件：

將矩陣組合則有如下線性方程成立：

可見(jiàn)，高斯基插值函數(shù)存在兩個(gè)重要參數(shù)需要確定，第一：插值系數(shù)矩陣λ；第二：局部支撐參數(shù)a。如果X、Y均為非奇異矩陣，則式(2)有唯一解，將插值系數(shù)矩陣λ代入式(1)便得到該局部重力異常基準(zhǔn)圖二維高斯基函數(shù)逼近解析式；局部支撐參數(shù)參數(shù)a決定了X、Y矩陣求逆的具體求逆誤差，應(yīng)賦予局部支撐參數(shù)a值，使矩陣X、Y求逆誤差較小；當(dāng)a過(guò)小時(shí)，由于局部支撐區(qū)間很小，使得插值曲線剛度小，容易出現(xiàn)鋸齒狀；當(dāng)a過(guò)大時(shí)，由于局部支撐區(qū)間很大，插值曲線剛度大，插值曲線過(guò)于平滑甚至使插值解析式在已知點(diǎn)的插值誤差很大，所以a太大或太小均失去插值意義，因此有必要對(duì)參數(shù)a進(jìn)行尋優(yōu)，結(jié)合插值模型精度評(píng)估方法進(jìn)行解算獲取參數(shù)a的最優(yōu)值。

1.2 插值系數(shù)矩陣可逆性和求逆誤差

1.2.1 插值系數(shù)矩陣可逆性研究

設(shè)上述式(2)中高斯對(duì)稱矩陣階數(shù)K，如下示：

上述系數(shù)矩陣記為Ψ，考慮Ψ在什么條件下滿足非奇異，由數(shù)值計(jì)算方法理論知只要該矩陣是嚴(yán)格對(duì)

角占優(yōu)即可。也就是對(duì)于第i行，：

嚴(yán)格對(duì)角占優(yōu)：

成立，上述式(5)左邊可放大為：

所以只要確保式(7)成立則可以保證系數(shù)矩陣Ψ非奇異，對(duì)曲線的高斯插值總是可以保證的。可見(jiàn)a值最大限度為。

1.2.2 求逆矩陣誤差研究

對(duì)于矩陣Ψ，設(shè)為定向量空間nR上的向量范數(shù)，矩陣范數(shù)由式(8)給出：

可逆矩陣Ψ的條件數(shù)由式(9)給出：

由于計(jì)算機(jī)計(jì)算系數(shù)矩陣Ψ時(shí)，每個(gè)元素存在計(jì)算誤差，于是設(shè)擾動(dòng)矩陣為E，則存在如下引理：

引理1：如果對(duì)任意滿足的擾動(dòng)矩陣E和一個(gè)與E無(wú)關(guān)的正常數(shù)δ，矩陣Ψ的條件數(shù)由：給出，對(duì)于方程Ψx=b和解x和y，式(10)成立，則。

又由式(11)知：2/(1-β)＞1，所以0＜β＜1。于是當(dāng)β單調(diào)增時(shí)，1-β單調(diào)減小，2(1-β)單調(diào)增，單調(diào)增。所以當(dāng)時(shí)，此時(shí)對(duì)應(yīng)的值就是方程的解x和y的最小的相對(duì)誤差上界。于是有以下定理結(jié)論：

定理1：如果對(duì)任意滿足的擾動(dòng)矩陣E，矩陣Ψ的條件數(shù)由：給出，方程的解x和y，于是關(guān)于矩陣Ψ、E、x、y存在如下不等式：

定理1說(shuō)明在誤差估計(jì)式(13)中，用求逆條件數(shù)可以得到最小的相對(duì)誤差上界。該等式給出了矩陣擾動(dòng)相對(duì)誤差和由擾動(dòng)引起的矩陣方程解相對(duì)誤差值間關(guān)系。

圖1 Ψ矩陣的條件數(shù)分析圖Fig.1 The condition number analysis of Ψ

1.3 插值模型精度評(píng)估

對(duì)局部支撐參數(shù)a進(jìn)行尋優(yōu)，為了保證插值函數(shù)與已知數(shù)據(jù)點(diǎn)保形性良好，同時(shí)考慮更多的函數(shù)非線性變化因素，研究利用分段二次拋物線代替所求插值函數(shù)對(duì)應(yīng)的準(zhǔn)確值，尋求分段二次拋物線與插值函數(shù)距離最小時(shí)對(duì)應(yīng)的局部支撐參數(shù)amin，同時(shí)將二維局部重力異常基準(zhǔn)圖問(wèn)題解耦為x方向一維高斯基函數(shù)和y方向一維高斯基函數(shù)分別構(gòu)造進(jìn)行尋優(yōu)建模。求解其對(duì)應(yīng)的x方向最優(yōu)局部支撐參數(shù)axmin，y方向最優(yōu)局部支撐參數(shù)aymin，由此選擇可使插值函數(shù)與已知點(diǎn)具有最好的保形性，差值誤差較小。

具體建模方法為：對(duì)n行格網(wǎng)點(diǎn)，每行分別進(jìn)行上述二維高斯插值局部支撐參數(shù)的尋優(yōu)建模，而后利用所有尋優(yōu)建模準(zhǔn)則函數(shù)之和作為二維高斯樣條函數(shù)插值x方向的尋優(yōu)模型，局部支撐優(yōu)化函數(shù)(Local Support Optimize Function)LSOF尋優(yōu)模型如式(15)示；對(duì)m列格網(wǎng)點(diǎn)，同x方向的尋優(yōu)處理方法，LSOF如式(17)示，于是將局部支撐參數(shù)a限定在[0.2，2.5]區(qū)間內(nèi)，分別按照式(15)(17)尋優(yōu)準(zhǔn)則進(jìn)行計(jì)算，即可得到最優(yōu)化局部支撐參數(shù)。

式中，ai，j、bi，j、ci，j為方程組(16)的解：

其中，ai，j、bi，j、ci，j為方程組(18)的解，

2 精準(zhǔn)重力輔助導(dǎo)航模型構(gòu)建算法

2.1 速度誤差精準(zhǔn)模型方程

實(shí)際系統(tǒng)中總存在各種誤差，所以實(shí)際系統(tǒng)中考慮各種誤差情況下，實(shí)際速度計(jì)算值應(yīng)由下述方程確定：

用式(20)減去式(19)，略去二階小量可得線性近似的捷聯(lián)慣導(dǎo)誤差方程為：

2.2 狀態(tài)方程

除速度誤差方程外導(dǎo)航系統(tǒng)其他誤差方程利用傳統(tǒng)思路方法參考相關(guān)文獻(xiàn)建立[13]，于是可以定義SINS的狀態(tài)變量如下，共15維狀態(tài)矢量：

于是得到局部連續(xù)重力輔助導(dǎo)航精準(zhǔn)模型方程:

式中，矩陣wSINS為6維系統(tǒng)干擾白噪聲，矩陣GSINS是15×6噪聲分配矩陣，矩陣FSINS是15×15狀態(tài)矩陣。

2.3 量測(cè)方程

本文系統(tǒng)量測(cè)值取為重力儀量測(cè)重力與慣導(dǎo)指示位置對(duì)應(yīng)重力圖上重力之差作為觀測(cè)方程。設(shè)定此時(shí)系統(tǒng)的量測(cè)方程Zg可由式(24)右邊計(jì)算得出：

由于繪制成的重力圖值，實(shí)質(zhì)上等于實(shí)際重力測(cè)量值經(jīng)過(guò)厄特弗斯修正而后減去標(biāo)準(zhǔn)重力值。所以對(duì)于實(shí)際位置點(diǎn)(L，)λ有：

代入方程(24)經(jīng)線性化處理后，可表示為方程(26)所示系統(tǒng)的量測(cè)方程為：

VGravity為重力傳感器、重力異常圖制作噪聲和線性化誤差總和，在實(shí)際計(jì)算中，重力異常梯度，可以采用九點(diǎn)擬合、全平面擬合等方法計(jì)算。

所以，根據(jù)式(32)確定其測(cè)量矩陣為

式中，Zgline為式(24)中計(jì)算得到的重力測(cè)量誤差值，為1×15維觀測(cè)矩陣。為重力傳感器的噪聲、重力異常圖制作噪聲和線形化誤差總和。

3 仿真試驗(yàn)與分析

本文采用第二炮兵工程大學(xué)慣性導(dǎo)航實(shí)驗(yàn)室基于Matlab軟件自行開(kāi)發(fā)研制的慣性導(dǎo)航軌跡生成軟件生成仿真需要的真實(shí)航跡和慣導(dǎo)指示航跡。重力異常數(shù)據(jù)使用渤海地區(qū)經(jīng)度范圍 118°～119°，緯度范圍 37°～38°的真實(shí)數(shù)據(jù)，該重力異常場(chǎng)的三維圖如圖2所示。重力圖分辨率為 1′×1′，重力儀的實(shí)時(shí)量測(cè)數(shù)據(jù)是利用真實(shí)航跡在重力數(shù)據(jù)庫(kù)中的采樣值疊加0.4 mGal的隨機(jī)量測(cè)噪聲獲取，陀螺常值漂移取 0.01°/h，隨機(jī)噪聲漂移為0.01/h，加速度計(jì)隨機(jī)常值零偏1×10-4g，隨機(jī)噪聲漂移為 1×10-5g；初始狀態(tài)：速度 10 m/s、經(jīng)度118.02°、緯度37.24°、慣導(dǎo)測(cè)量經(jīng)維度誤差0.02°、對(duì)準(zhǔn)水平姿態(tài)誤差角2′、對(duì)準(zhǔn)航向誤差角10′。運(yùn)載體保持俯仰、橫滾角 10°、航向角 30°的姿態(tài)，0.005 m/s2加速度每隔10 min加速一次，3 min采樣匹配濾波一次，行進(jìn)了70 min時(shí)間，得到24個(gè)測(cè)量點(diǎn)[12]。捷聯(lián)慣導(dǎo)的情況下進(jìn)行的仿真，對(duì)比分析了ICCP算法初次估計(jì)解算，局部連續(xù)場(chǎng)情況下傳統(tǒng)重力輔助導(dǎo)航模型和本文敘述的基于高斯插值的精準(zhǔn)重力輔助導(dǎo)航模型導(dǎo)航解算法導(dǎo)航效果（分別稱其為tradition-ICCP和gausscompensation-ICCP算法，簡(jiǎn)記為T(mén)RAI和GACOI算法）。

圖2 重力異常場(chǎng)三維圖Fig.2 3D map of gravity anomaly field

圖3(a)～圖3(b)為T(mén)RAI和本文GACOI新算法（最優(yōu)局部支撐參數(shù)）導(dǎo)航的經(jīng)、緯度位置誤差，航向、俯仰、橫滾姿態(tài)誤差角，東、北、天向速度誤差分布圖。如圖3中所示：TRAI算法導(dǎo)航位置誤差在 125～330 m、平均 210 m，水平姿態(tài)誤差在0.2′～0.75′、平均 0.50′，航向姿態(tài)誤差在 4′～9′、平均誤差7.5′，速度誤差在0.027～0.085 m/s、平均誤差0.045 m/s；本文 GACOI新算法導(dǎo)航能夠保持在100～200 m、平均110 m的位置誤差，0.01′～0.25′、平均 0.18′的水平姿態(tài)誤差，2′～6′、平均 2.5′的航向姿態(tài)誤差，0.002～0.05 m/s、平均0.016 m/s的速度誤差。可見(jiàn)，本文 GACOI導(dǎo)航新算法將傳統(tǒng)重力輔助慣性導(dǎo)航算法的位置導(dǎo)航精度提高了1倍左右，姿態(tài)角導(dǎo)航精度、速度導(dǎo)航精度提高了1～2倍?？梢?jiàn)新算法引入了先進(jìn)的局部支撐模型參數(shù)確定方法，同時(shí)建了考慮重力擾動(dòng)矢量、標(biāo)準(zhǔn)重力值誤差、厄特弗斯修正計(jì)算值的精準(zhǔn)重力輔助導(dǎo)航模型進(jìn)行導(dǎo)航解算提高了導(dǎo)航精度，達(dá)到了接近于高精度導(dǎo)航的水平。

圖3 導(dǎo)航仿真誤差分析圖Fig.3 The analysis of navigation simulation errors

4 結(jié) 論

本文從分析重力輔助導(dǎo)航所需的地球重力場(chǎng)數(shù)學(xué)建模開(kāi)始，通過(guò)對(duì)高斯樣條優(yōu)化插值所涉及到的系數(shù)矩陣可逆性進(jìn)行研究，得出了局部支撐參數(shù)的最大限度值，緊接著通過(guò)分析與之相關(guān)的求逆誤差進(jìn)一步印證了結(jié)論的準(zhǔn)確性，給出了局部支撐參數(shù)的數(shù)值區(qū)間，接著對(duì)插值模型精度評(píng)估理論進(jìn)行研究，設(shè)計(jì)了給定區(qū)間上的最優(yōu)局部支撐參數(shù)求解方法，而后建了考慮重力擾動(dòng)矢量、標(biāo)準(zhǔn)重力值誤差、厄特弗斯修正計(jì)算值的精準(zhǔn)重力輔助導(dǎo)航模型。最終實(shí)現(xiàn)了重力輔助導(dǎo)航算法的精確導(dǎo)航定位計(jì)算。本文研究?jī)?nèi)容不但可以為重力輔助慣性導(dǎo)航算法的實(shí)踐化做一定的理論鋪墊，而且相關(guān)理論研究對(duì)推進(jìn)數(shù)學(xué)高斯插值相關(guān)研究也有一定的意義，當(dāng)慣性導(dǎo)航系統(tǒng)導(dǎo)航觀測(cè)模型為非線性情況下，如何在高斯樣條逼近模型的基礎(chǔ)上展開(kāi)相應(yīng)的濾波算法有待進(jìn)一步研究。

References)

[1]Mcphail S.Autosub6000:A deep diving long range AUV[J].Journal of Bionic Engineering，2009(6):55-62.

[2]李?yuàn)檴?，吳曉平，馬彪.水下重力異常相關(guān)極值匹配算法[J].測(cè)繪學(xué)報(bào)，2011，40(4)：464-470.LI Shan-shan，WU Xiao-ping，MA Biao.Correlative extremum matching algorithm using underwater gravity anomalies[J].Acta Geodaetica et Cartographica Sinica，2011，40(4):464-470.

[3]袁贛南，張紅偉，袁克非，吳簡(jiǎn)彤.基于重力梯度輔助定位的概率神經(jīng)網(wǎng)絡(luò)改進(jìn)方法[J].中國(guó)慣性技術(shù)學(xué)報(bào)，2013，21(3)：369-374.YUAN Gan-nan，ZHANG Hong-wei，YUAN Ke-fei，WU Jian-tong.Improved probabilistic neural network method based on gravity gradient aided location[J].Journal of Chinese Inertial Technology，2013，21(3):369-374.

[4]Jiang F，Wu Y M，Zhang Z S，et al.Combinational seabed terrain matching algorithm basing on probability data associate filtering and iterative closest contour point[C]//2009 Second International Conference on Intelligent Computation Technology and Automation.IEEE Computer Society，Piscataway，2009:245-249.

[5]Yuan G N，Zhang H W，Yuan K F.A Combinational underwater aided navigation algorithm based on TERCOM/ICCP and Kalman filter[J].IEEE Journal of Sciences and Optimization，2011，23(1):952-955.

[6]Wang Z G，Bian S F，Xiao S H.A local geopotential model for implementation of underwater passive navigation[J].Progress in Natural Science，2008，18(9):1139-1145.

[7]童余德，邊少鋒，將東方.一種新的基于局部重力圖逼近的組合匹配算法[J].地球物理學(xué)報(bào)，2012，55(9)：2917-2925.TONG Yu-de，BIAN Shao-feng，JIANG Dong-fang.A new integrated gravity matching algorithm based on approximated local gravity map[J].Chinese Journal of Geophysics，2012，55(9):2917-2925.

[8]Dransfield M，ZENG Y.Airbone Gravity gradiometry:Terrian corrections and elevation error[J].Geophysics，2009，74(5):137-142.

[9]李勝全，歐陽(yáng)永忠，常國(guó)賓.慣性導(dǎo)航系統(tǒng)重力擾動(dòng)量補(bǔ)償技術(shù)[J].中國(guó)慣性技術(shù)學(xué)報(bào)，2012，20(4)：410-414.LI Sheng-quan，OUYANG Yong-zhong，CHANG Guo-bin.Compensation technology of gravity disturbance vector in inertial navigation system[J].Journal of Chinese Inertial Technology，2012，20(4):410-414.

[10]DeGregoria A.Gravity gradiometry and map matching an aid to aircraft inertial navigation systems[D].Ohio:Air Force Institute of Technology，2010.

[11]Senobari M.New results in airborne vector gravimetry using strapdown INS/DGPS[J].Journal of Geodesy，2010，84:277-291.

[12]Neumeyer J，Sch?ferb U，Kremerc J，et al.Derivation of gravity anomalies from airborne gravimeter and IMU recordings-Validation with regional analytic models using ground and satellite gravity data[J].Journal of Geodynamics，2009，47:191-200.

Constructive model algorithm for gravity aided navigation in local anomaly field

WEN Chao-bin1，2，WANG Yue-gang1，TIAN Qi1，GUO Zhi-bin1，ZUO Zhao-yang1，YANG Jia-sheng1
(1.304 Unit，The Second Artillery Engineering University，Xi’an 710025，China;2.Unit 96124 of the Chinese people’s liberation army，Tonghua 134000，China)

In studying high-precision gravity aided navigation，it is indispensable to build a local continuous gravity anomaly field model with high precision and good analytic property，and accurately compensate the gravity in error equation.To solve the choice problem of local support optimal parameters in local continuous gravity anomaly field model with 2D gauss spline interpolation，an local support optimal calculation method is proposed by studying the related problems such as，coefficient of matrix，solution errors，and interpolation model evaluation，etc.Then a high-precision construct algorithm based on the gravity aided navigation model is put forward，which has compensated the gravity disturbance vector，standard gravity value error and eotvos effect.The experiment results show that the positioning accuracy is increased by about 1 times，the attitude and the velocity are both increased by 1～2 times.Its latitude and longitude errors are kept within 100～300 m.

gravity aided navigation;gauss spline interpolation;gravity error compensation;navigation model

U666.1

A

1005-6734(2014)03-0368-06

10.13695/j.cnki.12-1222/o3.2014.03.017

2013-12-25；

2014-04-11

國(guó)防預(yù)研（103030203）

文超斌（1986—），男，博士研究生，從事導(dǎo)航、制導(dǎo)與控制研究。E-mail：weijing123@126.com

聯(lián) 系 人：王躍鋼（1958—），男，教授，博士生導(dǎo)師。E-mail：wangyueg@163.com