張?chǎng)瓮? 吳秀芹

北京林業(yè)大學(xué)水土保持學(xué)院, 北京 100083;

北京林業(yè)大學(xué)水土保持國(guó)家林業(yè)局重點(diǎn)實(shí)驗(yàn)室, 北京 100083;

北京林業(yè)大學(xué)云南建水荒漠生態(tài)系統(tǒng)國(guó)家定位研究站, 云南建水 654399

植被凈初級(jí)生產(chǎn)力(NPP), 是指綠色植物經(jīng)過(guò)光合作用合成的有機(jī)物總量減去自身呼吸作用消耗的有機(jī)物之后, 剩余的可用于植物生長(zhǎng)、發(fā)育和繁殖的能量值, 用來(lái)表現(xiàn)綠色植物的固碳量(尹鍇等,2015)。利用NPP對(duì)地表植被固碳量進(jìn)行反演模擬,既可以反映植物的生長(zhǎng)狀況、完成對(duì)植被數(shù)量的監(jiān)測(cè), 又可以對(duì)碳循環(huán)以及整個(gè)全球氣候的變化研究起到重要輔助作用, 有很強(qiáng)的現(xiàn)實(shí)可用性。NPP估算模型主要有: 氣候生產(chǎn)力模型、生態(tài)學(xué)過(guò)程模型、光能利用率模型。其中氣候生產(chǎn)力模型只考慮氣候因子而忽略了許多影響植物干物質(zhì)積累的其他因素和植物對(duì)環(huán)境的反饋?zhàn)饔? 因此誤差較大(朱文泉等, 2005a; 生態(tài)學(xué)過(guò)程模型涉及多種復(fù)雜參數(shù)的采集, 難度很大。相比之下, 光能利用率模型數(shù)據(jù)獲取難度較小且精度較高, 另外其可通過(guò)遙感手段輔助計(jì)算, 節(jié)省了大量繁雜的野外試驗(yàn)測(cè)定步驟, 在近年來(lái)遙感技術(shù)的迅速發(fā)展推動(dòng)下, 逐漸成為 NPP測(cè)算研究方法的主流(朱文泉等, 2005b)。其中,CASA模型通過(guò)遙感數(shù)據(jù)、氣象數(shù)據(jù)對(duì)NPP進(jìn)行估算(Potter et al., 1993), 既避免了復(fù)雜參數(shù)的采集,又能較為準(zhǔn)確地進(jìn)行大范圍的NPP模擬, 在國(guó)內(nèi)外廣泛應(yīng)用于不同地區(qū)不同地表植被類(lèi)型固碳量的反演模擬。

我國(guó)曾對(duì)不同區(qū)域進(jìn)行了 NPP的測(cè)算并展開(kāi)研究。如基于CEVSA模型分析了青藏高原2000—2014年植被凈初級(jí)生產(chǎn)力時(shí)空變化(許潔等, 2020);基于CASA模型研究了中國(guó)西北地區(qū)植被NPP多時(shí)間尺度的變化, 并分析了其對(duì)氣候變化的響應(yīng)(賈俊鶴等, 2019); 結(jié)合 MODIS數(shù)據(jù)模擬了福建省植被NPP(余弘泳和余會(huì)康, 2018); 運(yùn)用CASA模型模擬了新疆近十年草地生態(tài)系統(tǒng) NPP并分析了其時(shí)空格局變化(楊紅飛等, 2014)。通過(guò)這些研究, 對(duì)區(qū)域生態(tài)系統(tǒng)變化情況有了進(jìn)一步的了解, 對(duì)于采取下一步的生態(tài)治理措施、研究區(qū)域發(fā)展起到了重要作用。在我國(guó)西南喀斯特地區(qū), 特殊的地質(zhì)構(gòu)造、雨熱同期的氣候條件、廣泛分布的碳酸鹽巖使其形成了獨(dú)特的碳循環(huán)和巨大的碳匯潛力, 在區(qū)域和全球碳循環(huán)過(guò)程中發(fā)揮著重要作用(王世杰等, 2017),巖溶地區(qū)植被 NPP的研究是植被建設(shè)和生態(tài)系統(tǒng)恢復(fù)的重要基礎(chǔ), 掌握其動(dòng)態(tài)變化情況對(duì)于該區(qū)域生態(tài)恢復(fù)與可持續(xù)發(fā)展具有重要意義。巖溶斷陷盆地位于云貴高原, 共有43個(gè)縣, 其中, 有32個(gè)位于云南, 由于其特殊的地形、地貌、氣候條件, 地表極易受到不可恢復(fù)的破壞和發(fā)生嚴(yán)重的水土流失現(xiàn)象, 使得云南省成為除貴州之外第二大石漠化分布地區(qū)(王宇等, 2017), 全省巖溶面積占總面積的29%。由于特殊的地質(zhì)構(gòu)造體系、強(qiáng)烈的構(gòu)造運(yùn)動(dòng)及各類(lèi)人為干擾(張?jiān)频? 2010), 云南石漠化面積不斷擴(kuò)大發(fā)展, 治理工作迫在眉睫。因此, 從2008年開(kāi)始, 為推進(jìn)石漠化防治, 全省開(kāi)展了石漠化綜合治理試點(diǎn)工程建設(shè), 至2012年, 65個(gè)石漠化重點(diǎn)縣全面開(kāi)展了綜合石漠化治理, 將生物措施、工程措施、農(nóng)業(yè)措施相結(jié)合, 突出重點(diǎn)、集中治理, 取得了明顯成效, 區(qū)域植被蓋度明顯增加、區(qū)域環(huán)境逐步改善。

