秦巴山區MODIS LST時序數據重建及特征分析

韓 曉 勇，韓 玲，戚 鵬 程

(1．長安大學地質工程與測繪學院，陜西 西安 710054；2.南陽師范學院環境科學與旅游學院，河南 南陽 473061)

秦巴山區MODIS LST時序數據重建及特征分析

韓 曉 勇1，韓 玲1，戚 鵬 程2

(1．長安大學地質工程與測繪學院，陜西 西安 710054；2.南陽師范學院環境科學與旅游學院，河南 南陽 473061)

地表溫度(LST)是大氣—地表水熱平衡的重要參數，但受云等大氣狀況影響，MODIS LST 時間序列產品存在大量缺失數據，限制了LST 數據使用。該文以陜南秦巴山區為研究區，從LST與高程、歸一化植被指數(NDVI)、經度和緯度的相關性特征出發，首先，針對不同的相關性特征自適應確定重建因子；然后，設計了以時點重建模型為主算法、總重建模型為備用算法的LST時間序列重建流程；最后，對研究區的地溫特征進行分析。結果表明，高程與LST呈較為穩定的負相關關系，NDVI隨植被的物候期變化與LST的相關性呈周期性變化，緯度與LST白天時點以正相關為主，夜晚則呈負相關特征；白天重建數據的平均絕對誤差(MAE)變化范圍較大，夜晚的變化范圍小且平均值較低;71.8%的白天重建數據和78.2%的夜晚重建數據的誤差可控制在3℃以內；研究區內平均地溫隨著高程的升高而下降，海拔1～2 km區域的年均地溫比海拔低于1km區域低2.0℃，海拔大于2 km區域的年均地溫比海拔1～2 km區域低2.6℃。

MODIS；地表溫度；相關性分析；數據重建；秦巴山區

0 引言

地表溫度(LST)是衡量大氣—地表水熱平衡極其重要的參數，在氣候變化、水文循環、生態環境等研究中具有重要意義[1]。傳統的LST數據來源于氣象站點觀測，并利用地統計插值獲得空間連續數據，我國中、西部地區氣象站點的密度低，分布不均勻，地形狀況復雜，地物類型多樣，使得地統計插值方法很難達到最優效果[2]。隨著遙感技術的不斷發展，可以利用熱紅外遙感獲取區域范圍的LST數據[3]。其中MODIS傳感器因為其優化的時空分辨率成為比較理想的數據源，在農業、生態環境以及全球變化等領域得到了廣泛應用[4，5]。

目前，MODIS LST產品在天氣晴朗條件下平均精度可達1 K[6]，但由于傳感器對地觀測時，受云、氣溶膠、觀測角度和太陽光照角度等影響，產生部分觀測值失真數據。MODIS LST標準產品已采用MODIS全球云掩膜產品進行失真像元的剔除，這又導致頻繁出現無效像元，從而制約了監測應用效果。為消除數據缺失和數據質量不高對數據使用帶來的影響，插補缺失觀測值，優化時間序列數據為相關研究提供更完備的基礎數據，國內外學者已開展對MODIS LST產品重建的研究。有學者從時間序列數據濾波處理的角度進行LST數據重建，如Julie等[7]使用HANTS變換對歐洲地區1982-1999年的LST數據進行去云重建，并對LST的年際變化特征進行分析；徐永明等[8]利用HANTS方法對長江三角洲地區2005年8天LST合成數據進行諧波擬合，得到重建的LST時間序列數據；蘇紅等[9]綜合考慮LST年變化和日變化的周期性，利用諧波分析方法對河南省2011年的LST時間序列進行擬合。此外，有學者從LST與其他地理因子的相關性出發,通過建立回歸模型實現LST數據重建。如柯靈紅等[10]利用高程與溫度的回歸關系重建青藏高原MODIS 8 天LST產品，并用氣象站0 cm LST進行驗證；劉梅等[11]利用無云區像元NDVI與LST 的負相關關系，基于短期NDVI 穩定的特點，重建臨近時段云區LST。

