基于曲率統(tǒng)計(jì)的LiDAR點(diǎn)云二次濾波方法

萬(wàn)劍華,黃榮剛,周 行,曾 喆

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580；2.中國(guó)石油大學(xué)地球科學(xué)學(xué)院,北京 102249)

基于曲率統(tǒng)計(jì)的LiDAR點(diǎn)云二次濾波方法

萬(wàn)劍華1,黃榮剛1,周 行2,曾 喆1

(1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580；2.中國(guó)石油大學(xué)地球科學(xué)學(xué)院,北京 102249)

針對(duì)傳統(tǒng)偏度平衡方法濾波結(jié)果中存在低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云問(wèn)題,在傳統(tǒng)偏度平衡方法的點(diǎn)云一次濾波算法的基礎(chǔ)上,提出一種基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法。對(duì)該方法進(jìn)行試驗(yàn),并將試驗(yàn)結(jié)果與傳統(tǒng)偏度平衡方法濾波結(jié)果進(jìn)行對(duì)比分析。結(jié)果表明:基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法比傳統(tǒng)偏度平衡法能夠多濾除83%的植被點(diǎn)云、5%的建筑物點(diǎn)云,能夠有效地濾除傳統(tǒng)偏度平衡方法濾波結(jié)果中的低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云。

激光雷達(dá)；點(diǎn)云；濾波方法；曲率統(tǒng)計(jì)

機(jī)載激光雷達(dá)(LiDAR)點(diǎn)云濾波是獲取高精度數(shù)字高程模型的關(guān)鍵問(wèn)題。目前,已研究出了許多濾波算法,有數(shù)學(xué)形態(tài)學(xué)濾波方法[1-2]、基于坡度的濾波方法[3]、漸進(jìn)三角網(wǎng)加密方法[4-5]、線性預(yù)測(cè)方法[6]等。但是,這些算法都需要根據(jù)經(jīng)驗(yàn)進(jìn)行閾值的設(shè)定,而自然世界的地形、地物變化復(fù)雜,難以確定合適的閾值。Bartels等[7]提出了一種非監(jiān)督、免閾值的濾波算法——偏度平衡方法(skewness balancing),Bartels[8]、Bao等[9-10]對(duì)算法進(jìn)行了改進(jìn),但該算法仍存在濾波結(jié)果中遺留有低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云的問(wèn)題。為此,筆者在傳統(tǒng)偏度平衡方法點(diǎn)云一次濾波的基礎(chǔ)上,提出一種基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法,以濾除傳統(tǒng)偏度平衡方法濾波結(jié)果中剩余的非地面點(diǎn)云。

1 濾波方法

對(duì)LiDAR點(diǎn)云濾波采用兩次濾波處理,首先對(duì)點(diǎn)云進(jìn)行傳統(tǒng)偏度平衡方法點(diǎn)云一次濾波處理,然后對(duì)剩余點(diǎn)云進(jìn)行基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法處理,減少非地面點(diǎn)云被誤判為地面點(diǎn)云的錯(cuò)誤,提高提取地形數(shù)據(jù)的精度,如圖1所示。

圖1 LiDAR點(diǎn)云濾波方法流程Fig.1 Flow chart of LiDAR point cloud filter method

在現(xiàn)實(shí)中,LiDAR點(diǎn)云一般都含有無(wú)法引起曲面局部突變的非地面點(diǎn)云區(qū)塊,如建筑物頂面和側(cè)面等。這部分點(diǎn)云的曲率信息與地面點(diǎn)云的曲率信息具有相似的特性,容易造成這部分點(diǎn)云的曲率信息與地面點(diǎn)云曲率信息混淆。如果直接采用基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波處理,容易出現(xiàn)大量非地面點(diǎn)云被誤判為地面點(diǎn)云的情況。因此,需要先采用基于傳統(tǒng)偏度平衡方法的點(diǎn)云一次濾波來(lái)濾除高程值較大且無(wú)法引起曲面局部突變的點(diǎn)云區(qū)塊。