綜上, 為了解近年石漠化治理工程實(shí)施以來(lái)云南斷陷盆地的生態(tài)環(huán)境變化, 本文用CASA模型在縣域尺度上模擬估算了其 2005—2019年間的 NPP值, 并分析了其時(shí)間、空間分異及動(dòng)態(tài)變化, 旨在提示NPP對(duì)石漠化治理效果、生態(tài)環(huán)境恢復(fù)狀況的響應(yīng), 以期為斷陷盆地下一步生態(tài)環(huán)境保護(hù)建設(shè)提供科學(xué)參考。

1 研究區(qū)概況

巖溶斷陷盆地位于我國(guó)云貴高原, 處于100°3′—105°10′E、22°39′—29°10′N(xiāo) 之間, 包括滇東至四川攀西鹽源及貴州西部的43個(gè)縣(史晨璐和吳秀芹, 2020), 其中云南省境內(nèi)有32個(gè)縣(圖 1)。斷陷盆地盆山共存、地形變化劇烈, 季節(jié)性干旱缺水且地下河深埋等自然條件, 導(dǎo)致其水土資源不匹配,生產(chǎn)條件落后。近些年, 城鎮(zhèn)化和經(jīng)濟(jì)開(kāi)發(fā)建設(shè)過(guò)程中斷陷盆地土地被過(guò)度占用, 導(dǎo)致人地關(guān)系矛盾日益突出, 加速了石漠化發(fā)展, 導(dǎo)致地區(qū)經(jīng)濟(jì)貧困。其中, 滇東為該區(qū)域巖溶發(fā)育最嚴(yán)重和重度石漠化地區(qū)分布最多的區(qū)域(曹建華等, 2016)。滇東區(qū)域終年太陽(yáng)高度角大、是西南巖溶斷陷盆地光照強(qiáng)度最大的區(qū)域。屬于熱帶季風(fēng)氣候, 干濕兩季分明,光照時(shí)間長(zhǎng)且季節(jié)性干旱嚴(yán)重, 且由于斷陷盆地氣候反差大、較之周?chē)降亟邓俣舭l(fā)快, 導(dǎo)致缺水極為嚴(yán)重、巖溶地貌發(fā)育顯著(李丹, 2019)。云南省石漠化分布廣, 尤其以東部、南部最為嚴(yán)重。按照地質(zhì)環(huán)境及石漠化發(fā)育狀況, 可將滇東劃分為:滇東北中山峽谷巖溶盆地區(qū)石漠化區(qū)、滇東巖溶斷陷盆地區(qū)石漠化區(qū)、滇東南巖溶峰叢洼谷地區(qū)石漠化區(qū)三個(gè)巖溶分布區(qū)(呂文麗, 2014)。自2000年以來(lái), 云南省積極進(jìn)行生態(tài)環(huán)境建設(shè)、治理工作, 積極應(yīng)對(duì)水土流失、石漠化等環(huán)境難題。從2008年起,政府投入大量資金對(duì) 12個(gè)重點(diǎn)縣實(shí)施石漠化綜合治理試點(diǎn)工程, 全省持續(xù)采用石漠化綜合治理、天然林生態(tài)林公益林保護(hù)、退耕還林還草等生態(tài)工程,采取人工造林、封山育林(草)、陡坡改梯等方式, 加快巖溶地區(qū)林草植被保護(hù)與石漠化土地修復(fù)。2012年, 石漠化綜合治理從 12個(gè)縣增加到 65個(gè)縣, 基本涵蓋了斷陷盆地所在縣域。

