產品幾何規(guī)范中非理想表面的多尺度表征

石照耀，姜海洋，張　敏

(北京工業(yè)大學 北京市精密測控技術與儀器工程技術研究中心，北京 100124)

?

產品幾何規(guī)范中非理想表面的多尺度表征

石照耀，姜海洋，張敏

(北京工業(yè)大學 北京市精密測控技術與儀器工程技術研究中心，北京 100124)

考慮目前多數(shù)計算機輔助公差工具僅能針對具有理想幾何表面的CAD模型，無法從物理幾何角度真正反映制造誤差，本文研究了非理想表面的多尺度表征，提出了一種能夠表征產品表面及其截面輪廓宏觀及微觀多尺度形貌誤差的非理想表面模型。首先，提出了一種利用離散小波實現(xiàn)形貌誤差多尺度仿真的方法；其次，針對實際工件，利用離散小波對其表面及截面輪廓采樣數(shù)據進行形貌誤差多尺度仿真；最后，對形貌誤差各尺度成分進行合成，得到具有多尺度形貌誤差成分的工件三維表面及其二維截面輪廓的非理想表面模型。仿真及實驗結果表明：利用提出方法可以實現(xiàn)具有多尺度形貌誤差的非理想表面模型表征；仿真所得粗糙度Ra值與白光干涉儀測量所得值的平均相對誤差不超過4%。得到的結果證明了提出方法的正確性和可行性，為更加全面地表征產品的非理想表面模型提供了有效途徑。

產品幾何技術規(guī)范；非理想表面模型；離散幾何；小波分析；多尺度表征

1　引　言

幾何公差用于控制零件幾何要素的形狀、方向、位置與其設計值的偏差[1-2]。目前，大多數(shù)計算機輔助公差工具可以幫助設計者實現(xiàn)公差規(guī)范，但是其研究對象僅針對具有理想幾何表面的CAD模型，不能從物理幾何角度真正反映制造誤差，且缺少對于產品幾何量變化的控制能力，不利于公差信息在整個產品生命周期的傳遞，也不利于計算機集成制造技術在現(xiàn)代制造業(yè)中的有效利用。同時，由于通用商業(yè)CAD軟件的局限性，公差信息僅以文本信息形式存在于計算機中，不利于產品信息的數(shù)字化傳遞，這些都阻礙了公差信息在產品設計、制造及檢驗過程中的有效應用。

針對上述問題，Mathieu和Ballu提出了用于產品規(guī)范的GeoSpelling模型，該模型通過基于非理想表面模型的一系列概念和操作算子可實現(xiàn)公差在計算機中的表達及計算機輔助公差設計[3-4]。法國LURPA課題組對非理想表面模型[5]進行了一系列研究，張敏等[6-7]率先提出了利用離散幾何算法實現(xiàn)非理想表面模型的建模仿真方法，以及針對基本零件曲面類型的形位公差的仿真算法，實現(xiàn)了GD&T與非理想幾何形狀模型的信息集成，對非理想表面模型的表征及仿真研究進行了初步探索。Schleich等[8]在張敏研究工作的基礎上，提出了非理想表面模型表征的綜合建模框架，對非理想表面模型的仿真與表征進行了完善，并成功應用于具有高副轉動機構的公差分析中。上述關于非理想表面模型表征的研究均將制造誤差劃分為系統(tǒng)誤差和隨機誤差來進行仿真分析，雖然能夠表征出工件宏觀的尺寸及形位公差變化信息，但無法描繪其微觀表面形貌誤差信息。因此，本研究提出一種能夠反映產品加工表面宏觀及微觀多尺度形貌誤差的非理想表面模型的表征和仿真方法，對于公差的數(shù)學模型以數(shù)據結構的存儲形式有機集成到計算機系統(tǒng)中具有重要作用。

本文以新一代產品幾何技術規(guī)范(GPS)為基礎，以GeoSpelling模型為對象，以真實測量所獲得的離散點云數(shù)據為切入點，通過離散幾何及小波分析等現(xiàn)代工具，研究了一種可表征產品表面及其截面輪廓宏觀及微觀多尺度形貌誤差的非理想表面模型的離散數(shù)據模型建模方法，可解決產品在設計、制造、檢驗過程中規(guī)范表達統(tǒng)一的難題，從而實現(xiàn)新一代GPS標準中規(guī)范與認證的對偶性要求。