雖然對LST重建開展了大量研究，但由于LST在空間上差異較大，時間上變化較快，單從濾波處理的角度將LST數據分解成不同周期的諧波變量，會丟失LST的突變信息。若僅考慮高程或NDVI與LST的回歸關系，又將LST的回歸模型限制為單因子模型，難以保證長時間序列重建模型的穩健性。為此本文從長時間序列LST與高程、NDVI、經度和緯度4個因子的相關性出發，針對不同的相關性特征自適應確定回歸因子和回歸模型,以保證LST重建模型的穩健性，在此基礎上對所選研究區的地溫特征進行分析。

1 研究區域與數據預處理

1.1 研究區概況

本文選擇陜西省內秦巴山區(31°43′～35°6′N，105°29′～111°0′E)為研究區，主要包括商洛、漢中和安康，此外由于西安、寶雞南部也屬于秦巴山區,故本文將兩市也納入研究區，土地總面積98 713 km2。秦嶺是黃河、長江兩大水系的分水嶺和我國南北氣候分界線，其北坡為暖溫帶氣候，南坡為北亞熱帶氣候，在陜西境內連綿約500 km，海拔多在1 000 m以上，主峰太白山海拔3 767 m。巴山是嘉陵江和漢江的分水嶺，其北麓為北亞熱帶氣候，南麓為中亞熱帶，海拔平均1 500～2 500 m。該區是陜西省最大最主要的林區，屬于暖溫帶落葉闊葉林地帶和北亞熱帶常綠落葉闊葉混交林地。秦嶺與巴山之間，分布著漢中盆地、安康盆地、商丹盆地和洛南盆地，盆地內耕地集中，農業生產水平較高，是陜南重要的農業生產區和人口聚居區。此外研究區還包括了部分渭河平原和渭北旱塬丘陵溝壑區。

1.2 研究數據及預處理

