孫時(shí)雨，宋承運(yùn)，周 露

(1.安徽理工大學(xué) 空間信息與測繪工程學(xué)院，安徽 淮南 232001；2.安徽理工大學(xué) 礦山采動(dòng)災(zāi)害空天地協(xié)同監(jiān)測與預(yù)警安徽普通高校重點(diǎn)實(shí)驗(yàn)室，安徽 淮南 232001)

0 引言

土壤水分是能量、水和碳循環(huán)的重要組成部分,為陸地-大氣模型、水文模型和降水預(yù)測模型提供了重要信息[1]。在多尺度相互作用中,土壤水分在控制地表能量分配、調(diào)節(jié)地表徑流和土壤排水、調(diào)節(jié)冠層蒸騰和碳排放等方面發(fā)揮著重要作用[2]。長時(shí)序的土壤水分產(chǎn)品可用于更準(zhǔn)確和可靠地估計(jì)深層土壤水分和蒸散發(fā)等,也可用于數(shù)據(jù)同化、天氣預(yù)報(bào)和干旱監(jiān)測等。因此,迫切需要可靠、準(zhǔn)確和完整的土壤水分信息[3]。

由于土壤水分空間異質(zhì)性強(qiáng),用來觀測土壤水分的設(shè)備成本高,不僅耗費(fèi)人力,還耗費(fèi)物力,因此很難建立區(qū)域和全球范圍的高分辨率觀測網(wǎng)絡(luò)。傳統(tǒng)的土壤水分觀測方法一般只適合范圍比較小的土壤水分的信息采集,無法應(yīng)用到大范圍的地表水分的測量。隨著遙感技術(shù)的快速發(fā)展,衛(wèi)星遙感是獲取地表土壤水分(通常在0～5 cm深度)最有效的工具[4]。特別是微波遙感被認(rèn)為是監(jiān)測土壤水分最主流的方法,因?yàn)樗軌蛉旃ぷ?不受天氣的限制,可穿透云層、植被和土壤[5]。各種衛(wèi)星研究和應(yīng)用機(jī)構(gòu)的微波傳感器可以生成全球遙感產(chǎn)品,如高級微波掃描輻射計(jì)EOS(AMSR-E)、土壤水分主-被動(dòng)遙感衛(wèi)星(SMAP)、歐洲航天局(ESA)的土壤濕度和海洋鹽度觀測衛(wèi)星(SMOS)、日本的先進(jìn)微波掃描輻射計(jì)2(AMSR2)和中國的微波輻射成像儀(MWRI)[6]。然而上述衛(wèi)星產(chǎn)品的時(shí)間覆蓋范圍太短,無法滿足水文氣象的應(yīng)用需求。因此, ESA CCI發(fā)布了一個(gè)完整和一致的全球土壤水分?jǐn)?shù)據(jù)集。但由于衛(wèi)星軌道變化、射頻干擾(RFI)和微型傳感器的物理限制,該數(shù)據(jù)集存在數(shù)據(jù)缺失[7]。數(shù)據(jù)缺失限制了ESA CCI 土壤水分產(chǎn)品的實(shí)際應(yīng)用,特別是在需要空間完整和時(shí)間連續(xù)的土壤水分?jǐn)?shù)據(jù)的領(lǐng)域,如氣候模擬、干旱監(jiān)測和水資源管理[8]等。