2　形貌誤差多尺度表征方法

非理想表面模型能夠在物理幾何角度真正反映制造誤差信息，是對于物體真實形狀的表征[3]，所以在構造非理想表面模型時要充分考慮研究物體真實形狀的多尺度形貌誤差。對于任意工件的加工表面，其表面形貌均可劃分為宏觀和微觀，包括形狀誤差、波紋度以及粗糙度等多尺度形貌誤差成分[9-12]。所以，為了獲得含有宏觀和微觀形貌誤差成分的非理想表面模型，針對上述3種尺度的形貌誤差成分進行仿真分析以實現(xiàn)其誤差表征。

2.1離散小波變換原理

為了實現(xiàn)3種尺度形貌誤差成分的表征，首先要選擇一種能夠準確提取出多尺度誤差成分的方法。目前，傳統(tǒng)的濾波技術大多基于傅里葉變換，這種方法是將表面信號看做不同頻率正弦波的線性疊加[13]。在幾何產品表面成分分析過程中，對于高頻信號和低頻信號的分析宜分別采用小時間窗和大時間窗。然而，這種變時間窗的要求是時間窗固定的傅里葉變換無法滿足的，其時域定位精度有限，無法實現(xiàn)時頻同時定位。

基于傅里葉變換的上述局限性，20世紀80年代小波分析理論逐漸興起。小波分析具有多尺度分析的特點，在時域和頻域都有表征信號局部信息的能力，時間窗和頻率窗都可以根據信號的具體形態(tài)進行動態(tài)調整?；谶@些特點小波分析可以探測正常信號中的瞬態(tài)成分并展示其頻率成分，具有對信號同時進行空間和頻率定位的能力，能夠實現(xiàn)表面信號多尺度頻率成分的提取。因此，采用小波分析作為提取并表征產品表面及其截面輪廓多尺度形貌誤差的方法。離散后的小波基函數(shù)在理想情況下滿足正交完備性條件，此時計算出的小波變換系數(shù)無任何冗余度，可以最大限度地壓縮數(shù)據并減小計算量。因此，采用離散小波變換對表面形貌采樣數(shù)據進行多尺度濾波處理，可以解決同時表達微觀和宏觀表面數(shù)據量大的問題。

2.2輪廓誤差小波模型

假設產品表面的截面輪廓為p(x)，將p(x)看作一維信號的離散采樣點，且假設p(x)為所能獲得的粗糙表面最精細的綜合形貌，是包括表面形狀誤差、波紋度以及粗糙度的多尺度頻率成分，頻率為ω0～ω3。則有如圖1所示形狀誤差f(x)對應低頻成分[ω0，ω1]，波紋度w(x)對應帶通成分[ω1，ω2]，粗糙度r(x)對應高頻成分[ω2，ω3]，相應數(shù)學模型為：

p(x)=f(x)=w(x)+r(x).

(1)

圖1　形貌誤差成分頻帶分布

根據多尺度分析的Mallat算法，原始截面輪廓信號p(x)經過在尺度空間的j次分解后被分解為子集aj和dj，分別代表輪廓信號中的低頻和高頻成分在尺度2-j上的信息。p(x)分解后得到的小波分析結果為[15]：

Wj[p(x)]=(a1，d1)=(a2,d1,d2)=…=

(aj，d1，…，dk，…，dj).

(2)

2.3小波基選取

采用小波分析方法對截面輪廓形貌誤差進行分析，首先要對其在不同尺度上進行分解，然后通過重構獲得所需要的3種不同尺度的誤差成分信息。由于小波分析的基函數(shù)具有不唯一性，采用不同小波基對同一產品表面進行分析時獲得的誤差成分信息結果相差甚遠[15]，因此對于小波基的選取至關重要。

典型的正交小波基有Haar小波、Symlets小波、Coiflets小波及Daubechies小波，正交小波具有相同的分解和重構尺度及小波函數(shù)。然而，當這二者對應相同時，線性相位和完全重構條件不可兼得[13]。雙正交小波分解與重構的小波尺度函數(shù)并不相同，因此這里選擇雙正交小波基，在可以實現(xiàn)完全重構的條件下獲取線性相位，從而實現(xiàn)對于表面信號在時域及頻域上的準確提取，進而實現(xiàn)表面特征的多尺度表征。此外，在傳輸特性上雙正交小波基具有平滑的截止特性，且低通幅頻特性接近水平，不存在過沖現(xiàn)象，故適合分析表面形貌信號[16]。

