張志遠,樓偉,蔡夕方,李斌,尹朝暉,王毅,高姍
(1.海軍海洋水文氣象中心,北京100161;2.清華大學計算機科學與技術(shù)系,北京100084;3.國家海洋環(huán)境預報中心,北京100081)
北印度洋風浪流數(shù)值預報系統(tǒng):II-檢驗分析
張志遠1,2,樓偉1,蔡夕方1,李斌1,尹朝暉3,王毅3,高姍3
(1.海軍海洋水文氣象中心,北京100161;2.清華大學計算機科學與技術(shù)系,北京100084;3.國家海洋環(huán)境預報中心,北京100081)
摘要:基于系統(tǒng)構(gòu)建工作[1],開展北印度洋風浪流數(shù)值預報系統(tǒng)后報和準業(yè)務(wù)化預報,并利用2013年9月—2014年3月共6個月的資料對預報結(jié)果進行了統(tǒng)計檢驗。結(jié)果顯示北印度洋風浪流數(shù)值預報業(yè)務(wù)運行穩(wěn)定可靠,大氣模式(WRF)72 h預報的500 hPa位勢高度距平相關(guān)系數(shù)達到89%,海浪模式(SWAN)的72 h有效波高預報的相對誤差低于20%,海流模式(ROMS)的72 h海表溫度預報的均方根誤差在0.5℃左右;同時對2013年10月期間孟加拉灣的超級氣旋風暴“PHAILIN”的預報結(jié)果進行了分析。該風、浪、流預報系統(tǒng)能夠較好地預報“PHAILIN”的移動路徑、最低氣壓及相應(yīng)的海浪和海流過程。該系統(tǒng)的試運行和檢驗分析結(jié)果,對建立新一代海洋環(huán)境數(shù)值預報系統(tǒng)具有一定借鑒意義。
關(guān)鍵詞:北印度洋;數(shù)值預報;模式檢驗;誤差分析
目前,數(shù)值預報產(chǎn)品已經(jīng)成為各級水文氣象部門發(fā)布天氣和海洋預報的主要依據(jù)。盡管數(shù)值預報模式日益完善,計算精度不斷提高,但現(xiàn)在數(shù)值預報模式還不能完全達到模擬真實海洋環(huán)境變化程度。而準確的業(yè)務(wù)預報、高質(zhì)量的服務(wù)都必須建立在盡可能準確的數(shù)值預報產(chǎn)品的基礎(chǔ)上的,為此,要進一步用好北印度洋風浪流數(shù)值預報系統(tǒng)的預報產(chǎn)品,必須對模式誤差進行統(tǒng)計分析,掌握預報規(guī)律,同時利用該系統(tǒng)對影響該地區(qū)的主要天氣系統(tǒng)活動的關(guān)鍵區(qū)的預報誤差進行研究。很多學者對北印度洋區(qū)域的數(shù)值預報產(chǎn)品的檢驗做了大量研究[2-3]。
一個數(shù)值預報系統(tǒng)從研制建立到正式運行階段中間,必須要有業(yè)務(wù)化試運行階段。在試運行階段,要從數(shù)值預報系統(tǒng)完整性、準確性、穩(wěn)定性、時效性、自動化和可視化程度等方面考核模式的成熟度,以便通過驗收進入正式業(yè)務(wù)化運行階段。預報系統(tǒng)的主要預報要素的準確性是否滿足業(yè)務(wù)化需求,主要根據(jù)后報檢驗報告和業(yè)務(wù)化試運行評估報告進行評價,因此業(yè)務(wù)運行單位必須開展預報結(jié)果精度檢驗評估工作。
本文從評估北印度洋風浪流數(shù)值預報系統(tǒng)模式性能出發(fā),對預報結(jié)果進行分析,為用戶提供預報產(chǎn)品準確度范圍和可信度,供其在使用數(shù)值預報產(chǎn)品時作為參考依據(jù)。本文一方面利用歷史觀測數(shù)據(jù)對系統(tǒng)進行多月份的預報結(jié)果統(tǒng)計檢驗,另一方面選取重要天氣系統(tǒng)開展典型個例檢驗。
2013年7月2日投人準業(yè)務(wù)運行以來,北印度洋風浪流數(shù)值預報系統(tǒng)已穩(wěn)定運行約1年時間,運行故障率低,出圖率達到100%。該系統(tǒng)一直不間斷地提供每日2次北印度洋區(qū)域大氣、海浪、海流要素場和形勢場預報,同時根據(jù)水文氣象保障部門的需要,提供有針對性的二次開發(fā)定制產(chǎn)品。系統(tǒng)可輸出0—72 h的各種預報產(chǎn)品,對特殊保障任務(wù),還可以提供有較好參考價值的72 h以上預報產(chǎn)品。由于本系統(tǒng)及時、準確、客觀定量的預報和高度的自動化運行服務(wù),使預報保障人員能夠得到更多、更詳細的預報信息,對日常水文氣象保障工作起到了較好的指導作用。
該系統(tǒng)由大氣模式WRF、海浪模式SWAN和海流模式ROMS構(gòu)成,具體模式設(shè)計結(jié)構(gòu)、分辨率、嵌套關(guān)系、驅(qū)動場、強迫場設(shè)置等設(shè)計實現(xiàn)內(nèi)容,請參見《北印度洋風浪流數(shù)值預報系統(tǒng):I-設(shè)計與實現(xiàn)》[1]。
2.1統(tǒng)計檢驗方法
數(shù)值預報產(chǎn)品檢驗一直是本領(lǐng)域研究人員非常觀測的研究熱點問題,本文給出較為常用的統(tǒng)計檢驗方法。
2.1.1大氣模式檢驗方法
本文主要選擇形勢場檢驗來評估WRF模式在指定統(tǒng)計范圍內(nèi)的系統(tǒng)性誤差、預報穩(wěn)定性和可用預報時效等。這里采用的是世界氣象組織WMO (World Meteorological Organization)基本系統(tǒng)委員會推薦的數(shù)值預報產(chǎn)品標準化檢驗方法即該區(qū)域500 hPa位勢高度預報距平相關(guān)系數(shù)(ACC)[4],如式(1):