國內(nèi)外已經(jīng)提出了一些方法來填補(bǔ)衛(wèi)星數(shù)據(jù)集的空白。這些方法主要可分為統(tǒng)計(jì)、插值和機(jī)器學(xué)習(xí)等幾類。統(tǒng)計(jì)方法(例如多元線性回歸、非線性回歸)主要依賴輸入和目標(biāo)數(shù)據(jù)的線性和非線性關(guān)系。Wang等[9]測試了一種基于三維離散余弦變換(DCT-PLS)的最小二乘回歸的統(tǒng)計(jì)建模方法來填補(bǔ)衛(wèi)星產(chǎn)品空白。Yan等[10]運(yùn)用了反距離加權(quán)(IDW)根據(jù)目標(biāo)變量的空間模式填充缺失數(shù)據(jù)。其他方法已用于重建,例如Turlapaty等[11]使用改進(jìn)的奇異光譜分析來填補(bǔ)AMSR-E土壤水分?jǐn)?shù)據(jù)集中的缺失。近年來,支持向量機(jī)(SVM)和人工神經(jīng)網(wǎng)絡(luò)(ANN)被廣泛用于填補(bǔ)遙感數(shù)據(jù)的缺失[12]。然而,基于機(jī)器學(xué)習(xí)算法模型的優(yōu)劣很大程度上取決于輸入變量的選擇以及輸入變量和目標(biāo)變量的關(guān)系,模型的結(jié)構(gòu)與輸入變量、目標(biāo)變量和周期高度相關(guān)[13]。收斂速度慢、易過擬合是傳統(tǒng)的機(jī)器學(xué)習(xí)模型及BP神經(jīng)網(wǎng)絡(luò)的缺點(diǎn),而且得到數(shù)據(jù)的精度誤差大[14]。

基于以上分析,以青藏高原那曲為研究區(qū)域,利用ESA CCI V04.2組合土壤水分產(chǎn)品、數(shù)字高程數(shù)據(jù)(DEM)以及基于MODIS數(shù)據(jù)重建的地表溫度數(shù)據(jù)、歸一化植被指數(shù)(NDVI)和反照率數(shù)據(jù),采用原子搜索優(yōu)化(Atom Search Optimization,ASO)算法去優(yōu)化BP神經(jīng)網(wǎng)絡(luò)模型的權(quán)值、閾值,從而建立ASO-BP土壤水分重建模型,得到空間連續(xù)土壤水分,并利用青藏高原中部土壤水分多尺度觀測網(wǎng)(Soil Moisture and Temperature Monitoring Network on the Central TP,TP-SMTMN)數(shù)據(jù)對結(jié)果進(jìn)行驗(yàn)證分析。

1 數(shù)據(jù)來源

1.1 研究區(qū)域概況

研究區(qū)域?yàn)榍嗖馗咴乔貐^(qū)(83°55′E～95°05′E,29°55′N～36°30′N)。那曲約有40萬km2的面積,地勢開闊,干濕季分明,屬于高原亞寒帶氣候區(qū),月平均降水量在7月最大,在12月或1月最小,年平均降水量為476.5 mm[15]。

那曲地區(qū)在3月下旬凍土開始消融,土壤水分增加。從 5月開始表層土壤水分會繼續(xù)增加,受季風(fēng)降水的影響,7月達(dá)到峰值,8月以后降水減小,土壤水分逐漸減小,10月溫度降低,土壤水分逐漸由液態(tài)變?yōu)楣虘B(tài),土壤水分減小,1月達(dá)到谷值,一直持續(xù)到 3月,完成一年的循環(huán)周期[16]。研究區(qū)高程及地面觀測站點(diǎn)示意如圖1所示。

(a)研究區(qū)高程

(b)地面觀測站點(diǎn)

1.2 遙感數(shù)據(jù)來源與處理

(1)ESA CCI土壤水分?jǐn)?shù)據(jù)。ESA CCI是一種長時(shí)間序列、多衛(wèi)星融合的土壤水分產(chǎn)品,它基于主動(dòng)和被動(dòng)微波傳感器,生產(chǎn)的數(shù)據(jù)集包含主動(dòng)、被動(dòng)和融合數(shù)據(jù)集3種,空間分辨率為0.25°,時(shí)間分辨率為1 d[17]。它也是第一個(gè)基于遙感的長時(shí)間序列土壤水分產(chǎn)品,具有很好的一致性,而衛(wèi)星軌道、植被、凍土、積雪和冰川等因素影響了反演后土壤水分的空間連續(xù)性與精度[18]。本文采用ECV V04.2組合土壤水分產(chǎn)品(每日,0.25°),作為主要數(shù)據(jù)源 (http:∥www.esa-soilmosture-cci.org/)[19]。

