韓曉磊 張慶君 劉杰 張潤寧 袁智

(1 北京空間飛行器總體設(shè)計部，北京 100094)(2 航天東方紅衛(wèi)星有限公司，北京 100094)

敏捷SAR衛(wèi)星聚束模式姿態(tài)機動策略研究

韓曉磊1張慶君1劉杰1張潤寧2袁智1

(1 北京空間飛行器總體設(shè)計部，北京 100094)(2 航天東方紅衛(wèi)星有限公司，北京 100094)

提出了敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法，充分考慮了相位中心位置和離軸角的影響，通過矢量計算方法得到敏捷SAR衛(wèi)星聚束模式成像所需的初始姿態(tài)角，通過迭代計算消除相位中心位置影響，得到準(zhǔn)確姿態(tài)角。文章通過Matlab和STK軟件聯(lián)合仿真，驗證了上述計算方法的有效性。

合成孔徑雷達；敏捷SAR衛(wèi)星； 聚束模式；姿態(tài)機動

1 引言

合成孔徑雷達(SAR)是一種微波成像雷達，可以對地面目標(biāo)進行高分辨率的探測，具有全天候、全天時的優(yōu)勢。分辨率是星載SAR 系統(tǒng)的關(guān)鍵性指標(biāo)之一，高分辨率可以提高SAR對目標(biāo)的檢測和識別性能，擴展SAR系統(tǒng)的應(yīng)用領(lǐng)域。根據(jù)SAR工作原理，距離向分辨率可以通過增加發(fā)射信號帶寬提高；方位向分辨率可以通過減小方位向天線孔徑長度提高[1]。但隨著方位孔徑長度減小，天線的增益也隨之下降，受限于系統(tǒng)靈敏度決定的固定功率孔徑積，不能無限制地通過減小天線孔徑提高星載SAR系統(tǒng)方位向分辨率。而另一種實現(xiàn)方位向高分辨率的方式是通過雷達波束掃描增加合成孔徑時間，實現(xiàn)方位向高分辨率成像，由此形成了一種常用的高分辨率星載SAR工作模式：聚束模式。它通過雷達波束掃描使波束中心始終瞄準(zhǔn)地面固定點，從而使成像區(qū)域始終處于雷達波束的照射之下，延長了合成孔徑時間，實現(xiàn)了方位向高分辨率成像[2]。德國的TerraSAR-X衛(wèi)星在設(shè)計壽命末期進行了聚束模式試驗，它通過將方位向波束掃描范圍擴展到±2.2°，在軌實現(xiàn)了聚束模式，將TerraSAR-X衛(wèi)星的最高方位向分辨率從1.1 m提升到約0.2 m[3]。

敏捷SAR衛(wèi)星有效載荷固定在衛(wèi)星上，依靠姿態(tài)控制系統(tǒng)控制衛(wèi)星整體繞俯仰、橫滾、偏航3個軸向擺動，實現(xiàn)SAR成像所需的波束掃描。敏捷SAR衛(wèi)星具有較高的靈活性和成像能力，能高效地實現(xiàn)高分辨率聚束模式成像[4-6]。文獻[6]給出了敏捷SAR衛(wèi)星滑動聚束模式對俯仰姿態(tài)機動控制的要求，文獻[7]給出了輕型SAR衛(wèi)星聚束模式姿態(tài)快速機動和高精度追蹤控制策略及實現(xiàn)方式。但上述文獻的分析結(jié)果都是基于直線三角幾何關(guān)系得到，未考慮衛(wèi)星軌道和地球表面的彎曲特性影響，而且未考慮雷達天線相位中心和波束指向偏離星本體坐標(biāo)系原點和坐標(biāo)軸的影響，不適用于高分辨率聚束模式SAR衛(wèi)星。本文從精確的軌道和地球模型出發(fā)，提出了一種基于迭代計算的敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法，該方法能適應(yīng)SAR天線在星體上各種不同的安裝位置，并充分考慮了可能存在的波束離軸角。

2 衛(wèi)星坐標(biāo)系及姿態(tài)描述方法

2.1 坐標(biāo)系

