高 勇,劉 家 駿,郭 瀟,鄔 倫
(北京大學遙感與地理信息系統研究所,北京 100871)
面向大規模動態地形可視化的LOD組織與調度技術
高 勇,劉 家 駿,郭 瀟,鄔 倫
(北京大學遙感與地理信息系統研究所,北京 100871)
地貌演化等地表過程的分析模擬實時產生系列地形數據,該動態過程的可視化對性能要求較高,傳統基于預處理的方法計算量大、繪制延遲明顯。面向大規模動態地形可視化,提出一種LOD預處理與實時更新相結合的混合調度技術,基于場景圖對動態地形LOD建模與組織,建立多線程模型實現地形瓦片的動態更新和繪制。對于每次迭代計算產生的新地形,該方法以LOD瓦片的動態調度替代傳統的靜態預處理方法,與地表過程計算同步更新瓦片,且僅重新計算和渲染當前視域和分辨率下地形發生變化的局部瓦片,通過LOD的局部動態更新避免可視化數據的全局重新生成,顯著降低過程模型計算和結果更新對渲染過程的影響,提高整體的可視化性能。
地形可視化;細節層次模型;場景圖;多線程;動態調度
地形可視化是三維地理信息系統的基礎,是地形分析應用中重要的技術環節,而在諸如地貌演化、水土流失等地表過程分析和模擬中,更是對地形可視化產生了新的需求。這類基于過程的地形數值分析和模擬,表現為地形因子,特別是高程h隨時間t的動態變化過程[1],即h=f(t)或Δh=f(Δt)。將該過程通過三維可視化的方法實時動態展現出來,可以更直觀有效地探查、分析和檢驗地形演變過程的規律和特征,特別是在大規模地形分析中尤為如此。
當前大規模地形的三維可視化技術,主要以細節層次模型(Level of Detail,LOD)為基礎[2,3],以到觀察點的距離評價地形塊的重要程度,距離觀察點較遠則以粗糙模型顯示,反之則顯示精細模型。在實現中,首先通過LOD算法從原始的地形數據生成地形模型,然后在可視化中調用生成的模型完成顯示。其中比較有代表性的算法有ROAM算法[4]、限制性四叉樹算法[5]、Chunked LOD算法[6]等。從原始地形到地形模型,需要一個處理過程,這個過程不能在一幀的時間內完成。因此,大規模地形可視化所用的地形模型,通常通過預處理生成。
從地形可視化到面向過程的動態地形可視化,其技術挑戰主要來自于LOD生成。預處理生成的LOD模型是靜態的,且生成效率通常較低,盡管目前已經提出了一些LOD動態調度[6,7]、基于GPU的并行模型生成[8-11]等加速技術,但這些方法仍然是先整體預處理再渲染的兩階段過程。而地貌演化等地表過程是高程隨時間變化的函數,其最大的特點是隨時間演替不斷生成新的地形數據,即本文所稱的動態地形。當前的可視化方法對每次產生的新地形數據都要重新進行預處理并重新構建LOD,然后交給渲染過程,才能顯示新地形的變化。但由于在大數據量下預處理耗時長,很難在地表過程的一個時間步長內完成,這就必然會導致地形實時可視化的停頓、更新不及時等性能問題,難以滿足動態地形流暢漫游的需求。因此,對于動態地形,需要結合可視化過程實時生成LOD模型,當前已提出了一些局部地形的實時生成和更新方法,用以支持交互式地形編輯和渲染[1],如結合體素的混合地形表達[12]、基于時空立方體的地形變形可視化[13]等,但這些方法主要面向局部地形更新,還不能很好處理諸如地貌過程模擬等大規模全局動態地形的可視化。
因此,本文提出一種LOD預處理與實時生成更新相結合的混合調度技術,提高動態地形過程的可視化性能。先通過預處理生成基礎的地形模型數據,保證可視化時存在可調用的地形模型。在地形過程計算中,對于新產生的地形數據,根據視點特征,動態查找需要更新的LOD瓦片,實時生成局部的新的地形模型,并完成該局部地形的動態繪制。該方法避免了全局地形模型的重新生成,而只進行少量的局部的必要更新,可以較好地解決大規模動態地形可視化的性能問題。
結合預處理與實時生成更新方法進行LOD的混合調度,在完成地表過程動態計算的同時,實現動態地形的實時可視化和流暢漫游。具體的,先通過傳統的預處理過程生成基礎的地形瓦片,之后的可視化過程與地表過程計算實時結合,地表過程模型在原始地形數據的基礎上不斷迭代產生新的地形數據,而在同步的地形可視化漫游中,首先采用與靜態模型相同的方法加載預處理生成的地形瓦片,同時動態判斷在當前視域和分辨率下新動態地形導致變化的瓦片,系統根據變化后的新地形,實時的為需要更新的地形生成新的瓦片,并替代之前的瓦片進行渲染。在該過程中,僅重新計算和渲染當前視域和LOD層次上受新地形影響的局部瓦片,避免了全局的、所有層次上的瓦片更新,因而可以大幅度降低因數據準備而導致的渲染延遲。
LOD技術控制了同屏渲染的三角形的數量,也使得屏幕中同時出現的地形瓦片不會特別多,即使在大范圍漫游的情況下,實際更新的瓦片數量也不會很大。因此,新地形模型的生成并不會消耗太多時間,瓦片更新的延遲也不會很大,而且沒有變化的部分仍以基礎瓦片顯示而無需計算,因此完全可以滿足動態地形的可視化需求。
該LOD混合調度技術的實現分為3個步驟:1)動態地形的計算。讀取當前地形,應用過程模型Δh=f(Δt),迭代計算經過一定時間步長后的新地形。這部分工作比較常規,單純使用原始的地形數據即可完成,不涉及可視化部分。2)瓦片的生成以及可視化前的準備工作。經過過程計算后的新地形,要先更新LOD瓦片才能實施渲染,這可以看作可視化前的準備。準備工作包括但不限于:動態判定并實時生成受新地形影響而需更新的瓦片,從原始數據提取未受影響的瓦片、提前對顯示列表進行編譯、對紋理作提前編譯等。3)可視化。將準備好的數據加入渲染流程。
要實現流暢的可視化效果,一幀的渲染需要在至少1/24 s內完成,而上述3個步驟不可能在一幀之內完成。因此本文使用多線程技術,將LOD實時生成的第1步和第2步在其他線程完成,只有第3步才在渲染線程中完成。具體的,本文的調度技術涉及3個線程:計算線程、渲染準備線程(簡稱“準備線程”)、渲染線程;3組數據:地形數據、渲染準備任務列表(簡稱“任務列表”)、渲染準備結果列表(簡稱“結果列表”)。3個線程與3組數據之間的基本調度關系如圖1所示。