點(diǎn)云一次濾波處理后剩余的點(diǎn)云中存在遺留有低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云問(wèn)題。綜合分析后發(fā)現(xiàn),濾波后剩余的點(diǎn)云數(shù)據(jù)具有以下特征:地面點(diǎn)云占主體且數(shù)量較大,而遺留的非地面點(diǎn)云只占少數(shù)；地面點(diǎn)云高程變化比較緩慢,表現(xiàn)為曲率較小,而非地面點(diǎn)云數(shù)量小、分布范圍小,能夠引起曲面的局部突變,導(dǎo)致非地面點(diǎn)云具有較大的曲率值。據(jù)此,根據(jù)中心極限定理可以假設(shè):地面點(diǎn)云的曲率信息符合正態(tài)分布,非地面點(diǎn)的曲率看作是對(duì)地面點(diǎn)云曲率正態(tài)分布的一種擾動(dòng)。因此,當(dāng)Li-DAR點(diǎn)云存在非地面點(diǎn)云時(shí),其點(diǎn)云曲率信息的分布表現(xiàn)為右尾拉長(zhǎng)現(xiàn)象(偏度值大于零),可以采用如圖2所示方法從曲率方面來(lái)區(qū)分地面和非地面點(diǎn)云高程間的不同,以對(duì)非地面點(diǎn)云進(jìn)行濾除。

(1)點(diǎn)云離散曲率計(jì)算。對(duì)剩余點(diǎn)云構(gòu)建不規(guī)則三角網(wǎng)(TIN),計(jì)算點(diǎn)云的離散曲率。在基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波中,曲率一般選擇最大主曲率較為合適。因?yàn)橐栏皆诘孛娴牡桶堑孛纥c(diǎn)云,都能引起局部曲面的尖銳突變,表現(xiàn)為曲率較大；而地形一般具有連續(xù)、平緩變化的特性,曲率相對(duì)于非地面點(diǎn)云較小。選用最大主曲率能夠使得非地面點(diǎn)云曲率和地面點(diǎn)云曲率信息在頻率分布圖中的距離變大,有利于非地面點(diǎn)云的濾除。采用Taubin方法[11-12]進(jìn)行離散主曲率估算。首先計(jì)算連續(xù)曲面上點(diǎn)的主曲率,然后對(duì)連續(xù)曲面公式進(jìn)行離散化,使其適合三角網(wǎng)各頂點(diǎn)的最大主曲率、最小主曲率等的計(jì)算。

圖2 基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波流程Fig.2 Flow chart of the secondary filter method of LiDAR point cloud based on curvature statistics

設(shè)T為曲面上一點(diǎn)p處的某一單位切向量,N為p處的法向量,T1、T2為頂點(diǎn)p處切平面內(nèi)的主方向向量,則T可以表示為nN+t1T1+t2T2(n、t1、t2為N、T1、T2的系數(shù))。

p點(diǎn)處沿T方向的法曲率kp(T)按下式計(jì)算:

根據(jù)Euler公式,對(duì)沿T方向的法曲率進(jìn)行變換表示:

式中,θ為T與主方向T1、T2的夾角,且-π≤θ≤π。

定義對(duì)稱矩陣Mp:

根據(jù)式(5)、(6),可以得出主曲率是對(duì)稱矩陣Mp非零特征值的函數(shù):

對(duì)于TIN的主曲率估算,關(guān)鍵是對(duì)連續(xù)曲面主曲率計(jì)算過(guò)程中的Mp、法曲率、法向量等進(jìn)行離散化。設(shè)TIN的某頂點(diǎn)為vi,頂點(diǎn)vi對(duì)應(yīng)的法向量為Nvi,對(duì)應(yīng)的1環(huán)鄰點(diǎn)集為V。頂點(diǎn)vi的1環(huán)鄰域(鄰域內(nèi)的三角形集記為F)內(nèi)某一三角形(fk)的單位法向量為Nfk,三角形(fk)的面積為Afk。