為研究敏捷SAR衛(wèi)星聚束模式的姿態(tài)機動參數(shù)需求，須定義如下坐標(biāo)系[8]：地心慣性坐標(biāo)系OE-XIYIZI，地心固連坐標(biāo)系OE-XEYEZE，衛(wèi)星軌道坐標(biāo)系OS-XOYOZO，衛(wèi)星本體坐標(biāo)系OS-XBYBZB和天線坐標(biāo)系OP-XAYAZA等，各坐標(biāo)系之間的關(guān)系如圖1所示。

圖1 坐標(biāo)系定義示意圖Fig.1 Schematic diagram of definition of the coordinate systems

2.2 姿態(tài)描述方法

衛(wèi)星本體坐標(biāo)系相對于參考坐標(biāo)系的姿態(tài)描述有多種方法，如方向余弦矩陣、歐拉角、歐拉軸/角參數(shù)、歐拉四元數(shù)、羅德里格參數(shù)等。本文用到其中3種常用的衛(wèi)星姿態(tài)描述方法，即歐拉軸/角參數(shù)、歐拉角和歐拉四元數(shù)。實際應(yīng)用中需要在上述姿態(tài)描述方法之間進行轉(zhuǎn)換，歐拉軸/角參數(shù)到歐拉四元數(shù)的轉(zhuǎn)換關(guān)系如下[9]：

(1)

歐拉四元數(shù)到3-1-2轉(zhuǎn)序下歐拉角的轉(zhuǎn)換關(guān)系如下[10]：

(2)

3-1-2轉(zhuǎn)序下，歐拉角到歐拉四元數(shù)的轉(zhuǎn)換關(guān)系如下[10]：

(3)

式中：ex、ey、ez分別為轉(zhuǎn)軸方向矢量e在參考坐標(biāo)系中的三個方向余弦，φ為繞轉(zhuǎn)軸的轉(zhuǎn)動角度，q1、q2、q3和q4為歐拉四元數(shù)，φ、θ和ψ為3-1-2轉(zhuǎn)序下繞X軸、Y軸和Z軸轉(zhuǎn)動的歐拉角。而其他轉(zhuǎn)序下的轉(zhuǎn)換關(guān)系可以類推得到。

3 敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法

敏捷SAR衛(wèi)星一般采用整星機動實現(xiàn)雷達波束方位向掃描，不具備方位向電掃描能力。但距離向一般具備一定的電掃描能力，成像過程中存在波束指向偏離星本體坐標(biāo)系Z軸的離軸角。此外，敏捷SAR衛(wèi)星天線安裝方式一般不能保證天線相位中心與衛(wèi)星質(zhì)心重合，該中心的軌道系坐標(biāo)隨整星姿態(tài)機動發(fā)生變化，傳統(tǒng)衛(wèi)星姿態(tài)機動參數(shù)計算方法沒有考慮這個變化的影響，將造成指向偏差，如圖2所示。因此，有必要針對敏捷SAR衛(wèi)星聚束模式自身特點，開發(fā)適合它的姿態(tài)機動參數(shù)計算方法，以減小甚至消除姿態(tài)指向偏差。

圖2 算法必要性說明示意圖Fig.2 Schematic diagram of necessity of the algorithm

本節(jié)給出一種新的敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法，該方法包含3個主要步驟，分別為瞄準(zhǔn)點和時序規(guī)劃、需求姿態(tài)角粗算和需求姿態(tài)角精算，如圖3所示。

圖3 計算流程圖Fig.3 Block diagram of the computation algorithm

1)第1步，瞄準(zhǔn)點和時序規(guī)劃

聚束模式工作時雷達波束中心始終瞄準(zhǔn)地面待觀測場景的幾何中心，因此選擇該幾何中心作為地面瞄準(zhǔn)點。首先將瞄準(zhǔn)點經(jīng)緯度坐標(biāo)變換成地心慣性系坐標(biāo)，根據(jù)此坐標(biāo)及星歷數(shù)據(jù)，尋找瞄準(zhǔn)點回波多普勒中心為零的星歷時刻，多普勒中心計算方法如下[11]：

(4)