該檢驗利用NCEP FNL(Final Global Analysis)全球大氣分析資料作為檢驗數(shù)據(jù),對預報產(chǎn)品進行檢驗。
2.1.2海洋模式檢驗方法
海洋模式的統(tǒng)計檢驗包括時空配準和統(tǒng)計計算兩部分。
首先,進行預報數(shù)據(jù)與觀測數(shù)據(jù)(比如高度計資料有效波高數(shù)據(jù)、XBT、ADCP數(shù)據(jù)等)的時間空間配準。(1)時間配準:因為觀測數(shù)據(jù)的觀測時間分布不規(guī)則,無法與整點時刻產(chǎn)生一一對應(yīng)關(guān)系,而海洋是慢變過程,在1個小時內(nèi),觀測數(shù)據(jù)所含物理要素的變化幅度不大。因此,本文將任意整點時刻前后0.5 h內(nèi)的觀測數(shù)據(jù)作為該整點時刻的數(shù)據(jù)。這樣,在時間維上,將觀測數(shù)據(jù)時間整點化。而預報數(shù)據(jù)是每小時整點時刻輸出的,這樣就實現(xiàn)了預報數(shù)據(jù)與觀測數(shù)據(jù)的時間配準;(2)空間配準:本文不再討論其本身的誤差。因此本文的空間配準方式是將格點化的預報數(shù)據(jù)插值到觀測數(shù)據(jù)點上。讀取0—24 h預報時段內(nèi)整點時刻的預報場,按照觀測數(shù)據(jù)的經(jīng)緯度坐標,將預報數(shù)據(jù)利用雙線性插值方法插到該點上,作為1對(觀測和預報)檢驗數(shù)據(jù)。所有1個月內(nèi)的預報檢驗數(shù)據(jù)作為一個基本統(tǒng)計數(shù)據(jù)集。
其次,按照上述步驟形成的0—24 h、24—48 h、48—72 h預報檢驗數(shù)據(jù)集。對每一個數(shù)據(jù)集,按照張志遠等[5]給出的相對誤差RE(Relative Error)和均方根誤差(RMS)定義予以統(tǒng)計計算。
2.2大氣預報效果檢驗
圖1是針對2013年9月—2014年3月做的500 hPa位勢高度距平相關(guān)系數(shù)統(tǒng)計。從結(jié)果上看,72 h內(nèi)的距平相關(guān)系數(shù)均在89%以上,該相關(guān)系數(shù)較高表明預報準確率較高。尤其是在冬季的預報效果更好,這與北印度洋的天氣變化規(guī)律基本相符。北印度洋冬季盛行東北季風,這段時期印度洋的風向以東北向為主,風力、風向穩(wěn)定,相應(yīng)距平相關(guān)系數(shù)值較高,9—10月份和3—4月份為夏季風向冬季風或冬季風向夏季風轉(zhuǎn)換的季節(jié),風力較小,風向多變不穩(wěn)定,相應(yīng)距平相關(guān)系數(shù)稍低。
從表1看出,除2014年9月份72 h的500 hPa位勢高度距平相關(guān)系數(shù)外,其他數(shù)據(jù)均在90%以上,這個數(shù)據(jù)可以證明大氣模式對該區(qū)域短期大氣環(huán)流的預報結(jié)果可信。