MODIS地溫數據： 從EOSDIS(http://reverb.echo.nasa.gov/)下載2005-2014年的MODIS LST 8天合成產品，每期數據包括白天/夜晚4個時點(1:30、10:30、13:30、22:30)，該產品是8天內每日對應時點有效LST的算術平均，空間分辨率1 km，每年4個時點各46期數據。產品中還包括與每個像元對應的質量控制標記數據，使用MRT工具將每期HDF文件中的LST和LST質量控制數據分別進行拼接、裁剪，并統一轉換到WGS84 UTM N49投影帶，利用質量控制文件保留質量標記最優的LST像元值(標記值為0)。

MODIS植被指數數據：從EOSDIS獲取研究區對應年度的MODIS植被指數產品，空間分辨率1 km，每16天一期數據，10年共230期數據，利用質量控制文件剔除質量差的像元。使用較常用的Savitzky-Golay(S-G)濾波方法對MODIS NDVI產品進行數據重建，S-G濾波能有效抑制噪聲，并保證NDVI的變化趨勢[12]。為了使植被指數數據與LST數據的時間尺度一致，將S-G濾波后的相鄰兩期NDVI數據的均值作為二者中間期的NDVI數據，使得NDVI的時間分辨率與LST數據一致。S-G濾波處理后首尾幾期數據無法重建,故本文最終統一設置為每年37期數據。

DEM數據：本文從地理空間數據云網站(http://www.gscloud.cn/)下載了研究區的SRTM DEM數據，空間分辨率為90 m，利用雙線性插值法將DEM降尺度至1 km空間分辨率。

氣象數據：從中國氣象數據網(http://data.cma.gov.cn/)下載了國際交換站氣候資料日值數據集，使用數據集中西安站和漢中站的0 cm地溫數據對MODIS LST數據重建效果進行評價。

Z-Score標準化：由于研究數據間的量綱不同，為了使數據間具有可比性，保證LST重建模型的穩定性，將LST、DEM、經度和緯度數據進行Z-Score標準化如下:

(1)

式中：x為原始數據，μ和σ分別為有效LST抽樣數據均值與標準差。

由于2009年度NDVI數據質量整體較差，本文未對該年度進行LST重建。

2 研究方法

MODISLST數據進行重建的基本思想：由于LST在空間上差異大，時間上變化快，單從LST數據自身進行數據重建，難以保證數據精度。通過LST數據與其他地理因子的相關性分析，確定與當前時點LST數據相關性較強的因子，利用多元回歸方法構建重建模型。換言之，利用相關性分析確定能夠在某一時點反映LST數據特征的地理因子，進而建立多元回歸的LST重建模型，具體流程如圖1。

圖1 LST數據重建流程

Fig.1 Flow chart of LST reconstruction

2.1 隨機抽樣

為避免樣本值過度集中，根據研究區的地形特征，利用高程分段進行隨機抽樣，高程分段及各段抽樣比例見表1，隨機抽樣后只保留有效LST樣點數據。本文每次抽取660個樣點數據進行LST重建模型構建，另按此抽取方案抽取500個樣點的LST數據用于精度評價與控制。

表1 隨機抽樣點各高程分段比例

Table 1 Proportion of random sampling points in different elevation section

高程(km)<0.50.5～11～1.51.5～22～2.5>2.5比例(%)8304010102

2.2 重建模型

首先，將抽樣點海拔、NDVI、經度和緯度 4個因子與LST進行相關性分析，將相關性顯著的因子作為重建模型的構建因子。本文根據LST時序數據特點進行了總體相關性分析和時點相關性分析。總體相關性是將研究區內不同年度同期同時點全部LST樣點數據進行LST與4個因子的相關性分析，該參數可從宏觀上判斷不同時間LST與各因子的相關性，計算公式如式(2)。時點相關性是將LST樣點數據逐期同時點進行LST與4個因子的相關性分析。然后利用相關系數閾值判斷條件(白天時點絕對值大于0.2、夜晚時點大于0.15)自適應確定重建模型的構建因子。

(2)

對時點相關性分析確定的重建因子，采用逐步線性回歸的方法建立時點重建模型如式(3)，然后利用式(4)計算得到真實LST數據。時點重建模型用于有效抽樣數據占比大于25%，此時有效樣點數量較多，可保證回歸模型的可信度，是LST數據重建的主算法。對于研究區數據缺失面積過大(理論上大于75%的數據缺失)，則使用對應時點總體相關性分析確定的構建因子，基于逐步線性回歸的方法構建此時點的總重建模型以實現該類數據的重建，該方法是LST重建的備用算法。

(3)

LST=σLST*+μ

(4)

式中：αi為第i個影響因子的權重系數，xi為第i個影響因子，b為常數項，m為相關分析確定的影響因子數；μ和σ分別為有效LST抽樣數據均值與標準差，研究區內取值范圍分別為-12.2～32.9℃、1.8～5.52。

2.3 質量控制

本文引入MAE作為質量控制條件，根據蘇紅等[9]對LST數據修復精度，將MAE的閾值設置為3.0℃。利用重建LST與原始LST計算MAE，判斷是否超過閾值，若小于閾值則此模型可用，否則重新抽樣再次構建重建模型，并加入循環次數控制條件，本文設定為10，若10次循環均無法滿足MAE的閾值，則此時點數據無法完成重建。

(5)

式中：LST與LST′分別為原始數據與重建數據，n為驗證點數量。

2.4 地溫數據特征分析

數據重建完成后，本文按照小于1km、1～2km、大于3km3個高程區間對研究區內LST特征進行分類研究。首先利用重建后的LST數據平均計算月平均地溫，然后對月平均地溫數據按年度求出年平均地溫，并利用回歸趨勢線的方法對地溫數據的年際變化率進行了計算,公式如下[13]:

(6)

式中:n為時間序列長度，本文取n=9；LSTYi為第i年年均地溫值；θslope為趨勢線年際變化率；θslope>0，說明研究區的LST在該時段內呈增加趨勢；θsiope<0，表示呈減少趨勢；θslope=0，則表示未發生變化。

3 結果分析

3.1LST與影響因子相關性

3.1.1 總相關性特征 總相關系數可反映LST與4種影響因子的總體規律，圖2是研究區內4個時點按期序的LST與4種影響因子的相關系數折線圖。白天兩個時點(圖2a、圖2c)的相關系數時序變化較相似，夜晚兩個時點(圖2b、圖2d)的變化也較為相似，白天與夜晚差異較大。從各個因子看，經度在4個時點與LST的相關系數表現為波動較小的正相關，但相關性較低，白天兩個時點相關系數在0.2附近波動，晚上在0.15附近波動，為此本文將0.2和0.15分別設置為白天與夜晚相關系數的判斷閾值，相關系數的絕對值大于該閾值的影響因子才可進入重建模型的回歸分析。高程與LST的相關系數表現為較為穩定的負相關關系，白天兩個時點的相關系數均值為-0.43，夜晚兩個時點的均值為-0.39，說明高程在白天和夜晚對LST的影響較穩定，夜晚兩個時點1-5期的相關性較低，基本在-0.2左右，在9-34期內高程是與LST的相關性最高的因子，此時段高程對LST的影響最大。NDVI與LST的相關系數變化較大，白天兩個時點在1-10期(2月-4月中旬)與LST的相關性不顯著，此時段植被從休眠期到萌芽期，11-37期NDVI與LST表現為顯著的負相關特征，其中11-22期(4月下旬到7月底)二者間的相關性大于高程與LST相關性，此時植被從快速生長過渡到穩定生長階段， 23-33期(8月-10月中旬)NDVI與LST的相關系數變化較小,在-0.40附近波動，相關性略低于高程，此時植被處于穩定生長期，34-37期(10月下旬-11月)NDVI與LST的相關性急劇下降，此時段植被從衰落期過渡到休眠期；夜晚NDVI與LST的相關性整體較低，其中1-10期以正相關為主，除5-7期外均大于0.15，這應與此時段NDVI高值區主要分布在緯度較低的常綠林區有關，12-32期二者之間為負相關關系，其中15-27期的相關系數小于-0.15，33-37期又表現為正相關，夜晚的相關性變化也與植被的物候特征有關，但不如白天顯著。緯度與LST的相關性在白天兩個時點9-31期二者基本為正相關關系，其他期則表現為負相關，高值區主要集中在17-22期，該時段LST隨緯度增加呈增加的特征，但在夜晚二者之間均為負相關，除18-25期(夏季)的相關系數大于-0.15外，其他期的相關系數均小于-0.15，這說明區內除夏季之外夜晚的LST緯向逆分布明顯。

3.1.2 時點相關性特征 同總相關系數相比各影響因子與LST的時點相關系數時序變化隨機性較強，這與地溫短時間變化隨機性強的特點一致，海拔與LST的相關系數為較穩定的負相關特征，白天時點NDVI與LST的相關系數也能反映植被不同生長階段與LST的關系。利用時點相關系數最終確定參與重建模型的構建因子，表2中為最終各因子能參與重建模型構建的數量，其中高程占比達86%，其次為NDVI因子，且NDVI在白天時點的參與率比夜晚模型高，經度參與率最低，這說明了高程因子在LST數據重建中的重要性，但其他3個因子均擁有一定的參與度，因此LST重建時需綜合考慮幾個因子的影響。

圖2 LST與影響因子總相關性(時點：a.10:30;b.22:30;c.13:30;d.1:30 )

Fig.2 Total-correlation of LST and impact factors(Time point:a.10:30;b.22:30;c.13:30;d.1:30)

表2 各因子參與重建模型的期數Table 2 Periods of each impact factor participating reconstruction model

3.2 LST數據重建質量

圖3為利用有效樣點對9年4個時點的LST數據重建前后逐期計算的MAE，圖中黑點為平均的MAE，細線表示不同年份同期MAE的變化區間。白天LST修復數據的MAE(圖3a、圖3b)變化區間較大，其中2-16期的平均MAE均大于2℃，從17期開始逐漸降至2℃以下。夜晚兩個時點的平均MAE(圖3c、圖3d)較低，基本分布在1.5～2℃，且MAE的變化區間較小。第30期數據9年共36景影像，僅有4景可用主算法重建，其他景均為備用算法實現，導致該期的MAE整體偏高，這與備用算法較主算法的精度低有關。從表3中可知,71.8%的白天重建數據和78.2%的夜晚重建數據的數據偏差可控制在3℃以內，總體上90%以上數據誤差可控制在5℃以內，表明本文方法對原始影像中的缺失數據的重建是可行的。圖4是研究區代表性的一景(2005年15期 13:30)LST影像重建前后對比，左圖中的黑色區域為無值區，可見LST修復數據空間分布特征合理，并能抑制原始數據中噪聲數據的影響。

表3 重建LST與原始LST數據偏差頻率統計

Table 3 The frequency statistics of origin and reconstructed LST difference

時點<2[2,3)[3,5]>510:3051.119.819.29.913:3053.519.319.67.622:3058.018.615.38.11:3062.517.311.98.3

圖3 LST重建數據平均MAE(時點：a.10:30;b.13:30;c.22:30;d.1:30)

Fig.3 Average MAE of reconstructed LST(Time point:a.10:30;b.13:30;c.22:30;d.1:30)

圖4 原始影像( 左) 與重建影像( 右) 對比

Fig.4 Comparison of original(left)and reconstructed images(right)

利用研究區內的兩個基本氣象站點觀測的日均LST數據對重建數據進行驗證。先將氣象站點日均LST重新計算為與重建LST數據時間對應的均值數據，重建的LST數據則按期計算4個時點的均值。圖5為兩者之間的絕對誤差分布，重建數據與站點數據的絕對誤差較小，其中西安站的最大誤差為2.8℃，漢中站的最大誤差為-2.3℃，重建數據與站點數據絕對誤差在9-27期以低估為主，其他時段以高估為主，存在地溫較高的時段以低估為主、地溫較低的時段以高估為主的特點。

圖5 重建地溫與站點數據的絕對誤差Fig.5 Absolute errors of reconstructed LST and the meteorological data

3.3 地溫變化分析

圖6左圖為9年內不同高程的月均LST，月均LST年度周期性變化特征明顯，隨著高程的升高呈下降趨勢；右圖為年均LST，不同海拔的年均LST呈平行分布的特點，其中海拔1～2 km的區域年均LST比海拔小于1 km的區域平均低2.0℃，海拔大于2 km的區域年均LST比海拔1～2 km區域平均低2.6℃。利用回歸趨勢線的方法計算三者的年際LST變化率分別為：0.04、0.02和0.06，這表明不同海拔區域的年均LST均呈緩慢增長趨勢，且海拔大于2 km的區域年均LST增長速率最快，海拔1～2 km的區域年均LST增長較慢。

4 結論與討論

本文以秦巴山區為研究區域，在LST與影響因子相關性分析的基礎上對9年的MODIS LST數據產品進行重建，實現了研究區LST數據的高質量重建。研究取得了如下結論：1)LST數據與各影響因子的相關性不盡相同，高程與LST呈較為穩定的負相關關系；NDVI則隨植被的物候期變化與LST的相關性也呈周期性變化，白天兩個時點尤為明顯；緯度與LST的相關性白天與夜晚的差異較大，白天以正相關為主，而在夜晚全年呈負相關。2)白天LST重建數據的MAE變化區間較大，夜晚MAE的變化區間小且均值較低，基本分布在1.5～2℃，71.8%的白天重建數據和78.2%的夜晚重建數據的誤差可控制在3℃以內，總體上90%的數據誤差可控制在5℃以內。3)研究區內隨著高程的升高平均LST呈下降趨勢，海拔1～2 km的區域年均LST比海拔小于1 km的區域低2.0℃，海拔大于2 km的區域年均LST比海拔1～2 km的區域低2.6℃。不同海拔的年際LST變化率分別為：0.04、0.02和0.06，表明不同海拔區域的年均地溫均呈緩慢增長趨勢。