式中：Rs為地心慣性坐標(biāo)系下衛(wèi)星的位置矢量，Rt為地心慣性坐標(biāo)系下瞄準(zhǔn)點的位置矢量，Vs為地心慣性坐標(biāo)系下衛(wèi)星的速度矢量，Vt為地心慣性坐標(biāo)系下瞄準(zhǔn)點的速度矢量，λ為載波波長，Rst為衛(wèi)星與瞄準(zhǔn)點之間的距離。以瞄準(zhǔn)點回波多普勒中心為零的星歷時刻t0為整個成像時序的中心。此步驟將傳統(tǒng)SAR衛(wèi)星的姿態(tài)導(dǎo)引包含到成像姿態(tài)機動參數(shù)計算中，保證了回波信號多普勒中心為0Hz，因此，敏捷SAR衛(wèi)星聚束模式無須單獨進行姿態(tài)導(dǎo)引。

根據(jù)需求方位向分辨率ρa和多普勒調(diào)頻率fa計算成像時間為

(5)

式中：Vg為零多普勒線掃過地面的速度，kwa為多普勒信號處理加權(quán)擴展因子，fa為多普勒調(diào)頻率。

(6)

最后，將地面瞄準(zhǔn)點和地心坐標(biāo)變換到軌道坐標(biāo)系，為后續(xù)姿態(tài)角計算提供輸入。

2)第2步，需求姿態(tài)角粗算

計算需求姿態(tài)角時，應(yīng)按照時序規(guī)劃得到的成像星歷時刻逐點進行計算。天線相位中心測定后，可以計算出天線相位中心在衛(wèi)星本體坐標(biāo)系下的坐標(biāo)，但因姿態(tài)角未知，還不能得到天線相位中心在軌道坐標(biāo)系下的精確坐標(biāo)。粗算階段，以其在衛(wèi)星本體坐標(biāo)系下的坐標(biāo)代替軌道坐標(biāo)系坐標(biāo)，由于上述兩個坐標(biāo)系原點同為衛(wèi)星質(zhì)心，且星體尺寸一般較小，兩坐標(biāo)系下的相位中心坐標(biāo)差一般為米級，可以用于姿態(tài)角粗算，這種近似引入的誤差將在精算階段予以消除。

如圖4所示，根據(jù)天線相位中心OP和地面瞄準(zhǔn)點T在軌道系下的坐標(biāo)，可得到期望的波束中心指向向量，以VPT表示。根據(jù)衛(wèi)星質(zhì)心OS和天線相位中心OP在軌道系下的坐標(biāo)，可得到衛(wèi)星質(zhì)心到相位中心的向量，以VSP表示。由于姿態(tài)機動時是星體圍繞衛(wèi)星質(zhì)心OS轉(zhuǎn)動，應(yīng)以衛(wèi)星質(zhì)心為起點計算期望指向向量，以VST表示。

(7)

圖4 需求姿態(tài)角計算方法示意圖Fig.4 Schematic diagram of algorithm for computing the required attitude angle

根據(jù)地心OE在軌道系下的坐標(biāo)，可得到起始Z軸指向向量，以VSE表示。由于雷達波束存在離軸角，根據(jù)旋轉(zhuǎn)變換關(guān)系，可得到波束起始指向向量VSL在軌道系下的投影矩陣VSL_T，具體如下：

(8)

式中：C(θ)為旋轉(zhuǎn)轉(zhuǎn)換矩陣，θ為離軸角，VSE_T為VSE在軌道系下的投影矩陣。

根據(jù)歐拉定理，衛(wèi)星繞質(zhì)心的姿態(tài)機動可以由繞通過質(zhì)心的某一固定軸轉(zhuǎn)動一定角度得到。該轉(zhuǎn)軸方向矢量用e表示，轉(zhuǎn)動角度用φ表示，這即是歐拉軸/角參數(shù)式。根據(jù)最短路徑旋轉(zhuǎn)準(zhǔn)則，可由VSL和VST計算得到[12]

(9)

(10)

根據(jù)式(1)可得到歐拉四元數(shù)，再利用式(2)得到3-1-2轉(zhuǎn)序下的初始歐拉角[9]。由于本階段計算過程中，使用天線相位中心在衛(wèi)星本體坐標(biāo)系的坐標(biāo)代替軌道坐標(biāo)系坐標(biāo)，所以初始歐拉角中包含一定的誤差，下一步將重點剔除這個誤差。

3)第3步，需求姿態(tài)角精算

