李志才，張 鵬，孫占義，王 凡

(1. 國(guó)家基礎(chǔ)地理信息中心大地測(cè)量部，北京 100830； 2. 導(dǎo)航與位置服務(wù)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室，北京 100830)

大規(guī)模GNSS基準(zhǔn)站網(wǎng)快速同步處理方法研究

李志才1，2，張 鵬1，2，孫占義1，2，王 凡1，2

(1. 國(guó)家基礎(chǔ)地理信息中心大地測(cè)量部，北京 100830； 2. 導(dǎo)航與位置服務(wù)國(guó)家測(cè)繪地理信息局重點(diǎn)實(shí)驗(yàn)室，北京 100830)

目前我國(guó)GNSS連續(xù)運(yùn)行基準(zhǔn)站網(wǎng)已基本建設(shè)完畢，全國(guó)范圍內(nèi)建設(shè)完成約2000個(gè)GNSS連續(xù)運(yùn)行基準(zhǔn)站。隨著站點(diǎn)規(guī)模的加大，數(shù)據(jù)計(jì)算的效率也迫切需要提高，采用傳統(tǒng)的高精度數(shù)據(jù)處理軟件已適應(yīng)不了大規(guī)模GNSS網(wǎng)的數(shù)據(jù)解算要求。本文基于BERNESE5.2軟件研究了我國(guó)“陸態(tài)網(wǎng)絡(luò)工程”260個(gè)基準(zhǔn)站的大規(guī)模GNSS網(wǎng)同步數(shù)據(jù)處理方法，通過修改源程序及利用并行計(jì)算技術(shù)，成功實(shí)現(xiàn)了陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站快速、高效、高精度的數(shù)據(jù)計(jì)算能力。實(shí)例驗(yàn)證表明，陸態(tài)網(wǎng)絡(luò)單天260個(gè)站的數(shù)據(jù)在無需分區(qū)的情況下，可在1 h內(nèi)獲得全球框架下的約束解，解算的框架點(diǎn)坐標(biāo)精度在毫米量級(jí)，大大提高了國(guó)家數(shù)據(jù)中心的大規(guī)模GNSS網(wǎng)數(shù)據(jù)處理能力。

參考框架；大規(guī)模GNSS網(wǎng)；并行計(jì)算；處理方法