2.4頻率成分劃分

因為二維輪廓各尺度的形貌誤差成分的主要區(qū)別在于所處的頻率帶不同，所以可以通過帶通濾波器將其區(qū)分開。其中，尺度系數(shù)aj代表低頻帶[ω0，ω1]的信息，即形狀誤差；dk+1，…，dj，代表頻率帶[ω1，ω2]的信息，即波紋度；小波系數(shù)d1，…，dk，代表高頻帶[ω2，ω3]的信息，即粗糙度。其中，ω1=1/λf，λf為波紋度成分的頻率下限；ω2=1/λc，λc為粗糙度成分的頻率下限(粗糙度截止頻率)；ω3=1/λs，λs為采樣間隔決定的頻率上限[17]，如圖2所示。

圖2　形貌誤差成分傳輸特性

形狀誤差：f(x)=IW(aj)=aj(x)-μ[dj(x)]

(3)

3　形貌誤差多尺度表征方法驗證

3.1仿真驗證

(4)

利用離散小波對p(n)進行多尺度分解及重構，得到如圖3所示仿真效果圖。其中，采樣點數(shù)為10 000個，采樣間隔為1，其中橫坐標表示采樣點數(shù)，縱坐標表示模擬信號的幅值，單位為微米。由圖中的仿真結果可以看出，離散小波變換可以成功地分離出不同尺度的正弦波信號和輪廓誤差信號。

圖3　輪廓形貌誤差多尺度仿真

3.2實驗驗證

為了驗證提出的方法可以應用在基于實測離散點云數(shù)據的形貌誤差成分表征中，采用如圖4所示的Talysurf CCI系列白光干涉儀(WLI)對實驗工件表面形貌進行測量。

利用Talymap軟件對測量區(qū)域沿y軸方向進行8次不重復隨機輪廓截取，獲得8組不同的產品截面輪廓離散采樣數(shù)據及Ra值。其中，白光干涉儀每次取樣區(qū)域為0.847 76 mm×0.847 84 mm，截取后的表面輪廓沿x軸方向取樣長度為lr=0.847 76 mm，采樣點數(shù)為1 024個，采樣間距為0.828 7 μm。由于加工表面不均勻，故在實際評定表面粗糙度時取評定長度為 5lr[18]，獲得4.27 mm×4.27 mm的評定區(qū)域。如圖5所示為利用Matlab軟件對實測離散點云數(shù)據進行處理所得到的評定區(qū)域表面形貌仿真結果，其中箭頭所示方向為輪廓截取方向，為了使仿真結果能夠更加清晰地表現(xiàn)，采用紅、藍兩色反映高度極值，紅色表示極大值，藍色表示極小值。

圖5　基于實測數(shù)據的表面形貌仿真

Fig.5Surface topography simulation based on actual measurement data

表1所示為8組輪廓數(shù)據經Talymap軟件進行高斯濾波后獲得的粗糙度值與本方法所獲得的粗糙度值的比較。通過分析表1中數(shù)據發(fā)現(xiàn)各組誤差均在允許范圍內，平均相對誤差值小于4%。因此，所提出的方法可以有效實現(xiàn)基于實測離散點云數(shù)據的多尺度輪廓形貌誤差成分表征。

表1白光干涉儀實測與仿真所得粗糙度值比對

Tab.1Comparison for WLI measurement roughness and simulation roughness

序數(shù)測量值/μm仿真值/μm相對誤差值/%10.501140.49191.84420.584690.59722.14030.475030.47990.33140.438970.41615.21050.548990.49859.19760.467360.50006.98470.541930.54230.06880.464380.44105.035

4　非理想表面模型建模實例

新一代GPS標準中的非理想表面模型是設計者頭腦中想象的能夠反映產品功能要求且考慮了加工誤差的幾何模型[3]。圖6所示為實現(xiàn)具有宏觀及微觀多尺度形貌誤差成分的非理想表面模型表征的整體設計思路。

圖6　生成非理想表面模型的過程