利用上一步得到的初始歐拉角，計算星本體坐標(biāo)系到軌道坐標(biāo)系的轉(zhuǎn)換矩陣，并將天線相位中心星本體坐標(biāo)系坐標(biāo)變換成軌道坐標(biāo)系坐標(biāo)。然后以此坐標(biāo)為輸入，重復(fù)第2步的計算過程，得到新的歐拉角。比較新歐拉角與前一次計算得到的歐拉角之間的差值，如大于預(yù)先設(shè)定的門限值，則用新歐拉角代替前一次計算得到的歐拉角，重復(fù)上述過程，直到新歐拉角與前一次計算得到的歐拉角之間的差值小于預(yù)先設(shè)定的門限值。

由于計算新歐拉角所用的轉(zhuǎn)換矩陣總是比前一次計算歐拉角時的轉(zhuǎn)換矩陣更精確，因此，新歐拉角逐漸向準(zhǔn)確的需求姿態(tài)角逼近，而新歐拉角與前一次計算得到的歐拉角之差逐漸向零收斂。通過選擇合適的門限值，能得到滿足精度需求的姿態(tài)角。這個迭代計算過程有效剔除了天線相位中心與質(zhì)心不重合導(dǎo)致的指向誤差，得到滿足精度要求的聚束模式成像需求姿態(tài)角。利用姿態(tài)運動學(xué)方程，可進一步得到需求的角速度三軸分量。

4 設(shè)計實例及仿真驗證

通過Matlab和STK軟件聯(lián)合仿真，驗證上述計算方法的有效性，驗證流程如圖5所示。首先根據(jù)輸入條件，按本文給出的計算方法，在Matlab軟件中計算姿態(tài)機動參數(shù)，然后將姿態(tài)機動參數(shù)計算結(jié)果導(dǎo)入STK軟件，顯示姿態(tài)機動和波束指向變換，并將地面瞄準(zhǔn)點位置信息輸回Matlab軟件，分析STK軟件中得到的瞄準(zhǔn)點位置信息與期望的瞄準(zhǔn)點位置信息之間的差異，并給出相應(yīng)性能分析結(jié)果。由于STK軟件的正確性已經(jīng)經(jīng)過工程驗證，它的分析結(jié)果能真實地反映衛(wèi)星運動規(guī)律，故上述仿真實驗具有較高的可信度。

圖5 仿真驗證流程圖Fig.5 Block diagram of the simulation and validation

本次仿真實驗的SAR衛(wèi)星聚束模式輸入?yún)?shù)如表1所示，根據(jù)本文給出的計算方法，得到實現(xiàn)聚束模式的姿態(tài)機動參數(shù)，圖6為以軌道系為參考系的需求歐拉四元數(shù)，圖7為以地心慣性系為參考系的需求歐拉四元數(shù)，圖8為以軌道系為參考系，3-1-2轉(zhuǎn)序下的需求歐拉角?？梢钥吹骄凼Ｊ匠上襁^程中俯仰機動范圍最大，達到1.883 5°，偏航機動次之，達到0.476 8°，橫滾機動最小，為0.004 8°，實際工程中此方向可以保持不動。姿態(tài)角速度在三軸的角速度分量如圖9所示。

表1 輸入?yún)?shù)表Table 1 Input parameters

圖6 以軌道系為參考系的需求歐拉四元數(shù)Fig.6 Required Euler quaternion referring to the orbit coordinate system

圖7 以地心慣性系為參考系的需求歐拉四元數(shù)Fig.7 Required Euler quaternion referring to the earth inertial coordinate system

圖8 以軌道系為參考系3-1-2轉(zhuǎn)序下的需求歐拉角Fig.8 Required Euler angle referring to the orbit coordinate system with rotation sequence 3-1-2

圖9 姿態(tài)角速度三軸分量Fig.9 Angle velocity components on the three axes