vi處法向量Nvi計(jì)算公式為

vj為vi的1環(huán)鄰點(diǎn),由vi到vj的向量記為Vij。Vij在點(diǎn)vi切平面上投影為Tij,則沿Tij方向的法曲率為

頂點(diǎn)vi的Mvi估算,實(shí)際上就是對(duì)Mp進(jìn)行離散化,表示為

“新時(shí)代”催生“新任務(wù)”，在新時(shí)代背景下，為主動(dòng)適應(yīng)市場(chǎng)競(jìng)爭(zhēng)，傳統(tǒng)專業(yè)在人才培養(yǎng)目標(biāo)、核心課程設(shè)置、實(shí)驗(yàn)實(shí)訓(xùn)課程等方面必然有內(nèi)在的“新要求”和改革內(nèi)涵。為進(jìn)一步說(shuō)明專業(yè)綜合改革和發(fā)展內(nèi)涵，本文以安徽財(cái)經(jīng)大學(xué)財(cái)政學(xué)專業(yè)為例展開分析。

其中,ωij為面積加權(quán)；E為單位權(quán)矩陣。計(jì)算Mvi的特征向量和特征值,然后根據(jù)式(7)和(8)計(jì)算出頂點(diǎn)vi的主曲率、。

高斯曲率(KGauss)和平均曲率(KMean)的計(jì)算公式為

(2)點(diǎn)云二次濾波處理。以曲率信息作為研究對(duì)象,進(jìn)行基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波處理。利用下式計(jì)算點(diǎn)云曲率信息分布的偏度值:

式中,S為所剩點(diǎn)云曲率信息分布的偏度值；n為所剩點(diǎn)云Xi的總數(shù)量；σ為所剩點(diǎn)云Xi曲率的標(biāo)準(zhǔn)方差；μ為所剩點(diǎn)云Xi曲率的算術(shù)平均值。

如果偏度值為0,則濾波處理結(jié)束；如果偏度值大于0,表明所剩點(diǎn)云中存在非地面點(diǎn),則將曲率信息最大的點(diǎn)濾除,再次計(jì)算剩余點(diǎn)云曲率的偏度值。按這種方式循環(huán)剔除點(diǎn)云中的非地面點(diǎn)云。

基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法是一個(gè)循環(huán)濾波的過(guò)程(圖2)。由于使用基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波處理后,點(diǎn)云重構(gòu)三角網(wǎng)的離散點(diǎn)曲率發(fā)生了變化,導(dǎo)致三角網(wǎng)重構(gòu)后的曲率信息不再服從正態(tài)分布。因此,需要對(duì)點(diǎn)云進(jìn)行循環(huán)濾波,直到基于曲率統(tǒng)計(jì)的LiDAR點(diǎn)云二次濾波濾除的點(diǎn)云數(shù)為零為止。

2 試驗(yàn)結(jié)果分析

2.1 試驗(yàn)數(shù)據(jù)

圖3 區(qū)域A點(diǎn)云分類參考數(shù)據(jù)Fig.3 Reference data of classification of LiDAR point cloud in region A

2.2 試驗(yàn)過(guò)程及其結(jié)果分析

試驗(yàn)主要是在對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行基于傳統(tǒng)偏度平衡方法的點(diǎn)云一次濾波后,進(jìn)行基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波。