為了準確獲得二維輪廓的非理想表面模型，選擇表1中相對誤差值最小的第7組數(shù)據進行多尺度輪廓形貌誤差成分的提取及表征。圖7所示為基于第7組數(shù)據利用Matlab仿真得到的截面輪廓，圖中橫坐標為采樣點序數(shù)，縱坐標為采樣點高度值，單位為μm。

圖7　截面輪廓仿真

假設n為測量工件表面形貌時測量區(qū)域內某一截面輪廓上的采樣點數(shù)，則在進行小波變換時依據最大理論尺度有：

j=[lbn].

(5)

則表面最多能進行j次小波分解。為了使采樣結果更加穩(wěn)定、準確，故取原始輪廓數(shù)據上第2 049至3 072個采樣點，總計1 024個采樣點進行多尺度誤差成分的分離及提取。根據式(5)可知，采樣信號最多可進行10次小波分解。然而，當分解到第10層時僅剩下1個采樣點，以1個點代替一條輪廓曲線是不合理的，所以在此僅對采樣信號進行9層分解。如圖8所示為經過小波分解、重構、頻帶劃分及成分合成等處理后得到的表面輪廓的3種不同尺度形貌誤差成分，其中橫坐標為采樣點序數(shù)，縱坐標為不同成分所對應的高度值，單位為μm。

圖8　輪廓誤差成分多尺度表征

通過仿真分析得到了產品截面輪廓多尺度形貌誤差成分后，將各成分再進行合成，可最終得到如圖9所示的具有多尺度形貌誤差成分的產品的一個截面的二維輪廓的非理想表面模型。

圖9　非理想表面模型(二維輪廓)

輪廓是客觀評定表面的基礎，表面的功能特性在很大程度上取決于表面的多尺度形貌誤差[19]，因此可將此研究方法推廣至三維離散表面。由于組成圖5評定區(qū)域的每塊取樣區(qū)域的實測數(shù)據為1 024×1 024個采樣點的三坐標信息，這會造成進行三維表面形貌誤差成分表征時的計算量較大，因此首先對三維采樣數(shù)據進行數(shù)據抽稀，保留128×128個采樣點進行后續(xù)分析處理。為了使采樣數(shù)據更為準確，取評定區(qū)域中心一塊取樣區(qū)域作為研究對象，利用二維離散小波變換，通過分解、重構、頻帶劃分、成分合成等處理，得到如圖10所示的三維表面多尺度形貌誤差成分。其中，采樣點X及Y坐標表示沿x軸方向及y軸方向對應的采樣序數(shù)，縱坐標為采樣點高度值z，單位為μm。

(a)形狀誤差(a) Form error　(b)波紋度(b) Waveness　(c)粗糙度(c) Roughness

將表面誤差各尺度成分進行合成后，得到了如圖11所示的能夠表征具有宏觀及微觀多尺度形貌誤差成分的三維表面的非理想表面模型。

圖11　非理想表面模型(三維表面)

4　結　論

本文以實際工件表面為研究對象，利用白光干涉儀獲取了其測量區(qū)域內的三維離散點云采樣數(shù)據。針對對三維表面進行輪廓提取獲得的二維離散采樣數(shù)據，利用離散小波對二維采樣數(shù)據進行多尺度分解與重構，獲得了不同頻譜范圍內的輪廓形貌誤差成分信息。將Talymap軟件進行高斯濾波后獲得的粗糙度值與利用提出方法所獲得的粗糙度值比對，發(fā)現(xiàn)平均相對誤差值不超過4%，從而證明了研究方法對于二維輪廓形貌誤差成分的多尺度表征行之有效。通過將多尺度形貌誤差成分進行合成，得到了能夠表征工件截面輪廓宏觀及微觀多尺度形貌誤差的非理想表面模型，并成功推廣至工件三維表面。

為表征和仿真能夠反映加工表面多尺度宏觀及微觀形貌誤差信息的非理想表面模型提供了一種可行可靠的新方法，對于完善GeoSpelling模型、新一代GPS標準體系及推動計算機輔助公差技術的發(fā)展均具有重要作用。

[1]蔣向前. 新一代GPS標準理論與應用[M]. 北京：高等教育出版社，2007： 148-150.

JIANG X Q.TheoryandApplicationsofNew-GenerationGeometricalProductSpecifications[M]. Beijing： Higher Education Press， 2007： 148-150.(in Chinese)

[2]郭崇穎，劉檢華，唐承統(tǒng)， 等. 基于三維凸包的公差基準軸線擬合[J]. 光學 精密工程， 2014， 22(12)： 3247-3257.

