馮曉萌,吳玲達,2,于榮歡,楊 超
(1. 裝備學院 復雜電子系統(tǒng)仿真實驗室, 北京 101416;2. 國防科技大學 信息系統(tǒng)工程重點實驗室, 湖南 長沙 410073)
?
單硬件限制下的電磁環(huán)境加速繪制*
馮曉萌1,吳玲達1,2,于榮歡1,楊超1
(1. 裝備學院 復雜電子系統(tǒng)仿真實驗室, 北京101416;2. 國防科技大學 信息系統(tǒng)工程重點實驗室, 湖南 長沙410073)
單硬件實現的高效電磁環(huán)境繪制適用范圍更廣;但是,并行光線投射繪制電磁環(huán)境時,其效率受硬件性能制約。在研究硬件限制并行光線投射效率的基礎上,提出一種面向硬件制約的像素插值方法。當硬件限制并行光線投射繪制不能實時完成時,減少并行投射的光線數量,即部分圖像像素由光線投射生成,其余像素插值生成。像素插值以圖像質量換取執(zhí)行效率,當圖像更新停頓時重新使用光線投射生成插值獲得的像素,以恢復圖像內容。實驗結果表明,低硬件配置條件下,像素插值能夠大幅度提高繪制圖像的生成效率。同時,對比多個體數據的繪制效果和誤差統(tǒng)計得出:電磁環(huán)境數據場最適合使用像素插值方法。
電磁環(huán)境;光線投射;硬件限制;像素插值;統(tǒng)一計算設備架構
電磁環(huán)境已經成為現代信息化戰(zhàn)場中決定勝負的關鍵因素,對其進行可視化是虛擬戰(zhàn)場環(huán)境的重要研究內容。在現有技術中,效果最好的電磁環(huán)境可視化方法是體繪制方法[1]。
體繪制方法將三維數據場直接映射成二維圖像[2],以展示三維數據場中的相關信息,方便用戶觀察和理解數據場。體繪制方法生成的繪制結果展示了復雜電磁環(huán)境的內部細節(jié)[3],為用戶觀察和理解電磁環(huán)境態(tài)勢提供了支持。光線投射算法[4]是繪制效果最好的體繪制算法。但是,其計算量大,達到實時繪制通常需要基于圖形硬件GPU的并行實現[5]。統(tǒng)一計算設備架構(Compute Unified Device Architecture,CUDA)[6]出現后,GPU進行通用計算更加方便,使用其加速體繪制的算法[7]越來越普遍。
文獻[5]總結了基于GPU的體繪制研究現狀,體繪制算法在GPU中并行實現后加速效果明顯,能夠交互繪制大規(guī)模體數據。文獻[3]中使用GPU對電磁環(huán)境體繪制進行加速,提高了繪制效率。文獻[8]基于CUDA架構加速繪制動態(tài)電磁環(huán)境,繪制效率相比文獻[3]有所提升。
但是,CUDA加速依賴于硬件的自身性能,當單個硬件限制體繪制執(zhí)行效率時,使用多個硬件可以提高效率[9]。然而,單硬件環(huán)境比多硬件環(huán)境更容易獲得,使用單硬件完成電磁環(huán)境繪制的應用范圍將更為廣泛,特別是針對指揮員的使用,低硬件配置更容易滿足需求。因此,研究單硬件條件下的電磁環(huán)境加速繪制問題具有很強的實用意義。
光線投射算法中光線相互獨立,易于并行實現,從而提高其執(zhí)行效率。但是,同一并行光線投射算法在不同型號的顯卡上執(zhí)行時間不同,如圖1所示。
圖1 并行光線投射效率Fig.1 Efficiency of parallel ray-casting
圖1中,不僅硬件不同對應的執(zhí)行時間不同,當圖像分辨率超過一定規(guī)模后,兩個硬件的執(zhí)行時間均與分辨率成正比。而并行光線投射算法中,每個像素對應一條光線,說明并行投射光線越多,執(zhí)行時間越長。因此,在硬件保持不變且無法修改光線投射算法時,若需要提高執(zhí)行效率,可適當減少投射光線的數量。
減少投射光線數量即減少光線投射生成的像素,此時能夠提高算法執(zhí)行效率,但是應該同時保持圖像分辨率才有意義。因此,使用光線投射算法生成圖像中的部分像素,其余像素使用比光線投射計算量小的方法生成。
2.1像素插值可行性
光線投射算法是將三維數據場映射為二維圖像效果最好的方法,因此,在提高算法執(zhí)行效率的同時應該保證光線投射對圖像生成的作用。圖像中部分像素使用光線投射算法生成,然后利用這些像素插值生成其余的像素,通過插值間接獲得光線投射的繪制效果。
在光線投射算法中,體數據為三維離散數據場,光線穿越數據場時,對數據的采樣點并非全部命中數據點,通常需要使用采樣點附近的多個數據點進行插值獲取采樣點值。說明,插值過程在光線投射算法中已有使用,插值得到的數據值轉換為像素值的過程由傳遞函數[10]決定,可以看作是相鄰的可插值數據經過傳遞函數變成了相鄰的像素。因此,相鄰像素間插值避開了傳遞函數,節(jié)省光線投射操作時間的同時引入了插值與傳遞函數不同所致的誤差。
電磁環(huán)境的可視化技術中[3,11],其數據場由連續(xù)的電磁信號值離散得到,數據點與數據點之間相關且無突變,較適合進行插值。因此,電磁環(huán)境數據場相對一般體數據而言,其數據場內部較適宜插值,從而經過傳遞函數后得到的像素也會相對比較適宜插值。這一推論將在后面的實驗結果中得到驗證。
2.2保持圖像分辨率
為便于后續(xù)討論,定義下述4個用語。
1)源像素:光線投射生成的像素。
2)插值像素:通過插值生成的像素。
3)行插值:某像素利用其同一行左右相鄰的兩個像素進行插值生成。
4)列插值:某像素利用其同一列上下相鄰的兩個像素進行插值生成。
對某像素,插值生成的值與光線投射生成的值之間存在差異。因此,像素插值影響生成圖像的質量。為盡量保證圖像質量,插值像素數量應該盡量少,以減少像素值的誤差。
圖2 源像素與插值像素分布Fig.2 Distribution of source and interpolation pixels
圖2為源像素與插值像素數量比例為1 ∶1時在生成圖像中的分布示意圖,其中顯示了兩種分布情況:圖2(a)至圖2(d)為第一種;圖2(e)和圖2(f)為第二種。這兩種分布中,插值像素與源像素的相對位置有規(guī)律性,利于并行程序的編寫。
觀察圖2中的6幅子圖:當插值像素位于圖像內部時可以統(tǒng)一進行行插值;當插值像素位于圖像邊緣時,需要特殊考慮。據此可知,圖2的第一種分布情況中,需要特殊考慮左右兩個邊緣的像素,并區(qū)分奇偶行;第二種分布情況則只需特殊考慮右側邊緣的像素,且不區(qū)分奇偶行。因此,圖2中的第二種分布情況插值類型少,更易于實現并行編程。
2.3提高插值像素比例
圖2中,源像素與插值像素數量比為1 ∶1,由圖1可知,減少一半的投射光線數量可以節(jié)省不到一半的執(zhí)行時間。當需要節(jié)省更多的執(zhí)行時間時,可以繼續(xù)減少投射光線的數量,即提高插值像素的比例。在圖2的基礎上繼續(xù)減少一半的源像素后,像素分布如圖3所示,源像素與插值像素的比例約為1 ∶3。此時,相對圖2所示的分布,在執(zhí)行過程中光線投射時間減少、插值操作時間增加,總時間將會減少,從而進一步提高生成圖像的效率。
圖3 源像素與插值像素近似比例1 ∶3Fig.3 Ratio of source pixel to interpolation pixel is 1 ∶3
圖3中的插值像素可以分兩步生成,首先生成圖中源像素下方相鄰的插值像素,將此時插值生成的像素看作源像素即與圖2中的第二種分布相同,然后再利用圖2中第二種分布的插值方法生成其余像素。在上述兩個插值步驟中,除邊緣像素需特殊考慮外,其他插值像素在第一步中進行列插值,在第二步中進行行插值,操作統(tǒng)一利于并行處理。
繪制電磁環(huán)境數據場時,圖像邊緣通常對應數據場的邊緣,而數據場邊緣的數據信息量很小,可以不精確顯示。因此,繪制電磁環(huán)境數據場時,將位于圖像邊緣的插值像素設置為與其相鄰的源像素值,以盡量節(jié)省執(zhí)行時間。
圖3所示的分布中,源像素所占的比例已經比較小,若再降低,將使誤差繼續(xù)擴大,嚴重影響生成圖像的質量。因此,非特殊需要將不再繼續(xù)提高插值像素比例。
像素插值在犧牲圖像質量的前提下提高圖像生成速度,以支持圖像的高更新頻率。但是,當圖像更新停頓時,應使用光線投射重新生成插值像素,以恢復圖像內容。
3.1插值像素的恢復流程
從像素插值到恢復插值像素的流程如圖4所示。像素插值是一個循環(huán)過程,保證實時更新圖像,通過判斷“是否空閑”決定繼續(xù)循環(huán)還是跳出循環(huán)。跳出循環(huán)后,對插值生成的像素使用光線投射算法進行恢復。
圖4 像素插值及恢復流程圖Fig.4 Flow of pixel interpolation and recovery
圖4中,上下兩個虛線框分別表示設備端(GPU)和主機端(CPU)中的操作,由主機端確定“是否空閑”,“是”表示暫時不需要更新圖像,“否”表示需要立即更新圖像。當用戶通過主機端對體繪制進行交互控制或者電磁環(huán)境數據場發(fā)生變化時,均需要重新生成并顯示圖像,此時采用像素插值的方法提高圖像生成效率??臻e時,設備端恢復圖像并顯示。此時,若主機端需要繼續(xù)更新圖像,則從“開始”操作重新執(zhí)行設備端的處理流程。圖4中表示GPU操作的虛線框中明確標識了“開始”未標識“結束”,這是因為是否結束由CPU控制,可以在任何操作時結束。
3.2是否空閑的判斷
圖4中,是否空閑的含義為是否需要更新圖像,由主機端決定,其因素主要有用戶交互、數據場變化等。用戶交互調整傳遞函數或者觀察視角時,圖像需要根據用戶的調整及時更新,以反饋交互效果。當硬件限制而不能實時更新圖像時,可采用像素插值的方法提高圖像生成效率。然而,用戶的交互操作會有停頓即“空閑”,此時進行像素恢復。同理,數據場變化有停頓時也進行像素恢復。使用時間閾值ε判斷“空閑”狀態(tài):在ε時間內既沒有用戶交互操作,數據場又未變化時,認定此時為“空閑”。
“空閑”狀態(tài)下恢復圖像是為了消除像素插值造成的生成圖像誤差,恢復圖像質量,便于用戶更好地觀察電磁環(huán)境的繪制結果。
依據圖4所示的流程進行編碼實現,程序主要分為兩個部分:主機端代碼和設備端代碼。主機端除負責將繪制電磁環(huán)境需要的數據和控制參數傳遞到設備端外,最主要的任務就是實時監(jiān)控“空閑”狀態(tài),并將結果反饋給設備端。設備端使用CUDA架構并行實現,光線投射、像素插值、像素替換等操作均并行實現。
算法1中顯示了主機端控制像素插值流程的關鍵代碼,其中插值像素的分布如圖3所示。其中,第3行獲取圖像分辨率;第4和第5行取圖像分辨率的1/4;第6至8行設置開啟的線程數,每個源像素對應1個線程;第9行進行光線投射生成源像素;第10行進行列插值;第11至13行更新開啟的線程數,為圖3中的行插值做準備;第14行進行行插值。
算法1 像素插值流程關鍵代碼
由于圖3中行插值像素個數大約為列插值像素個數的2倍,因此在11至13行對線程個數進行了重置,表1中行插值函數g_InsertPix2()比列插值函數g_InsertPix1()開啟的線程數多1倍。光線投射函數g_CastRay()的具體執(zhí)行過程與文獻[8]中相同。函數g_InsertPix1()的偽代碼見算法2。
算法2 像素插值內核線程
算法2中,前四行根據線程的ID號確定其對應插值像素的坐標,當插值像素位于圖像邊緣時直接將相鄰源像素的值賦予它(第7行),否則進行線性列插值(第9行)。函數g_InsertPix2()的執(zhí)行代碼與算法2中相似,插值時使用行插值。
基于CUDA架構實現本算法,編程開發(fā)環(huán)境為集成了NVIDIA GPU Computing SDK 4.2的Microsoft Visual Studio 2005。硬件支持環(huán)境為:Intel(R) Core(TM) i5-2400 CPU 3.10 GHz、4 GB內存、NVIDIA GeForce GT 530顯卡。由圖1可知,此顯卡并行計算能力相對較低,對并行光線投射算法的硬件限制較大,當其并行能力不能支持實時光線投射時,可以使用像素插值的方法進行加速。
5.1執(zhí)行效率
使用體繪制中的DVR(直接體繪制)算法和MIDA算法[12]進行了兩組實驗,生成圖像分辨率為1200×900,連續(xù)生成100幀圖像的時間統(tǒng)計如圖5所示。每組實驗分別統(tǒng)計了三種情況下的生成時間:①全部像素使用光線投射生成;②源像素與插值像素數量比例為1 ∶1;③源像素與插值像素數量比例為1 ∶3。圖5的兩幅圖中,上中下三條曲線分別對應第1、第2、第3種情況。
(a) DVR生成圖像時間(a) Implementation time of DVR
(b) MIDA生成圖像時間(b) Implementation time of MIDA圖5 算法執(zhí)行時間Fig.5 Time for algorithm implementation
通過圖5兩幅子圖中三條曲線的具體時間可以看出,像素插值能夠將圖像的生成時間從40 ms以上降到40 ms以下,即達到實時生成圖像。說明,在由于硬件能力限制而不能實時生成圖像的情況下,像素插值能夠加速圖像生成,使其能夠實時完成。同時,將各子圖中的三條曲線進行對比可以得出如下結論:生成圖像的分辨率不變時,插值像素越多圖像的生成時間越短。因此,適當提高插值像素在圖像中所占的比例能夠獲得更好的加速效果。
圖5(a)中三條曲線的平均每幀時間為51.84 ms,29.68 ms,18.96 ms,中下兩條曲線的平均時間分別是最上面曲線的57.25%,36.57%;上述數據值對應到圖5(b)中分別是59.06 ms,33.79 ms,21.08 ms,57.21%,35.69%。對比可知,雖然MIDA算法比DVR算法用時要長,但是像素插值對兩個算法的加速效率相差不多,甚至MIDA算法的加速效率高些。說明,像素插值對光線投射生成圖像的加速效率比較穩(wěn)定,可適用算法范圍較廣。
5.2繪制效果與誤差統(tǒng)計
記電磁環(huán)境數據場為EME,使用文獻[8]中的方法計算得到。實驗中,使用4個體數據場(EME:200×200×100,Sphere:64×64×64,Engine:256×256×110,stagbeetle:277×277×164)生成分辨率為1200×900的圖像時,4個體數據對應的圖像中人眼均無法辨別出使用了像素插值的圖像。而將生成圖像分辨率調整為與體數據場規(guī)模的最大兩維相等后,生成圖像中部分區(qū)域如圖6所示。圖6中,每行對應一個體數據,從上至下依次為EME,Sphere,Engine,stagbeetle;3列從左到右依次為全部光線投射生成、圖2中第二種分布的像素插值生成、圖3所示分布的像素插值生成。
圖6 繪制效果對比Fig.6 Contrast of rendering results
仔細觀察對比圖6中每個體數據的不同繪制結果可知,EME和Sphere的繪制結果幾乎看不出差別,而Engine和stagbeetle的繪制結果則插值像素越多越模糊、效果越差。而將生成圖像分辨率改為1200×900后,Engine和stagbeetle的繪制結果也幾乎看不出差別,說明小規(guī)模的數據場生成大分辨率的圖像時,插值像素對繪制結果的影響減小。而當數據場規(guī)模與圖像分辨率相當時,EME和Sphere這種計算獲得的數據場要比CT掃描獲得的數據場(Engine和stagbeetle)更適合使用像素插值方法進行加速繪制。
從數據的角度對像素插值生成的圖像中的誤差進行統(tǒng)計分析:保持數據場規(guī)模不變,生成圖像分辨率為1200×900時的統(tǒng)計結果見表1;生成圖像分辨率與體數據場規(guī)模的最大兩個維度相等時的統(tǒng)計結果見表2。其中,所有數據均使用相同的傳遞函數繪制生成。表1和表2中數據均為10幀圖像數據的統(tǒng)計均值,“最大誤差”列中數據只保留了整數。兩個表的數據中“B”通道統(tǒng)計數值均為0是因為傳遞函數設置中未使用藍色通道。
表1 插值像素誤差統(tǒng)計1
表2 插值像素誤差統(tǒng)計2
對于某個像素點,稱其由插值生成的值與光線投射生成的值之間的差為誤差,最大誤差指圖像中誤差最大值。平均誤差計算方法為所有誤差的和除以存在誤差的像素總數。誤差像素比是存在誤差的像素總數除以圖像中的非零(4個通道不全為零)像素總數。每個體數據均對應兩行數據,上面一行為插值像素與源像素比例為1 ∶1的情況,下面一行為3 ∶1的情況。
分別對比表1和表2中每個數據對應的兩行數據可知,插值像素個數較少時,誤差情況相對較小。對比所有數據的“最大誤差”列,前兩個體數據明顯優(yōu)于后兩個,而EME和Sphere均為數值計算得到的數據,后兩者則是實物掃描數據,說明數值計算得到的數據場較適于使用像素插值的方法進行加速。“平均誤差”與“誤差像素比”兩列的數據同樣是EME的最好,說明電磁環(huán)境數據場更適合使用像素插值的方法進行繪制時的加速。
對比兩個表中數據,印證了圖像分辨率對像素插值結果的影響:數據場規(guī)模一定時,生成的圖像分辨率越大,像素插值的誤差越小。同時,表2與圖6相對應,其中EME對應的數據說明當使用更大規(guī)模電磁環(huán)境數據場生成高分辨率的圖像時,仍然可以使用像素插值的方法加速體繪制。而其他掃描得到的實體數據則將不適用像素插值方法。
5.3像素恢復
表1、表2顯示繪制電磁環(huán)境時使用像素插值方法存在誤差,為保證圖像質量恢復插值像素是有必要的。實驗中,根據插值像素分布設置空閑判斷閾值ε:源像素與插值像素為1 ∶1時設置為使用像素插值時生成第一幀圖像的執(zhí)行時間;1 ∶3 時設置為生成第一幀圖像執(zhí)行時間的1.5倍。如此,基本保證閾值ε大約為不進行像素插值的一半,當已空閑1個ε時使用ε或1.33ε的時間進行像素恢復。如此設置閾值ε在實驗中未影響用戶的交互操作,是可行的。
上述實驗和分析說明,像素插值方法適用于繪制電磁環(huán)境數據場,不僅繪制誤差小,同時還能夠在繪制效率受硬件性能制約時進行加速繪制,從而得到更高的繪制效率。因此,本算法能夠應用于文獻[3]和文獻[8]中,進一步加速其中的電磁環(huán)境繪制效率。同時,若接受一定的繪制效果損失,可以將本算法應用于繪制其他體數據,繼續(xù)加速文獻[5]中提到的基于GPU的體繪制算法。
提出一種適用于電磁環(huán)境數據場繪制的像素插值算法,其可以在硬件性能限制繪制效率時對繪制進行加速,獲得更高的繪制效率為其實時交互和應用提供支持。實驗結果表明:像素插值算法的加速效果明顯,且生成圖像的誤差明顯優(yōu)于其他體數據。算法的不足是,犧牲生成圖像的質量換取繪制效率的提升,對其應用具有一定的限制。提高插值像素的精度可以改善這一限制,是下一步的研究工作。
References)
[1]吳迎年, 張霖, 張利芳, 等. 電磁環(huán)境仿真與可視化研究綜述[J]. 系統(tǒng)仿真學報, 2009, 21(20): 6332-6338.WU Yingnian, ZHANG Lin, ZHANG Lifang, et al. Survey on electromagnetic environment simulation and visualization[J]. Journal of System Simulation, 2009, 21(20): 6332-6338. (in Chinese)
[2]Arie K, Klaus M. Overview of volume rendering[M]//Hansen C D, Johnson C R. The Visualization Handbook, USA:Academic Press, 2005.
[3]楊超, 徐江斌, 吳玲達. 硬件加速的虛擬電磁環(huán)境體可視化[J]. 北京郵電大學學報, 2011, 34(1): 55-59.
YANG Chao, XU Jiangbin, WU Lingda. Hardware accelerated volume visualization in virtual electromagnetic environment[J]. Journal of Beijing University of Posts and Telecommunications, 2011, 34(1): 55-59. (in Chinese)
[4]馬千里, 李思昆, 白曉征, 等. CFD非結構化網格格心格式數據高質量體繪制方法[J]. 計算機學報, 2011, 34(3): 508-516.
MA Qianli, LI Sikun, BAI Xiaozheng, et al. High-quality volume rendering of unstructured-grid cell-centered data in CFD[J]. Chinese Journal of Computers, 2011, 34(3): 508-516. (in Chinese)
[5]Beyer J, Hadwiger M, Pfister H. A survey of GPU-based large-scale volume visualization[C]// Proceedings of the Eurographics Conference on Visualization 2014, Swansea Wales, UK, 2014: 105-123.
[6]Rosen P. A visual approach to investigating shared and global memory behavior of CUDA kernels[C]// Proceedings of the Eurographics Conference on Visualization, Leipzig, Germany, 2013: 161-170.
[7]Zhang Y B, Dong Z, Ma K L. Real-time volume rendering in dynamic lighting environments using precomputed photon mapping[J]. IEEE Transactions on Visualization and Computer Graphics, 2013, 19(8): 1317-1330.
[8]馮曉萌, 吳玲達, 董士偉. CUDA加速的動態(tài)電磁環(huán)境數據場實時繪制[J]. 系統(tǒng)仿真學報, 2014, 26(9): 2044-2049.
FENG Xiaomeng, WU Lingda, DONG Shiwei. CUDA accelerated real-time rendering for dynamic electromagnetic environment volume data [J]. Journal of System Simulation, 2014, 26(9): 2044-2049. (in Chinese)
[9]Stuart J A, Chen C K, Ma K L, et al. Multi-GPU volume rendering using MapReduce[C]//Proceedings of the 19th ACM International Symposium on High Performance Distributed Computing, 2010: 841-848.
[10]Arens S, Domik G. A survey of transfer functions suitable for volume rendering[C]//Proceedings of the 8th IEEE/EG International Conference on Volume Graphics, 2010: 77-83.
[11]陳鵬, 楊超, 吳玲達. 硬件加速的三維雷達作用范圍表現[J]. 國防科技大學學報, 2007, 29(6): 49-53.
CHEN Peng, YANG Chao, WU Lingda. Hardware accelerated 3D radar detection range visualization[J]. Journal of National University of Defense Technology, 2007, 29(6): 49-53. (in Chinese)
[12]Stefan B, Eduard M G. Instant volume visualization using maximum intensity difference accumulation [J]. Computer Graphics Forum, 2009, 28(3): 775-782.
Accelerated rendering for electromagnetic environment under single device restriction
FENG Xiaomeng1, WU Lingda1,2, YU Ronghuan1, YANG Chao1
(1. Science and Technology on Complex Electronic System Simulation Laboratory, Equipment Academy, Beijing 101416, China;2. Key Laboratory of Information System Engineering, National University of Defense Technology, Changsha 410073, China)
Electromagnetic environment with a high efficiency based on single device supports has wide range of applications. But the efficiency of parallel ray-casting for rendering electromagnetic environment is restricted by the device. Based on researching the restriction of device for parallel ray-casting, a pixel interpolation method focusing on the restriction was presented. The number of rays was reduced when the parallel ray-casting rendering under device restriction couldn’t be completed immediately, namely, a part of pixels in the rendering image were generated by ray-casting and the rest pixels through interpolation. Pixel interpolation got rendering efficiency at the cost of image quality, so when image update paused, the interpolated pixels were regenerated to recover image quality. The experiments show that pixel interpolation obviously improves rendering efficiency when implemented on a low device. Compared with the rendering images of some volume data and the errors in these images, the electromagnetic environment data has the best rendering result, which proves that pixels interpolation is useable especially for rendering electromagnetic environment on a low device.
electromagnetic environment; ray-casting; device restriction; pixel interpolation; compute unified device architecture
10.11887/j.cn.201604011http://journal.nudt.edu.cn
2015-04-28
國家自然科學基金資助項目(61202129)
馮曉萌(1986—),男,河北石家莊人,博士研究生,E-mail:130123feng@163.com;吳玲達(通信作者),女,研究員,博士,博士生導師,E-mail:wld@nudt.edu.cn
TP391.9
A
1001-2486(2016)04-069-07