隨著GNSS衛(wèi)星導(dǎo)航定位技術(shù)的快速發(fā)展，利用地面GNSS基準(zhǔn)站實(shí)現(xiàn)高精度實(shí)時(shí)導(dǎo)航定位的技術(shù)已經(jīng)在中國(guó)獲得大范圍的應(yīng)用，從國(guó)家到地方建設(shè)了近2000個(gè)GNSS連續(xù)運(yùn)行基準(zhǔn)站，其中國(guó)家級(jí)基準(zhǔn)站網(wǎng)的規(guī)模在200～300站之間(如“陸態(tài)網(wǎng)絡(luò)工程”260站、國(guó)家現(xiàn)代測(cè)繪基準(zhǔn)一期工程360站)，而省級(jí)CORS站網(wǎng)一般在100站以下，也有一些省級(jí)站規(guī)模超過100站以上(如山東省、甘肅省、黑龍江省CORS網(wǎng)等)，為我國(guó)高精度實(shí)時(shí)定位服務(wù)提供了起算基準(zhǔn)，也為維持我國(guó)2000國(guó)家坐標(biāo)系提供了更多的框架點(diǎn)[1-5]。這些國(guó)家和省級(jí)基準(zhǔn)站必須進(jìn)行精密解算，以實(shí)現(xiàn)我國(guó)坐標(biāo)框架的動(dòng)態(tài)維護(hù)，建立國(guó)家坐標(biāo)框架與國(guó)際坐標(biāo)框架、省級(jí)坐標(biāo)框架的聯(lián)系，而隨著基準(zhǔn)站網(wǎng)的規(guī)模越來越大，相應(yīng)的高精度數(shù)據(jù)處理的效率問題就迫切需要解決。當(dāng)前對(duì)大規(guī)模網(wǎng)的處理策略主流觀點(diǎn)是分區(qū)處理和分布式處理[6-9]。如對(duì)于100站網(wǎng)以內(nèi)的省級(jí)CORS中心，利用GAMIT的分區(qū)技術(shù)可以解決計(jì)算效率的問題[7]，但是隨著基準(zhǔn)站點(diǎn)的大幅增加，如采用GAMIT軟件進(jìn)行日常數(shù)據(jù)處理，即使分區(qū)解算效率仍難以有效解決，必須增加硬件來進(jìn)行彌補(bǔ)。根據(jù)國(guó)家數(shù)據(jù)處理中心的分析報(bào)告，如果采用GAMIT軟件進(jìn)行單日解的日常數(shù)據(jù)處理，260個(gè)國(guó)家站的單天數(shù)據(jù)處理大約需分6～7個(gè)區(qū)[10]，如果用單臺(tái)服務(wù)器總耗時(shí)大約24 h。分布式方法可以極大地提高計(jì)算效率，但是其方法的提出仍然是基于試驗(yàn)，并沒有在國(guó)內(nèi)大規(guī)模GNSS網(wǎng)中得到應(yīng)用，如何應(yīng)用分布式處理策略，還需盡快從研究轉(zhuǎn)入到實(shí)用中來[8]。從當(dāng)前國(guó)內(nèi)外的處理策略來看，對(duì)于大規(guī)模的GNSS基準(zhǔn)站網(wǎng)，主流仍然是采用了分區(qū)的處理策略，不管是各大分析中心還是國(guó)內(nèi)的數(shù)據(jù)處理中心，大都以分區(qū)處理策略為主。那么對(duì)于大規(guī)模的GNSS基準(zhǔn)站網(wǎng)，是否可以提出一種非分區(qū)的處理方法？本文正是基于這一想法，研究利用BERNESE 5.2軟件來實(shí)現(xiàn)大規(guī)模GNSS基準(zhǔn)站的數(shù)據(jù)處理能力[11]，解決我國(guó)陸態(tài)網(wǎng)絡(luò)大規(guī)模GNSS基準(zhǔn)站網(wǎng)日常處理需要分區(qū)的問題，為下一步測(cè)繪基準(zhǔn)工程360個(gè)站同步處理打下良好基礎(chǔ)。

1 大規(guī)模GNSS網(wǎng)快速同步處理方法

(2) 利用并行處理算法實(shí)現(xiàn)了計(jì)算效率的大幅度提升。BERNESE5.2內(nèi)核采用并行處理算法實(shí)現(xiàn)了數(shù)據(jù)的讀取、轉(zhuǎn)換、單點(diǎn)定位、基線模糊度固定、法方程計(jì)算等過程的并行計(jì)算，根據(jù)服務(wù)器的硬件性能，在分配任務(wù)時(shí)可進(jìn)行多時(shí)段、多任務(wù)并行處理，使得計(jì)算效率大大提高，大幅提高了數(shù)據(jù)處理速度。在實(shí)際解算過程中，采用了8 CPU、4核的計(jì)算服務(wù)器，可同時(shí)進(jìn)行32個(gè)并行任務(wù)，大大提高了計(jì)算效率。

(3) 采用了與GAMIT不同的平差處理策略。GAMIT軟件僅提供各單天解的自由網(wǎng)平差結(jié)果，而不是一個(gè)精密的數(shù)據(jù)處理成果，其單天解成果仍需要利用GLOBK軟件進(jìn)行約束來獲得多天解的網(wǎng)平差結(jié)果，這一步驟也非常消耗時(shí)間。根據(jù)BERNESE5.2的平差模塊，完善了BERNESE的單天解平差策略，解決了該軟件在平差過程中只能實(shí)現(xiàn)區(qū)域小網(wǎng)的平差，而不能進(jìn)行全球解的聯(lián)合平差問題。將其單天解的SINEX文件與全球IGS分析中心的SINEX進(jìn)行合并，以全球框架站為約束，實(shí)現(xiàn)了全球框架下單天解與全球網(wǎng)的平差，從而大大節(jié)約了網(wǎng)平差的計(jì)算時(shí)間，獲得了相對(duì)精確的單天解的精密成果。

2 陸態(tài)網(wǎng)絡(luò)工程GNSS網(wǎng)同步處理方法驗(yàn)證