然而本文研究還存在如下問題需要改善：1)本文的精度評價以絕對誤差為主，在LST較高的情況下比較合適，對于LST較低的時段雖絕對誤差較低，但相對誤差可能較高，今后的研究中考慮LST溫度較低的情況下重建問題。2)本文的LST重建未考慮云陰影的影響，雖有學者開展了相關研究，認為將重建LST數據乘一系數轉換為真實LST，該系數的合理確定仍需研究。3)本文僅引入了4個影響因子，探討了空間分辨率較低的LST重建對于分辨率較高的遙感影像適用性問題。4)引入變量中僅包含一個NDVI時間變量，此因子是否削弱了LST的時間變化特征也需進一步研究。此外，本文中的算法是在秦巴山區建立、驗證和應用的，若應用于與研究區地貌、氣候差異較大的區域，則可能需要引入其他地理因子。

圖6 不同高程月平均LST與年平均LST

Fig.6 Monthly mean and annual mean LST in different elevation section

[1] QIN Z H， KARNIELI A．A mono-window algorithm for retrieving land surface temperature from Landsat TM data and its application to the Israel-Egypt Border region[J]．International Journal of Remote Sensing，2001，22(18)：3719-3746．

[2] 田輝，馬耀明，胡曉，等．使用MODIS 陸地產品LST 和NDVI監測中國中、西部干旱[J]．高原氣象,2007,26(5):1086-1095．