圖1 500 hPa位勢高度距平相關(guān)系數(shù)

表1 大氣模式WRF統(tǒng)計檢驗結(jié)果

表2 有效波高檢驗結(jié)果(檢驗數(shù)據(jù)源:Jason-2衛(wèi)星高度計數(shù)據(jù))
從表1看,本文還完成了對該區(qū)域大氣模式500 hPa風速、溫度的距平相關(guān)系數(shù)檢驗。對于海洋環(huán)境保障來講,不僅要關(guān)注中尺度天氣的溫壓風濕等要素情況,更重要是關(guān)注海面風場和潛熱、感熱等通量對海浪、海流模式的影響。因此,在后續(xù)的工作中,我們還將進行海面風場等指標的檢驗。
2.3海浪預報效果檢驗
該系統(tǒng)將HY-2A高度計資料中的有效波高數(shù)據(jù)同化到海浪模式中,這里主要利用Jason-2有效波高觀測值作為檢驗數(shù)據(jù)。統(tǒng)計結(jié)果見表2。秋冬季的印度洋24 h海浪有效波高的各月平均相對誤差為13.3%,各月平均均方根誤差均在0.40 m;48 h海浪有效波高的各月平均相對誤差為14.1%,各月平均均方根誤差均在0.41 m;72 h海浪有效波高的各月平均相對誤差為15.1%,各月平均均方根誤差均在0.44 m。
圖2給出了4條Jason-2高度計觀測軌跡上北印度洋有效波高的分布以及對應(yīng)時刻沿衛(wèi)星軌道模式值和觀測值的比較,圖中(a, c, e, g)分布是對應(yīng)時刻Jason-2高度計觀測軌跡,圖中(b, d, f, h)是對應(yīng)軌跡附近模式預報結(jié)果與觀測值分布示意圖。從圖中看出,模式較好地預報了北印度洋海區(qū)海浪的變化狀況。
2.4海流預報效果檢驗
本文選取2013年11月—2014年2月期間,在北印度洋海區(qū)獲得的XBT和ADCP資料對海流模式進行檢驗。檢驗要素主要包括印度洋海區(qū)的海溫剖面和0 m、50 m、75 m、100 m、150 m、200 m和500 m標準層的海流流速、流向。海溫剖面檢驗結(jié)果顯示,海表溫度的均方根誤差在0.5℃左右,絕對平均誤差在0.3℃。
如圖3所示,隨著水深增加,在通常的50—200 m左右的溫度躍層附近,均方根誤差和絕對平均誤差都有明顯的增大。在溫躍層以深水體中,海溫的均方根誤差或絕對平均誤差變化也較小。這顯示出在波動較大的溫躍層附近,預報誤差較大,海流模式在躍層附近的模擬精度和技巧還需要進一步提高。

圖2 不同時刻模式預報值與高度計觀測值的比較