將上述姿態(tài)角導(dǎo)入STK軟件，通過STK中虛擬SAR衛(wèi)星的雷達波束指向，驗證上述計算結(jié)果的正確性。STK軟件中虛擬衛(wèi)星同樣采用表1中的輸入?yún)?shù)，得到成像起始、中間和結(jié)束時刻波束指向，以及地面瞄準(zhǔn)點位置，如圖10所示。其中黃色叉形區(qū)域為期望的成像場景中心，可見雷達波束在成像起始、中間和結(jié)束時刻都準(zhǔn)確指向期望的場景中心，證明了上述姿態(tài)機動參數(shù)計算結(jié)果的正確性。再將STK軟件中得到的實際瞄準(zhǔn)點位置信息導(dǎo)回到Matlab軟件中進行分析，得到瞄準(zhǔn)誤差，如圖11所示?？梢娒闇?zhǔn)誤差在10-8m量級，可以忽略不計。

根據(jù)得到的姿態(tài)信息，計算成像場景中心點的距離徙動曲線，如圖12所示。成像場景中心點的瞬時多普勒曲線如圖13所示，可見場景中心點的多普勒帶寬為12 625 Hz，零多普勒線掃過地面的速度為6894 m/s，計算得到場景中心點分辨率為0.6 m，滿足設(shè)計要求，證明了上述計算過程的有效性。

為了進行對比分析，本文在相同輸入?yún)?shù)條件下，采用傳統(tǒng)方法計算得到聚束模式姿態(tài)機動參數(shù)。然后，同樣利用STK軟件進行分析，最終得到的瞄準(zhǔn)誤差如圖14所示，其中圖14(a)為沒有考慮離軸角和天線安裝位置時，得到的瞄準(zhǔn)誤差；圖14(b)為沒有考慮天線安裝位置時，得到的瞄準(zhǔn)誤差。由圖14可見，敏捷SAR衛(wèi)星實現(xiàn)聚束模式時必須考慮離軸角，而當(dāng)精度要求較高時，還需要考慮天線安裝位置。上述仿真實驗證明了本文提出的敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法的必要性和有效性。

圖10 STK軟件中虛擬SAR衛(wèi)星的波束指向Fig.10 Beam direction of virtual satellite in STK

圖11 地面瞄準(zhǔn)誤差Fig.11 Aiming error on the ground

圖12 場景中心點的距離徙動曲線Fig.12 Range migration graph of central point of scene

圖13 場景中心點的瞬時多普勒頻率Fig.13 Instantaneous Doppler frequency graph of central point of scene

圖14 傳統(tǒng)姿態(tài)機動需求計算方法得到的瞄準(zhǔn)誤差Fig.14 Aiming error of conventional algorithms computing attitude maneuver requirement

5 結(jié)論

本文提出了一種基于迭代計算的敏捷SAR衛(wèi)星聚束模式姿態(tài)機動參數(shù)計算方法，該方法能適應(yīng)SAR天線在星體上各種不同的安裝位置，并充分考慮了可能存在的波束離軸角，得到了準(zhǔn)確的敏捷SAR衛(wèi)星聚束模式成像所需的姿態(tài)參數(shù)，經(jīng)試驗驗證，地面瞄準(zhǔn)誤差可達到10-8m量級，可為基于衛(wèi)星平臺姿態(tài)機動實現(xiàn)高分辨率星載SAR成像提供支撐。此外，本文提出的方法還可以應(yīng)用于相控陣體制高分辨率SAR衛(wèi)星，使其可以通過整星機動減小成像所需的天線電掃描范圍，降低對相控陣天線T/R組件的需求，實現(xiàn)整星質(zhì)量和成本的下降。

References)

[1] J C Curlander，R N McDonough. Synthetic aperture radar: systems and signal processing[M]. Washington D.C.: John Wiley & Sons,Inc,1991

[2]W Carrara,R Goodman，R Majewski. Spotlight synthetic aperture radar signal processing algorithm[M]. Boston: Artech House,1995

[3]J Mittermayer,S Wollstadt,P Prats Iraola，et al. The TerraSAR-X staring spotlight mode concept[J]. IEEE Transactions on Geoscience and Remote Sensing,2014,52(6): 3695-3706

[4]韓曉磊,李世強,王宇,等. 基于敏捷衛(wèi)星平臺的星載SAR Mosaic模式研究[J]. 宇航學(xué)報，2013,34(7): 971-979

Han Xiaolei,Li Shiqiang,Wang Yu,et al. A study on space-borne SAR Mosaic mode for an agile satellite [J]. Journal of Astronautics,2013,34(7): 971-979 (in Chinese)