(2)LST、NDVI、Albedo與DEM遙感數(shù)據(jù)。LST、NDVI和Albedo數(shù)據(jù)是由MODIS MOD11C1 LST產(chǎn)品(0.05°)、MOD13C1 NDVI產(chǎn)品(0.05°)與MOD43C1反照率產(chǎn)品(0.05°)分別采用多時(shí)相穩(wěn)健回歸方法[20]和基于時(shí)間濾波統(tǒng)計(jì)方法[21]進(jìn)行數(shù)據(jù)重建,得到空間連續(xù)的LST、NDVI和Albedo數(shù)據(jù),分辨率為0.05°(數(shù)據(jù)來自 http:∥reverb.echo.nasa.gov/)。美國國家航空航天局(National Aeronautics and Space Administration,NASA)航天飛機(jī)雷達(dá)地形任務(wù)(SRTM)公開的DEM(30 m)(數(shù)據(jù)來自 http:∥reverb.echo.nasa.gov/DEM)。

(3)地面觀測數(shù)據(jù)。由研究區(qū)內(nèi)青藏高原中部土壤水分多尺度觀測網(wǎng)的數(shù)據(jù)集提供。此數(shù)據(jù)網(wǎng)提供3種尺度和4種深度的土壤溫度與土壤水分的觀測數(shù)據(jù),3種尺度分別為大尺度1°、中尺度0.3°與小尺度0.1°,4種深度分別為 0～5、10、20 、40 cm,觀測間隔為0.5 h。由于被動(dòng)微波土壤水分產(chǎn)品的空間分辨率為0.25°,選取中尺度(如圖1(b)所示) 范圍內(nèi)深度為 0～5 cm的土壤水分觀測數(shù)據(jù)去驗(yàn)證分析結(jié)果[22]。

2 研究方法與原理

2.1 BP神經(jīng)網(wǎng)絡(luò)原理

BP神經(jīng)網(wǎng)絡(luò)是一種多層前饋網(wǎng)絡(luò),其核心算法為誤差反向傳播,在函數(shù)逼近、模式識別、分類等領(lǐng)域應(yīng)用較為廣泛,具有自組織、自學(xué)習(xí)能力等特點(diǎn)[23]。

基本原理:基本BP算法計(jì)算誤差輸出時(shí)是按從輸入到輸出的方向執(zhí)行的,而調(diào)整權(quán)值和閾值則相反。通過調(diào)整輸入節(jié)點(diǎn)與隱層節(jié)點(diǎn)、隱層節(jié)點(diǎn)與輸出節(jié)點(diǎn)的聯(lián)接強(qiáng)度以及閾值,從而使BP網(wǎng)絡(luò)的預(yù)測輸入不斷接近期望輸出。BP神經(jīng)網(wǎng)絡(luò)還可以用來建立多元非線性關(guān)系,因此,在遙感反演估算中被廣泛應(yīng)用[24]。BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖2所示。

圖2 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.2 BP neural network structure

2.2 ASO算法原理

ASO算法是一種基于分子動(dòng)力學(xué)模型的新穎智能算法,它通過種群間各個(gè)原子間的相互作用力來指導(dǎo)群體,從而進(jìn)行智能優(yōu)化搜索[25]。

基本原理:根據(jù)分子動(dòng)力學(xué)中原子的物理運(yùn)動(dòng)規(guī)律建立算法模型。由牛頓第二定律可知,如果Fi是作用在第i個(gè)原子上的相互作用力,Gi為作用在第i個(gè)原子上的約束力,原子質(zhì)量為mi,第i個(gè)原子的加速度可表示為:

(1)

式中:t為當(dāng)前迭代的次數(shù),d為原子所在的維數(shù)。下面對Fi、Gi、mi三個(gè)值的求解進(jìn)行闡述。

(2)

式中:Kbest為適應(yīng)度函數(shù)值比較好的k個(gè)原子的集合。

(3)

(4)

(5)