圖1 研究區(qū)位置示意圖Fig. 1 Location map of the study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來(lái)源及數(shù)據(jù)處理

2.1.1 NDVI數(shù)據(jù)來(lái)源與處理

NDVI數(shù)據(jù)來(lái)源于MODIS數(shù)據(jù), 采用NASA官網(wǎng)MOD13Q1系列產(chǎn)品, 該產(chǎn)品是Terra衛(wèi)星采集的以16天為周期的全球250 m分辨率NDVI數(shù)據(jù), 根據(jù)研究需要和影像采集規(guī)則, 選取研究區(qū)云南所在的h26v06、h27v06兩景影像, 并進(jìn)行影像的批量拼接、轉(zhuǎn)投影、最大值合成處理。借助 ArcGIS工具提取出云南省斷陷盆地 2005、2010、2015、2019年各月份 NDVI數(shù)據(jù), 并統(tǒng)一定義為 Albers投影,將一年12個(gè)月數(shù)據(jù)進(jìn)行波段組合, 形成年均NDVI分布圖, 保存?zhèn)溆谩?/p>

2.1.2 氣象數(shù)據(jù)來(lái)源及處理

氣象數(shù)據(jù)來(lái)源于中國(guó)地面氣候資料月值數(shù)據(jù)集, 經(jīng)過(guò)數(shù)據(jù)整理, 在 Arcmap中加載云南省各氣象站點(diǎn)點(diǎn)位信息, 并通過(guò)克里金插值法插值得到2005年、2010年、2015、2019年1—12月云南氣溫、降水分布圖, 設(shè)置像元大小為250 m×250 m。再運(yùn)用波段組合工具將每年12個(gè)月的氣溫、降水圖分別組合, 得到年平均氣溫圖、年總降水量圖, 最后根據(jù)云南斷陷盆地邊界范圍進(jìn)行裁剪。同理, 利用中國(guó)地面日氣象資料數(shù)據(jù)集, 根據(jù)其日照信息文件(SSD)及溫度信息, 計(jì)算出各站點(diǎn)太陽(yáng)輻射值,后利用距離倒數(shù)插值法進(jìn)行插值, 再用與制作氣溫、降水圖相同的方法, 得到太陽(yáng)輻射分布圖。

2.1.3 CASA模型靜態(tài)參數(shù)配置

在一定范圍內(nèi), FPAR與NDVI、SR存在較好的線(xiàn)性關(guān)系(周偉等, 2017), 故通過(guò)NDVI、SR加權(quán)平均或平均值估算CASA模型中的FPAR參數(shù)。靜態(tài)參數(shù)文件用于配置各植被類(lèi)型的最大光能利用率及NDVI最大、最小值, 本實(shí)驗(yàn)參考朱文泉研究成果(朱文泉等, 2006)配置靜態(tài)參數(shù)文件。植被類(lèi)型及對(duì)應(yīng)的參數(shù)值如表1所示。

表1 各植被類(lèi)型靜態(tài)參數(shù)Table 1 Static parameters of various vegetation types

2.1.4 植被分布圖處理

植被分布圖來(lái)源于 ESA CCI的全球植被覆蓋數(shù)據(jù), 分別下載2005、2010、2015、2019年植被覆蓋數(shù)據(jù), 掩膜提取出云南省斷陷盆地植被覆蓋數(shù)據(jù),并轉(zhuǎn)化為與其他數(shù)據(jù)相同的Albers投影類(lèi)型。同時(shí),對(duì)其進(jìn)行重分類(lèi)操作, 分成與靜態(tài)參數(shù)表相對(duì)應(yīng)的22種類(lèi)型以進(jìn)行CASA模型運(yùn)算。

2.2 研究方法