首先對(duì)區(qū)域A點(diǎn)云數(shù)據(jù)進(jìn)行基于傳統(tǒng)偏度平衡方法的點(diǎn)云一次濾波處理,濾波結(jié)果如圖4(a)所示,灰色點(diǎn)云為地面點(diǎn)云、黑色點(diǎn)云為非地面點(diǎn)云。由于區(qū)域A地形略有起伏,出現(xiàn)地面點(diǎn)云高于某些非地面點(diǎn)云的情況,導(dǎo)致部分非地面點(diǎn)云混淆于地面點(diǎn)云高程分布中。并且該部分點(diǎn)云占據(jù)了較大比例,成為地面點(diǎn)云高程正態(tài)分布中必不可少的一部分。因此,濾波結(jié)果中存在低矮植被以及建筑物側(cè)面墻腳等非地面點(diǎn)云被誤判為地面點(diǎn)的情況。部分被誤判的非地面點(diǎn)云如圖4(a)中A、B、C、D等處黑色點(diǎn)云所示,其中D處為建筑物側(cè)面墻腳點(diǎn)云,其余3處為低矮植被點(diǎn)云。

經(jīng)過(guò)基于傳統(tǒng)偏度平衡方法的點(diǎn)云一次濾波處理后,對(duì)剩余點(diǎn)云的曲率進(jìn)行綜合分析發(fā)現(xiàn):地面點(diǎn)云的曲率信息占主體,遺留的非地面點(diǎn)云占少數(shù)；遺留的非地面點(diǎn)都能引起擬合曲面(試驗(yàn)采用不規(guī)則三角網(wǎng)TIN)局部突變,導(dǎo)致局部曲率值突變。這些特征保證了地面點(diǎn)云曲率符合正態(tài)分布,非地面點(diǎn)云曲率是對(duì)其正態(tài)分布的一種擾動(dòng),并且非地面點(diǎn)云曲率值偏大,使得曲率分布呈右尾拉長(zhǎng)狀。因此,可以采用基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波進(jìn)行處理。在曲率的選擇上,經(jīng)過(guò)對(duì)地面點(diǎn)云和非地面點(diǎn)云曲率特性的分析,以及多組濾波試驗(yàn)結(jié)果分析,發(fā)現(xiàn)基于最大主曲率統(tǒng)計(jì)的點(diǎn)云二次濾波對(duì)點(diǎn)云的濾波效果最好,而基于平均曲率、高斯曲率和最小曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法對(duì)點(diǎn)云基本沒(méi)有濾波效果,所以在試驗(yàn)中選擇最大主曲率來(lái)進(jìn)行基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波,濾波結(jié)果如圖4(b)所示,結(jié)果中剩余較少的非地面點(diǎn)云。

圖4 基于傳統(tǒng)偏度平衡方法的點(diǎn)云濾波與基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波結(jié)果對(duì)比Fig.4 Result comparison of conventional skewness balancing and the secondary filter method of LiDAR point cloud based on curvature statistics

由圖4可以發(fā)現(xiàn):圖4(a)中A、B兩處植被點(diǎn)云,在基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波后基本被濾除；圖4(a)中C處的點(diǎn)云是分布在道路旁邊的密集且高度大致的人造植被,導(dǎo)致部分點(diǎn)云曲率值較小,從而在基于曲率統(tǒng)計(jì)的二次濾波后還有部分植被點(diǎn)云被保留；圖4(a)中D處的建筑物側(cè)面墻腳點(diǎn)云,在基于曲率統(tǒng)計(jì)的二次濾波后基本被濾除,只剩余零星的建筑物側(cè)面墻腳點(diǎn)云。由此,基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波能夠有效的濾除基于傳統(tǒng)偏度平衡方法濾波結(jié)果中的大部分低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云。

將基于傳統(tǒng)偏度平衡方法的點(diǎn)云濾波結(jié)果、基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波結(jié)果與區(qū)域A的分類參考數(shù)據(jù)進(jìn)行定量統(tǒng)計(jì)分析,統(tǒng)計(jì)結(jié)果如表1、圖5所示。其中圖5中柱狀之上的百分率為非地面點(diǎn)云與參考點(diǎn)云中相應(yīng)非地面點(diǎn)云數(shù)量的比值。表1和圖5都反映了各類別點(diǎn)云數(shù)量隨點(diǎn)云濾波處理的變化。