上述方法需采用相應(yīng)的數(shù)據(jù)解算進(jìn)行驗(yàn)證，當(dāng)前國(guó)內(nèi)已建成的最大規(guī)模的國(guó)家級(jí)基準(zhǔn)站網(wǎng)是由地震局牽頭組織6個(gè)部委建設(shè)而成的陸態(tài)網(wǎng)絡(luò)工程，總設(shè)計(jì)點(diǎn)數(shù)260個(gè)，一般日常運(yùn)行的點(diǎn)數(shù)約為255個(gè)左右。主要應(yīng)用于地震監(jiān)測(cè)研究，可滿足測(cè)繪、導(dǎo)航、氣象、科研等需要，所有數(shù)據(jù)處理均為事后精密處理[4]。采用陸態(tài)網(wǎng)絡(luò)2015年整年的數(shù)據(jù)驗(yàn)證本文提出的整網(wǎng)快速同步處理計(jì)算方案。選擇中國(guó)大陸周邊6個(gè)IGS站進(jìn)行區(qū)域框架約束，在單天解平差過程中，下載了CODE分析中心的每天全球解SINEX文件進(jìn)行聯(lián)合平差，以全球框架站為約束，獲得單天解在全球基準(zhǔn)框架的平差結(jié)果。

2.1 數(shù)據(jù)準(zhǔn)備

必要的數(shù)據(jù)準(zhǔn)備是進(jìn)行GNSS精密數(shù)據(jù)處理的必要步驟，陸態(tài)網(wǎng)絡(luò)工程的數(shù)據(jù)準(zhǔn)備主要包括如下幾項(xiàng)工作：

(1) 準(zhǔn)備相關(guān)測(cè)站文件，包括測(cè)站信息表、測(cè)站代碼縮略表、測(cè)站所屬板塊表；測(cè)站信息表主要為測(cè)站代碼、測(cè)站接收機(jī)、天線型號(hào)、天線高。

(2) 原始數(shù)據(jù)，為節(jié)省空間，選擇壓縮的原始觀測(cè)D文件，利用CRX2RNX模塊將原始數(shù)據(jù)解壓成標(biāo)準(zhǔn)的RINEX O文件。

(3) 精密星歷文件，主要包括精密星歷文件EPH、精密鐘文件CLK、地球自轉(zhuǎn)參數(shù)文件ERP、全球電離層模型文件ION及CODE全球解的SNX文件。

(4) 基準(zhǔn)控制文件，根據(jù)不同年代的全球框架下載不同的框架文件，當(dāng)前國(guó)際主要框架為ITRF2008框架[12]，下載IGB08_R.CRD，IGB08_R.VEL文件，分別為全球框架下IGS站的坐標(biāo)和速度場(chǎng)，IGB08.FIX文件包括全球框架下的基準(zhǔn)點(diǎn)，這些文件可以將陸態(tài)網(wǎng)絡(luò)點(diǎn)的單天解由區(qū)域框架轉(zhuǎn)換到全球框架下。

(5) 各類地球物理模型及其他輔助文件，如大地水準(zhǔn)面模型文件EGM2008、章動(dòng)模型文件IAU2000R06.nut、IERS常量文件IERS2010XY.sub、海潮模型FES2004、地球固體潮模型TIDE2000、行星星歷等文件，以及接收機(jī)文件、天線文件、衛(wèi)星狀態(tài)文件。其中衛(wèi)星狀態(tài)文件需及時(shí)更新。

2.2 單天解數(shù)據(jù)計(jì)算方案

在拷貝必要的文件后，開始進(jìn)行單天解的數(shù)據(jù)處理，主要包括軌道文件生成、數(shù)據(jù)預(yù)處理、轉(zhuǎn)換和同步觀測(cè)數(shù)據(jù)、基線組合及相位數(shù)據(jù)預(yù)處理；整周模糊度消除、計(jì)算模糊度固定后的網(wǎng)解等。在基線解算時(shí)，主要采用的參數(shù)設(shè)置見表1。

表1 陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站單天解數(shù)據(jù)解算方案