圖1 動態地形可視化的多線程模型
Fig.1 Multithread model for dynamic terrain visualization
每個線程的具體工作如下:1)計算線程完成第1步動態地形的迭代計算,并負責寫入新產生的地形數據。2)渲染線程完成第3步可視化的工作,并負責寫入任務列表和讀取結果列表。渲染線程在每次渲染中檢查是否存在新的地形,當計算線程完成一次迭代計算并生成新的地形后,渲染線程會檢測到地形發生了變化,進而根據當前視域和分辨率判斷需要更新的瓦片,并生成渲染準備任務,并放入任務列表。渲染線程還將查看提交的任務是否完成,如果任務被處理完,在結果列表就有了新的地形瓦片,渲染線程將用新的瓦片替代原有的瓦片完成渲染工作。這兩項工作通過揀選回調和更新回調完成。3)準備線程完成第2步的瓦片生成和準備工作,負責讀取地形數據以及任務列表,實時生成需要更新的局部瓦片,并將結果寫入結果列表。準備線程會不斷地查詢任務列表,當任務列表中有待完成的任務時,就取出任務,并將任務生成的新瓦片放入結果列表。
該動態地形可視化的LOD混合調度技術的具體實現細節將在下文詳述。
2.1 動態地形LOD的場景圖模型
在動態地形可視化的混合調度技術中,地形LOD基于場景圖進行組織管理。場景圖是用于組織場景信息的圖或樹結構。一棵場景樹包含一個根節點、多級內部的組節點、多個末端的葉子節點,其中葉子節點記錄幾何體信息,主要是地形瓦片數據。通過修改組節點上的變換函數,可直接作用于子節點,而無需重復更改子節點的屬性,因此可以極大簡化修改操作,有利于大規模場景的變換和管理。
為了支持動態地形,將其LOD場景圖組織為四叉樹結構,如圖2所示(以3層為例)。其中LOD節點用于組織外存數據,如果其下一層模型的子節點多于1個(最多4個),則用一個組節點組織在一起。與傳統場景圖最大的不同是,該模型擴展了自定義的動態瓦片節點代替原有的幾何體節點,負責管理動態地形的LOD瓦片的信息,包括該節點所對應的地形瓦片的地理范圍、分辨率等參數,以及上次可視化時的幀數、上次可視化時所用數據的地形迭代計算計數、是否需要更新等過程狀態,還需記錄地形數據和準備數據的指針等。這部分新增的節點信息用于可視化過程中動態判斷該節點瓦片是否需要更新。同時,在動態瓦片節點上添加更新回調和揀選回調,負責響應線程模型的調度,具體完成對應節點在必要條件下的瓦片更新工作。瓦片更新的判斷和重新生成工作都在單個節點上完成,可以有效保證其局部性,并且基于場景圖的結構,可以實時的適應并匹配當前的視域和分辨率。