結(jié)合表1和圖5,對(duì)試驗(yàn)結(jié)果數(shù)據(jù)進(jìn)行定量分析,可以得出:

(1)基于傳統(tǒng)偏度平衡方法的點(diǎn)云濾波結(jié)果中植被點(diǎn)云剩余1 196個(gè)點(diǎn),占原始植被點(diǎn)云的94.3%,建筑物點(diǎn)云剩余163個(gè)點(diǎn),占原始建筑物點(diǎn)云的6.355%,表明算法在建筑物點(diǎn)云濾波方面具有較好效果,在植被點(diǎn)云濾波方面效果非常差；而基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波,濾除的植被點(diǎn)云占原始植被點(diǎn)云的89.117%、建筑物點(diǎn)云占原始建筑物點(diǎn)云的98.752%,較傳統(tǒng)偏度平衡法對(duì)植被點(diǎn)云的濾除效果提高了83%,對(duì)建筑物點(diǎn)云的濾波效果提高了5%。因此,基于曲率統(tǒng)計(jì)的二次濾波方法較傳統(tǒng)偏度平衡法更有效地克服了其對(duì)植被濾波效果不好的缺點(diǎn),有效地提高了濾波結(jié)果的精度,減少了自動(dòng)濾波后人工編輯工作量。

(2)基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波后仍存在少量的非地面點(diǎn)云。其原因是:區(qū)域A地形略有起伏,導(dǎo)致少量非地面點(diǎn)云與部分地面點(diǎn)云的曲率信息出現(xiàn)混淆,使得部分非地面點(diǎn)云曲率無(wú)法擾動(dòng)地面點(diǎn)云曲率的正態(tài)分布。因此,基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法目前僅適用于較平坦區(qū)域。

綜上,基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波,在較平坦區(qū)域?qū)趥鹘y(tǒng)偏度平衡方法點(diǎn)云濾波結(jié)果中的低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云的濾除是有效的。

表1 各類別點(diǎn)云數(shù)量統(tǒng)計(jì)Table 1 Quantity statistics of LiDAR point cloud for each class

圖5 非地面點(diǎn)云數(shù)量變化柱狀圖Fig.5 Bar chart of quantity change in non-ground points

3 結(jié)束語(yǔ)

針對(duì)傳統(tǒng)偏度平衡方法在地形略有起伏的情況下濾波結(jié)果中遺留有低矮植被、建筑物側(cè)面墻腳等非地面點(diǎn)云問(wèn)題,發(fā)展了一種基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波對(duì)傳統(tǒng)偏度平衡方法濾波結(jié)果進(jìn)行二次濾波處理,提高了點(diǎn)云自動(dòng)濾波的正確率。通過(guò)對(duì)地形略有起伏的某居民區(qū)LiDAR點(diǎn)云進(jìn)行濾波處理試驗(yàn),并與傳統(tǒng)偏度平衡方法濾波試驗(yàn)結(jié)果對(duì)比分析,發(fā)現(xiàn)基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波方法能夠有效地濾除較平坦地區(qū)的LiDAR點(diǎn)云經(jīng)過(guò)傳統(tǒng)偏度平衡方法濾波處理后遺留的大部分非地面點(diǎn)云。基于曲率統(tǒng)計(jì)的點(diǎn)云二次濾波采用的曲率為最大主曲率,對(duì)某些地形突變細(xì)節(jié)特征具有腐蝕影響,而采用其他曲率信息,則濾波結(jié)果中仍過(guò)多地保留非地面點(diǎn)云。

[1] LINDENBERGER J.Laser-profilmenssungen zur Topographischen Gelandeaufnahme[D].Stuttgart:Universitat Stuttgart,1993.