[5]Lemaitre M,Verfaillie G,Jouhaud F. Selecting and scheduling observations of agile satellites[J]. Aerospace Science and Technology,2002,6(5): 367-381

[6]侯小瑾，楊芳，劉勝利，等. 敏捷SAR衛(wèi)星滑動聚束模式下的控制要求分析[C]//2011年小衛(wèi)星技術(shù)交流會. 北京: 航天東方紅衛(wèi)星有限公司,2011

Hou Xiaojin,Yang Fang,Liu Shengli,et al. Analysis of control demand of sliding spotlight mode based on the agile SAR[C]// 2011 Small Satellite Technique Symposium. Beijing: DFH Co．Ltd.,2011 (in Chinese)

[7]邵曉巍,陶久亮,龔德仁,等. 輕型SAR二維姿態(tài)牽引規(guī)劃與高精度追蹤控制[C]//2011年小衛(wèi)星技術(shù)交流會. 北京:航天東方紅衛(wèi)星有限公司,2011

Shao Xiaowei,Tao Jiuliang,Gong Deren,et al. Two dimension attitude planning and high precision tracking control of light SAR[C]// 2011 Small Satellite Technique Symposium. Beijing: DFH Co.,Ltd.,2011 (in Chinese)

[8]張守信. 外彈道測量與衛(wèi)星軌道測量基礎(chǔ)[M]. 北京：國防工業(yè)出版社,1992

Zhang Shouxin.Exterior ballistic measurement and satellite orbit measurement[M]. Beijing: National Defence Industry Press,1992 (in Chinese)

[9]魏靜波. 視頻小衛(wèi)星姿態(tài)控制技術(shù)研究[D]. 長沙：國防科學(xué)技術(shù)大學(xué)研究生院,2011

Wei Jingbo. Attitude control technology of a small video satellite [D]. Changsha: Graduate School of National University of Defense Technology,2011 (in Chinese)

[10] 章仁為．衛(wèi)星軌道姿態(tài)動力學(xué)與控制 [M]. 北京：北京航空航天大學(xué)出版社，1998

Zhang Renwei. Satellite orbit attitude dynamics and control[M].Beijing：Beihang University Press,1998

[11]王睿,楊汝良. 橢圓軌道下雷達多普勒特性估計[J]. 電子與信息學(xué)報,2004,26(1): 107-111

Wang Rui,Yang Ruliang. Radar Dopplerpro perties evaluation in elliptic orbits[J]. Journal of Electronics&Information Technology,2004,26(1): 107-111 (in Chinese)

[12]孫志遠,張劉,金光,等. 視頻小衛(wèi)星凝視姿態(tài)跟蹤的仿真與實驗[J]. 光學(xué)精密工程,2011,19(11): 2715-2723

Sun Zhiyuan,Zhang Liu,Jin Guang,et al. Simulation and experiment on attitude tracking control of small TV satellite[J]. Optics and Precision Engineering,2011,19(11): 2715-2723 (in Chinese)

(編輯：張小琳)

Study on Attitude Maneuver Strategy of Spotlight Mode of Agile SAR Satellite

HAN Xiaolei1ZHANG Qingjun1LIU Jie1ZHANG Running2YUAN Zhi1

(1 Beijing Institute of Spacecraft System Engineering,Beijing 100094,China) (2 DFH Satellite Co., Ltd.，Beijing 100094,China)

In this paper,a new algorithm computing attitude maneuver parameters of spotlight mode is presented for agile SAR satellite. It adequately considers the effects of the position of antenna phase center and the angle departing from coordinate axis. Through vector computation,it acquires initial attitude angles required by spotlight mode of agile SAR satellite. Through iterative computation the errors caused by the position of antenna phase center are removed,and a precise attitude angle is achieved. Finally,the combining simulation results by Matlab and STK confirm the validity of the proposed algorithm.

synthetic aperture radar; agile SAR satellite; spotlight mode; attitude maneuver

2016-03-09；

2016-07-12

國家重大科技專項工程

韓曉磊，男，博士，研究方向為合成孔徑雷達衛(wèi)星工作體制與信號處理。Email：hanxiaolei23@163.com。

TN958

A

10.3969/j.issn.1673-8748.2016.04.003