典型喀斯特地區(qū)土壤水分遙感反演
——以漣江流域?yàn)槔?/h1>
2019-10-17 02:36:14梁建方周秋文韋小茶馬龍生顏紅
人民珠江 2019年9期
梁建方，周秋文，韋小茶，馬龍生，顏紅
(貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院，貴州貴陽(yáng)550001)
土壤是重要的自然資源，土壤水分作為陸地和大氣能量交換過(guò)程中的重要因子，對(duì)陸表蒸散、水的運(yùn)動(dòng)和碳循環(huán)具有控制作用，對(duì)自然環(huán)境和人類社會(huì)的可持續(xù)發(fā)展有著重要影響[1-3]。獲取區(qū)域乃至全球尺度的土壤水分空間分布信息，是研究全球范圍內(nèi)水循環(huán)、能量平衡和氣象預(yù)報(bào)的先決條件，對(duì)改善區(qū)域及全球氣候模式預(yù)報(bào)結(jié)果、進(jìn)行農(nóng)作物旱情監(jiān)測(cè)及估產(chǎn)、自然和生態(tài)環(huán)境問(wèn)題的研究起重要作用[4-5]。傳統(tǒng)土壤水分測(cè)量方法只能得到有限測(cè)量點(diǎn)數(shù)據(jù)，僅能反應(yīng)測(cè)量點(diǎn)周圍有限區(qū)域的土壤水分情況。而通過(guò)星載遙感衛(wèi)星傳感器提供的監(jiān)測(cè)手段，能獲取大區(qū)域的土壤水分信息，同時(shí)具有時(shí)效快、動(dòng)態(tài)性較強(qiáng)及時(shí)間空間分辨率高等優(yōu)點(diǎn)，能為區(qū)域尺度土壤水分信息的獲取提供良好手段[6-8]。
光學(xué)遙感監(jiān)測(cè)土壤水分的方法主要有熱慣量法、植被指數(shù)法等，但熱慣量法僅適用于裸地及低值被覆蓋區(qū)域，而植被指數(shù)法僅適用于植被覆蓋較高區(qū)域[1,8]。Price和Carlson研究發(fā)現(xiàn)植被指數(shù)和地表溫度之間存在負(fù)相關(guān)關(guān)系，且以植被指數(shù)和地表溫度為橫縱坐標(biāo)生成的散點(diǎn)圖呈三角形狀[9-10]。Sandholt等[11]進(jìn)一步研究植被指數(shù)與地表溫度的關(guān)系，提出利用溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index，TVDI)來(lái)監(jiān)測(cè)和估算土壤水分的方法。國(guó)內(nèi)外學(xué)者基于該方法進(jìn)行區(qū)域地表水分狀況等開(kāi)展大量研究，并取得了豐碩成果[12]。王秀君等[13]基于Landsat TM數(shù)據(jù)反演了河北省黃驊市土壤水分，結(jié)果表明TVDI與土壤濕度顯著相關(guān)；曹雷等[14]利用TVDI在艾比湖地區(qū)進(jìn)行試驗(yàn)，發(fā)現(xiàn)TVDI能有效反演區(qū)域土壤水分，且精度較高；陳斌等[15]以內(nèi)蒙古錫林郭勒盟地區(qū)為研究對(duì)象，發(fā)現(xiàn)TVDI可以用來(lái)監(jiān)測(cè)旱情變化及評(píng)價(jià)草原干旱狀況；沙莎等[16]利用歷史遙感數(shù)據(jù)建立的特征空間，發(fā)現(xiàn)TVDI能夠較好地指示每年土壤相對(duì)濕度RSM的時(shí)空分布特征；吳孟泉、張喆等[17-18]研究均發(fā)現(xiàn)TVDI能夠很好地用于山區(qū)干旱預(yù)警與監(jiān)測(cè)。
喀斯特地區(qū)生態(tài)系統(tǒng)敏感度高, 災(zāi)變承受能力低, 環(huán)境容量小, 因而成為典型的生態(tài)脆弱區(qū)[19]。其土層淺薄、土壤貧瘠、成土速率低、水文過(guò)程變化迅速、地表土壤水受特殊地質(zhì)地貌及水文條件等自然因素影響較大[20-24]。以貴州省為中心的中國(guó)西南喀斯特地區(qū)是世界上面積最大的喀斯特連續(xù)帶，該地區(qū)地表破碎、地形復(fù)雜，由于森林砍伐、過(guò)度放牧、農(nóng)業(yè)耕種和地下水過(guò)度開(kāi)采等行為，地表植被的生長(zhǎng)和人類生產(chǎn)生活深受土壤水分貧瘠影響[25-30]。郭巍等[26]采用TVDI對(duì)云南省2010年春季旱情進(jìn)行監(jiān)測(cè)，揭示了云南省入春以來(lái)的旱情動(dòng)態(tài)變化及旱情分布特點(diǎn)；王海等[27]運(yùn)用云南省境內(nèi)的TVDI與昆明3個(gè)站點(diǎn)的土壤濕度做線性相關(guān)分析表明，TVDI對(duì)云南境內(nèi)的干旱監(jiān)測(cè)具有較好的適用性和可靠性；康為民等[28]利用同一時(shí)期TVDI與當(dāng)?shù)貧庀笳军c(diǎn)的土壤濕度數(shù)據(jù)進(jìn)行定量驗(yàn)證，發(fā)現(xiàn)TVDI與土壤濕度顯著相關(guān)，適用于大區(qū)域復(fù)雜地形的干旱檢測(cè)與預(yù)警。
然而已有研究大都基于MODIS數(shù)據(jù)與Landsat TM/ETM+數(shù)據(jù)來(lái)進(jìn)行干旱監(jiān)測(cè)研究，而利用Landsat8 OLI_TIRS數(shù)據(jù)進(jìn)行研究干旱監(jiān)測(cè)研究較少。喀斯特地貌的典型特征是地表破碎、地形復(fù)雜，土壤水分分布情況隨降雨及地表情況的不同而表現(xiàn)出不同區(qū)域特征。因此，本文基于Landsat8 OLI_TIRS數(shù)據(jù)，選取貴州省典型喀斯特地區(qū)的漣江流域?yàn)槔瑢?duì)流域的土壤水分進(jìn)行遙感反演，并實(shí)現(xiàn)干旱時(shí)空分布情況的動(dòng)態(tài)監(jiān)測(cè)，同時(shí)結(jié)合該區(qū)野外實(shí)測(cè)土壤含水量數(shù)據(jù)進(jìn)行驗(yàn)證分析。研究結(jié)果可為相關(guān)部門制定合理有效的抗旱、防旱措施提供科學(xué)依據(jù)，對(duì)喀斯特地區(qū)水土保持政策的制定、石漠化治理措施實(shí)施具有非常重要的指導(dǎo)作用，是喀斯特地區(qū)自然、經(jīng)濟(jì)、社會(huì)可持續(xù)發(fā)展的重要需求。
1 研究區(qū)概況
漣江流域(圖1)位于貴州省西南部(106°14′～106°51′E，25°32′～26°26′N)，面積約3 006.67 km2，流域海拔高度在517～1 650 m之間，地勢(shì)北高南低，中部平緩，其境內(nèi)主要分布碳酸鹽巖，喀斯特地貌典型發(fā)育。屬亞熱帶季風(fēng)氣候，降水豐富，但降水量年際、季節(jié)性分布不均，而地表和地下水系構(gòu)成的二元結(jié)構(gòu)往往會(huì)導(dǎo)致降水滲漏，土壤儲(chǔ)水能力低，加上不合理的人類活動(dòng)造成植被不斷遭受破壞，全區(qū)一年四季均存在輕重程度不相同的旱情。
圖1 漣江流域地理位置及高程
2 數(shù)據(jù)來(lái)源與方法
2.1 數(shù)據(jù)來(lái)源及預(yù)處理
遙感數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)(http://www.gscloud.cn)，獲取2013年9月29日、2014年1月19日與2017年4月1日的Landsat8 OLI_TIRS數(shù)據(jù)，以及ASTER GDEM高程數(shù)據(jù)，空間分辨率為30 m。利用ENVI 5.1軟件對(duì)遙感影像進(jìn)行預(yù)處理，主要包括輻射定標(biāo)、裁剪和大氣校正等。選取漣江子流域(揚(yáng)眉河流域)中具有代表性的林地、灌木林地、草地、耕地，于2017年4月2日使用TDR-300土壤水分儀測(cè)量得到24個(gè)樣點(diǎn)表層(0～10 cm)土壤體積含水量，并用手持GPS獲取采樣點(diǎn)的經(jīng)緯度信息。
2.2 溫度植被干旱指數(shù)估算
Price和Carlson等研究發(fā)現(xiàn)以NDVI與Ts為橫縱坐標(biāo)生成散點(diǎn)圖構(gòu)成的空間關(guān)系為三角形關(guān)系。Sandholt等對(duì)NDVI-LST特征空間進(jìn)行簡(jiǎn)化，提出溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index，TVDI)概念，計(jì)算公式如下：
(1)
式中TVDI——溫度植被干旱指數(shù)；Ts——地表真實(shí)溫度；Tsmax——某一NDVI對(duì)應(yīng)最大地表溫度，即干邊；Tsmin——地表最低溫度，即濕邊。干邊與濕邊方程分別為：Tsmax=a·NDVI+b,Tsmin=c·NDVI+d，其中a、b、c、d為干、濕邊擬合方程系數(shù)。
模型中地表真實(shí)溫度Ts由輻射傳導(dǎo)方程法獲得，具體步驟如下。
a) 計(jì)算植被指數(shù)NDVI。NDVI由遙感影像中，近紅外波段的反射值與紅波段的反射值之差比上兩者之和，其計(jì)算公式為：
(2)
式中NIR——近紅外波段像元亮度值；R——紅波段像元亮度值。
b) 計(jì)算地表比輻射率ε。本文采用Van經(jīng)驗(yàn)公式[19]計(jì)算地表比輻射率：
ε=1.0094 + 0.047·ln·(NDVI)
(3)
式中，NDVI值在0.157～0.727之間，NDVI小于0的像元只要是水體，其體表比輻射率近似0.995；城鎮(zhèn)很多建筑物像元NDVI位于0～0.157之間，其地表比輻射率可設(shè)置為0.923；而對(duì)于NDVI大于0.727像元可看作植被完全覆蓋，地表比輻射率為0.986。
c) 計(jì)算相同溫度下黑體的輻射亮度值B(TS)。衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ由三部分組成：大氣向上輻射亮度L↑；地面的真實(shí)輻射亮度經(jīng)過(guò)大氣層之后到達(dá)衛(wèi)星傳感器的能量；大氣向下輻射到達(dá)地面后反射的能量。其輻射傳輸方程為：
Lλ=[ε·B(TS) + (1-ε)L↓]·τ+L↑
(4)
根據(jù)普朗克定律得到的黑體熱輻射亮度B(TS)，其公式為：
B(TS) = [Lλ-L↑-τ·(1-ε)L↓]/(τ·ε)
(5)
式中Lλ——TIRS影像第十波段大氣校正后的輻射定標(biāo)值；ε——地表輻射率；TS——地表真實(shí)溫度；τ——大氣在熱紅外波段透過(guò)率。
d) 計(jì)算地表真實(shí)溫度TS。在獲取B(TS)后，根據(jù)普朗克公式的反函數(shù)，求得地表真實(shí)溫度TS，其公式為：
TS=K2/ln(K1/B(TS)+ 1)
(6)
式中B(TS)——溫度為TS黑體熱輻射亮度；K1、K2——為傳感器定標(biāo)常數(shù)。
2.3 表層土壤水計(jì)算
利用野外表層土壤水分(0～10 cm)采樣數(shù)據(jù)，建立土壤水分與溫度植被干旱指數(shù)之間的線性相關(guān)模型：
W=a·TVDI+b
(7)
式中W——表層土壤水分；TVDI——溫度植被干旱指數(shù)；a、b——線性系數(shù)。
3 結(jié)果與分析
3.1 干、濕邊擬合
由于NDVI值小于0的像元主要是水體和云，濕度接近100%，NDVI值在0～0.157時(shí)的像元主要是城市建筑物與水泥地等。為避免水體、云及建筑物對(duì)特征空間的影響，在干濕邊擬合時(shí)不考慮此類像元。因此，選取NDVI值大于0.157的像元，利用IDL編程實(shí)現(xiàn)相同NDVI值對(duì)應(yīng)的TS的最大值和最小值提取，步長(zhǎng)值為0.002，構(gòu)成NDVI-TS特征空間的形狀近似呈三角形(圖2)。從3個(gè)季節(jié)的擬合結(jié)果可以看出，干邊與濕邊的決定系數(shù)R2分別在0.9和0.5以上，擬合精度較高。且NDVI與地表溫度LST近似呈線性關(guān)系，隨著NDVI值的增大，地表溫度最大值LSTmax逐漸減小，最小值LSTmin逐漸增大。
a) 秋季
b) 冬季
c) 春季圖2 Ts-NDVI特征空間
3.2 TVDI干旱等級(jí)劃分
為從宏觀了解漣江流域干旱情況，根據(jù)已擬合的干、濕邊方程，計(jì)算研究區(qū)內(nèi)各個(gè)像元TVDI，根據(jù)TVDI灰度值范圍進(jìn)行干旱等級(jí)劃分指標(biāo)[21]，將干旱情況劃分為五級(jí)：濕潤(rùn)(0～0.2)、正常(0.2～0.4)、輕旱(0.4～0.6)、干旱(0.6～0.8)、重旱(0.8～1.0)。研究結(jié)果表明：研究區(qū)3個(gè)季節(jié)干旱情況在時(shí)間分布上總體一致，都為正常>輕旱>濕潤(rùn)>干旱>重旱。3個(gè)季節(jié)中干旱等級(jí)主要以濕潤(rùn)和正常為主，面積分別為1 876.911、1 983.149、1 905.917 km2，均占流域面積的62%以上(表1、圖3)；全區(qū)3個(gè)季節(jié)均有輕旱分布，面積所占比例達(dá)30%以上，其中春季>秋季>冬季；干旱分布較少，具體分布情況為秋季>冬季>春季，研究區(qū)幾乎無(wú)重旱。而研究區(qū)3個(gè)季節(jié)的旱情在空間分布上存在一定差異，全區(qū)3個(gè)季節(jié)輕旱分布較廣，分布區(qū)域基本一致；秋季干旱主要分布在流域中部地勢(shì)較為平坦的耕作區(qū)域，南部有較少分布；在冬季和春季，干旱則主要分布在流域中惠水縣城至花溪大學(xué)城一帶，該地區(qū)城鎮(zhèn)化程度高，周圍有大量裸地，可能是導(dǎo)致其干旱的原因。
3.3 土壤含水量反演
將GPS記錄實(shí)測(cè)土壤水分采樣點(diǎn)經(jīng)緯度信息導(dǎo)入ArcGIS，提取采樣點(diǎn)對(duì)應(yīng)2017年4月1日TVDI值，與實(shí)測(cè)土壤體積含水量進(jìn)行最小二乘法回歸分析，得到兩者關(guān)系(圖4)及回歸方程：y= -149.47TVDI+97.255(R2= 0.7572)，且擬合結(jié)果通過(guò)t檢驗(yàn)，表明兩者之間具有良好的線性關(guān)系。把該線性模型輸入ENVI 5.1軟件中進(jìn)行反演，即可得到土壤水分分布(圖5)。由圖3、5可以看出，土壤水分與干旱情況分布基本一致，且土壤水分與TVDI呈現(xiàn)出明顯負(fù)相關(guān)關(guān)系，土壤水分隨著TVDI的增大而逐漸減小，土壤含水量越低，干旱情況則越嚴(yán)重，這說(shuō)明了TVDI能夠較好地反映研究區(qū)地表土壤水分狀況。
表1 土壤干旱等級(jí)、面積及比例
a) 秋季
b) 冬季
c) 春季圖3 不同季節(jié)土壤干旱等級(jí)空間分布
圖4 土壤表層含水量和TVDI關(guān)系
a) 秋季圖5 不同季節(jié)下土壤水分含量分布
b) 冬季
c) 春季續(xù)圖5 不同季節(jié)下土壤水分含量分布
3.4 不同坡度土壤水分季節(jié)變化分析
基于DEM高程數(shù)據(jù)，利用ArcGIS空間分析得到研究區(qū)的坡度因子，根據(jù)1984年中國(guó)農(nóng)業(yè)區(qū)劃委員會(huì)頒發(fā)的《土地利用現(xiàn)狀調(diào)查技術(shù)規(guī)程》把坡度分為五級(jí)，即Ⅰ(≤2°)、Ⅱ(2°～6°)、Ⅲ(6°～15°)、Ⅳ(15°～25°)、V(>25°)，通過(guò)分區(qū)統(tǒng)計(jì)，得到不同坡度等級(jí)下土壤水分柱狀圖(圖6)。從圖6可知，冬季和春季土壤水分均值隨著坡度增大而增加，秋季則是先隨著坡度增加，再到25°以上陡坡區(qū)域有所降低。這是由于秋季降水較少，而地表溫度較高，且在坡度較陡、植被覆蓋較差或地表裸露率高的區(qū)域，土壤水分相對(duì)較低。在第Ⅰ、Ⅱ、Ⅲ坡度等級(jí)下，3個(gè)季節(jié)的土壤水分均值中，冬季相對(duì)秋季和春季較高，并在第Ⅰ、Ⅱ坡度等級(jí)差異較為明顯，分析其原因可能是該坡度等級(jí)下區(qū)域主要為耕地和建設(shè)用地，秋季和春季地表溫度相對(duì)冬季較高，土壤水分則較少。
圖6 不同季節(jié)坡度等級(jí)下土壤水分含量
3.5 不同坡向土壤水分季節(jié)變化分析
通過(guò)研究區(qū)DEM數(shù)據(jù)進(jìn)行坡向提取，選取0°為正北方向，沿順時(shí)針旋轉(zhuǎn)，將坡向劃分為9種類型，即FLAT(-1)、N(0°～22.5°，337.5°～360°)、NE(22.5°～67.5°)、E(67.5°～112.5°)、SE(112.5°～157.5°)、S(157.5°～202.5°)、SW(202.5°～247.5°)、W(247.5°～292.5°)、NW(292.5°～337.5°)。通過(guò)分區(qū)統(tǒng)計(jì)得到不同坡向等級(jí)土壤水分折線(圖7)。從圖7可以看出，坡向?qū)?個(gè)季節(jié)的土壤水分影響基本一致，土壤水分均值在不同坡向上存在著差異性變化規(guī)律。土壤水分均值按坡向等級(jí)劃分，從左到右先增加再減少，然后又先減少再增大，其在北坡和西北坡相對(duì)較高，出現(xiàn)峰值，而在無(wú)坡向和東南坡相對(duì)較低，出現(xiàn)低谷，可以看出研究區(qū)土壤水分含量陰坡比陽(yáng)坡和無(wú)坡向多。相同坡向土壤水分均值相對(duì)冬季來(lái)看，春季與秋季差值相對(duì)較小，冬季則較大，產(chǎn)生此現(xiàn)象原因可能是冬季氣溫較低，地表蒸散發(fā)較?。磺锛倦m植被覆蓋較好，但地表溫度較高，蒸散發(fā)相對(duì)較大；春季與冬季相比，春季地表溫度上升，蒸散發(fā)也隨之增大。
圖7 不同季節(jié)下坡向等級(jí)土壤水分含量
3.6 不同土地利用類型土壤水分季節(jié)變化分析
利用ENVI 5.1軟件對(duì)2017年4月的Landsat8 OLI_TIRS數(shù)據(jù)進(jìn)行假彩色合成，然后采取監(jiān)督分類法對(duì)研究區(qū)土地利用類型進(jìn)行分類，主要分為耕地、林地、灌木林地、草地、建設(shè)用地、水域和其他用地7個(gè)類型(圖8)。選取前5個(gè)土地利用類型對(duì)3個(gè)季節(jié)土壤水分進(jìn)行統(tǒng)計(jì)分析，得到不同土地利用類型土壤水分均值折線圖(圖9)。從圖9可以看出，林地在3個(gè)季節(jié)中的土壤水分含量最高，灌木林地次之，其他3個(gè)地類則存在差異。而3個(gè)季節(jié)中不同土地利用類型的土壤水分分布在春季和秋季是一致的，主要表現(xiàn)為林地>灌木林地>草地>耕地>其他用地，冬季則為林地>灌木林地>耕地>草地>其他用地。5個(gè)土地利用類型中，林地、灌木林地和草地在3個(gè)季節(jié)土壤水分差異較小，其他用地和耕地差異則較大，以冬季>春季>秋季。其中，耕地在秋季地表溫度相對(duì)春季和冬季較高，而秋季又是農(nóng)作物豐收時(shí)期，地表植被覆蓋減少，導(dǎo)致蒸散發(fā)增大，所以表層土壤水分最低；冬季則是一年中溫度最低季節(jié)，地表蒸散發(fā)較少，所以土壤水分含量可能變高；春季是耕作季節(jié)，表層土壤被翻轉(zhuǎn)，加之溫度上升，地表蒸散發(fā)相對(duì)較大，所以土壤水分比冬季低。
圖8 土地利用類型
圖9 不同季節(jié)下土地類型土壤水分含量
4 結(jié)論
利用Landsat8 OLI_TIRS數(shù)據(jù)，結(jié)合溫度植被干旱指數(shù)TVDI，對(duì)貴州省喀斯特漣江流域土壤水分反演評(píng)價(jià)，得出以下結(jié)論。
a) 研究區(qū)內(nèi)3個(gè)季節(jié)旱情分布較廣，以輕旱旱情為主，面積所占比例達(dá)30%以上，春季>秋季>冬季。干旱次之，幾乎無(wú)重旱。該區(qū)旱情主要發(fā)生在流域北部、中部和南部地勢(shì)較為平坦的耕作區(qū)域和城鎮(zhèn)周圍，TVDI適用于流域的土壤水分含量反演及干早監(jiān)測(cè)評(píng)價(jià)。
b) 冬季和春季土壤水分均值隨著坡度增大而增加，秋季則是先隨著坡度增加，再到25°以上陡坡區(qū)域減少。相同坡度條件下，冬季土壤水分總體比秋季和春季多，春季與秋季差值相對(duì)較小，冬季則較大。
c) 土壤水分均值在不同坡向有較大差異，由無(wú)坡向、北坡按順時(shí)針旋轉(zhuǎn)至西北坡，土壤水分先增加再減少，然后又先減少再曾加。總體上同一坡向條件下土壤水分含量陰坡比陽(yáng)坡和無(wú)坡向多，冬季比秋季和春季多，春季與秋季的差值相對(duì)較小，冬季則較大。
d) 在耕地、林地、灌木林地、草地、其他用地五個(gè)地類中，林地在不同季節(jié)土壤水分含量最高，灌木林地次之，其他3個(gè)地類則存在差異。相同地類在不同季節(jié)、林地、灌木林地和草地的土壤水分差異較小，其他地類和耕地的差異則較大，冬季>春季>秋季。
影響植被指數(shù)與地表溫度空間特征關(guān)系的因素還有很多，而在貴州喀斯特山地地區(qū)，海拔落差、云霧、空氣水分含量、降雨及研究區(qū)域大小和地表覆蓋物類型都會(huì)對(duì)植被指數(shù)與地表溫度空間特征關(guān)系產(chǎn)生影響。若在該經(jīng)驗(yàn)?zāi)Ｐ椭胁荒芫C合考慮這些因素，會(huì)引起一定偏差。
猜你喜歡
FMS與YBT相關(guān)性的實(shí)證研究體育科技文獻(xiàn)通報(bào)(2022年3期)2022-05-23 13:46:54 永久基本農(nóng)田集中區(qū)域“禁廢”今日農(nóng)業(yè)(2021年9期)2021-11-26 07:41:24 2020年國(guó)內(nèi)翻譯研究述評(píng)天津外國(guó)語(yǔ)大學(xué)學(xué)報(bào)(2021年3期)2021-08-13 08:32:18 遼代千人邑研究述論遼金歷史與考古(2021年0期)2021-07-29 01:06:54 分割區(qū)域發(fā)明與創(chuàng)新·小學(xué)生(2021年3期)2021-03-25 11:48:49 視錯(cuò)覺(jué)在平面設(shè)計(jì)中的應(yīng)用與研究科技傳播(2019年22期)2020-01-14 03:06:54 EMA伺服控制系統(tǒng)研究民用飛機(jī)設(shè)計(jì)與研究(2019年4期)2019-05-21 07:21:24 新版C-NCAP側(cè)面碰撞假人損傷研究汽車工程學(xué)報(bào)(2017年2期)2017-07-05 08:13:02 關(guān)于四色猜想中國(guó)科技博覽(2016年2期)2016-04-25 20:32:39 分區(qū)域小學(xué)生導(dǎo)刊(2016年34期)2016-04-11 00:49:44
梁建方，周秋文，韋小茶，馬龍生，顏紅
(貴州師范大學(xué)地理與環(huán)境科學(xué)學(xué)院，貴州貴陽(yáng)550001)
土壤是重要的自然資源，土壤水分作為陸地和大氣能量交換過(guò)程中的重要因子，對(duì)陸表蒸散、水的運(yùn)動(dòng)和碳循環(huán)具有控制作用，對(duì)自然環(huán)境和人類社會(huì)的可持續(xù)發(fā)展有著重要影響[1-3]。獲取區(qū)域乃至全球尺度的土壤水分空間分布信息，是研究全球范圍內(nèi)水循環(huán)、能量平衡和氣象預(yù)報(bào)的先決條件，對(duì)改善區(qū)域及全球氣候模式預(yù)報(bào)結(jié)果、進(jìn)行農(nóng)作物旱情監(jiān)測(cè)及估產(chǎn)、自然和生態(tài)環(huán)境問(wèn)題的研究起重要作用[4-5]。傳統(tǒng)土壤水分測(cè)量方法只能得到有限測(cè)量點(diǎn)數(shù)據(jù)，僅能反應(yīng)測(cè)量點(diǎn)周圍有限區(qū)域的土壤水分情況。而通過(guò)星載遙感衛(wèi)星傳感器提供的監(jiān)測(cè)手段，能獲取大區(qū)域的土壤水分信息，同時(shí)具有時(shí)效快、動(dòng)態(tài)性較強(qiáng)及時(shí)間空間分辨率高等優(yōu)點(diǎn)，能為區(qū)域尺度土壤水分信息的獲取提供良好手段[6-8]。
光學(xué)遙感監(jiān)測(cè)土壤水分的方法主要有熱慣量法、植被指數(shù)法等，但熱慣量法僅適用于裸地及低值被覆蓋區(qū)域，而植被指數(shù)法僅適用于植被覆蓋較高區(qū)域[1,8]。Price和Carlson研究發(fā)現(xiàn)植被指數(shù)和地表溫度之間存在負(fù)相關(guān)關(guān)系，且以植被指數(shù)和地表溫度為橫縱坐標(biāo)生成的散點(diǎn)圖呈三角形狀[9-10]。Sandholt等[11]進(jìn)一步研究植被指數(shù)與地表溫度的關(guān)系，提出利用溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index，TVDI)來(lái)監(jiān)測(cè)和估算土壤水分的方法。國(guó)內(nèi)外學(xué)者基于該方法進(jìn)行區(qū)域地表水分狀況等開(kāi)展大量研究，并取得了豐碩成果[12]。王秀君等[13]基于Landsat TM數(shù)據(jù)反演了河北省黃驊市土壤水分，結(jié)果表明TVDI與土壤濕度顯著相關(guān)；曹雷等[14]利用TVDI在艾比湖地區(qū)進(jìn)行試驗(yàn)，發(fā)現(xiàn)TVDI能有效反演區(qū)域土壤水分，且精度較高；陳斌等[15]以內(nèi)蒙古錫林郭勒盟地區(qū)為研究對(duì)象，發(fā)現(xiàn)TVDI可以用來(lái)監(jiān)測(cè)旱情變化及評(píng)價(jià)草原干旱狀況；沙莎等[16]利用歷史遙感數(shù)據(jù)建立的特征空間，發(fā)現(xiàn)TVDI能夠較好地指示每年土壤相對(duì)濕度RSM的時(shí)空分布特征；吳孟泉、張喆等[17-18]研究均發(fā)現(xiàn)TVDI能夠很好地用于山區(qū)干旱預(yù)警與監(jiān)測(cè)。
喀斯特地區(qū)生態(tài)系統(tǒng)敏感度高, 災(zāi)變承受能力低, 環(huán)境容量小, 因而成為典型的生態(tài)脆弱區(qū)[19]。其土層淺薄、土壤貧瘠、成土速率低、水文過(guò)程變化迅速、地表土壤水受特殊地質(zhì)地貌及水文條件等自然因素影響較大[20-24]。以貴州省為中心的中國(guó)西南喀斯特地區(qū)是世界上面積最大的喀斯特連續(xù)帶，該地區(qū)地表破碎、地形復(fù)雜，由于森林砍伐、過(guò)度放牧、農(nóng)業(yè)耕種和地下水過(guò)度開(kāi)采等行為，地表植被的生長(zhǎng)和人類生產(chǎn)生活深受土壤水分貧瘠影響[25-30]。郭巍等[26]采用TVDI對(duì)云南省2010年春季旱情進(jìn)行監(jiān)測(cè)，揭示了云南省入春以來(lái)的旱情動(dòng)態(tài)變化及旱情分布特點(diǎn)；王海等[27]運(yùn)用云南省境內(nèi)的TVDI與昆明3個(gè)站點(diǎn)的土壤濕度做線性相關(guān)分析表明，TVDI對(duì)云南境內(nèi)的干旱監(jiān)測(cè)具有較好的適用性和可靠性；康為民等[28]利用同一時(shí)期TVDI與當(dāng)?shù)貧庀笳军c(diǎn)的土壤濕度數(shù)據(jù)進(jìn)行定量驗(yàn)證，發(fā)現(xiàn)TVDI與土壤濕度顯著相關(guān)，適用于大區(qū)域復(fù)雜地形的干旱檢測(cè)與預(yù)警。
然而已有研究大都基于MODIS數(shù)據(jù)與Landsat TM/ETM+數(shù)據(jù)來(lái)進(jìn)行干旱監(jiān)測(cè)研究，而利用Landsat8 OLI_TIRS數(shù)據(jù)進(jìn)行研究干旱監(jiān)測(cè)研究較少。喀斯特地貌的典型特征是地表破碎、地形復(fù)雜，土壤水分分布情況隨降雨及地表情況的不同而表現(xiàn)出不同區(qū)域特征。因此，本文基于Landsat8 OLI_TIRS數(shù)據(jù)，選取貴州省典型喀斯特地區(qū)的漣江流域?yàn)槔瑢?duì)流域的土壤水分進(jìn)行遙感反演，并實(shí)現(xiàn)干旱時(shí)空分布情況的動(dòng)態(tài)監(jiān)測(cè)，同時(shí)結(jié)合該區(qū)野外實(shí)測(cè)土壤含水量數(shù)據(jù)進(jìn)行驗(yàn)證分析。研究結(jié)果可為相關(guān)部門制定合理有效的抗旱、防旱措施提供科學(xué)依據(jù)，對(duì)喀斯特地區(qū)水土保持政策的制定、石漠化治理措施實(shí)施具有非常重要的指導(dǎo)作用，是喀斯特地區(qū)自然、經(jīng)濟(jì)、社會(huì)可持續(xù)發(fā)展的重要需求。
1 研究區(qū)概況
漣江流域(圖1)位于貴州省西南部(106°14′～106°51′E，25°32′～26°26′N)，面積約3 006.67 km2，流域海拔高度在517～1 650 m之間，地勢(shì)北高南低，中部平緩，其境內(nèi)主要分布碳酸鹽巖，喀斯特地貌典型發(fā)育。屬亞熱帶季風(fēng)氣候，降水豐富，但降水量年際、季節(jié)性分布不均，而地表和地下水系構(gòu)成的二元結(jié)構(gòu)往往會(huì)導(dǎo)致降水滲漏，土壤儲(chǔ)水能力低，加上不合理的人類活動(dòng)造成植被不斷遭受破壞，全區(qū)一年四季均存在輕重程度不相同的旱情。
圖1 漣江流域地理位置及高程
2 數(shù)據(jù)來(lái)源與方法
2.1 數(shù)據(jù)來(lái)源及預(yù)處理
遙感數(shù)據(jù)來(lái)源于中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心地理空間數(shù)據(jù)云平臺(tái)(http://www.gscloud.cn)，獲取2013年9月29日、2014年1月19日與2017年4月1日的Landsat8 OLI_TIRS數(shù)據(jù)，以及ASTER GDEM高程數(shù)據(jù)，空間分辨率為30 m。利用ENVI 5.1軟件對(duì)遙感影像進(jìn)行預(yù)處理，主要包括輻射定標(biāo)、裁剪和大氣校正等。選取漣江子流域(揚(yáng)眉河流域)中具有代表性的林地、灌木林地、草地、耕地，于2017年4月2日使用TDR-300土壤水分儀測(cè)量得到24個(gè)樣點(diǎn)表層(0～10 cm)土壤體積含水量，并用手持GPS獲取采樣點(diǎn)的經(jīng)緯度信息。
2.2 溫度植被干旱指數(shù)估算
Price和Carlson等研究發(fā)現(xiàn)以NDVI與Ts為橫縱坐標(biāo)生成散點(diǎn)圖構(gòu)成的空間關(guān)系為三角形關(guān)系。Sandholt等對(duì)NDVI-LST特征空間進(jìn)行簡(jiǎn)化，提出溫度植被干旱指數(shù)(Temperature Vegetation Dryness Index，TVDI)概念，計(jì)算公式如下：
(1)
式中TVDI——溫度植被干旱指數(shù)；Ts——地表真實(shí)溫度；Tsmax——某一NDVI對(duì)應(yīng)最大地表溫度，即干邊；Tsmin——地表最低溫度，即濕邊。干邊與濕邊方程分別為：Tsmax=a·NDVI+b,Tsmin=c·NDVI+d，其中a、b、c、d為干、濕邊擬合方程系數(shù)。
模型中地表真實(shí)溫度Ts由輻射傳導(dǎo)方程法獲得，具體步驟如下。
a) 計(jì)算植被指數(shù)NDVI。NDVI由遙感影像中，近紅外波段的反射值與紅波段的反射值之差比上兩者之和，其計(jì)算公式為：
(2)
式中NIR——近紅外波段像元亮度值；R——紅波段像元亮度值。
b) 計(jì)算地表比輻射率ε。本文采用Van經(jīng)驗(yàn)公式[19]計(jì)算地表比輻射率：
ε=1.0094 + 0.047·ln·(NDVI)
(3)
式中，NDVI值在0.157～0.727之間，NDVI小于0的像元只要是水體，其體表比輻射率近似0.995；城鎮(zhèn)很多建筑物像元NDVI位于0～0.157之間，其地表比輻射率可設(shè)置為0.923；而對(duì)于NDVI大于0.727像元可看作植被完全覆蓋，地表比輻射率為0.986。
c) 計(jì)算相同溫度下黑體的輻射亮度值B(TS)。衛(wèi)星傳感器接收到的熱紅外輻射亮度值Lλ由三部分組成：大氣向上輻射亮度L↑；地面的真實(shí)輻射亮度經(jīng)過(guò)大氣層之后到達(dá)衛(wèi)星傳感器的能量；大氣向下輻射到達(dá)地面后反射的能量。其輻射傳輸方程為：
Lλ=[ε·B(TS) + (1-ε)L↓]·τ+L↑
(4)
根據(jù)普朗克定律得到的黑體熱輻射亮度B(TS)，其公式為：
B(TS) = [Lλ-L↑-τ·(1-ε)L↓]/(τ·ε)
(5)
式中Lλ——TIRS影像第十波段大氣校正后的輻射定標(biāo)值；ε——地表輻射率；TS——地表真實(shí)溫度；τ——大氣在熱紅外波段透過(guò)率。
d) 計(jì)算地表真實(shí)溫度TS。在獲取B(TS)后，根據(jù)普朗克公式的反函數(shù)，求得地表真實(shí)溫度TS，其公式為：
TS=K2/ln(K1/B(TS)+ 1)
(6)
式中B(TS)——溫度為TS黑體熱輻射亮度；K1、K2——為傳感器定標(biāo)常數(shù)。
2.3 表層土壤水計(jì)算
利用野外表層土壤水分(0～10 cm)采樣數(shù)據(jù)，建立土壤水分與溫度植被干旱指數(shù)之間的線性相關(guān)模型：
W=a·TVDI+b
(7)
式中W——表層土壤水分；TVDI——溫度植被干旱指數(shù)；a、b——線性系數(shù)。
3 結(jié)果與分析
3.1 干、濕邊擬合
由于NDVI值小于0的像元主要是水體和云，濕度接近100%，NDVI值在0～0.157時(shí)的像元主要是城市建筑物與水泥地等。為避免水體、云及建筑物對(duì)特征空間的影響，在干濕邊擬合時(shí)不考慮此類像元。因此，選取NDVI值大于0.157的像元，利用IDL編程實(shí)現(xiàn)相同NDVI值對(duì)應(yīng)的TS的最大值和最小值提取，步長(zhǎng)值為0.002，構(gòu)成NDVI-TS特征空間的形狀近似呈三角形(圖2)。從3個(gè)季節(jié)的擬合結(jié)果可以看出，干邊與濕邊的決定系數(shù)R2分別在0.9和0.5以上，擬合精度較高。且NDVI與地表溫度LST近似呈線性關(guān)系，隨著NDVI值的增大，地表溫度最大值LSTmax逐漸減小，最小值LSTmin逐漸增大。
a) 秋季
b) 冬季
c) 春季圖2 Ts-NDVI特征空間
3.2 TVDI干旱等級(jí)劃分
為從宏觀了解漣江流域干旱情況，根據(jù)已擬合的干、濕邊方程，計(jì)算研究區(qū)內(nèi)各個(gè)像元TVDI，根據(jù)TVDI灰度值范圍進(jìn)行干旱等級(jí)劃分指標(biāo)[21]，將干旱情況劃分為五級(jí)：濕潤(rùn)(0～0.2)、正常(0.2～0.4)、輕旱(0.4～0.6)、干旱(0.6～0.8)、重旱(0.8～1.0)。研究結(jié)果表明：研究區(qū)3個(gè)季節(jié)干旱情況在時(shí)間分布上總體一致，都為正常>輕旱>濕潤(rùn)>干旱>重旱。3個(gè)季節(jié)中干旱等級(jí)主要以濕潤(rùn)和正常為主，面積分別為1 876.911、1 983.149、1 905.917 km2，均占流域面積的62%以上(表1、圖3)；全區(qū)3個(gè)季節(jié)均有輕旱分布，面積所占比例達(dá)30%以上，其中春季>秋季>冬季；干旱分布較少，具體分布情況為秋季>冬季>春季，研究區(qū)幾乎無(wú)重旱。而研究區(qū)3個(gè)季節(jié)的旱情在空間分布上存在一定差異，全區(qū)3個(gè)季節(jié)輕旱分布較廣，分布區(qū)域基本一致；秋季干旱主要分布在流域中部地勢(shì)較為平坦的耕作區(qū)域，南部有較少分布；在冬季和春季，干旱則主要分布在流域中惠水縣城至花溪大學(xué)城一帶，該地區(qū)城鎮(zhèn)化程度高，周圍有大量裸地，可能是導(dǎo)致其干旱的原因。
3.3 土壤含水量反演
將GPS記錄實(shí)測(cè)土壤水分采樣點(diǎn)經(jīng)緯度信息導(dǎo)入ArcGIS，提取采樣點(diǎn)對(duì)應(yīng)2017年4月1日TVDI值，與實(shí)測(cè)土壤體積含水量進(jìn)行最小二乘法回歸分析，得到兩者關(guān)系(圖4)及回歸方程：y= -149.47TVDI+97.255(R2= 0.7572)，且擬合結(jié)果通過(guò)t檢驗(yàn)，表明兩者之間具有良好的線性關(guān)系。把該線性模型輸入ENVI 5.1軟件中進(jìn)行反演，即可得到土壤水分分布(圖5)。由圖3、5可以看出，土壤水分與干旱情況分布基本一致，且土壤水分與TVDI呈現(xiàn)出明顯負(fù)相關(guān)關(guān)系，土壤水分隨著TVDI的增大而逐漸減小，土壤含水量越低，干旱情況則越嚴(yán)重，這說(shuō)明了TVDI能夠較好地反映研究區(qū)地表土壤水分狀況。
表1 土壤干旱等級(jí)、面積及比例
a) 秋季
b) 冬季
c) 春季圖3 不同季節(jié)土壤干旱等級(jí)空間分布
圖4 土壤表層含水量和TVDI關(guān)系
a) 秋季圖5 不同季節(jié)下土壤水分含量分布
b) 冬季
c) 春季續(xù)圖5 不同季節(jié)下土壤水分含量分布
3.4 不同坡度土壤水分季節(jié)變化分析
基于DEM高程數(shù)據(jù)，利用ArcGIS空間分析得到研究區(qū)的坡度因子，根據(jù)1984年中國(guó)農(nóng)業(yè)區(qū)劃委員會(huì)頒發(fā)的《土地利用現(xiàn)狀調(diào)查技術(shù)規(guī)程》把坡度分為五級(jí)，即Ⅰ(≤2°)、Ⅱ(2°～6°)、Ⅲ(6°～15°)、Ⅳ(15°～25°)、V(>25°)，通過(guò)分區(qū)統(tǒng)計(jì)，得到不同坡度等級(jí)下土壤水分柱狀圖(圖6)。從圖6可知，冬季和春季土壤水分均值隨著坡度增大而增加，秋季則是先隨著坡度增加，再到25°以上陡坡區(qū)域有所降低。這是由于秋季降水較少，而地表溫度較高，且在坡度較陡、植被覆蓋較差或地表裸露率高的區(qū)域，土壤水分相對(duì)較低。在第Ⅰ、Ⅱ、Ⅲ坡度等級(jí)下，3個(gè)季節(jié)的土壤水分均值中，冬季相對(duì)秋季和春季較高，并在第Ⅰ、Ⅱ坡度等級(jí)差異較為明顯，分析其原因可能是該坡度等級(jí)下區(qū)域主要為耕地和建設(shè)用地，秋季和春季地表溫度相對(duì)冬季較高，土壤水分則較少。
圖6 不同季節(jié)坡度等級(jí)下土壤水分含量
3.5 不同坡向土壤水分季節(jié)變化分析
通過(guò)研究區(qū)DEM數(shù)據(jù)進(jìn)行坡向提取，選取0°為正北方向，沿順時(shí)針旋轉(zhuǎn)，將坡向劃分為9種類型，即FLAT(-1)、N(0°～22.5°，337.5°～360°)、NE(22.5°～67.5°)、E(67.5°～112.5°)、SE(112.5°～157.5°)、S(157.5°～202.5°)、SW(202.5°～247.5°)、W(247.5°～292.5°)、NW(292.5°～337.5°)。通過(guò)分區(qū)統(tǒng)計(jì)得到不同坡向等級(jí)土壤水分折線(圖7)。從圖7可以看出，坡向?qū)?個(gè)季節(jié)的土壤水分影響基本一致，土壤水分均值在不同坡向上存在著差異性變化規(guī)律。土壤水分均值按坡向等級(jí)劃分，從左到右先增加再減少，然后又先減少再增大，其在北坡和西北坡相對(duì)較高，出現(xiàn)峰值，而在無(wú)坡向和東南坡相對(duì)較低，出現(xiàn)低谷，可以看出研究區(qū)土壤水分含量陰坡比陽(yáng)坡和無(wú)坡向多。相同坡向土壤水分均值相對(duì)冬季來(lái)看，春季與秋季差值相對(duì)較小，冬季則較大，產(chǎn)生此現(xiàn)象原因可能是冬季氣溫較低，地表蒸散發(fā)較?。磺锛倦m植被覆蓋較好，但地表溫度較高，蒸散發(fā)相對(duì)較大；春季與冬季相比，春季地表溫度上升，蒸散發(fā)也隨之增大。
圖7 不同季節(jié)下坡向等級(jí)土壤水分含量
3.6 不同土地利用類型土壤水分季節(jié)變化分析
利用ENVI 5.1軟件對(duì)2017年4月的Landsat8 OLI_TIRS數(shù)據(jù)進(jìn)行假彩色合成，然后采取監(jiān)督分類法對(duì)研究區(qū)土地利用類型進(jìn)行分類，主要分為耕地、林地、灌木林地、草地、建設(shè)用地、水域和其他用地7個(gè)類型(圖8)。選取前5個(gè)土地利用類型對(duì)3個(gè)季節(jié)土壤水分進(jìn)行統(tǒng)計(jì)分析，得到不同土地利用類型土壤水分均值折線圖(圖9)。從圖9可以看出，林地在3個(gè)季節(jié)中的土壤水分含量最高，灌木林地次之，其他3個(gè)地類則存在差異。而3個(gè)季節(jié)中不同土地利用類型的土壤水分分布在春季和秋季是一致的，主要表現(xiàn)為林地>灌木林地>草地>耕地>其他用地，冬季則為林地>灌木林地>耕地>草地>其他用地。5個(gè)土地利用類型中，林地、灌木林地和草地在3個(gè)季節(jié)土壤水分差異較小，其他用地和耕地差異則較大，以冬季>春季>秋季。其中，耕地在秋季地表溫度相對(duì)春季和冬季較高，而秋季又是農(nóng)作物豐收時(shí)期，地表植被覆蓋減少，導(dǎo)致蒸散發(fā)增大，所以表層土壤水分最低；冬季則是一年中溫度最低季節(jié)，地表蒸散發(fā)較少，所以土壤水分含量可能變高；春季是耕作季節(jié)，表層土壤被翻轉(zhuǎn)，加之溫度上升，地表蒸散發(fā)相對(duì)較大，所以土壤水分比冬季低。
圖8 土地利用類型
圖9 不同季節(jié)下土地類型土壤水分含量
4 結(jié)論
利用Landsat8 OLI_TIRS數(shù)據(jù)，結(jié)合溫度植被干旱指數(shù)TVDI，對(duì)貴州省喀斯特漣江流域土壤水分反演評(píng)價(jià)，得出以下結(jié)論。
a) 研究區(qū)內(nèi)3個(gè)季節(jié)旱情分布較廣，以輕旱旱情為主，面積所占比例達(dá)30%以上，春季>秋季>冬季。干旱次之，幾乎無(wú)重旱。該區(qū)旱情主要發(fā)生在流域北部、中部和南部地勢(shì)較為平坦的耕作區(qū)域和城鎮(zhèn)周圍，TVDI適用于流域的土壤水分含量反演及干早監(jiān)測(cè)評(píng)價(jià)。
b) 冬季和春季土壤水分均值隨著坡度增大而增加，秋季則是先隨著坡度增加，再到25°以上陡坡區(qū)域減少。相同坡度條件下，冬季土壤水分總體比秋季和春季多，春季與秋季差值相對(duì)較小，冬季則較大。
c) 土壤水分均值在不同坡向有較大差異，由無(wú)坡向、北坡按順時(shí)針旋轉(zhuǎn)至西北坡，土壤水分先增加再減少，然后又先減少再曾加。總體上同一坡向條件下土壤水分含量陰坡比陽(yáng)坡和無(wú)坡向多，冬季比秋季和春季多，春季與秋季的差值相對(duì)較小，冬季則較大。
d) 在耕地、林地、灌木林地、草地、其他用地五個(gè)地類中，林地在不同季節(jié)土壤水分含量最高，灌木林地次之，其他3個(gè)地類則存在差異。相同地類在不同季節(jié)、林地、灌木林地和草地的土壤水分差異較小，其他地類和耕地的差異則較大，冬季>春季>秋季。
影響植被指數(shù)與地表溫度空間特征關(guān)系的因素還有很多，而在貴州喀斯特山地地區(qū)，海拔落差、云霧、空氣水分含量、降雨及研究區(qū)域大小和地表覆蓋物類型都會(huì)對(duì)植被指數(shù)與地表溫度空間特征關(guān)系產(chǎn)生影響。若在該經(jīng)驗(yàn)?zāi)Ｐ椭胁荒芫C合考慮這些因素，會(huì)引起一定偏差。
猜你喜歡
體育科技文獻(xiàn)通報(bào)(2022年3期)2022-05-23 13:46:54
今日農(nóng)業(yè)(2021年9期)2021-11-26 07:41:24
天津外國(guó)語(yǔ)大學(xué)學(xué)報(bào)(2021年3期)2021-08-13 08:32:18
遼金歷史與考古(2021年0期)2021-07-29 01:06:54
發(fā)明與創(chuàng)新·小學(xué)生(2021年3期)2021-03-25 11:48:49
科技傳播(2019年22期)2020-01-14 03:06:54
民用飛機(jī)設(shè)計(jì)與研究(2019年4期)2019-05-21 07:21:24
汽車工程學(xué)報(bào)(2017年2期)2017-07-05 08:13:02
中國(guó)科技博覽(2016年2期)2016-04-25 20:32:39
小學(xué)生導(dǎo)刊(2016年34期)2016-04-11 00:49:44