[3] 趙少華，邱國玉，秦其明，等． MODIS數據反演地表溫度的應用研究[J]．干旱區資源與環境，2010，24(3)：51-54．

[4] PARIDA B R，OINAM B，PATEL N R，et al． Land surface temperature variation in relation to vegetation type using MODIS satellite data in Gujarat state of India[J]．International Journal of Remote Sensing，2008，29(14)：4219-4235．

[5] 于曉靜，劉煥軍，曲長祥，等． 基于MODIS數據的黑龍江省黑土帶地溫空間格局分析[J]．水土保持研究，2012，19(2)：66-70．

[6] WAN Z M． New refinements and validation of the MODIS land-surface temperature/emissivity products[J]．Remote Sensing of Environment，2008，112(1)： 59-74．

[7] JULIEN Y，SOBRINO J，VERHOEF W． Changes in land surface temperatures and NDVI values over Europe between 1982 and 1999[J]．Remote Sensing of Environment，2006，103(1)：43-55．

[8] 徐永明，覃志豪，沈燕． 長江三角洲地區地表溫度年內變化規律與氣候因子的關系分析[J]．國土資源遙感，2010，83(1)：60-64．

[9] 蘇紅，劉峻明，王春艷,等． 基于時間序列MODIS LST產品的重構研究[J]．中國農業科技導報，2014，16(5)：99-107．