GUO C Y， LIU J H， TANG CH T，etal.. Datum axis fitting for tolerances based on 3D convex hulls [J].Opt.PrecisionEng.， 2014， 22(12)： 3247-3257.(in Chinese)

[3]BALLU A， MATHIEU L， DANTAN J Y. Formal language for GeoSpelling [J].JournalofComputingandInformationScienceinEngineering， 2015， 15(2)： 1-6.

[4]ZHANG M， SHI ZH Y， MATHIEU L，etal.. Geometric product specification of gears： the GeoSpelling perspective [C].Procedingsofthe13thCIRPConferenceonComputerAidedTolerancing，Zhejiang，China， 2015， 27： 90-96.

[5]ISO 17450-1： 2011Geometric product specification(GPS)—General concepts—Part 1： Model for geometrical specification and verification [S].Geneva：InternationalOrganizationforStandardiza-tion， 2011.

[6]ZHANG M， ANWER N， MATHIEU L，etal.. A discrete geometry framework for geometrical product specifications [C].InProcedingsofthe21stCIRPDesignConference，Korea， 2011： 142-148.

[7]張敏，姜海洋，石照耀. 新一代GPS標準中非理想表面模型的建模方法[J]. 儀器儀表學報, 2016, 37(1): 91-98.

ZHANG M, JIANG H Y, SHI ZH Y. Modeling method for the skin model in new-generation geometrical product specifications [J].ChineseJournalofScientificInstrument, 2016, 37(1): 91-98. (in Chinese)

[8]SCHLEICH B， WARTZACK S. Tolerance analysis of rotating Mmchanism based on skin model shapes in discrete Ggometry [C].InProcedingsofthe13thCIRPConferenceonComputerAidedTolerancing，Zhejiang， 2014： 1-6.

[9]李海燕， 劉國棟， 劉炳國， 等. 雙密度小波在表面形貌信號分離中的應用[J]. 光學 精密工程， 2008， 16(6)： 1093-1097.

LI H Y， LIU G D， LIU B G，etal.. Application of double density wavelet transform to surface topographic signal separation [J].Opt.PrecisionEng.， 2008， 16(6)： 1093-1097.(in Chinese)

[10]張浩，袁怡寶，張峰. 巴特沃斯小波在表面形貌信號分離中的應用[J]. 光學 精密工程， 2010， 18(7)： 1661-1667.

ZHANG H， YUAN Y B， ZHANG F. Application of Butterworth wavelet to surface topographic signal separation [J].Opt.PrecisionEng.， 2010， 18(7)： 1661-1667.(in Chinese)

[11]任志英，高誠輝，申丁， 等. 雙樹復小波穩(wěn)健濾波在工程表面粗糙度評定中的應用[J]. 光學 精密工程， 2014， 22(7)： 1820-1827.

REN ZH Y， GAO CH H， SHEN D，etal.. Application of DT-XWT robust filtering to evaluation of engineering surface roughness [J].Opt.PrecisionEng.， 2014， 22(7)： 1820-1827.(in Chinese)

[12]LU W L， ZHANG G P， LIU X J. Prediction of Surface topography at the end of sliding running-in wear based on areal surface parameters[J].TribologyTransactions， 2014， 57(3)：553-560.

[13]曾文涵. 雙樹復小波表面分析模型及加工過程形貌辨識方法研究[D]. 武漢：華中科技大學， 2005.

ZENG W H.DualTreeComplexWaveletBasedSurfaceAnalysisModelandTopographyRecognitionforMachiningProcess[D]. Wuhan：Huazhong University of Science&Technology， 2005.(in Chinese)

[14]王明祥， 寧宇蓉， 王晉國. 基于Mallat算法的一維離散小波變換的實現(xiàn)[J].西北大學學報， 2006， 36(3)： 364-368.

WANG M X， NING Y R， WANG J G. Realization of one-dimentional discrete wavelet transform based on Mallat arithmetic [J].JournalofNorthwesternPolytechnicalUniversity，2006， 36(3)： 364-368.(in Chinese)

[15]張志航， 崔海， 丁海娟， 等. MWEDM表面三維粗糙度的小波評定基準面[J]. 哈爾濱工程大學學報， 2011， 32(9)： 1185-1189.