[2] KILIAN J,HAALA N,ENGLICH M.Capture and evaluation of airborne laser data[J].International Archives of Photogrammetry and Remote Sensing,1996,31(3):383-388.

[3] VOSSELMAN G.Slope based filtering of laser altimetry data[J].International Archives of Photogrammetry and Remote Sensing,2000,33:958-964.

[4] AXELSSON P.DEM generation from laser scanner data using adaptive TIN models[J].International Archives of Photogrammetry and Remote Sensing,2000,33:85-92.

[5] 趙影.基于機(jī)載LIDAR數(shù)據(jù)的城市建筑物提取研究[D].長(zhǎng)春:吉林大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,2011.

ZHAO Ying.Extraction study of urban buildings based on airborne LIDAR data[D].Changchun:College of Computer Science and Technology in Jilin University,2011.

[6] KRAUS K,PFEIFER N.Determination of terrain models in wooded areas with airborne laser scanner data[J].ISPRS Journal of Photogrammetry and Remote Sensing, 1998,53(4):193-203.

[7] BARTELS M,HONG W.Segmentation of LIDAR data using measures of distribution[J].International Archives of Photogrammetry,Remote Sensing and Spatial Information Sciences,2006,36(7):426-31.

[8] BARTELS M,HONG W.Threshold-free object and ground point separation in LIDAR data[J].Pattern Recognition Letters,2010,31:1089-1099.

[9] BAO Yun-fei,CAO Chun-xiang,CHANG Chao-yi,et al. Segmentation to the point clouds of LIDAR data based on change of kurtosis[C]//International Symposium on Photoelectronic Detection and Imaging Technology and Application(ISPDI)2007,Beijing,China,2007.Bellingham, USA:SPIE,c2008.

[10] BAO Yun-fei,LI Guo-ping,CAO Chun-xiang,et al. Classification of LIDAR point cloud and generation of DTM from LIDAR height and intensity data in forested area[J].International Archives of the Photogrammetry, Remote SensingandSpatialInformationSciences, 2008,37(B3b):313-318.

[11] TAUBIN G.A signal processing approach to fair surface design[C]//Proceeding of the 22nd annual conference on Computer Graphics and Interactive Techniques, 1995.New York,USA:ACM,1995:351-358.

[12] 齊寶明.三角網(wǎng)格離散曲率估計(jì)和Taubin方法改進(jìn)[D].大連:大連理工大學(xué)數(shù)學(xué)科學(xué)學(xué)院,2008.

QI Bao-ming.Curvatures estimation and the improvement of Taubin?s method on triangular mesh[D].Dalian: School of Mathmatical Sciences in Dalian University of Technology,2008.

(編輯 修榮榮)

A secondary filter method of LiDAR point cloud based on curvature statistics

WAN Jian-hua1,HUANG Rong-gang1,ZHOU Hang2,ZENG Zhe1

(1.School of Geosciences in China University of Petroleum,Qingdao 266580,China； 2.College of Geosciences in China University of Petroleum,Beijing 102249,China)

The rusults of conventional skewness balancing exist non-ground points of low vegetation and side of buildings. Aimed at this,a secondary filter method of LiDAR point cloud was developed based on curvature statistics to filter the remained non-ground points in the result of filter based on conventional skewness balancing.The results of the secondary filter method proposed and conventional skewness balancing were compared.The results show that the former can remove 83% more pionts of vegetation and 5%more pionts of building than the latter.The secondary filter method of LiDAR point cloud based on curvature statistics can efficiently filter the non-ground points remained in the conventional skewness balancing method.

LiDAR；point cloud；filter method；curvature statistics

P 237

A

1673-5005(2013)01-0056-05

10.3969/j.issn.1673-5005.2013.01.009

2011-11-30

中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)(10CX05007A)

萬(wàn)劍華(1966-),男,教授,博士,主要從事數(shù)字油田、數(shù)字海洋等方面的研究工作。E-mail:wjh66310@163.com。