圖3 利用XBT數(shù)據(jù)檢驗海流模式海溫誤差垂向分布圖
同時,從表3海流流速流向的檢驗結(jié)果看,海流流速的平均誤差呈現(xiàn)較穩(wěn)定的負偏差,而較大值也出現(xiàn)在50—200 m的躍層范圍內(nèi);海流流向的平均誤差均在15°以下,可以看出在100 m以淺的水深范圍內(nèi),海水的流動方向不穩(wěn)定,有各種上混合層和躍層物理過程相互影響流向平均誤差均在10°以上,而100 m以深的水深范圍內(nèi),海流的流動方向主要受大洋洋流趨勢影響,相對穩(wěn)定,因此誤差有明顯的降低。
除了完成統(tǒng)計檢驗外,對于數(shù)值預報效果的評估,更為直觀和有比對效果的是選取具有典型海區(qū)特點的重要天氣過程,比如臺風、熱帶氣旋等進行短期預報效果檢驗。本文選取北印度洋風浪流數(shù)值預報系統(tǒng)準業(yè)務(wù)化運行過程中2013年10月的超級氣旋風暴“PHAILIN”(02B)進行分析,檢驗大氣、海浪和海流模式對超級氣旋風暴的預報效果。

表3 海流流速流向檢驗結(jié)果(檢驗數(shù)據(jù)源:ADCP數(shù)據(jù))
3.1超級氣旋風暴“PHAILIN”
2013年10月出現(xiàn)的超級氣旋風暴“PHAILIN”是北印度洋海區(qū)有完整記載以來出現(xiàn)的最強熱帶氣旋。10月5日,在泰國灣海面上因強擾動生成熱帶低壓,并于7日進入安達曼海。9日聯(lián)合臺風警報中心JTWC(Joint Typhoon Warning Center)將其升格為熱帶風暴。在進入孟加拉灣后,“PHAILIN”沿著副熱帶高壓的南部邊緣向西北方向移動,并逐漸加強發(fā)展。并于9日下午并形成穩(wěn)定的風眼。10日下午2時升格其為三級熱帶氣旋,6 h后升格為四級熱帶氣旋。在11日中午12時達到最大強度,中心最低氣壓為940 hPa,風速達59 m/s,升格為五級熱帶氣旋。在12日17時“PHAILIN”以超級氣旋風暴的強度在印度奧里薩邦戈巴爾布爾附近登陸。登陸時中心最低氣壓為942 hPa,最大風速56m/s。由于孟加拉灣水汽充沛,導致“PHAILIN”引起的暴雨規(guī)模很大。13日減弱為熱帶風暴。
3.2“PHAILIN”的路徑預報
圖4給出大氣模式WRF模擬的超級氣旋風暴“PHAILIN”路徑和印度氣象局臺風預警中心的客觀定位路徑[6],模式超級氣旋風暴中心位置由整層平均的流場氣旋中心確定,從預報結(jié)果與客觀定位比較看,WRF模式預報的登陸地點比實際稍偏北,但與客觀定位非常接近。

圖4 超級氣旋風暴“PHAILIN”的路徑與印度氣象局臺風預警中心客觀定位路徑
本系統(tǒng)做出的24 h預報登陸點與印度氣象局臺風預警中心的客觀定位登陸點之間只差14.11 km,48 h預報登陸點差44.91 km,72 h預報登陸點與客觀定位登陸點差98.17 km。從結(jié)果看,24 h、48 h的登陸路徑預報效果較好,72 h稍差,整體看本系統(tǒng)WRF模式對超級氣旋風暴“PHAILIN”路徑模擬是成功的。
3.3“PHAILIN”的氣壓強度預報
氣旋風暴強度包括氣旋中心氣壓強度和氣旋風力預報。圖5是WRF模式模擬的氣旋中心氣壓強度變化與印度氣象局臺風預警中心發(fā)布的實況結(jié)果比較示意圖。結(jié)果表明,登陸前的強度變化比較小,登陸后中心氣壓強度迅速減弱,這與實況較為一致。24 h預報的絕對誤差為8.11 hPa,48 h預報的絕對誤差為8.45 hPa,72 h預報的絕對誤差為10.72 hPa。但WRF對氣旋登陸前后的模擬均偏弱,這與本文試驗中沒有人工干預采用Bogus技術(shù)[7]有關(guān)系。但WRF模擬的氣旋強度強弱的變化過程還是比較成功的。
從模擬臺風登陸前海平面氣壓和風矢量場,如圖6可見,風速在臺風中心的西北象限最大,呈NW-SE向非對稱分布。
加密觀測資料數(shù)值試驗研究[8]表明,呈NW-SE向非對稱分布的氣旋風暴往往使氣旋呈現(xiàn)偏北方向運動。而印度氣象局臺風預警中心發(fā)布的后期客觀定位路徑,如圖4所示,也證明了這一點。

