陳婷婷 潘耀忠 孫 林
(1 山東科技大學測繪科學與工程學院,山東 青島 266590)
(2 遙感科學國家重點實驗室(北京師范大學地理科學學部),北京 100875)
(3 環境演變與自然災害教育部重點實驗室,北京師范大學地理科學學部,北京 100875)
土壤水作為地表水存儲的重要組成部分,不僅直接影響著陸地和大氣間的物質和能量交換,而且是水文模型、氣候模型、生態模型以及陸面過程模型的關鍵輸入參數[1]。土壤水分監測方法從數據獲取方式上分為基于站點的土壤水分觀測、基于氣象數據和基礎地理數據的土壤水分模擬與同化,以及基于遙感數據的土壤水分反演[2]。站點測量土壤水分的測量深度和精度均比較高[3],但是由于土壤水分的時空異質性,站點觀測并不能完全表征區域內的土壤水分變化[4]。基于氣象數據和基礎地理數據驅動、利用模型模擬的土壤水分盡管在時空尺度上有連續性,但模擬精度極大程度上依賴于參數化方案和參數化過程的選擇,參數繁多較難實用。遙感技術直接以面的形式對地表進行大面積同步觀測,其監測范圍不受站點分布位置的限制[5],為土壤水分的動態監測提供了有效的方法[6]。
遙感反演土壤水分分為光學和微波遙感兩種。其中可見光—近紅外方法主要依賴于地物的反射光譜特性,由于受云及太陽光照射條件的限制,可見光、近紅外等光學遙感方法對土壤水分的研究存在一定的局限性[7]。微波遙感通過探測土壤的介電屬性變化監測土壤水分信息,而土壤介電常數與土壤水分密切相關;且微波遙感不受太陽輻射的影響,大氣的影響微乎其微,具有全天候、全天時的優勢,從而逐漸成為監測土壤水分的一種有效手段。
微波遙感反演土壤水分主要分為主動和被動微波遙感兩種,被動微波遙感重復觀測頻率高,且被動微波輻射信號對介電屬性較表面粗糙度更為敏感,因而對地表特征要求不高,但其監測空間分辨率低,適于大尺度調查。與之相比,合成孔徑雷達(Synthetic Aperture Radar, SAR)具有更高的空間分辨率(1 m~1 km),可以獲取較大范圍更為精細尺度上的土壤水分信息[8]。主動微波探測土壤水分的主要難點在于,土壤表面的后向散射系數強度與土壤水分之間并非呈簡單的線性關系,它除了與土壤水分含量有關外,還與土壤物理特性(結構、成分)、植被(數量、結構)以及雷達系統參數相關。因此,土壤水分的反演在本質上屬于“病態”問題,雷達后向散射系數和土壤水分之間的關系必然存在著不確定性。為降低這種不確定性,可借助于多種觀測模式(多極化[9-11]、多角度[9,12]、多波段[13-14])的數據或者多源數據(被動微波[15-16]、光學遙感[14,17])來反演土壤水分,但是多種觀測模式或者多源的觀測數據通常難以獲取。
已有的研究表明,利用重復觀測的SAR數據進行土壤水分反演具有較大的潛力。多時相土壤水分反演方法基于以下假設,在觀測期內,地表粗糙度和植被相較于土壤水分變化緩慢(不考慮耕作的發生),從而可將其忽略,并假定時間序列的雷達后向散射系數的變化僅由土壤水分的變化造成,該方法適于裸土及低植被覆蓋區域[18-20]。Wickel等[21]使用一個月內獲取的10景Radarsat影像對收割后的小麥地進行土壤水分監測,結果表明后向散射系數σ0的變化與土壤水分的變化存在較好的相關性(R2=0.89);Kim等[22]利用多期觀測數據建立成本函數,通過最小化成本函數在預先建立的查找表中得到多期土壤水分的解,反演結果與地面實測數據有較好的相關性;Balenzano等[23]基于上述假設提出了Alpha近似模型,采用兩個時相雷達后向散射系數的比值將粗糙度、植被覆蓋和土壤水分對后向散射系數的影響進行解耦,并構建觀測方程組反演得到每期的土壤水分值,研究表明,利用HH極化數據反演的土壤水分精度可以達到0.05 m3? m-3。該方法為利用時間序列SAR數據反演土壤水分提供了一種新思路,He等[24]對Alpha算法進行了擴展,在求解觀測方程組的過程中,通過使用Juxtaposition方法自適應確定土壤水分的上下界,不再依賴于地面實測數據;Ouellette等[25]使用L波段機載雷達數據將Alpha模型應用于植被覆蓋區域的土壤水分反演,并將觀測方程組擴展到多極化(HH和VV),結果表現出較好的適用性;Zhang等[26]針對觀測方程組欠定問題,通過融合多源時序SAR數據增加了有效觀測方程,從而實現觀測方程數大于未知數個數,提高了土壤水分反演的可靠性。
影響雷達回波強度的因素可分為雷達系統參數和地表特性,前者包括頻率、極化方式、入射角等,后者包含地表的物理特性以及幾何特性。在雷達系統參數確定的情況下,雷達回波強度主要與地表特性相關[7]。Alpha近似模型將兩個時相雷達后向散射系數的變化歸結為土壤水分的變化,這要求前后期的雷達觀測系統參數一致或者變化較小。然而在實際情況中,當雷達觀測系統確定后,頻率、極化方式隨之確定,入射角會有發生變化的情況,特別是在不同軌道過境時獲取地面影像。由于自然界大多數地表并不完全是粗糙的朗伯表面,其反射并非各向同性,因此同一地物在不同入射角下得到的后向散射會表現出差異。在應用Alpha近似模型時,如果兩個時期雷達觀測的入射角不同,那么后向散射系數的變化也包含來自入射角變化的影響。針對上述問題,本文提出改進的Alpha模型,將時間序列SAR數據的入射角的變化考慮進來,選取了兩個典型實驗區:我國黑河流域的農田區域和加拿大曼尼托巴省農田區域,利用歐洲太空局 (European Space Agency, ESA) “哨兵系列”的雷達衛星Sentinel-1數據及對應區域的土壤水分觀測網絡數據對算法進行驗證,并對結果進行分析與討論。
本文所選實驗區均位于農田區域,地表在作物收成之后一段時間處于裸露狀態,土壤水分反演可不考慮耕作對土壤粗糙度的影響及植被覆蓋的影響;地形較為平坦,因此地表起伏對雷達后向散射系數的影響也可忽略;此外均有長時間序列的站點觀測數據對土壤水分的反演結果做驗證支持。實驗區1位于我國黑河流域(100°10′~100°30′E,38°45′~38°55′N),屬大陸性溫帶干旱氣候;其中游屬于農業經濟區,也是西北地區重要的商品糧基地。實驗區2位于加拿大的曼尼托巴省溫尼伯市(98°30′~97°W,49°~50°30′N),地處北紅河與阿西尼泊因河交匯的廣闊平原,北部緊鄰溫尼伯湖,氣候屬溫帶大陸性氣候。
Sentinel-1是歐洲航天局哥白尼計劃(GMES)中的地球觀測衛星,由A、B兩顆星組成,主要服務于陸地和海洋監測,由歐空局(European Space Agency,ESA)組織研發。本文所用數據為干涉寬幅模式(Interferometric wide swath,IW)下經過多視處理和地距轉換的GRD格式產品。利用SNAP軟件對雷達影像進行預處理,主要包括熱噪聲去除、軌道文件校正、輻射定標、濾波校正、多普勒地形校正,最終獲得入射角和后向散射信息。表1給出了本文所用黑河和溫尼伯實驗區各期時序Sentinel-1影像的具體信息。以黑河實驗區為例,當衛星按相同的方向且以同一軌道過境時,雷達波束的入射角不變,如2016-11-07和2016-12-01獲取的影像,站點平均入射角均為43.7°;當衛星過境方向或者軌道發生變化,則入射角隨之改變,有的變化較小,如2016-11-07和2016-11-13獲取的影像,有的變化較大如2016-12-07和2016-12-14的影像,相對于整個時間序列的Sentinel-1影像入射角發生了10°左右的變化,當角度變化較小時其對后向散射的影響可忽略不計,而當角度明顯變化時其影響將不能忽略,Alpha近似模型的應用會受到影響。