整周模糊度消除策略中采用了偽距、相位兩種基線處理模式，分別采用了基于偽距的寬巷(WL)、基于相位的寬巷(L5)、準(zhǔn)無電離層影響(quasi-ionosphere-free，QIF)等模糊度消除方法，經(jīng)過多種模糊度處理策略，可大幅降低模糊度解，實(shí)現(xiàn)高精度數(shù)據(jù)處理結(jié)果。將這些模糊度消除后，單天解的RMS大部分在1 mm左右。在單天解計(jì)算過程中，為確保中國(guó)大陸的區(qū)域框架，并建立與全球框架的聯(lián)系，采用了中國(guó)大陸B(tài)JFS及周邊的DAEJ、IISC、POL2、TCMS、ULAB等6個(gè)國(guó)際IGS站進(jìn)行約束。

2.3 單天網(wǎng)解平差方案

陸態(tài)網(wǎng)絡(luò)工程作為中國(guó)目前一個(gè)較大的區(qū)域網(wǎng)，必須將其單天解結(jié)果轉(zhuǎn)換到全球框架下。下載了CODE數(shù)據(jù)處理分析中心的全球單日解的SINEX文件，通過相應(yīng)模塊將其轉(zhuǎn)換成法方程文件，然后與陸態(tài)網(wǎng)絡(luò)單天解的法方程疊加，通過固定前述6個(gè)IGS站及全球框架站，進(jìn)行區(qū)域單天解與全球框架單天解的融合來獲得單天網(wǎng)解。在基準(zhǔn)約束時(shí)，采用了IGB08框架下所有站作為基準(zhǔn)框架點(diǎn)，所有基準(zhǔn)點(diǎn)的先驗(yàn)σ0設(shè)為0.001 m，以進(jìn)行基準(zhǔn)約束，最小二乘約束采用全球七參數(shù)轉(zhuǎn)換方法，而不是區(qū)域網(wǎng)平差的三參數(shù)轉(zhuǎn)換方法。區(qū)域網(wǎng)的單天解通過6個(gè)固定的IGS站與IGB08框架進(jìn)行關(guān)聯(lián)并約束，最終得到陸態(tài)網(wǎng)絡(luò)全球IGB08框架下的單天網(wǎng)解及單天解法方程，其中單天解法方程可作為多天解的輸入值，最后進(jìn)行多天解的聯(lián)合平差，以獲得速度場(chǎng)結(jié)果。

2.4 單天網(wǎng)解質(zhì)量檢核

基線解算是GNSS數(shù)據(jù)處理中占用時(shí)間最長(zhǎng)、工作量最大的一部分，也是進(jìn)行網(wǎng)平差的基礎(chǔ)，基線解算質(zhì)量的好壞將直接影響GNSS網(wǎng)的定位精度和工作效率。而基線的解算質(zhì)量影響數(shù)據(jù)解算工作的結(jié)果和使用。NRMS值是衡量數(shù)據(jù)處理解算結(jié)果的一個(gè)重要指標(biāo)，用來表示單時(shí)段解算出的基線值偏離其加權(quán)平均值的程度，是從歷元的模糊度解算中得出的殘差。其計(jì)算公式如下

(1)

一般說來，NRMS值越小，基線估算精度越高；反之精度越低。根據(jù)國(guó)內(nèi)外GNSS數(shù)據(jù)處理經(jīng)驗(yàn)，其值一般應(yīng)小于3mm，若NRMS太大，則說明處理過程中周跳可能未得到完全修復(fù)。BERNESE的NRMS來自模糊度最終固定后的輸出文件驗(yàn)后單位權(quán)均方差(AposterioriRMSofunitweight)值。圖1給出了BERNESE解算陸態(tài)網(wǎng)絡(luò)觀測(cè)數(shù)據(jù)的NRMS值。其中NRMS的最大值為1.52mm)，最小值為1.37mm，平均值為1.42mm。從圖中可以很明顯地看出，2015年全年解的總體穩(wěn)定性較好，NRMS值全部都在1.5mm以內(nèi)，完全符合高精度數(shù)據(jù)處理的質(zhì)量要求。

圖1 陸態(tài)網(wǎng)絡(luò)基準(zhǔn)站單天解NRMS值

2.5 單天網(wǎng)解時(shí)間序列圖