ZHANG ZH H， CUI H， DING H J，etal.. The reference plane by wavelets for 3D roughness evaluation of micro wire electrical discharge machining(MWEDM)[J].JournalofHarbinEngineeringUniversity， 2011， 32(9)： 1185-1189.(in Chinese)

[16]FU SH Y， MURALIKRISHNAN B， RAJA J. Engineering surface analysis with different wavelet bases [J].JournalofManufacturingScienceandEngineering， 2003， 125(4)： 844-852.

[17]GB/T 3505-2009/ISO 4287：1997： 產品幾何技術規(guī)范(GPS)表面結構 輪廓法 術語、定義及表面結構參數(shù)[S]. 北京：中國標準出版社， 2009.GB/T 3505-2009/ISO4287：1997： Geometrical Product Specification(GPS)-Surface texture： Profile method-Terms， definitions and surface texture parameters [S].Standardization Administration of the People’s Republic of China. GB/T 3505-2009. Beijing： Standards Press of China， 2009.(in Chinese)

[18]GB/T 1031-2009：產品幾何技術規(guī)范(GPS)表面結構 輪廓法 表面粗糙度參數(shù)及其數(shù)值[S]. 北京：中國標準出版社， 2009.

GB/T 1031-2009: Geometrical Product Specification(GPS)-Surface texture: Profile method-Surface roughness parameters and their values [S]. Beijing: Standards Press of China, 2009.(in Chinese)

[19]張華， 石照耀， 張白. 三維幾何形狀判別及其誤差評定研究[J]. 儀器儀表學報， 2014， 35(6)： 1217-1222.

ZHANG H， SHI ZH Y， ZH B. Study on 3D geometric shape discriminant and form error evaluation [J].ChineseJournalofScientificInstrument， 2014， 35(6)： 1217-1222.(in Chinese)

姜海洋(1989-)，女，吉林吉林人，碩士研究生，2012年于上海師范大學獲得學士學位，主要研究方向為產品幾何技術規(guī)范與認證，精密測量技術，計算機輔助公差設計。E-mail： jiangspeed@126.com

導師簡介：

石照耀(1964-)，男，湖南岳陽人，博士，教育部長江學者特聘教授，博士生導師，1984年于合肥工業(yè)大學獲得學士學位，1988年于陜西機械學院獲得碩士學位，2001年于合肥工業(yè)大學獲得博士學位，主要研究方向為齒輪工程及精密測試技術與儀器。E-mail： shizhaoyao@bjut.edu.cn

(版權所有未經許可不得轉載)

Multi-scale representation for skin models in geometrical product specifications

SHI Zhao-yao， JIANG Hai-yang， ZHANG Min

(BeijingEngineeringResearchCenterofPrecisionMeasurementTechnology&Instruments，BeijingUniversityofTechnology，Beijing100124，China)

*Correspondingauthor，E-mail：shizhaoyao@bjut.edu.cn

As most of the Computer Aided Tolerancing tools can only deal with the CAD models with an ideal surface and can not reflect manufacturing errors in physics and geometry, this paper explores a multi-scale representation method for skin model shapes in the Geometrical Product Specification. A discrete data modeling method was proposed for the simulation of topographic errors of surfaces and section profiles of a product in macroscopic and microcosmic scales based on new-generation geometrical product specification. Firstly, a simulation method for multi-scale surface topographic errors based on discrete wavelet was presented. Then, discrete wavelet was used to simulate multi-scale surface topography errors for sampling data of the surfaces and section profiles of a part. Finally, multi-scale surface topography errors were composed and skin models for two-dimentional profile and three-dimentional surface were acquired. The simulation and experiment results show that the proposed method represents skin models with multi-scale surface topographic errors and the average relevant error between the results ofRaobtained by a white-light interferometer and the proposed simulation method is less than 4%. The results verify the correctness and applicability of the proposed method, and provide a valid way for more comprehensive representation of skin model in new-generation geometrical product specification.

geometrical product specification; skin model; discrete geometry; wavelet analysis; multi-scale representation

2016-03-15；

2016-04-20.

國家教育部博士點基金資助項目(No.20131103110001)；國家自然科學基金資助項目(No.51305006)；5北京市教委項目(No.JC001013201402)

1004-924X(2016)07-1647-08

TP391.7；TH701

Adoi:10.3788/OPE.20162407.1647