圖2 動態地形LOD的場景圖組織
Fig.2 Scene graph model for dynamic terrain LOD
2.2 更新瓦片的選擇
為了實現局部地形模型的實時生成和更新,需要根據視點特征,動態查找需要可視化的地形瓦片,僅僅更新該瓦片的地形模型數據。
依據場景圖的四叉樹結構,確定瓦片的可視判斷策略。具體的,從第0層地形模型開始判斷,當進入到第l層時,對該層第x行y列的瓦片tl,x,y,判斷視點到該瓦片的距離distancel,x,y。對于給定閾值thresholdl,如果該距離在[thresholdl,∞]內,則該瓦片可視,并停止判斷;否則該瓦片不可見,并進入該瓦片的所有子節點,繼續判斷。但需要注意,可視距離判斷的區間是[thresholdl,∞],而非[thresholdl,thresholdl-1],這是由于tl,x,y與tl-1,x,y的中心并不相同,使用后者進行判斷,會造成某些視點下瓦片顯示不完全。而當進入最底層瓦片時,則不需要判斷,直接顯示。
在瓦片的選擇中,需要決定如何設置距離閾值thresholdl。本文使用屏幕像素作為確定thresholdl的標準,即thresholdl應使得地形模型中的n個頂點投影在屏幕上時,正好對應屏幕上的n個像素,則:
(1)
式中:fovy為視場在垂直方向的夾角,一般為60°;p為屏幕在fovy相應維度的分辨率,一般為列分辨率;rl為l層地形模型的分辨率。
2.3 回調機制的選擇
對于場景圖的管理,除了執行繪制,系統還需要響應用戶操作,進而實現視角的調整、幾何體的更新,并通過揀選優化顯示效率。因此,一幀的顯示,執行的不只是繪制遍歷,而是更新、揀選、繪制3種遍歷。其中,更新遍歷允許修改場景中的幾何體、更改視點,以實現動態場景和場景漫游;揀選遍歷從待繪制的大量幾何體中挑選出當前視域可見的部分,之后的繪制遍歷只需執行這些幾何體的繪制,提高繪制效率;繪制遍歷使用揀選遍歷過程生成的渲染列表,實現幾何體的渲染。對于不同的應用場景,三種遍歷的行為也有所不同,可以通過圖形引擎提供的回調機制,編寫自己的回調函數,自定義更新、揀選、繪制遍歷的行為。
在渲染線程中,由于節點的可視判斷是在更新遍歷中完成,這表明在場景的更新遍歷完成之前還不能確定所有節點是否可視。所以發現需要更新的節點并生成任務寫入任務列表的工作,不能通過更新回調完成。在本文中,這部分工作通過揀選回調實現。而幾何體的更改必須在更新遍歷中完成,否則后續流程很有可能產生錯誤,甚至是線程沖突的。因此在本文中,任務完成后查詢結果列表并替換渲染幾何體的工作,通過更新回調實現。
同時還需注意一個任務組織的問題。由于準備線程會不斷查看任務列表,并取出待完成任務。如果任務的發起是由單個節點完成的,那么當單個節點發出任務后,就會馬上被準備線程取走,造成任務的碎片化,影響效率。所以本文中,任務的發出以幀為單位,將當前幀的所有任務收集后,再統一放入任務列表。具體方法為,動態瓦片葉節點在其裁剪回調中生成任務,并放入臨時任務列表,在下一幀根節點的更新回調中,查看臨時任務列表,并將其中的任務放入任務列表。
回調功能添加在LOD節點或葉子節點中都可以實現。考慮到葉子節點的裁減回調可以更方便明確地獲取需要可視化的節點,本文選擇在葉子節點,即動態地形瓦片節點上添加更新回調和揀選回調。
在該多線程模型的實現中,除渲染所用的數據外,還需要地形數據用于記錄動態變化的地形,需要準備數據用于記錄任務列表和結果列表。
(1)計算線程。計算線程不斷調用地形過程模型進行迭代計算。過程模型在外部定義,計算結果一般為新的DEM數據。每次迭代計算完成后,產生新的地形數據,并令迭代計數增加1。渲染線程可以查看對比自身的迭代計數,判斷是否需要更新。
(2)渲染線程。渲染線程的工作通過揀選回調和更新回調完成。揀選回調的工作是在當前視域和LOD層次的條件下,判斷因新產生的地形數據而需要更新的瓦片,將其生成任務并放入任務列表。具體的,首先記錄當前幀數,然后判斷當前瓦片節點是否需要更新。在更新判斷過程中,首先比較動態地形瓦片節點的范圍與當前動態區域的范圍是否重合,然后基于迭代計數判斷該動態地形瓦片的數據是否是最新的地形。如果既重合又不是最新地形,則該瓦片需要更新,這時將一個更新該瓦片的任務放入準備數據的任務列表,并且設置更新標記為true。揀選回調流程如圖3所示。