BERNESE單天網(wǎng)解計(jì)算給出的坐標(biāo)為全球地心坐標(biāo)框架下當(dāng)日歷元的X、Y、Z值，一般需轉(zhuǎn)換到站心坐標(biāo)系下進(jìn)行時(shí)間序列的展示，利用黃立人等的方法[15]，將所有單天解的結(jié)果轉(zhuǎn)換到相應(yīng)的站心坐標(biāo)系下，并用GMT畫出各測(cè)站的E、N、U方向時(shí)間序列圖。圖2給出中國(guó)區(qū)域東北(吉林延吉JLYJ)、西北(新疆庫(kù)爾勒XJKE)、華北(北京十三陵BJSH)、華南(廣東湛江GDZJ)、華中(湖南瀏陽HNLY)、西南(四川爐霍SCLH)共6個(gè)站作為典型性代表，從各圖時(shí)間序列結(jié)果來看，單天網(wǎng)解在2015年呈現(xiàn)自西向東、自北向南的運(yùn)動(dòng)，與其他數(shù)據(jù)處理軟件結(jié)果基本一致，說明單天網(wǎng)解的計(jì)算結(jié)果比較可靠。水平方向呈現(xiàn)明顯的趨勢(shì)項(xiàng)，垂向上能顯示一定季節(jié)性變化，但規(guī)律并不明顯，可能與數(shù)據(jù)解算時(shí)間較短有關(guān)系，需要根據(jù)長(zhǎng)時(shí)間段的數(shù)據(jù)進(jìn)行檢查。

圖2 中國(guó)大陸部分典型測(cè)站時(shí)間序列

2.6 時(shí)間效率

經(jīng)計(jì)算驗(yàn)證，采用單臺(tái)8 CPU、4核服務(wù)器同步并行5個(gè)時(shí)段的數(shù)據(jù)解算，CPU最大并行計(jì)算設(shè)置為32組，則一天約260站的數(shù)據(jù)處理過程大概需要1 h，并可獲得各單天在全球ITRF2008框架下的精密坐標(biāo)解，其解算精度在1 mm左右，而相同情況的GAMIT軟件如不分區(qū)處理至少需要24 h，本文方法時(shí)間效率大大提高，且完全避免了分區(qū)帶來的基準(zhǔn)不一致的問題。

3 結(jié)束語

目前我國(guó)GNSS連續(xù)運(yùn)行基準(zhǔn)站網(wǎng)已基本建設(shè)完畢，隨著站點(diǎn)規(guī)模的加大，數(shù)據(jù)計(jì)算的效率也迫切需要提高，采用傳統(tǒng)的高精度數(shù)據(jù)處理軟件已適應(yīng)不了大規(guī)模GNSS網(wǎng)的數(shù)據(jù)快速解算要求。本文基于BERNESE5.2軟件研究了陸態(tài)網(wǎng)絡(luò)工程260個(gè)站的大規(guī)模GNSS網(wǎng)同步數(shù)據(jù)快速處理方法，通過修改源程序及利用并行技術(shù)，實(shí)現(xiàn)國(guó)家級(jí)基準(zhǔn)站快速、高效、高精度的數(shù)據(jù)計(jì)算能力，通過實(shí)例驗(yàn)證，單天260個(gè)站的數(shù)據(jù)無需分區(qū)就可同步進(jìn)行解算，完全避免了分區(qū)帶來的問題；通過并行計(jì)算快速獲得全球框架下高精度的單天約束解，解算時(shí)間可控制在1 h以內(nèi)，極大地提高了國(guó)家數(shù)據(jù)中心的大規(guī)模GNSS網(wǎng)數(shù)據(jù)處理能力。

[1] 張鵬,武軍酈,孫占義. 國(guó)家測(cè)繪基準(zhǔn)體系基礎(chǔ)設(shè)施建設(shè)[J].測(cè)繪通報(bào),2015(10):9-11.DOI: 10.13474/j.cnki.11-2246.2015.0302.

[2] 陳俊勇.中國(guó)現(xiàn)代大地基準(zhǔn):中國(guó)大地坐標(biāo)系統(tǒng)2000(CGCS2000)及其框架[J].測(cè)繪學(xué)報(bào),2008,37(3): 269-271.

