周 瀟,楊功流,蔡慶中
(1. 北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院,北京 100191;2. 慣性技術(shù)國防重點(diǎn)實(shí)驗(yàn)室,北京 100191)
基于小波神經(jīng)網(wǎng)絡(luò)的高精度慣導(dǎo)重力擾動(dòng)補(bǔ)償方法
周 瀟,楊功流,蔡慶中
(1. 北京航空航天大學(xué) 儀器科學(xué)與光電工程學(xué)院,北京 100191;2. 慣性技術(shù)國防重點(diǎn)實(shí)驗(yàn)室,北京 100191)
重力擾動(dòng)矢量(空間同一點(diǎn)實(shí)際重力與正常重力之差,包括垂線偏差和重力異常兩部分) 一直是慣性導(dǎo)航系統(tǒng)的重要誤差源之一。針對(duì)重力擾動(dòng)誤差精確補(bǔ)償問題,推導(dǎo)并建立了考慮重力擾動(dòng)的慣導(dǎo)誤差方程,并提出了基于小波神經(jīng)網(wǎng)絡(luò)的重力擾動(dòng)補(bǔ)償方法。通過仿真驗(yàn)證了小波神經(jīng)網(wǎng)絡(luò)的重力擾動(dòng)補(bǔ)償方法對(duì)慣導(dǎo)導(dǎo)航精度的提高效果。24 h仿真結(jié)果表明:所提出的重力擾動(dòng)補(bǔ)償方法能有效減小慣導(dǎo)導(dǎo)航系統(tǒng)誤差,經(jīng)重力擾動(dòng)矢量補(bǔ)償后,速度誤差最大能減小約0.2 m/s,降低約30%,位置誤差最大能減小約3000 m,降低約25%。
捷聯(lián)慣導(dǎo)系統(tǒng);重力擾動(dòng)補(bǔ)償;慣導(dǎo)誤差方程;小波神經(jīng)網(wǎng)絡(luò)
慣性導(dǎo)航系統(tǒng)[1]不需要任何外來信息,也不會(huì)向外輻射任何信息,僅靠慣性導(dǎo)航系統(tǒng)本身就能在全天侯條件下,在全球范圍內(nèi)和任何介質(zhì)環(huán)境里自主地、隱蔽地進(jìn)行連續(xù)的三維定位和三維定向,具有諸多其他導(dǎo)航系統(tǒng)無法比擬的優(yōu)點(diǎn),故廣泛應(yīng)用于軍事,民用等諸多領(lǐng)域。隨著慣性器件精度的逐步提高以及高精度捷聯(lián)慣性導(dǎo)航系統(tǒng)應(yīng)用需求的提出,重力擾動(dòng)已經(jīng)成為高精度捷聯(lián)慣性導(dǎo)航系統(tǒng)最主要的剩余誤差源之一[2-3]。
為了解決這一問題,本文提出了一種基于小波神經(jīng)網(wǎng)絡(luò)的重力擾動(dòng)預(yù)測(cè)及補(bǔ)償方法,并通過仿真試驗(yàn)驗(yàn)證了該方法對(duì)重力擾動(dòng)的預(yù)測(cè)精度以及在高精度慣導(dǎo)中的補(bǔ)償效果。
在慣性導(dǎo)航系統(tǒng)中,陀螺儀和加速度計(jì)通過穩(wěn)定平臺(tái)安裝在載體上,分別測(cè)量載體相對(duì)慣性空間的角運(yùn)動(dòng)信息和線運(yùn)動(dòng)信息,在給定相關(guān)初值和其他必要信息(如重力信息)的情況下,經(jīng)過一系列運(yùn)算和控制可以得到載體相對(duì)于導(dǎo)航坐標(biāo)系的位置、速度和姿態(tài)等信息。圖1為本文所提出的考慮重力擾動(dòng)補(bǔ)償?shù)膽T導(dǎo)系統(tǒng)原理圖[4]。
圖1 慣性導(dǎo)航系統(tǒng)原理圖Fig.1 Principle of inertial navigation system
坐標(biāo)系定義:b系—載體坐標(biāo)系;n系—導(dǎo)航坐標(biāo)系(在本文中取為東-北-天當(dāng)?shù)厮街副弊鴺?biāo)系);e系—地固坐標(biāo)系;i系—慣性坐標(biāo)系。
重力擾動(dòng)矢量δg為實(shí)際重力矢量與正常重力矢量γ之間的差值[5]:
式中:0γ可以根據(jù)正常重力模型得到,通常采用1980年國際正常重力公式:
式(3)中:g0=9.780327,a=0.0053024,b=0.0000058。重力擾動(dòng)矢量可以表示為
式(4)中:η和ξ分別為東西向和南北向垂線偏差,一般在10″左右,個(gè)別地區(qū)可以達(dá)到30″;ΔgE、ΔgN、分別為重力擾動(dòng)矢量的東向、北向、垂向分量。
捷聯(lián)慣性導(dǎo)航誤差方程[6-8]:
姿態(tài)誤差方程為
式(5)中:
式(6)中:δKG為陀螺標(biāo)度系數(shù)誤差;δG為陀螺安裝誤差矩陣;φ為姿態(tài)角誤差;為坐標(biāo)系轉(zhuǎn)換矩陣;為載體系下的陀螺輸出;nε為導(dǎo)航坐標(biāo)系下的陀螺漂移;為指令角速度誤差,且可表示為
式中:L為載體所在位置的緯度;RM、RN、h分別為載體所在位置的子午圈半徑、卯酉圈半徑和海拔高度。
根據(jù)公式(5)可知,引起慣導(dǎo)姿態(tài)誤差的主要誤差源有陀螺的輸出誤差,除此之外還可以通過速度位置的耦合關(guān)系引起姿態(tài)誤差,體現(xiàn)于一項(xiàng)。
速度更新誤差方程為
式(9)中:
式中:δKA為加速度計(jì)標(biāo)度系數(shù)誤差;δA為加速度計(jì)安裝誤差矩陣;fn為導(dǎo)航坐標(biāo)系下的加速度計(jì)輸出;為導(dǎo)航坐標(biāo)系下加速度計(jì)零偏;δgn為導(dǎo)航坐標(biāo)系下的重力擾動(dòng)。
位置更新誤差方程為
式中:L、λ分別為地理緯度和經(jīng)度。
由公式(9)可知,引起慣導(dǎo)系統(tǒng)速度誤差的誤差源主要有加速度計(jì)輸出誤差慣導(dǎo)前一時(shí)刻積累的速度誤差在導(dǎo)航坐標(biāo)系下的投影、速度位置的耦合誤差項(xiàng)以及重力擾動(dòng)δgn。
由公式(11)可知,引起慣導(dǎo)系統(tǒng)位置誤差的主要誤差源為慣導(dǎo)前一時(shí)刻的速度誤差δVn以及速度位置的耦合誤差。
通過以上對(duì)重力擾動(dòng)在慣導(dǎo)系統(tǒng)中的誤差分析可以得出,重力擾動(dòng) δgn首先通過公式(9)引起慣導(dǎo)的速度誤差,速度誤差進(jìn)一步通過公式(11)引起慣導(dǎo)系統(tǒng)位置誤差,最終通過公式(5)引起慣導(dǎo)的姿態(tài)誤差,并且慣導(dǎo)誤差隨時(shí)間的增長(zhǎng)而發(fā)散。
在以往的慣導(dǎo)系統(tǒng)中,重力擾動(dòng)所引起的誤差遠(yuǎn)遠(yuǎn)小于由慣性器件所帶來的誤差,所以在慣導(dǎo)解算時(shí)通常對(duì)重力擾動(dòng)所帶來的系統(tǒng)誤差予以忽略。但是隨著慣性器件以及慣導(dǎo)標(biāo)定技術(shù)的快速發(fā)展,重力擾動(dòng)已成為慣導(dǎo)解算中影響系統(tǒng)誤差的主要來源,制約著慣導(dǎo)導(dǎo)航精度的進(jìn)一步提高,所以必須對(duì)慣導(dǎo)解算中的重力擾動(dòng)項(xiàng)進(jìn)行補(bǔ)償。
小波神經(jīng)網(wǎng)絡(luò)(Wavelet Neural Network,WNN)是將小波分析理論和神經(jīng)網(wǎng)絡(luò)相結(jié)合的一種神經(jīng)網(wǎng)絡(luò)模型,可以實(shí)現(xiàn)時(shí)域和頻域的同時(shí)分析,具有良好的局部特性,較強(qiáng)的學(xué)習(xí)能力和任意函數(shù)逼近能力[8-9]。WNN以傳統(tǒng)的BP神經(jīng)網(wǎng)絡(luò)構(gòu)架為支架,將神經(jīng)網(wǎng)絡(luò)隱含層的傳輸函數(shù)用小波函數(shù)來替代,實(shí)現(xiàn)信號(hào)的前向傳播和誤差的方向傳播。本文以3層神經(jīng)網(wǎng)絡(luò)為例,設(shè)計(jì)出基于局部區(qū)域重力擾動(dòng)數(shù)據(jù)的小波神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型,其網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。
圖2 小波神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.2 Structure of wavelet neural network
圖2中:X1,X2, …,Xk為小波神經(jīng)網(wǎng)絡(luò)的輸入?yún)?shù),本文選載體所在位置的經(jīng)緯度(λ,L)作為輸入?yún)?shù);Y1,Y2,…,Yk為小波神經(jīng)網(wǎng)絡(luò)的預(yù)測(cè)值輸出,這里為重力擾動(dòng)的水平分量,東西向和南北向的垂線偏差η、ξ和垂向重力擾動(dòng)Δg。ωij和ωjk為小波神經(jīng)網(wǎng)絡(luò)權(quán)值。在輸入信號(hào)為xi(i=1,2, …,k)時(shí),隱含層輸出公式為[10-11]
式中:h(j)為隱含層第j個(gè)節(jié)點(diǎn)輸出值;ωij為輸入層和隱含層之間的連接權(quán)值;bj為小波基函數(shù)hj的平移因子;aj為小波基函數(shù)hj的伸縮因子;hj為小波基函數(shù)。小波神經(jīng)網(wǎng)絡(luò)輸出層的輸出公式為
式中:ωik為隱含層到輸出層的連接權(quán)值;h(i)為第i個(gè)隱含層節(jié)點(diǎn)的輸出;l為隱含層節(jié)點(diǎn)數(shù);m為輸出層節(jié)點(diǎn)數(shù)。
在小波神經(jīng)網(wǎng)絡(luò)的應(yīng)用中,小波基函數(shù)的選擇對(duì)于網(wǎng)絡(luò)預(yù)測(cè)效果影響非常的大,根據(jù)不同的問題選擇不同的基函數(shù),但是目前尚沒有理論上的支撐。Morlet小波函數(shù)多用于分類、圖像識(shí)別和特征提?。桓咚挂浑A導(dǎo)數(shù)多用于函數(shù)估計(jì)和預(yù)測(cè);Mexican hat 小波函數(shù)多用于系統(tǒng)辨識(shí)等。本文通過反復(fù)計(jì)算得出,當(dāng)選擇高斯一階導(dǎo)數(shù)作為小波基函數(shù)時(shí),網(wǎng)絡(luò)模型有著較高的預(yù)測(cè)精度。
小波神經(jīng)網(wǎng)絡(luò)權(quán)值參數(shù)修正算法類似于BP神經(jīng)網(wǎng)絡(luò)權(quán)值得修正,采用梯度修正法修正網(wǎng)絡(luò)權(quán)值和小波基函數(shù)參數(shù),從而使小波神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)輸出不斷地逼近期望輸出[12]。
計(jì)算網(wǎng)絡(luò)預(yù)測(cè)誤差:
式中:yn(k)為期望誤差;y(k)為小波神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)輸出。
根據(jù)預(yù)測(cè)誤差e修正小波神經(jīng)網(wǎng)絡(luò)的權(quán)值和小波基函數(shù)系數(shù):
式中:η為小波神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)速率,取值在0到1之間。
本文提出的重力擾動(dòng)數(shù)據(jù)小波神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)算法的步驟如圖3所示。
圖3 小波神經(jīng)網(wǎng)絡(luò)算法流程圖Fig.3 Flow chart of wavelet neural network
1)網(wǎng)絡(luò)初始化。隨機(jī)初始化小波函數(shù)伸縮因子aj,平移因子bj以及網(wǎng)絡(luò)連接權(quán)值ωij、ωik,并設(shè)置網(wǎng)絡(luò)學(xué)習(xí)速率η。
2)樣本分類。將樣本分為訓(xùn)練樣本和測(cè)試樣本。
3)預(yù)測(cè)輸入。將訓(xùn)練樣本輸入網(wǎng)絡(luò),計(jì)算網(wǎng)絡(luò)輸出和網(wǎng)絡(luò)預(yù)測(cè)誤差e。
4)權(quán)值修正。根據(jù)網(wǎng)絡(luò)預(yù)測(cè)誤差e對(duì)網(wǎng)絡(luò)權(quán)值和小波函數(shù)參數(shù)進(jìn)行修正,使網(wǎng)絡(luò)預(yù)測(cè)值逼近期望值。
5)判斷e值。判斷網(wǎng)絡(luò)預(yù)測(cè)誤差e是否滿足給定誤差閾值,若滿足則結(jié)束訓(xùn)練,否則返回步驟3)重新訓(xùn)練。
仿真條件:仿真模擬數(shù)據(jù)來自美國國家大地?cái)?shù)據(jù)測(cè)量局提供的太平洋某片海域N40~43°、W160°~163°的重力擾動(dòng)格網(wǎng)數(shù)據(jù),格網(wǎng)數(shù)據(jù)的空間分辨率為1′×1′。網(wǎng)格內(nèi)重力擾動(dòng)數(shù)據(jù)的分布圖如圖4所示。
采用小波神經(jīng)網(wǎng)絡(luò)建立重力擾動(dòng)預(yù)測(cè)模型時(shí),將格網(wǎng)區(qū)域內(nèi)的重力擾動(dòng)數(shù)據(jù)分成兩個(gè)部分:一部分為訓(xùn)練樣本,占總點(diǎn)數(shù)的70%;另一部分為測(cè)試樣本,占總點(diǎn)數(shù)的30%。本文隨機(jī)選取10個(gè)點(diǎn)對(duì)小波神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)精度進(jìn)行評(píng)價(jià),并與反距離加權(quán)(IDW)插值結(jié)果相比較,計(jì)算結(jié)果如表1、表2所示。
圖4 格網(wǎng)區(qū)域內(nèi)重力擾動(dòng)分布圖Fig.4 Profile of gravity disturbance in the selected mesh region
由表1和表2的對(duì)比結(jié)果可以看出,利用小波神經(jīng)網(wǎng)絡(luò)對(duì)格網(wǎng)區(qū)域內(nèi)重力擾動(dòng)數(shù)據(jù)的預(yù)測(cè)精度明顯高于采用IDW插值方法的預(yù)測(cè)精度。當(dāng)采用小波神經(jīng)網(wǎng)絡(luò)進(jìn)行預(yù)測(cè)時(shí),重力擾動(dòng)各分量與真實(shí)值接近,求差后重力異常誤差的極值僅在1 mGal左右,垂線偏差誤差的極值也都在5″之內(nèi),從而說明小波神經(jīng)網(wǎng)絡(luò)對(duì)格網(wǎng)內(nèi)的重力擾動(dòng)數(shù)據(jù)有很好的預(yù)測(cè)效果,預(yù)測(cè)精度完全可以滿足高精度慣導(dǎo)對(duì)于重力擾動(dòng)精度的需求。
表1 IDW插值誤差評(píng)價(jià)表Tab.1 Performance of IDW interpolation
表2 小波神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)誤差評(píng)價(jià)表Tab.2 Performance of WNN estimation
仿真條件:三個(gè)陀螺的常值漂移為3×10-3(°)/h,隨機(jī)游走系數(shù)為2×10-4/h,加速度計(jì)的零漂均為10×10-6g。為了更好地分析重力擾動(dòng)對(duì)慣導(dǎo)誤差的影響,故忽略其他誤差的影響,設(shè)安裝誤差角為零,無初始對(duì)準(zhǔn)誤差。仿真過程中忽略高度通道的影響。仿真步長(zhǎng)設(shè)為1 s,仿真時(shí)間為24 h。模擬艦船實(shí)驗(yàn),故設(shè)載體的東向、北向速度為5 m/s,天向速度為0 m/s。仿真中利用提供的1′×1′的重力格網(wǎng)數(shù)據(jù)插值后得到網(wǎng)格區(qū)域內(nèi)的重力擾動(dòng)數(shù)據(jù)來模擬實(shí)際重力擾動(dòng),加入到慣導(dǎo)系統(tǒng)中可以得到無重力擾動(dòng)補(bǔ)償時(shí)捷聯(lián)慣導(dǎo)誤差,仿真結(jié)果如圖5所示。采用本文提出的小波神經(jīng)網(wǎng)絡(luò)重力擾動(dòng)補(bǔ)償方法,進(jìn)行重力擾動(dòng)補(bǔ)償后,得到的系統(tǒng)誤差如圖6所示。
圖5 無重力擾動(dòng)補(bǔ)償慣導(dǎo)系統(tǒng)誤差Fig.5 System errors of INS without gravity disturbance compensation
由圖5、圖6的仿真結(jié)果可以得出:經(jīng)過小波神經(jīng)網(wǎng)絡(luò)重力擾動(dòng)補(bǔ)償后的慣導(dǎo)速度誤差和位置誤差都有明顯減小。補(bǔ)償后的速度誤差最大能減小約0.2 m/s,降低約30%,位置誤差最大能減小約3000 m,降低約25%。以上分析表明,本文提出的重力擾動(dòng)補(bǔ)償方法可以有效地提高高精度慣導(dǎo)的導(dǎo)航精度。
圖6 重力擾動(dòng)補(bǔ)償后慣導(dǎo)系統(tǒng)誤差Fig.6 System errors of INS with gravity disturbance compensation
隨著慣性元件精度的不斷提高以及慣導(dǎo)技術(shù)的快速發(fā)展,重力擾動(dòng)已成為慣性導(dǎo)航系統(tǒng)的重要誤差源,必須對(duì)其進(jìn)行補(bǔ)償才能進(jìn)一步的提高慣導(dǎo)的導(dǎo)航精度。本文建立了考慮重力擾動(dòng)矢量的慣導(dǎo)系統(tǒng)誤差方程,并提出了基于小波神經(jīng)網(wǎng)絡(luò)的重力擾動(dòng)補(bǔ)償方法。通過仿真驗(yàn)證了小波神經(jīng)網(wǎng)絡(luò)的重力擾動(dòng)補(bǔ)償方法對(duì)慣導(dǎo)導(dǎo)航精度的提高效果。仿真結(jié)果表明:本文所提出的重力擾動(dòng)補(bǔ)償方法能有效地減小慣導(dǎo)導(dǎo)航系統(tǒng)誤差;經(jīng)重力擾動(dòng)矢量補(bǔ)償后,仿真24h速度誤差最大能減小約0.2 m/s,降低約30%,位置誤差最大能減小約3000 m,降低約25%。
(References):
[1]吳太旗, 歐陽永忠, 陸秀平, 等. 重力匹配導(dǎo)航的影響模式分析[J]. 中國慣性技術(shù)學(xué)報(bào), 2011, 19(5):559-564.Wu Tai-qi, Ouyang Yong-zhong, Lu Xiu-ping, et al. Analysis on effecting mode of several essential factors to gravity aided navigation[J]. Journal of Chinese Inertial Technology, 2011, 19(5): 559-564.
[2] 周瀟, 楊功流, 王晶, 等. 基于Kalman濾波原理對(duì)慣導(dǎo)中重力擾動(dòng)的估計(jì)及補(bǔ)償方法[J]. 中國慣性技術(shù)學(xué)報(bào),2015, 23(6): 721-726.Zhou Xiao, Yang Gong-liu, Wang Jing, et al. Estimation and compensation for gravity disturbance based on Kalman filtering in inertial navigation[J]. Journal of Chinese Inertial Technology, 2015, 23(6): 721-726.
[3]Man X, Fang J CH, Ning X L. An overview of autonomous navigation for a gravity-assist interplanetary spacecraft[J]. Progress in Aerospace Sciences, 2013, 63: 56-66.
[4]Rummel R, Yi W, Stummer C. GOCE gravitational gradiometry[J]. Journal of Geodesy, 2011, 85: 777-790.
[5]Cai Shao-kun, Zhang Kai-dong, Wu Mei-ping. Improving airborne strapdown vector gravimetry using stabilized horizontal components[J]. Journal of Applied Geophysics.2013: 79-89.
[6]Li X. Strapdown INS/DGPS airborne gravimetry tests in the Gulf of Mexico[J]. Journal of Geodesy, 2011, 85:597-605.
[7] 王虎彪, 王勇, 方劍. “最小均方差旋轉(zhuǎn)擬合法”重力輔助導(dǎo)航仿真研究[J]. 中國科學(xué): 地球科學(xué), 2012, 42(7):1055-1062.Wang Hu-biao, Wang Yong, Fang Jian. Simulation research on a minimum root-mean-square error rotationfitting algorithm for gravity matching navigation[J].Science China: Earth Sciences, 2012, 42(7): 1055-1062.
[8]戰(zhàn)德軍, 戴東凱, 張忠華, 等. 單軸旋轉(zhuǎn)INS/GPS組合導(dǎo)航中重力垂線偏差引起的姿態(tài)誤差分析[J]. 中國慣性技術(shù)學(xué)報(bào), 2014, 22(3): 301-305.Zhan De-jun, Dai Dong-kai, Zhang Zhong-hua et al.Analysis of gravity vertical deflection-induced attitude error in single-axis rotation INS/GPS integrated navigation system[J]. Journal of Chinese Inertial Technology, 2014, 22(3): 301-305.
[9]Falamarzi Y, Palizdan N, Huang Y F, et al. Estimating evapotranspiration from temperature and wind speed data using artificial and wavelet neural networks (WNNs)[J].Agric. Water Manag., 2014, 140: 26-36.
[10]趙忠, 王鵬. 高精度慣性導(dǎo)航系統(tǒng)垂線偏差影響與補(bǔ)償[J]. 中國慣性技術(shù)學(xué)報(bào), 2013, 21(6): 701-705.Zhao Zhong, Wang Peng. Analysis and compensation of vertical deflection effect on high accuracy inertial navigation system[J]. Journal of Chinese Inertial Technology,2013, 21(6): 701-705.
[11]Alexandridis A K, Zapranis A D. Wavelet neural networks:a practical guide[J]. Neural Networks, 2013, 42: 1-27.
[12]Guan C, Luh P B, Michel L D, et al. Very short-term load forecasting: wavelet neural networks with data pre-filtering[J]. IEEE Transactions on Power Systems, 2013,28(1): 30-41.
Compensation on gravity disturbance for high-precision INS based on wavelet neural network
ZHOU Xiao, YANG Gong-liu, CAI Qing-zhong
(1. School of Instrumentation Science and Opto-electronics Engineering, Beihang University, Beijing 100191, China;2. Science and Technology on Inertial Laboratory, Beijing 100191, China)
The gravity disturbance is one of the important error sources of inertial navigation system, which is the difference between the actual gravity and normal gravity on the same point, and is composed of two parts: deviation of plumb line and gravity anomaly. In this paper, the error equation of inertial navigation system considering gravity disturbance is derived and established, and the gravity disturbance compensation based on wavelet neural network is proposed. The 24 h simulation results show that the proposed gravity disturbance compensation method can effectively reduce the navigation errors of INS, and after the gravity disturbance compensation, the maximum velocity error and the maximum position error are reduced by about 0.2 m/s and 3000 m respectively, therefore the velocity accuracy and the position accuracy can be improved by about 30% and 25%, respectively.
strapdown inertial navigation system; gravity disturbance compensation; inertial navigation error equation; wavelet neural network
U666.1
A
1005-6734(2016)05-0571-06
10.13695/j.cnki.12-1222/o3.2016.05.003
2016-06-03;
2016-09-12
國家自然科學(xué)基金(61340044,11202010)
周瀟(1987—),男,博士研究生,從事慣性技術(shù)研究。E-mail: buaa17zhou@163.com
聯(lián) 系 人:楊功流(1967—),男,教授,博士生導(dǎo)師。E-mail: bhu17-yang@139.com