采用CASA模型對(duì)NPP進(jìn)行估算。CASA模型是一種基于過(guò)程的遙感模型, 基于CASA模型可對(duì)生態(tài)系統(tǒng)的生產(chǎn)力、土壤碳進(jìn)行模擬(朱文泉等,2005b, 2006, 2007)。其原理是通過(guò)植被光合有效輻射、實(shí)際光能利用率來(lái)估算 NPP。即 NPP(x,t)=APAR(x,y)*ε(x,y)(Potter et al., 1993; 朱文泉等,2007)。

(1)APAR(x,y)=SOL(x,t)×FPAR(x,t)×0.5

其中, SOL(x,t)表示像元x在t月的太陽(yáng)總輻射量, FPAR(x,t)表示植被層對(duì)入射光合有效輻射的吸收比例, 常數(shù) 0.5表示植被能利用的太陽(yáng)有效輻射占太陽(yáng)總輻射比例。

(2)ε(x,y)由最大光能利用率、溫度水分脅迫因子共同決定。

ε(x,y)=Tε1(x,y)×Tε2(x,y)×Wε(x,y)×εmax

由于CASA模型比較成熟, 具體公式計(jì)算方法與參數(shù)含義不再贅述。

3 結(jié)果分析

3.1 2005—2019年云南省斷陷盆地NPP時(shí)間變化

2005年到2019年, 云南省斷陷盆地NPP平均值總體上呈現(xiàn)出先減少后增加再減少的趨勢(shì), 為便于分析變化, 分別統(tǒng)計(jì)斷陷盆地各典型土地覆蓋類(lèi)型在該期間的變化情況。各植被類(lèi)型NPP值見(jiàn)表2與圖2。

圖2 各植被類(lèi)型NPP均值Fig. 2 Mean value of NPP of each vegetation type

表2 各植被類(lèi)型NPP值Table 2 NPP of various vegetation types

2005到2010期間, 林地類(lèi)型NPP平均值明顯下降, 其中, 常綠針葉林NPP平均值由626.92 gc/m2下降至218.49 gc/m2, 約降低408.43 gc/m2; 常綠闊葉林在 2005年至 2010年 NPP平均值由1 320.11 gc/m2下降至468.41 gc/m2; 落葉闊葉林由862.12 gc/m2降至305.75 gc/m2。草地NPP平均值也有明顯下降, 由614.95 gc/m2下降至218.96 gc/m2。此外, 灌叢、耕地NPP均值均有不同程度的下降。

2010到2015年間, NPP平均值略有上升, 且有持續(xù)上升趨勢(shì)。按照各土地覆蓋類(lèi)型來(lái)看, 常綠針葉林NPP均值約增加18.76 gc/m2, 常綠闊葉林NPP均值約增加 151.09 gc/m2, 落葉闊葉林均值約增加57.63 gc/m2, 草地NPP均值約增加53.15 gc/m2。

2015到2019年間, NPP平均值略有下降, 但相較于2010年NPP水平仍有一定提升。按各類(lèi)土地覆蓋類(lèi)型來(lái)看, 常綠闊葉林 NPP均值下降最明顯,約減少154.01 gc/m2; 落葉闊葉林NPP均值約減少37.6 gc/m2; 草地NPP均值約減少31.9 gc/m2。

從NPP總值變化來(lái)看, 2005—2010年, 斷陷盆地植被 NPP總量下降, 2010—2015年總量上升,2015—2019年 NPP總量略有下降。結(jié)合同期土地利用的轉(zhuǎn)移變化來(lái)看, 2005—2010年草地、林地、灌叢等植被面積均有減少, NPP總量也有一定下降。2010年以后, 隨著石漠化治理措施的實(shí)行以及大量植樹(shù)造林工程的推進(jìn), 林地、草地面積逐步增加,對(duì)提高 NPP水平起到了顯著作用, 但 2015—2019年, 各類(lèi)植被NPP水平又有一定下降。

2005—2019年NPP分布圖見(jiàn)圖3。

圖3 2005(a)、2010(b)、2015(c)、2019(d)年云南省斷陷盆地NPP分布圖Fig. 3 Distribution of NPP in Yunnan fault-depression basins in 2005 (a), 2010 (b) , 2015 (c) and 2019 (d)

3.2 云南省斷陷盆地NPP空間分異及變化