[10] 柯靈紅，王正興，宋春橋，等． 青藏高原東北部MODIS LST 時間序列重建及與臺站地溫比較[J]．地理科學進展，2011，30(7)：819-826．

[11] 劉梅，覃志豪，涂麗麗，等． 利用NDVI 估算云覆蓋地區的植被表面溫度研究[J]． 遙感技術與應用，2011，26(5)：689-697．

[12] 李杭燕，頡耀文，馬明國．時序NDVI數據集重建方法評價與實例研究[J]．遙感技術與應用， 2009，24(5)：596-602．

[13] STOW D D S，HOPE A． Variability of the seasonally integrated normalized difference vegetation index across the north slope of Alaska in the 1990s[J]．International Journal of Remote Sensing，2003，24(5)： 1111-1117．

Reconstruction of MODIS LST Time Series and Features Analysis in Qinling-Bashan Mountain Area

HAN Xiao-yong1,HAN Ling1,QI Peng-cheng2

(1.School of Geology Engineering and Geomatics,Chang′an University,Xi′an 710054；2.Environmental Science and Tourism College,Nanyang Normal University,Nanyang 473061,China)

Land surface temperature(LST)is an important parameter for measuring the atmosphere-surface water and heat balance.However,the LST data of MODIS time series product is commonly restricted by data missing owing to the influence of atmospheric conditions such as cloud.Qinling-Bashan Mountain area,south Shaanxi Province,was taken as the study area in this research.Based on the correlations among LST,elevation,NDVI,longitudinal and latitude,the reconstruction factor was determined self-adaptively.Then,a reconstruction process for LST time series data was designed,where the time-point-reconstruction model was considered as the master algorithm and total-reconstruction model was set as the alternate algorithm.Finally,the features of surface temperature were analyzed.The results suggest that the LST is negatively related with the elevation,while the correlation between NDVI and LST varies periodically with the phenology of vegetation.The relationship between latitude and LST data in daytime is mainly positive correlation,whereas the relationship is negative in nighttime.The mean absolute error of the daytime reconstruction data varies wider than that of the nighttime reconstruction data.Moreover,the mean value of the nighttime reconstruction data is relatively low.71.8% of data error about the daytime reconstruction data and 78.2% of data error about the nighttime reconstruction data was controlled within 3℃.The average temperature decreases with the increase of altitude in the research area.The annual average temperature of the area whose altitude is between 1 km and 2 km is 2.0℃ less than the area below 1 km,and the area above 2 km is 2.6℃ less than the area between 1 km and 2 km.

MODIS;LST;correlation analysis;data reconstruction;Qinling-Bashan Mountain area

2015-08-02；

2015-09-25

國家自然科學基金項目“基于空間化模擬方法修復遙感地表溫度圖像的數據缺失”(41201099)

韓曉勇(1981-)，男，博士研究生，主要從事遙感圖像處理與特征提取。E-mail:nytchxy@163.com

10.3969/j.issn.1672-0504.2016.01.014

TP79

A

1672-0504(2016)01-0071-07