圖5 模擬的超級氣旋風暴“PHAILIN”中心氣壓強度變化與實況比較示意圖
3.4海浪數(shù)值預報
氣旋能產(chǎn)生很強的海面風應(yīng)力,從而引發(fā)海上巨大的海浪,直接威脅海上作業(yè)和航行安全。圖7是模擬超級氣旋風暴“PHAILIN”過境時海浪響應(yīng)過程,預報結(jié)果顯示海浪浪高、浪向和浪高最大值的區(qū)域和變化趨勢與氣旋軌跡、中心強度等結(jié)構(gòu)特征相符。海浪的分布和演變受氣旋強度和移動的影響,浪高大值區(qū)與氣旋的風場大值區(qū)相對應(yīng),浪高的大小隨氣旋的增強而增大,氣旋減弱登陸后,海浪也相應(yīng)減小。氣旋風暴中心附近處,最大浪高超過14 m以上。而且可以較為明顯地看出海浪浪高大值區(qū)主要分布于氣旋行進方向的右側(cè)。
氣旋風暴的風場分布極大地影響了浪向的分布,在氣旋風暴中心的右側(cè),波浪一直處于強風速的持續(xù)作用下,風浪很大,波浪沿著氣旋風暴路徑的方向傳播,波向與風向基本一致;而在氣旋風暴中心的左側(cè),由于風速較右側(cè)小,涌浪占很大成分,波向由風浪和涌浪的方向共同決定,波浪的傳播方向與風向不一致。

圖6 2013101112時刻海平面氣壓和海面風場預報
3.5海流數(shù)值預報
氣旋風暴所帶來的強大風應(yīng)力,還給海洋上混合層和躍層的海水流動施加動量從而激起很強的海流運動,形成較強的風海流。圖7是模擬超級氣旋風暴“PHAILIN”過境過程中海表流場的響應(yīng)過程,預報結(jié)果顯示海流流速和流向,以及流速大值區(qū)域和變化趨勢與氣旋風暴軌跡、中心強度等結(jié)構(gòu)特征相符。隨著氣旋風暴的移動過程,我們發(fā)現(xiàn)氣旋風暴中心附近處出現(xiàn)了一個逆時針方向的較強的環(huán)流區(qū),逆時針方向的環(huán)流區(qū)隨著氣旋風暴中心的移動而移動,流速強度和范圍都隨臺風強度逐漸增大。氣旋風暴移動路徑的右側(cè)海表流速大于左側(cè),這與3.4節(jié)中有效波高的分布原因一致。

圖7 臺風“PHAILIN”過境時海浪浪高、浪向預報