借助Arcgis中分區(qū)統(tǒng)計(jì)的方法, 統(tǒng)計(jì)了云南省2005—2019年斷陷盆地所在縣域內(nèi)的NPP平均值,并繪制了柱狀圖進(jìn)行對(duì)比, 分析其空間分布特征。同時(shí), 結(jié)合各縣域石漠化分布狀況, 分析其石漠化治理工程實(shí)施以來(lái)NPP變化情況, 以探究NPP對(duì)石漠化的響應(yīng)。各縣域2005—2019年間NPP均值變化情況見(jiàn)表3與圖4。

表3 云南省2005—2019年斷陷盆地所在縣域NPP均值統(tǒng)計(jì)(單位: gc/m2)Table 3 Statistics of average NPP in Yunnan fault-depression basins in 2005–2019(unit: gc/m2)

圖4 2005(a)、2010(b)、2015(c)、2019(d)年云南省斷陷盆地縣域NPP分布均值圖Fig. 4 Bar graph of average NPP in Yunnan fault-depression basins in counties in 2005(a), 2010 (b), 2015 (c), 2019 (d)

由以上NPP分布圖和分區(qū)統(tǒng)計(jì)表來(lái)看, 云南省斷陷盆地NPP空間分布不均。從總體來(lái)看, 云南省斷陷盆地東部、南部地區(qū)NPP均值較高, 而西部、北部等地NPP值較低, 空間上呈現(xiàn)出自東向西、自南向北遞減的趨勢(shì)。以2015年NPP均值分布為例,其中, 西部五華區(qū)、盤(pán)龍區(qū)、呈貢縣、江川縣、澄江縣等NPP均值較低, 在101～213 gc/m2之間; 丘北縣、羅平縣、馬關(guān)縣、師宗縣NPP平均水平較高, 均大于300 gc/m2。

由NPP總值分布情況來(lái)看, 丘北縣、宣威縣、祿勸彝族苗族自治縣、彌勒縣、建水縣NPP總值較高, 均大于1.65×107gc; 通?？h、江川縣、 澄江縣、呈貢縣、盤(pán)龍區(qū)、五華區(qū) NPP總值最低, 均小于2.8×106gc。NPP總值在2005—2019年間也呈現(xiàn)出先下降后上升再下降的波動(dòng)趨勢(shì)。

由NPP動(dòng)態(tài)變化的空間分布來(lái)看, 2005—2010年, 云南省境內(nèi)斷陷盆地各縣NPP均值均有明顯下降。其中, 西北部、南部地區(qū)NPP均值下降幅度較大, 如馬關(guān)縣、個(gè)舊縣、師宗縣、富民縣、祿勸彝族苗族自治縣, NPP均值下降均超過(guò) 500 gc/m2。2010—2015年, 西南部通海縣、玉溪縣、江川縣等NPP均值變化不大, 而盤(pán)龍區(qū)、富民縣、馬龍縣、尋甸回族彝族自治縣等原 NPP均值較低的地區(qū),NPP水平均有較大提升。2015—2019年, 各縣NPP均值均有一定下降, 其中開(kāi)遠(yuǎn)縣、個(gè)舊縣、建水縣、昆明市各區(qū)下降幅度較大, 有些甚至低于 2010年NPP值。

4 結(jié)論與討論

4.1 結(jié)論

(1)從時(shí)間尺度上來(lái)看, 2005—2019年云南省內(nèi)斷陷盆地 NPP平均水平總體上呈現(xiàn)出先減少后增加再減少的趨勢(shì)。2005—2010年間, 林地、草地等植被NPP均值明顯下降、提供的生產(chǎn)力總值也有顯著下降, 但林地面積變化不顯著, 同時(shí), 城鎮(zhèn)面積有顯著提升。2010—2015年間, 林地、草地等植被NPP均值有一定上升, 且有持續(xù)上升趨勢(shì), 同時(shí),植被面積和提供的 NPP總量也有一定上升。2015—2019年, NPP均值又有小幅回落, 植被面積也有所減少。