[3] 陳俊勇, 張鵬, 武軍酈, 等.關(guān)于在中國(guó)構(gòu)建全球?qū)Ш叫l(wèi)星國(guó)家級(jí)連續(xù)運(yùn)行站系統(tǒng)的思考[J]. 測(cè)繪學(xué)報(bào)，2007,36(4):366-369.

[4] 游新兆,楊少敏,王琪. 中國(guó)大陸地殼運(yùn)動(dòng)觀測(cè)研究進(jìn)展[J].地震，2012，32(2)：31-39.

[5] 李強(qiáng)，游新兆，楊少敏，等. 中國(guó)大陸構(gòu)造變形高精度大密度GPS監(jiān)測(cè)現(xiàn)今速度場(chǎng)[J].中國(guó)科學(xué)：地球科學(xué)，2012，42(5)：629-632.

[6] GE M,DICK G，DICK G,et al.A New Data Procesing Strategy for Huge GNSS Global Network [J].Journal of Geodesy,2006,80(4):199-203.

[7] 姜衛(wèi)平，趙倩，劉鴻飛，等. 子網(wǎng)劃分在大規(guī)模GNSS基準(zhǔn)站網(wǎng)數(shù)據(jù)處理中的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2011，36(4)：389-391.

[8] 陳正生，呂志平，崔陽，等. 大規(guī)模GNSS數(shù)據(jù)的分布式處理與實(shí)現(xiàn)[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2015，40(3)：384-389.

[9] 程傳錄,蔣光偉,聶建亮,等.利用雙差的超大GNSS基準(zhǔn)站網(wǎng)解算方法改進(jìn)[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2014，39(5)：596-599.

[10] 張鵬，李志才，孫占義，等. 我國(guó)高精度地心坐標(biāo)框架確定方法研究[J].測(cè)繪通報(bào)，2015(12)：1-4.

[11] DACH R,HUGENTOBLER U,FRIDEZ P,et al. Bernese GPS Software Version 5.0[R].Bern:AIUB，2007.

[12] ALTAMIMI Z,COLLILIEUX X,MéTIVIER L.ITRF2008: An Improved Solution of the International Terrestrial Reference Frame[J]. Journal of Geodesy，2011,85(8):457-473.

[13] 黃立人，高硯龍，任立生. 關(guān)于NEU(ENU)坐標(biāo)系統(tǒng)[J].大地測(cè)量與地球動(dòng)力學(xué)，2006，26(1)：97-99.

Study on Fast Processing Method for Large Scale GNSS Network

LI Zhicai1,2，ZHANG Peng1,2，SUN Zhanyi1,2，WANG Fan1,2

(1. Geodesy Department, National Geomatics Center of China, Beijing 100830, China; 2. Key Laboratory of Navigation & Location Based Service,National Administration of Surveying, Mapping and Geoinformation, Beijing 100830, China)

At present, there are almost 2000 GNSS CORS stations across China as the construction of GNSS continuous operation reference station network nearly completed. With the increase of the amount of the sites, the efficiency of GNSS data processing urgently need to improve, while the traditional high precision data processing software represent the deficiency to compute so many more data. Based on BERNESE5.2 software, the data processing method of large scale GNSS network amount to 260 stations by CMONOC project is studied. The fast, efficient and high precision data computing capability had been got by modifying source program and utilizing the parallel computation technology. The experiment shows that the data of 260 stations in a single day can be processed within 1 hour without divided into sub-net to get the daily constrained coordinate in ITRF2008 frame with a millimeter scale accuracy, which highly improves the data processing capability of large scale GNSS network for national data center.

reference frame；large scale GNSS network; parallel computation; processing method

李志才，張鵬，孫占義，等.大規(guī)模GNSS基準(zhǔn)站網(wǎng)快速同步處理方法研究[J].測(cè)繪通報(bào)，2017(2)：65-69.

10.13474/j.cnki.11-2246.2017.0050.

2016-06-28

國(guó)家863計(jì)劃(2013AA122501)；國(guó)家測(cè)繪地理信息局公益性行業(yè)專項(xiàng)(201512001；201512004)；國(guó)家自然科學(xué)基金(41274030)

李志才(1976—)，男，博士，高級(jí)工程師，主要從事衛(wèi)星大地測(cè)量研究。E-mail：zcli@nsdi.gov.cn

P228.4

A

0494-0911(2017)02-0065-05