(6)

式中:xi和xj分別表示原子i和原子j的位置。

(7)

(3)原子質(zhì)量mi(t)。第t次迭代中第i個(gè)原子的質(zhì)量mi(t)由當(dāng)前種群個(gè)體的適應(yīng)度值大小決定,定義如下:

(8)

(9)

式中:fi(t)為第i個(gè)原子的適應(yīng)度值,fmax(t)、fmin(t)分別表示原子種群中最大適應(yīng)度值和最小適應(yīng)度值。

(10)

在每一次迭代過程中,根據(jù)得到的加速度更新原子i的速度和位置,更新如下:

(11)

(12)

2.3 ASO-BP算法模型

本文應(yīng)用一種基于ASO算法的反演模型,其關(guān)鍵在于通過優(yōu)化原子位置來確定BP神經(jīng)網(wǎng)絡(luò)的初始權(quán)值和閾值。首先利用ASO算法對原子位置進(jìn)行多次迭代搜索,以找到最優(yōu)的參數(shù)組合。再利用優(yōu)化后的BP神經(jīng)網(wǎng)絡(luò)模型進(jìn)行土壤水分的反演。算法的流程如圖3所示,其算法流程分為3步:

① 根據(jù)遙感影像選取所需的輸入數(shù)據(jù),將其進(jìn)行重建和均值升尺度處理后的數(shù)據(jù)作為整體的樣本集。在確定BP神經(jīng)網(wǎng)絡(luò)的拓?fù)浣Y(jié)構(gòu)之后對BP神經(jīng)網(wǎng)絡(luò)的權(quán)值與閾值進(jìn)行初始化。

② 計(jì)算ASO算法的決策變量長度,選取均方誤差作為優(yōu)化的目標(biāo)函數(shù),對BP神經(jīng)網(wǎng)絡(luò)的權(quán)值與閾值進(jìn)行優(yōu)化。ASO計(jì)算原子適應(yīng)度,更新目前找到的最優(yōu)解。根據(jù)式(8)、式(9)分別計(jì)算原子的個(gè)體加速度、速度及位置更新。設(shè)置當(dāng)達(dá)到最優(yōu)適應(yīng)度或者最大迭代次數(shù)時(shí)終止迭代的算法停止準(zhǔn)則,輸出原子個(gè)體的最優(yōu)位置,否則繼續(xù)重復(fù)ASO。

③ 使用優(yōu)化算法對神經(jīng)網(wǎng)絡(luò)的權(quán)值和閾值參數(shù)進(jìn)行優(yōu)化,將優(yōu)化后得到的權(quán)值和閾值參數(shù)賦給BP神經(jīng)網(wǎng)絡(luò)。輸入樣本集與測試集進(jìn)行反演,數(shù)據(jù)歸一化后最終獲得土壤水分。

圖3 ASO BP神經(jīng)網(wǎng)絡(luò)土壤水分重建流程Fig.3 Flowchart of ASO BP neural network soil moisture reconstruction

3 結(jié)果與討論

由NDVI、LST、Albedo與SRTM DEM數(shù)據(jù)采用均值升尺度后得到空間分辨率為0.25°的數(shù)據(jù),并與ECV V04.2作為輸入數(shù)據(jù),優(yōu)化后的土壤水分?jǐn)?shù)據(jù)作為輸出,建立BP神經(jīng)網(wǎng)絡(luò)。BP神經(jīng)網(wǎng)絡(luò)中采用了正切S型傳遞函數(shù)tansig作為傳輸函數(shù),最大的迭代次數(shù)為1 000,目標(biāo)0.000 001。采用ASO后的BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練及土壤水分反演。

為了驗(yàn)證土壤水分反演的有效性,評價(jià)傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)模型和ASO-BP神經(jīng)網(wǎng)絡(luò)模型的反演精度,選取了3個(gè)指標(biāo):相關(guān)系數(shù)(R)、均方根誤差(RMSE)和無偏均方根誤差(ubRMSE)。一般來說,R值越大,RMSE、ubRMSE越小,則反演結(jié)果越好、精度越高。