圖3 揀選回調流程
Fig.3 Culling callback flow diagram
更新回調是查看結果列表中是否存在已完成的瓦片更新任務的結果,若存在則用新的瓦片替換舊的瓦片。在判斷瓦片是否需要更新的過程中,先查看更新標記,然后比較上次幀數和當前幀數,判斷該節點上一次的渲染是否很接近。如果更新標記為true,同時幀數也在一定范圍內,則需要更新該節點,否則說明該節點已經不在視野內,則不必更新。當節點需要更新時,查看準備數據的結果列表是否存在當前節點的結果,如果存在則更新該動態地形瓦片節點。更新回調流程如圖4所示。

圖4 更新回調流程
Fig.4 Updating callback flow diagram
(3)準備線程。準備線程不斷查看準備數據的任務列表,如果不為空,則取出其中的任務進行處理,根據任務信息實時生成相應節點的新瓦片數據,并在任務處理完成后,將結果放入結果列表。處理任務時需要讀取地形數據,并利用LOD算法生成局部瓦片。在有些情況下,一些節點發出了準備任務,但是在任務完成前就離開了視域,而由此生成的結果則將在結果列表中滯留,不會被取走,因此準備線程還需要負責清除這部分滯留的結果。
4.1 實驗方法
針對本文提出的動態地形可視化的LOD調度技術,開展實驗測試并驗證其性能。實驗硬件環境使用PC機,主板芯片組Intel HM65,中央處理器Intel Core i7-4700MQ,顯示適配器NVIDIA GeForce GTX 765 M (2 GB),內存16 GB (DDR3-1600 SDRAM),硬盤Hitachi 750 GB 7200RPM。軟件環境為操作系統Microsoft Windows 7 Ultimate 64-bit(Service Pack 1),顯示設備驅動版本NVIDIA ForceWare 340.62。圖形開發環境為OpenGL 4.40,部分基礎功能在OpenSceneGraph 3.2.0基礎上擴展實現。OpenSceneGraph是一個基于OpenGL的由C實現的開源圖形引擎,以場景圖為核心組織三維場景(http://www.openscenegraph.org/)。
實驗測試用的原始地形數據是SRTM 90 m中國南方部分區域的格網DEM數據,經緯度范圍是102°~113°E,27°~32°N。將上述地形數據重采樣為173 MB、690 MB、1.5 GB、3.4 GB的4份DEM數據進行測試。
實驗選擇一個典型的流域地貌演化模型進行動態過程計算,該模型可簡化表示為[14]:
(2)式中:h為地形高程,t為時間,S為地形坡度;Q為單位寬度徑流量,基于降水量經過匯流累積計算獲得,為簡化模型復雜度,實驗中將降水量簡化為均勻持續降水;β、m、n、k是不隨高程和時間變化的模型參數,經率定得到β=0.01,m=1.1,n=0.3,k=0.02 m2/a。為增強測試效果,對模型迭代的時間步長進行了夸大。
4.2 實驗結果
使用上述測試數據和過程模型,采用遠近兩種視點、不同大小的動態地形區域進行測試。測試中關閉垂直同步,以便通過幀速率判斷性能。實驗的地貌演化動態效果如圖5所示。

圖5 動態地形可視化效果
Fig.5 Visualization results of dynamic terrain
動態可視化的性能測試結果如表1所示。另外,使用3.4 GB、動態范圍為20 000×10 000的DEM數據,測得一個典型可視化漫游過程的幀率變化如圖6所示。在各測試環境中都達到了大于60幀/s的顯示速率,滿足可視化需求。
表1 實驗性能
Table 1 Efficiency of the technology

數據量分辨率動態范圍視角幀率173MB10085×44702000×2000遠1329173MB10085×44702000×2000近324173MB10085×44705000×2000遠1329173MB10085×44705000×2000近329690MB20170×89402000×2000遠1360690MB20170×89402000×2000近273690MB20170×894010000×5000遠1325690MB20170×894010000×5000近2751.5GB30255×134102000×2000遠12671.5GB30255×134102000×2000近2521.5GB30255×1341015000×7000遠12911.5GB30255×1341015000×7000近2413.4GB45382×201152000×2000遠12083.4GB45382×201152000×2000近2103.4GB45382×2011520000×10000遠12813.4GB45382×2011520000×10000近208

圖6 典型漫游下的幀速率
Fig.6 Frame rate under classical roaming
動態地形的數據調度過程必然會降低渲染效率,因此有必要對本文提出的動態調度技術的性能進行評價。分別在開/關動態地形的情況下,對可視化幀率進行測試對比,其結果如表2所示,其中每行的測試數據和參數與表1對應。在測試環境下,動態過程造成的可視化性能損失均小于10%,證明本文方法具有很好的性能。
對比近年的已有研究,本文方法具有較好的性能。Koca等在2 048×2 048的DEM數據上疊加動態地形要素,其渲染性能平均可達235幀/s[12]。Kang等提出的多分辨率地形渲染方法,在16 000×16 000和32 000×32 000的DEM數據上,平均幀速率達423幀/s[11]。李欽等實現的基于分塊LOD的地形實時渲染算法,對于4 097×4 097的DEM數據,在中遠視角下漫游的平均幀速率達877幀/s[15]。而本文方法在考慮動態地形的情況下,多組測試數據上的渲染性能(表1、圖6)平均可達658幀/s,其中近視角下達225幀/s,中遠視角下達875幀/s,可視化性能均達到或高于已有研究成果。
表2 動態地形可視化的性能損失
Table 2 Efficiency loss of dynamic terrain visualization

動態幀率非動態幀率性能損失動態幀率非動態幀率性能損失132913572.06%126713002.54%3243301.82%2522571.95%132913572.06%126713002.54%3293300.30%2332579.34%136013640.29%120812745.18%2732792.15%1912109.05%132513642.86%126312740.86%2752791.43%2082100.95%
面向土壤侵蝕、地貌演化等地形過程模擬的應用需求,提出了一種預處理與實時更新相結合的LOD混合調度技術,實現大規模動態地形可視化。基于場景圖組織管理大規模地形的LOD模型,支持動態地形瓦片的實時查找、更新和繪制,提出并實現LOD瓦片動態調度和繪制的多線程模型,利用計算線程和準備線程分別負責動態地形的計算和地形瓦片的生成,利用渲染線程發現需更新的地形瓦片并及時更新。該技術以LOD瓦片的動態更新繪制替代靜態預處理方法,在地表過程模型迭代計算的同時,同步更新地形瓦片,并且僅更新在當前視域和分辨率下發生變化的局部瓦片,避免了全局更新的大量計算,使得瓦片更新過程不會對渲染速度產生明顯的影響,從而大幅度提高可視化性能。對于GB級以上的DEM數據,動態地形可視化幀速率達60幀/s以上,可有效滿足地形過程模擬可視化的性能需求。在本文的技術實現中還未加入對紋理的支持,另外在大視野場景下的漫游,會出現輕微的局部新地形的視覺抖動現象,這是下階段需要解決的問題。
[1] SMELIK R M,TUTENEL T,BIDARRA R,et al.A survey on procedural modeling for virtual worlds[J].Computer Graphics Forum,2014,33(6):31-50.
[2] 湯國安,劉學軍,閭國年.數字高程模型及地學分析的原理與方法[M].北京:科學出版社,2005.
[3] 李清泉,楊必勝,史文中,等.三維空間數據的實時獲取、建模與可視化[M].武漢:武漢大學出版社,2003.
[4] DUCHAINEAU M,WOLINSKYM,SIGETI D E,et al.ROAMing terrain:Real-time optimally adapting meshes[A].Proceedings of the 8th Conference on Visualization′97[C].Washington,DC:IEEE Computer Society Press,1997.81-88.
[5] R?TTGER S,HEIDRICH W,SEIDEL H P.Real-time generation of continuous levels of detail for height fields[A].Proceedings WSCG′98[C].Plzen:Union Agency,1998.315-322.
[6] URLICH T.Rendering massive terrains using chunked level of detail control[A].Proceedings of ACM SIGGRAPH′02[C].San Antonio,TX:ACM Press,2002.
[7] CIGNONI P.Planet-sized Batched Dynamic Adaptive Meshes (P-BDAM)[A].Proceedings of the 14th IEEE Conference on Visualization (VIS′03)[C].Washington,DC:IEEE Computer Society,2003.147-154.
[8] ASIRVATHAM A,HOPPE H.Terrain rendering using GPU-based geometry clipmaps[A].PHARR M,FERNANDO R.GPU Gems 2[C].Boston,MA:Addison-Wesley,2005.46-53.
[9] LIVNY Y,KOGAN Z,EL-SANA J.Seamless patches for GPU-based terrain rendering[J].The Visual Computer,2009,25 (3):197-208.
[10] HU L,SANDER P V,HOPPE H.Parallel view-dependent refinement of progressive meshes[A].Proceedings of the Symposium on Interactive 3D Graphics and Games (I3D′09)[C].Boston,MA.New York:ACM,2009.169-176.
[11] JANG H Y,CHO C S,HAN J H.Multi-resolution terrain rendering with GPU tessellation[J].The Visual Computer,2015,31(4):455-469.
[12] KOCA ,G D KBAY U.A hybrid representation for modeling,interactive editing,and real-time visualization of terrains with volumetric features[J].International Journal of Geographical Information Science,2014,28(9):1821-1847.
[13] TATEOSIAN L,MITASOVA H,THAKUR S,et al.Visualizations of coastal terrain time series[J].Information Visualization,2014,13(3):266-282.
[14] 劉光.基于GIS的現代黃土地貌演化過程動態仿真研究[D].北京:北京大學,2003.
[15] 李欽,戴樹嶺,趙永嘉,等.分塊LOD大規模地形實時渲染算法[J].計算機輔助設計與圖形學學報, 2013,25(5):708-713.
A LOD Scheduling Method for Massive Dynamic Terrain Visualization
GAO Yong,LIU Jia-jun,GUO Xiao,WU Lun
(Institute of Remote Sensing and Geographical Information System,Peking University,Beijing 100871,China)
The processes modeling and simulation,e.g.geomorphic processes and land processes,will iteratively generate a series of terrain data in real time.This dynamic process needs more efficient visualization technologies.But the traditional preprocessing-based methods cannot fulfill it for the rendering delay caused by huge computation.A hybrid LOD scheduling method is presented integrating preprocessing and real-time update to visualize massive dynamic terrain.The terrain LODs are modeled and managed based on the scene graph,and a multithread model is introduced to update and render the terrain tiles dynamically.Instead of reconstructing the whole visualization data,this method only update and render the necessary local LOD tiles for the new terrain data iteratively.So the visualization efficiency can be improved significantly.
terrain visualization;level of detail model;scene graph;multithread;dynamic scheduling
2015-09-10;
2015-10-16
國家863支持項目(2011AA120301、2011AA120303)
高勇(1974-),男,博士,副教授,研究方向為高性能地理計算與空間數據挖掘。E-mail:gaoyong@pku.edu.cn
10.3969/j.issn.1672-0504.2016.01.002
P208
A
1672-0504(2016)01-0006-06