表1 實驗區Sentinel-1 數據Table 1 List of Sentinel-1 SAR images of the experiment areas
黑河流域實驗區所用地面數據來自黑河生態水文遙感試驗(Heihe Watershed Allied Telemetry Experimental Research,HiWATER)。本文所用數據包括黑河流域中游張掖市的三個站點的數據[27-28],分別位于大滿、黨寨和花寨子,土地類型分別為農田、草地和荒漠。以上站點均提供全年逐日逐時(每10 min)不同深度(2 cm、4 cm、10 cm等)的土壤水分含量(%)、土壤溫度(℃)數據,本文使用其中4 cm深度的土壤數據,根據Sentinel-1影像的獲取時間選擇與其最接近時刻的觀測數據作為實測值;土壤質地、土壤容重數據來自2012年黑河流域土壤參數觀測數據集,其中位于黨寨的站點沒有相應土壤質地數據,以花寨子的數據代替。
加拿大曼尼托巴省實驗區的地面數據來自于RISMA(Real-time In-situ Soil Monitoring for Agriculture)土壤水分觀測網絡[29-31],由加拿大農業和農業食品部(Agriculture and Agri-Food Canada’s, AAFC)于2011年組織建立。本文所用12個監測站點的數據其中有9個位于溫尼伯市西南部的農田區域,3個位于溫尼伯市西北部的農田區域,每個站點均提供觀測期間不同深度(0~5 cm、5 cm、20 cm等)的土壤水分(%)、土壤溫度(℃)數據,并且每個深度均布設了3個探頭。本文使用其中5 cm深度的土壤數據,選擇與衛星過境時間最為接近時刻的觀測數據作為實測值,對于3個傳感器的觀測值本文取其平均;對于土壤質地、土壤容重,本文以世界土壤數據庫HWSD(V1.2)提供的1 km網格數據作為參考。
微波散射理論模型是通過建立后向散射系數與地表物理和幾何參數間的數學關系來研究粗糙表面的散射,以求解介電常數和土壤水分。小擾動模型SPM[32]是研究小尺度粗糙度表面的經典模型,其一階形式描述的是面散射,二階形式為體散射。由于大多數自然地表面散射占主要成分,因此常用其一階形式。以垂直極化VV為例,隨機粗糙表面的后向散射系數可以表示為[33]:

式中,k表示雷達波數,表達為:k=2π/λ,h為均方根高度,θ為入射角,εr、μr分別為介電常數和磁導率,W為表面粗糙度譜函數。當磁導率取值為1時,上式可簡化為:

式中,αVV為極化幅度,是雷達入射角度和介電常數的函數。因此式(2)可進一步表達為:

在低植被覆蓋情況下,假設植被參數對雷達后向散射的影響是乘性的,即在上式的基礎上乘植被衰減因子。假定地表粗糙度和植被狀況在一定時間內是不變的,只考慮土壤介電常數(主要受土壤水分影響)的變化,則將兩個時間觀測的后向散射系數作比值處理,不考慮入射角的變化,化簡后可得到式(5),即為Alpha近似模型[23]。

當時間序列SAR數據入射角變化較大時,跟入射角有關的各項將不能抵消,式(5)已不再成立。由于粗糙度譜函數主要是一個地表粗糙度相關的參數(表面相關長度),因此本文只考慮了與入射角直接相關的乘積項cos4θ,從而得到改進的Alpha模型如式(6)。

基于兩個時相的SAR數據,根據式(6)可得到觀測方程:

對于連續N期SAR影像(T1、T2、T3,…,TN),可以構建N-1 個有效觀測方程,組成如下方程組:

對于方程組(8),N-1方程中存在N個極化幅度未知數,是一個欠定方程組,存在無數多個解。可采用邊界約束的最小二乘法進行求解,根據實際的雷達入射角度和土壤水分范圍對αVV取值范圍進行限定。本文黑河實驗區所用站點分別位于不同的土地覆蓋類型區,曼尼托巴省的站點分布于不同區域的農田之中,因此每個站點各自設定αVV邊界。得到極化幅度αVV后,進而根據式(3)求得土壤的介電常數εr,最后利用介電混合模型[34]求得土壤水分mV。
本文使用黑河、溫尼伯實驗區各自獲取的6期VV極化時序Sentinel-1數據(T1、T2、T3、T4、T5、T6),分別基于原始Alpha模型及改進的Alpha模型構建觀測方程組,結合實測數據嚴格設定α系數邊界進行土壤水分反演,并利用相應時間的站點觀測數據對結果進行驗證。圖1給出了兩種方法反演的結果散點圖對比(所有站點所有時相的數據),其中橫軸為實測土壤水分(%),縱軸為模型反演結果(%),左邊的散點圖為原始模型的反演結果,右邊為相應的改進后模型的反演結果,明顯改進后模型的反演結果更靠近Y=X直線,即結果更接近真實值。
為定量評定基于時序Sentinel-1數據的改進Alpha模型的土壤水分反演精度,分別計算了反演結果與實測土壤水分之間的相關系數r和均方根誤差RMSE。由圖1可知,對于黑河實驗區,整體上土壤水分偏低(2%~16%),基于原始方法得到的土壤水分結果與實測值的相關系數為0.473(P>0.05),相關性不明顯(排除一個異常點,其反演得到的結果為復數);而在考慮入射角的變化之后,相關系數提高至0.851(P<0.01),均方根誤差也由0.053減小至0.023,結果更靠近Y=X線。對于溫尼伯實驗區,與黑河獲取的相同時間范圍的影像進行土壤水分反演,相對于黑河實驗區該區域的土壤更加濕潤,且變化范圍更大(9%~44%)。通過原始Alpha模型所得的土壤水分結果與真實值之間的相關系數r=0.601(P<0.01),表現出一定的相關性,但明顯有部分觀測點的反演結果偏高,遠離回歸直線;而通過改進的Alpha模型反演得到的結果如上圖所示,原來離散的點均被收聚至回歸線附近,更靠近Y=X線,相關系數由0.601提高至0.821(P<0.01),RMSE則由0.152降低至0.065。綜合兩個實驗區的所有站點數據如圖2所示,模型改進后得到的土壤水分與實測數據間的相關性由0.720提高至0.893,RMSE由0.138減小至0.059,改進后Alpha模型得到的結果優于原始模型。

圖1 實驗區土壤水分反演結果散點圖Fig.1 Scatter graph of the inversion of soil moisture in the experimental areas
為進一步評價模型質量本文引入了水文領域中常用的模型評價指標納什效率系數NSE(Nash Sutcliffe efficiency coefficient),該系數一般用以驗證水文模型模擬結果的好壞,也可用于其他模型。其計算公式如下:

為進一步說明改進Alpha模型土壤水分反演結果,分別在兩個實驗區選擇單個站點的反演結果進行對比分析。如圖3所示,橫軸代表觀測日期,縱軸代表土壤水分,方框實線代表實測土壤水分,空心圓虛線代表原始模型的反演結果,實心圓虛線代表改進模型的反演結果。由圖可知,兩個站點使用改進后模型的反演結果相較于原始方法更接近于實測情況。對于黑河實驗區的大滿站點,其2016-11-24的土壤水分相較于前一時期是減小的,但是原始模型得到的結果呈上升趨勢,由影像的入射角信息可知(表1),在該時期中大滿站點的入射角約為33°,較之前一時期減少了近10°,根據Ulaby等[35]的研究,在超過某一角度臨界值后,隨著觀測的入射角增大雷達接收的后向散射信號呈減小趨勢,由此推斷該時期雖然土壤水分降低使得后向散射系數減小,但由于入射角的減小使得后向散射系數顯著增高。本文使用cos4θ這一乘性因子將入射角對地物散射信號的影響考慮進來,使得反演結果捕捉到這一變化趨勢。對于2016-12-14的反演結果,其發生的角度變化同上述一致,該時期得到的結果并沒有捕捉到土壤水分減小趨勢,其可能原因是粗糙度變化造成的,但是相較于原始方法本文的方法仍在一定程度上減小了反演誤差。對于溫尼伯實驗區,入射角是在41°和32°之間交替變化的,因此原始方法得出的結果也呈高低交替的變化趨勢,應用改進后模型得到的結果相比于原始模型能較好地還原土壤水分變化的趨勢,且誤差也明顯減小。

圖3 單個站點土壤水分反演結果Fig.3 Inversion of soil moisture of a single site
以上從單個站點時間序列角度對反演結果進行了對比分析,接下來本文將從空間角度對反演的土壤水分進行分析與討論。表2列出了溫尼伯實驗區使用原始和改進的Alpha模型反演得到的各時期土壤水分結果精度對比,由表可知無論是原始模型還是改進后模型其結果的相關性均很顯著,相關系數r均大于0.8,改進后模型的反演結果相關性略有提高,r均在0.9以上。這是由于溫尼伯實驗區各個站點的南北縱向分布,使得衛星過境掃描時站土壤水分點間的入射角差異較小,最大角度差不超過3°,因此,雷達觀測的后向散射系數與土壤水分之間呈較好的相關性,受入射角影響較小。但是從回歸直線的斜率以及均方根誤差RMSE的角度來評價反演結果,原始模型則在整體上表現出較差的適用性,而改進后Alpha模型的反演結果則有明顯提高,回歸直線斜率更加接近于1,即更接近實測土壤水分值;RMSE(m3?m-3)變化范圍由原來的0.043~0.290減小至0.034~0.10。綜上所述,本文提出的改進Alpha近似模型可以很好地校正由于時序SAR數據入射角變化帶來的土壤水分反演結果誤差。

表2 溫尼伯實驗區各時期反演結果精度對比Table 2 Comparison of inversions of soil moisture in the Winnipeg experimental area in accuracy relative to time
值得注意的是,當入射角發生變化時不僅影響地物的后向散射,也會影響地表粗糙度的響應模式,即在不同入射角下粗糙度對后向散射信號的影響存在差異,因此利用本文的改進模型并不能完全消除入射角變化帶來的影響,粗糙度的影響仍存在其中。此外通過衛星觀測得到的是像元范圍內的土壤水分的整體情況,而站點觀測得到的是單個點的信息,土壤水分的空間異質性及測量過程的不確定性使得這種尺度轉換存在誤差,從而影響算法的結果驗證,上述問題均需要進一步的研究和分析。
基于Alpha模型,使用C波段的雷達數據進行土壤水分反演,不同時相數據不同的入射角影響結果的反演精度甚至使得模型不可用,通過本文提出的方法可以很好地校正角度變化帶來的影響,擴展了Alpha模型的應用。改進的Alpha模型是在小擾動模型SPM基礎上推導得到的,小擾動模型只適用局限的表面粗糙度范圍,通過本文實驗表明該模型在農田區域具有較好的適用性。
本文研究存在一些不足,實驗中使用Sentinel-1 C 波段雷達數據在農田區域進行模型驗證,未考慮在其他地表類型(如草地)及其他波段(如L、X波段)雷達數據的適用性,此外需要進一步探討改進后Alpha模型在高分三號衛星中的應用,其空間分辨率可達1 m,對于中小尺度的土壤水分反演研究具有重要意義。本研究也未考慮植被對于后向散射信號的影響,今后研究可結合相關的植被覆蓋區域的微波散射模型及光學影像對Alpha模型進行改進。