3.1 土壤水分反演結(jié)果與分析

由ASO-BP神經(jīng)網(wǎng)絡(luò)模型、BP神經(jīng)網(wǎng)絡(luò)模型重建的2015年5月1日土壤水分如圖4所示。由圖4可以看出,重建后的土壤水分與重建前ESA CCI土壤水分空間分布一致,中部和西部土壤水分值較小,東部值較大。通過重建,得到了空間連續(xù)的土壤水分。

(a)ESA CCI土壤水分圖

(b)BP重建土壤水分圖

(c)ASO-BP重建土壤水分圖

2015年1—12月ESA CCI土壤水分、重建的土壤水分和地觀測站土壤水分的時(shí)間序列如圖5所示。由圖5可以看出,重建后的土壤水分在對ESA CCI土壤水分中缺失的土壤水分進(jìn)行重建生成有效值的同時(shí),與實(shí)測土壤水分保持一致的變化趨勢。夏季土壤水分較高,冬季較低,在8月左右有一個(gè)明顯的低谷,ESA CCI土壤水分與重建土壤水分在此時(shí)間段也表現(xiàn)出明顯的低谷。但重建后的土壤水分值在4—6月以及10月低于實(shí)測值,主要是由于原始ESA CCI土壤水分在此期間低估了土壤水分,造成重建后的土壤水分值偏低。同時(shí),彌補(bǔ)了在ESA CCI土壤水分1—3月以及11—12月缺少的數(shù)據(jù),但與地面實(shí)測值相關(guān)較大,這主要是由于此時(shí)段為凍土,且在重建模型中,缺少相應(yīng)時(shí)期的土壤水分?jǐn)?shù)據(jù)作為輸入數(shù)據(jù)。因此,提高人工神經(jīng)網(wǎng)絡(luò)輸入數(shù)據(jù)的質(zhì)量是進(jìn)一步提高反演精度的關(guān)鍵步驟。

圖5 2015年土壤水分的時(shí)間序列比較Fig.5 Time series comparison of soil moisture in 2015

選取地面觀測網(wǎng)中尺度(0.3°)站點(diǎn)觀測平均值與反演得到的土壤水分進(jìn)行對比分析,圖6為BP重建土壤水分(圖6(a))和ASO-BP重建土壤水分(圖6(b))與地面觀測值的散點(diǎn)圖。 由圖6可以看出,在研究期內(nèi), ASO-BP重建土壤水分值不高于地面觀測值,土壤水分在大部分日期低于地面觀測值。而在土壤較為干燥時(shí),尤其是在土壤水分值小于 0.2 cm3·cm-3時(shí),土壤水分重建值則高于地面觀測值。這可能是由于土壤濕潤時(shí),多處于季風(fēng)期,降雨較多,植被生長較為茂盛,采用被動(dòng)微波反演土壤水分時(shí)植被透過率估算的誤差,造成重建得到的 ESA CCI土壤水分低于地面觀測0～5 cm土壤水分觀測值,而ESA CCI土壤水分在大部分日期偏低也是造成BP重建土壤水分與ASO-BP重建土壤水分值偏低的一個(gè)重要原因。相對于BP重建土壤水分,ASO-BP重建土壤水分在土壤較為干燥時(shí),與地面觀測值更為接近,主要原因是ASO對BP神經(jīng)網(wǎng)絡(luò)初始值的優(yōu)化,減小了干燥時(shí)反演土壤水分的誤差。

同時(shí),ASO-BP重建土壤水分與地面觀測值的RMSE值(0.029 cm3·cm-3),ubRMSE值為0.028 cm3·cm-3,略高于BP 重建土壤水分(0.034 cm3·cm-3),ubRMSE值為0.033 cm3·cm-3,整體精度得到提高。

整體而言,重建后的土壤水分的變化趨勢與站點(diǎn)實(shí)測數(shù)據(jù)基本保持一致,可以較好地監(jiān)測土壤水分變化情況。