圖8 2013年10月9日海表流速和流向預報
北印度洋風浪流數(shù)值預報系統(tǒng)自建立以來,保持準業(yè)務(wù)化運行,通過檢驗分析,本文得出如下結(jié)論:
(1)大氣模式WRF在北印度洋區(qū)域的模擬檢驗結(jié)果較好,臺站和浮標風的響應(yīng)分析反映出該模式在精細化預報方面具有較強的潛力。WRF在風場的模擬上具有較好的性能。所有檢驗月份的500 hPa的72 h內(nèi)的風速距平相關(guān)系數(shù)均大于0.89,滿足目前艦船航行對風速預測的精度要求。針對臺風等重要天氣過程,WRF模式能較好地預報出臺風的路徑、氣旋風暴中心最低氣壓和風力等情況,但數(shù)值預報結(jié)果相對實況整體偏弱,未來的業(yè)務(wù)化工作中,我們將在系統(tǒng)中加入人工Bogus方案[7],以便更好的完成“臺風”等重要天氣過程的模擬工作;
(2)海浪模式SWAN的模擬結(jié)果與觀測數(shù)據(jù)均基本吻合,均方根誤差時效及季節(jié)分布與絕對誤差分布基本一致,但同時次的值較絕對誤差略大,說明預報質(zhì)量隨預報時效增大而減弱。另外,風場精度的提高對海浪模式結(jié)果有重要影響。針對重要天氣系統(tǒng),可以較好地模擬臺風過程中的海浪變化情況,SWAN模擬臺風過程海浪的分布與臺風有較好的對應(yīng)關(guān)系,能較好地再現(xiàn)海浪的發(fā)展過程和合理地反映臺風浪的分布;
(3)海流模式ROMS的模擬結(jié)果與XBT和ADCP等觀測數(shù)據(jù)基本吻合,較大誤差主要出現(xiàn)在上混合層和躍層附近,海溫剖面的統(tǒng)計結(jié)果與美國
海軍的檢驗數(shù)據(jù)基本一致,說明我們的海流模式模擬效果較好。針對重要天氣系統(tǒng),風海流的響應(yīng)過程和程度都與氣旋風暴的發(fā)展吻合。因此,在未來工作中我們將結(jié)合收集到的各類觀測資料,對氣旋風暴引起的增水、海表溫度異常下降、垂直方向上各物理量的變化等情況進行分析檢驗。
參考文獻:
[1]蔡夕方,張志遠,樓偉,等.北印度洋風浪流數(shù)值預報系統(tǒng): I-設(shè)計與實現(xiàn)[J].海洋預報, 2014: 32(2): 7-13.
[2]齊鵬,范秀梅.高度計波高數(shù)據(jù)同化對印度洋海域海浪模式預報影響研究[J].海洋預報, 2013, 30(4): 70-78.
[3]楊永增,孫玉娟,王關(guān)鎖,等.基于MASNUM海浪預報系統(tǒng)的北印度洋波浪特征模擬與預報分析[J].海洋科學進展, 2011, 29(1): 1-9.
[4] Yang F L, Pan H L, Krueger S K, et al. Evaluation of the NCEP Global Forecast System at the ARM SGP Site[J]. Monthly Weather Review, 2006, 134(12): 3668-3690.
[5]張志遠,宋順強,劉利,等.浪流耦合模式數(shù)值模擬及檢驗分析[J].海洋技術(shù), 2011, 30(4): 87-92.
[6] Very Severe Cyclonic Storm, PHAILIN over the Bay of Bengal (08-14 October 2013): A Report[R]. New Delhi: Cyclone Warning Division, India Meteorological Department, 2013.
[7] Kurihara Y, Bender M A, Ross R J. An Initialization Scheme of Hurricane Models by Vortex Specification[J]. Monthly Weather Review, 1993, 121(7): 2030-2045.
[8]陳聯(lián)壽,羅哲賢.臺風科學、業(yè)務(wù)試驗和天氣動力學理論的研究(二)[M].北京:氣象出版社, 1996: 371-374.
North Indian Ocean wind-wave-circulation numerical forecast system: II-validation and analysis
ZHANG Zhi-yuan1,2,LOU Wei1,CAI Xi-fang1,LI Bin1,YIN Zhao-hui3,WANG Yi3,GAO Shan3
(1.Hydro-Meteorological Center of Navy, Beijing 100161 China; 2. Department of Computer Science and Technology, Tsinghua University, Beijing 100084 China; 3. National Marine Environmental Forecasting Center, Beijing 100081 China)
Abstract:Hindcasting and quasi-operating forecasting of the North Indian Ocean wind-wave-circulation numerical forecast system were implemented for system construction and some statistical tests and verification had been done using the 6 months(from September 2013 to March 2014)data in this paper. The results showed that the predictability and reliability of the system was perfect. The statistical results showed that the WRF simulated time series and trend analysis in 72hours of geopotential height anomaly correlation coefficient at 500 hPa was reached above 89%. The SWAN simulated result’s relative error of 72 hours of significant wave height (SWH) was less than 20%. The root mean square error (RMSE) of the sea surface temperature of 72 hour forecast results in ROMS was about 0.5℃. The verification of the case (the Very Severe Cyclonic Storm, PHAILIN) showed that the prediction of the track and the lowest central pressure of this storm, and the corresponding process of the wave and circulation were accurate. The validation and analysis of the wind-wave-circulation forecast system is expected to be a certain reference for the new generation of marine numerical prediction system.
Key words:NorthIndian Ocean;numerical forecast;model validation;error analysis
作者簡介:張志遠(1978-),男,工程師,博士,主要從事海洋環(huán)境信息化和數(shù)值預報研究。E-mail:generalzzy@139.com
基金項目:國家自然科學基金面上項目(41275098);國家海洋局海洋公益性行業(yè)科研專項(201005033)
收稿日期:2014-10-10
DOI:10.11737/j.issn.1003-0239.2015.03.007
中圖分類號:P731
文獻標識碼:A
文章編號:1003-0239(2015)03-0051-08