(2)從空間上來(lái)看, 首先, 云南省內(nèi)斷陷盆地NPP空間分布不均, 東部、南部地區(qū)NPP指數(shù)較高,而西部、北部等地NPP值較低, 空間上呈現(xiàn)出自東向西、自南向北遞減的趨勢(shì), 同時(shí), 由NPP動(dòng)態(tài)變化的空間分布來(lái)看, 2005—2010年, 除西部幾個(gè)NPP均值極低的縣之外, 其余大部分地區(qū)NPP均值有顯著下降; 2010—2015年, 大部分地區(qū)NPP均值均有提升, 其中, 斷陷盆地中部、西北部地區(qū) NPP變化較大, 其他如東部、南部地區(qū)NPP水平亦有上升趨勢(shì); 2015—2019年, 大部分地區(qū)NPP有一定下降, 其中南部、西部地區(qū)變化最為顯著。

4.2 討論

結(jié)合石漠化研究情況來(lái)看, 云南省石漠化面積居全國(guó)第二, 截止到2016年底, 全省石漠化面積達(dá)到 235.2公頃, 而由于云南斷陷盆地特殊的地質(zhì)構(gòu)造及脆弱的生態(tài)環(huán)境, 其所在的多個(gè)縣域均存在不同程度、不同面積的石漠化土地。其中, 馬關(guān)縣(32.69%)、丘北縣(23.91%)、羅平縣(22.15%)石漠化發(fā)生率最高, 均在 20%以上, 其他如師宗縣、蒙自縣、宣威縣、滬西縣石漠化發(fā)生率也在10%以上。2011—2016年, 云南省全省持續(xù)實(shí)施了石漠化綜合治理等生態(tài)工程, 植被覆蓋面積有顯著上升。經(jīng)過(guò)石漠化綜合治理, 截至2016年底, 云南全省石漠化土地面積減少48.8萬(wàn)公頃, 其中, 斷陷盆地地區(qū)羅平縣、馬龍縣、尋甸回族彝族自治縣、富民縣石漠化變化率較高, 分別為–20.92、–19.48、–33.13、–43.5。結(jié)合斷陷盆地所在縣域NPP變化來(lái)看, 2010—2015年, 羅平縣NPP均值由261.69 gc/m2上升至314.64 gc/m2, 約上升20.23%; 馬龍縣NPP均值由215 gc/m2上升至266.3 gc/m2, 約上升23.86%; 尋甸回族彝族自治縣 NPP均值由 228.24 gc/m2上升至275.64 gc/m2, 約上升 20.77%, 在一定程度上反映了石漠化治理效果。但2015—2019年間, NPP水平又出現(xiàn)了一定程度的回落, 說(shuō)明斷陷盆地區(qū)生態(tài)治理修復(fù)仍將是一項(xiàng)艱難的長(zhǎng)期工程。

另外, NPP變化受多種因素影響, 如氣候、地形等。黃曉云等研究得出NPP與氣溫、降水變化和分布情況顯著相關(guān), 氣溫上升對(duì)植被生長(zhǎng)有利, 因此氣溫與 NPP呈正相關(guān), 降水也與 NPP變化呈正相關(guān)。由于受到高太陽(yáng)輻射量影響, 溫度對(duì)斷陷盆地植被生長(zhǎng)的影響大于降水, 斷陷盆地氣候?qū)PP變化的影響在空間上呈現(xiàn) 4種模態(tài): 氣溫與降水量皆上升時(shí), NPP上升; 氣溫和降水量都下降時(shí), NPP上升; 氣溫上升但降水量下降時(shí), NPP上升; 氣溫下降、降水量上升時(shí), NPP下降(黃曉云, 2013)。研究表明, 2000年以來(lái), 喀斯特?cái)嘞菖璧貐^(qū)氣候存在逐漸變暖變干的趨勢(shì), 溫度的顯著升高和降水的顯著減少都對(duì)NPP變化產(chǎn)生了影響(莊義琳等, 2019)。左麗媛和高江波(2020)用地理探測(cè)器的方法證明了植被覆蓋度、土地利用類(lèi)型、溫度、海拔、坡度、降水等因子均可對(duì)NPP產(chǎn)生影響。因此, NPP的變化可能是各種因素共同作用的結(jié)果, 本文限于研究時(shí)間與研究方法, 尚未具體探究NPP變化與氣候、地形、巖性等影響因子的關(guān)系, NPP的變化原因還有待進(jìn)一步深入研究。

Acknowledgements:

This study was supported by National Key Research and Development Program of China (Nos.2016YFC0502500 and 2016YFC0502506), and National Natural Science Foundation of China (No.41671080).