3.2 不同植被覆蓋度下土壤水分精度分析

為了進(jìn)一步分析不同植被覆蓋度下ASO-BP土壤水分精度,選取NDVI值0.4作為不同植被覆蓋度閾值[26],分別分析NDVI≥0.4與NDVI<0.4時(shí)土壤水分精度。表1NDVI<0.4與NDVI≥0.4時(shí),土壤水分與地面觀測值的精度比較。

表1 不同土壤水分值與地面觀測值的精度比較

由表1可以看出,在高植被覆蓋度下(NDVI≥0.4),ASO-BP土壤水分與BP土壤水分的R值較接近,RMSE值分別為0.031、0.026 cm3·cm-3, ubRMSE從0.024 cm3·cm-3降到0.022 cm3·cm-3。在植被覆蓋區(qū)域,植被在一定程度上降低了微波信號對土壤水分的敏感性,造成微波反演土壤水分時(shí)產(chǎn)生誤差。而BP與ASO-BP算法能夠在一定程度上對植被覆蓋度較高區(qū)域的土壤水分進(jìn)行優(yōu)化,提高精度。當(dāng)NDVI<0.4時(shí), ASO-BP算法相關(guān)系數(shù)R值(0.58),高于BP重建土壤水分(0.39), ubRMSE從0.033 cm3·cm-3降到0.028 cm3·cm-3。

因此,ASO-BP優(yōu)化模型能夠提高土壤水分精度,并且在低植被覆蓋度區(qū)域(NDVI<0.4)時(shí)的效果比BP神經(jīng)網(wǎng)絡(luò)模型的效果要更好。

4 結(jié)論

本文使用ASO算法優(yōu)化BP神經(jīng)網(wǎng)絡(luò),融合了NDVI、LST、DEM、Albedo、ESA CCI土壤水分?jǐn)?shù)據(jù),建立了土壤水分重建模型。實(shí)驗(yàn)結(jié)果表明,本文算法能夠?qū)崿F(xiàn)土壤水分遙感產(chǎn)品的重建,同時(shí),提高了土壤水分反演精度,得出主要結(jié)論如下:

① 通過ASO算法對BP神經(jīng)網(wǎng)絡(luò)進(jìn)行優(yōu)化,改善了BP神經(jīng)網(wǎng)絡(luò)算法收斂速度慢、易過擬合的缺點(diǎn)。在算法中,采用ASO對BP神經(jīng)網(wǎng)絡(luò)的初始值進(jìn)行優(yōu)化,使得土壤水分重建精度較傳統(tǒng)BP算法獲得有效提高。

② 在不同植被覆蓋度下,體現(xiàn)了算法的適用性。高植被覆蓋度區(qū)域(NDVI≥0.4),ASO-BP算法優(yōu)于BP神經(jīng)網(wǎng)絡(luò)算法。ASO-BP土壤水分和BP土壤水分與地面觀測值的R值從0.79提高到0.80以上, ubRMSE從0.024 cm3·cm-3提高到0.022 cm3·cm-3。在低植被覆蓋度區(qū)域(NDVI<0.4),ASO-BP優(yōu)化模型土壤水分精度(R=0.58,RMSE=0.029 cm3·cm-3)比傳統(tǒng)BP神經(jīng)網(wǎng)絡(luò)精度(R=0.39,RMSE=0.034 cm3·cm-3)提高更為明顯。

本文中土壤水分空間分辨率低且時(shí)間序列較短,而影響土壤水分的因素復(fù)雜多樣,在一定程度上會影響土壤水分結(jié)果精度驗(yàn)證。同時(shí),在結(jié)果驗(yàn)證方面,采用像元內(nèi)地面觀測站點(diǎn)的平均值作為真實(shí)值,但易受站點(diǎn)的布設(shè)位置、區(qū)域環(huán)境等因素的影響。在進(jìn)一步研究中,將分析降水、土壤類型、土壤表面粗糙度等參數(shù)對模型的影響,以及土壤水分真實(shí)性驗(yàn)證等。