本發(fā)明涉及地理測繪技術(shù)領(lǐng)域，特別是涉及一種土壤含水量產(chǎn)品降尺度方法。

背景技術(shù)：

土壤含水量作為氣候與環(huán)境干旱化的重要指示因子之一，影響著土壤的理化性質(zhì)與植被的生長，進而影響我國的糧食產(chǎn)量。同時，土壤含水量是地表能量平衡和水循環(huán)的重要組成部分，是全球變化研究中的重要監(jiān)測因子。目前，各種基于全球觀測站點資料建立的土壤含水量數(shù)據(jù)集，由于觀測點的密度和空間代表性不足，模擬和預(yù)報的精度難以滿足應(yīng)用需求。

先進微波掃描輻射計AMSR-2(Advanced Microwave Scanning Radiometer for the Earth Observing System)，土壤水分和海洋鹽分傳感器SMOS(Soil Moisture and Ocean Salinity)，SMAP(Soil Moisture Active Passive)以及風(fēng)云三號氣象衛(wèi)星均具有全天時、全天候、觀測尺度大、重訪周期短等突出優(yōu)勢，能夠提供高覆蓋度高時效性的全球土壤含水量數(shù)據(jù)。但這些數(shù)據(jù)的空間分辨率(9-40km)較低，無法滿足流域尺度土壤含水量時空動態(tài)監(jiān)測的需求。光學(xué)遙感數(shù)據(jù)空間分辨率可以達到1km以下，高空間分辨率、低時間分辨率、易受天氣影響的特點，與被動微波遙感數(shù)據(jù)恰恰相反。

為了解決被動微波土壤含水量產(chǎn)品的空間分辨率低的問題，國內(nèi)外學(xué)者也提出了各種不同的降尺度方法：一類降尺度模型是以遙感技術(shù)獲得的土壤物理參數(shù)為基礎(chǔ)，如利用光學(xué)遙感數(shù)據(jù)反演土壤的蒸散量建立降尺度算法，但是這種算法的局限性在于沒有考慮微波遙感數(shù)據(jù)反演的土壤含水量和光學(xué)遙感數(shù)據(jù)反演的土壤有效蒸發(fā)量之間存在很強的非線性關(guān)系。還有一類模型是四維變分同化方法，就是在微波遙感數(shù)據(jù)的第四維尺度上分析土壤含水量。但是這種方法對地表資料的獲得要求較高。

由此可見，目前迫切需要發(fā)展一種新的土壤含水量產(chǎn)品降尺度方法，建立一種可行的空間降尺度模型，逐步提高土壤含水量數(shù)據(jù)的空間分辨率，推動流域尺度土壤含水量的時空動態(tài)監(jiān)測。

技術(shù)實現(xiàn)要素：

本發(fā)明要解決的技術(shù)問題是提供一種土壤含水量產(chǎn)品降尺度方法，使其提高被動微波土壤含水量數(shù)據(jù)產(chǎn)品的空間分辨率，滿足流域尺度水資源和農(nóng)業(yè)管理的應(yīng)用需求，克服目前土壤含水量產(chǎn)品數(shù)據(jù)集空間代表性不足、空間分辨率低等的不足。

為解決上述技術(shù)問題，本發(fā)明提供一種土壤含水量產(chǎn)品降尺度方法，所述方法包括以下步驟：

A.獲取待研究區(qū)的被動微波土壤含水量產(chǎn)品和同一時間的光學(xué)遙感影像數(shù)據(jù)；

B.基于多端元混合像元分解方法對所述光學(xué)遙感影像數(shù)據(jù)進行土壤光譜的提??；

C.利用GA-PLS建立所述土壤光譜反射特征與從所述被動微波土壤含水量產(chǎn)品中獲取的土壤含水量之間的定量關(guān)系模型；

D.基于步驟C建立的所述定量關(guān)系模型，利用泰勒級數(shù)展開形式構(gòu)建土壤含水量降尺度模型，獲得高空間分辨率的土壤含水量數(shù)據(jù)。

作為本發(fā)明的一種改進，所述步驟A中被動微波土壤含水量產(chǎn)品采用SMAP土壤含水量數(shù)據(jù)；所述光學(xué)遙感影像數(shù)據(jù)采用MODIS影像數(shù)據(jù)。

進一步改進，所述步驟B中提取所述土壤光譜的方法為：將高空間分辨率的MODIS圖像重采樣到與SMAP數(shù)據(jù)同樣的低空間分辨率，分別對該高空間分辨率的MODIS影像和重采樣后低空間分辨率的MODIS應(yīng)用MESMA方法進行土壤光譜的提取。

進一步改進，所述MESMA方法包括光譜庫創(chuàng)建、最優(yōu)光譜庫選取和多端元混合像元分析步驟，

所述光譜庫創(chuàng)建包括基于ROI創(chuàng)建光譜庫、光譜庫元數(shù)據(jù)制作和光譜庫管理；

所述最優(yōu)光譜庫的選取包括創(chuàng)建方形陣列和光譜庫優(yōu)選；

所述多端元混合像元分析采用植被-不透水面-土壤模型，將優(yōu)選的植被、不透水面、土壤光譜集組合構(gòu)成2EM、3EM、4EM混合像元分析模型，基于所述最優(yōu)光譜庫的優(yōu)選結(jié)果對MESMA結(jié)果進行陰影歸一化處理，得到各端元豐度值和表示結(jié)果精度的均方根誤差，再利用所述各端元豐度值和下列公式得到研究區(qū)土壤端元光譜，

其中，R_s(λ)為土壤光譜在波段λ的反射率，R(λ)為像元在波段λ上的反射率，R(i，λ)為第i個端元在波段λ上的反射率，f_i為第i個端元豐度值，N為端元個數(shù)，ε_λ是殘差，所有端元組分的豐度值之和定義為1。

進一步改進，所述步驟C建立所述土壤光譜反射特征與所述土壤含水量定量關(guān)系模型是基于所述土壤光譜計算得到的每個像元的土壤光譜中各波段的反射率、波段比值、曲率與所述土壤含水量的定量關(guān)系。

進一步改進，所述步驟D中泰勒級數(shù)展開形式的表達式為：

其中，θ_n-1和θ_n分別代表低空間分辨率和高空間分辨率土壤含水量，R_n-1(λ_i)_s和R_n(λ_i)_s分別表示低空間分辨率和高空間分辨率的土壤光譜在第i波段的反射率，Ratio_n-1(j)_s)和Ratio_n(j)_s分別表示低空間分辨率和高空間分辨率的土壤光譜波段比值，Curv_n-1(k)_s和Curv_n(k)_s分別表示高空間分辨率和低空間分辨率的土壤光譜曲率，M、N和L分別代表反射率i、波段比值j和曲率k的變量總數(shù)。

進一步改進，所述步驟D中降尺度模型的建立方法采用逐步遞減方式進行降尺度運算，逐步達到高空間分辨率的要求。

進一步改進，所述步驟D中高空間分辨率的土壤含水量數(shù)據(jù)是指500m空間分辨率的土壤含水量數(shù)據(jù)。

進一步改進，所述步驟A中還同時獲取所述待研究區(qū)的土壤輔助數(shù)據(jù)，所述土壤輔助數(shù)據(jù)包括：土地利用或土地覆蓋分類數(shù)據(jù)，以及DEM數(shù)據(jù)。

進一步改進，所述步驟C建立的GA-PLS模型還可以增加每個像元的DEM數(shù)據(jù)或坡度數(shù)據(jù)作為GA-PLS模型的輸入?yún)?shù)。

采用上述技術(shù)方案，本發(fā)明至少具有以下優(yōu)點：

(1)本發(fā)明基于光學(xué)遙感和被動微波遙感數(shù)據(jù)，獲取高空間分辨率的土壤含水量方法，可滿足大范圍流域尺度區(qū)域研究，準(zhǔn)確度高，易于建立，省時省力。

(2)本發(fā)明可擴展性高，在應(yīng)用的過程中，可根據(jù)實際情況，進行輔助數(shù)據(jù)或土壤含水量定量反演模型關(guān)系項的增減，還可以采用逐級回歸式的降尺度方法，以逐步推進的方式提高被動微波土壤含水量產(chǎn)品的空間分辨率，不斷提高模型的計算精度。

附圖說明

上述僅是本發(fā)明技術(shù)方案的概述，為了能夠更清楚了解本發(fā)明的技術(shù)手段，以下結(jié)合附圖與具體實施方式對本發(fā)明作進一步的詳細說明。

圖1是本發(fā)明土壤含水量產(chǎn)品降尺度方法的原理流程圖。

圖2是本發(fā)明中土壤最優(yōu)光譜選取結(jié)果圖。

圖3是本發(fā)明中植被最優(yōu)光譜選取結(jié)果圖。

圖4是本發(fā)明中9km空間分辨率MODIS圖像MESMA結(jié)果圖。

圖5是本發(fā)明中樣本點分布圖。

圖6是本發(fā)明中9km尺度GA-PLS模型精度評價結(jié)果。

圖7是本發(fā)明中5km空間分辨率MODIS圖像MESMA結(jié)果圖。

圖8是本發(fā)明中5km土壤含水量降尺度精度評價結(jié)果。

圖9是本發(fā)明中500m空間分辨率MODIS圖像MESMA結(jié)果圖。

圖10是本發(fā)明中5km尺度GA-PLS模型精度評價結(jié)果。

圖11是本發(fā)明中500m土壤含水量降尺度精度評價結(jié)果。

具體實施方式

本發(fā)明在現(xiàn)有技術(shù)的基礎(chǔ)上，利用土壤的反射光譜特性，即土壤含水量是影響土壤光譜反射率的一個重要因素，又由于土壤含水量具有非常大的時空變異性，高空間分辨率的土壤光譜信息能夠更好的反映土壤含水量在時空范圍上的變化特征，則本發(fā)明利用土壤的光譜特征與土壤含水量之間的定量關(guān)系來完成土壤含水量的降尺度研究，為土壤含水量的高時空分辨率監(jiān)測提供了一個全新的思路。

又鑒于被動微波遙感在全球土壤含水量數(shù)據(jù)獲取方面的優(yōu)勢，創(chuàng)設(shè)一種可行的空間降尺度模型，逐步提高土壤含水量數(shù)據(jù)的空間分辨率，通過光學(xué)遙感和被動微波遙感的綜合應(yīng)用，提高被動微波土壤含水量數(shù)據(jù)產(chǎn)品的實用性，最終推動流域尺度土壤含水量的時空動態(tài)監(jiān)測。其具體的土壤含水量產(chǎn)品降尺度方法如下。

參照附圖1所示，本發(fā)明土壤含水量產(chǎn)品降尺度方法，主要包括以下步驟：

A.獲取被動微波土壤含水量產(chǎn)品和同一時間的光學(xué)遙感影像數(shù)據(jù)；

其具體的數(shù)據(jù)獲取方法為：

獲得待研究區(qū)SMAP土壤含水量產(chǎn)品。SMAP衛(wèi)星是于2015年1月由美國NASA發(fā)射，其上搭在一個L波段的雷達和一個L波段的輻射計。下述實施例選取的是SMAP L4空間分辨率為9km的土壤水日平均數(shù)據(jù)。下載地址為：https://ns idc.org/。

目前高空間分辨率的光學(xué)圖像時間分辨率往往較低，受氣象條件的影響，影像獲取成功率低，相比之下，中低空間分辨率的光學(xué)影像因時間分辨率高更易獲取。下述實施例光學(xué)遙感影像數(shù)據(jù)選取的是與SMAP土壤含水量數(shù)據(jù)同步的空間分辨率為500m的MODIS陸地產(chǎn)品中的MOD09地面反射率數(shù)據(jù)。其中采用該產(chǎn)品的1-7波段(620～670，841～876，459～479，545～565，1230～1250，1628～1652，2105～2155nm)數(shù)據(jù)建立土壤含水量反演模型。

該步驟中同時還獲取待研究區(qū)土壤輔助數(shù)據(jù)，其包括：(1)土地利用/土地覆蓋分類數(shù)據(jù)，可以獲得MODIS土地覆蓋類型產(chǎn)品，用于混合像元分解過程中端元組分的選取和驗證，最大程度去除植被等其它地物對農(nóng)業(yè)土壤光譜的影響；(2)DEM數(shù)據(jù)，在受到地形影響的地區(qū)，該數(shù)據(jù)可用于本發(fā)明降尺度模型的修正和驗證。

B.基于多端元混合像元分解方法進行土壤光譜的提取；

多端元混合像元分解法(Multiple Endmember Spectral Mixture Analysis，MESMA)，首先以端元的物理意義為理論基礎(chǔ)，為每類地物選取多條光譜，并以此生成多個端元組合(每個端元組合由不同地物中的某一條光譜組成)，接著對每個像元尋找最小二乘法誤差最小的端元組合，進而求出每個像元的端元比例。

本發(fā)明中，將500m空間分辨率的MODIS圖像重采樣到與SMAP數(shù)據(jù)同樣的空間分辨率9km，分別對500m的MODIS影像和重采樣后9km的MODIS應(yīng)用MESMA方法進行土壤光譜的提取，獲得9km空間分辨率和500m空間分辨率的土壤光譜。

其中，MESMA處理分析過程主要包括光譜庫創(chuàng)建、最優(yōu)光譜庫選取、多端元混合像元分析三部分。其光譜庫創(chuàng)建包括基于ROI創(chuàng)建光譜庫、光譜庫元數(shù)據(jù)制作和光譜庫管理；其最優(yōu)光譜庫的選取是多端元混合像元分解成功的關(guān)鍵，該最優(yōu)光譜庫選取過程包括創(chuàng)建方形陣列和光譜庫優(yōu)選，其優(yōu)選的方法是結(jié)合COB(Count-based Endmember Selection)、EAR(Endmember Average RMES)、MASA(Minimum Average Spectral Angle)計算，得到最能代表各類地物的光譜庫；其混合像元分解過程主要采用V-I-S(植被-不透水面-土壤)模型，將優(yōu)選的植被、不透水面、土壤光譜集組合構(gòu)成2EM、3EM、4EM混合像元分析模型，基于優(yōu)選結(jié)果對MESMA結(jié)果進行陰影歸一化處理，最終得到不同端元組合的混合像元分解結(jié)果，即各端元豐度值和表示結(jié)果精度的均方根誤差(RMSE)。

再利用上述MESMA得到的端元豐度值和下列公式(1)得到研究區(qū)土壤端元光譜。

其中，R_s(λ)為土壤光譜在波段λ的反射率，R(λ)為像元在波段λ上的反射率，R(i，λ)為第i個端元在波段λ上的反射率，f_i為第i個端元豐度值。N為端元個數(shù)，ε_λ是殘差。所有端元組分的豐度值之和定義為1。

本實施例還結(jié)合ENVI和Matlab設(shè)計了純凈像元光譜的自動化提取方法。利用ENVI軟件提取各個采樣點端元組合模型的豐度值、RMSE和殘差，確定采樣點選取的端元組合模型并輸出采樣點組合豐度值f_i。

首先在隨機選取的采樣點中根據(jù)土地利用分類類型剔除掉包含不透水面的采樣點，對影像只建立S-V模型，對S-V光譜的組合模型賦屬性值，提取各組合模型中土壤端元和植被端元的坐標(biāo)值，提取光譜庫各條光譜坐標(biāo)值，并對兩者進行匹配，得到S-V光譜組合模型的組合光譜；提取各個采樣組合模型屬性值，與S-V光譜組合模型庫中屬性值匹配得到各個采樣點的S-V端元組合光譜R(i,λ)。根據(jù)上述混合像元分解公式(1)計算得到每一個采樣點土壤光譜。

C.利用GA-PLS建立土壤光譜和被動微波獲取的土壤含水量之間的定量關(guān)系模型；

在上述9km空間分辨率的土壤光譜獲取的基礎(chǔ)上，基于土壤光譜隨土壤含水量變化的特征規(guī)律，引入偏最小二乘－遺傳算法(GA-PLS)，構(gòu)建土壤含水量與光譜的定量關(guān)系模型，作為降尺度的基礎(chǔ)。

其中PLS是多因變量對多自變量的回歸建模方法，可以應(yīng)用在全譜數(shù)據(jù)或部分譜數(shù)據(jù)的分析。它將數(shù)據(jù)矩陣的分解和回歸相結(jié)合，得出與預(yù)測組分相關(guān)的特征值向量，這使其可以應(yīng)用于復(fù)雜的分析體系，得到更穩(wěn)健的結(jié)果。遺傳算法(GA)的基本思想是在演化過程中進行自然選擇，該演化過程由基因重組和變異來實現(xiàn)。該算法最顯著的優(yōu)點是避免了初始值選擇的問題。本發(fā)明為了減小計算量，提高PLS模型的精度，使用遺傳算法(GA)對光譜參數(shù)進行篩選，剔除對含水量變化反應(yīng)不明顯的變量。

本發(fā)明使用GA-PLS模型建立土壤含水量反射光譜特征模型，如將9km的MODIS反射率、波段比值、曲率及SMAP L4土壤水產(chǎn)品作為GA-PLS模型的輸入變量，得到9km尺度上的MODIS反射率、波段比值、曲率與土壤含水量之間的GA-PLS模型，并進行精度評價。

D.基于步驟C建立的定量關(guān)系模型，利用泰勒級數(shù)展開的形式構(gòu)建土壤含水量降尺度模型，獲得空間分辨率為500m的土壤含水量數(shù)據(jù)。

本發(fā)明降尺度的依據(jù)是土壤含水量與土壤光譜反射率之間的定量關(guān)系，即在步驟C中闡述的基于被動微波土壤水產(chǎn)品與MODIS土壤反射特征如各波段的反射率、波段比值(如MODIS3/MODIS1)、曲率(如：MODIS3×MODIS1/(MODIS2)²)的GA-PLS模型。在空間尺度的變化過程中，土壤含水量與反射率、波段比值、曲率等的關(guān)系變化是非線性的，所以本發(fā)明利用泰勒級數(shù)展開的形式，如公式(2)，來構(gòu)建不同尺度下的土壤含水量之間的關(guān)系。

其中，θ_n-1和θ_n分別代表低空間分辨率(如：9km)和高空間分辨率(如：500m)土壤含水量，R_n-1(λ_i)_s和R_n(λ_i)_s分別表示低空間分辨率和高空間分辨率的土壤光譜在第i波段的反射率，Ratio_n-1(j)_s)和Ratio_n(j)_s分別表示低空間分辨率和高空間分辨率的土壤光譜波段比值，Curv_n-1(k)_s和Curv_n(k)_s分別表示高空間分辨率和低空間分辨率的土壤光譜曲率。M、N和L分別代表反射率(i)、波段比值(j)和曲率(k)的變量總數(shù)。關(guān)系式f表示由GA-PLS模型分別得到的關(guān)于低空間分辨率土壤光譜反射率、波段比值和曲率與土壤含水量的關(guān)系式。

該步驟D可根據(jù)情況，采用逐步遞減的方式進行降尺度運算，如先選擇一個尺度間隔t(如4km)將9km土壤含水量數(shù)據(jù)降到(9-t)尺度，通過反復(fù)使用公式(2)，將前一步計算出的高空間分辨率數(shù)據(jù)作為下一次公式低空間分辨率數(shù)據(jù)再次計算得到更高空間分辨率的土壤含水量，最終逐步逼近達到最后500m空間分辨率的要求。其中最優(yōu)尺度間隔t，可以根據(jù)實際數(shù)據(jù)獲取情況和結(jié)果精度來確定。

該分步降尺度方法的另一個目的是在降尺度計算過程中，不斷控制土壤含水量結(jié)果的反演精度。每降一級尺度反演的土壤含水量結(jié)果，將與同一空間分辨率土壤含水量實測數(shù)據(jù)建立相關(guān)關(guān)系，若驗證數(shù)據(jù)相關(guān)性較高，降尺度計算繼續(xù)；若驗證數(shù)據(jù)相關(guān)性較差，降尺度運算適時停止，并進行影響因素分析。如果降尺度的精度受到地形、植被等因素的影響，可為降尺度模型增加新的關(guān)系項，及時進行調(diào)整和優(yōu)化。最后，模型的降尺度結(jié)果將通過同步野外土壤含水量試驗數(shù)據(jù)得到驗證。

下述以對美國中部平原地區(qū)的SMAP土壤含水量數(shù)據(jù)進行降尺度的具體實例來說明本發(fā)明基于土壤光學(xué)特性的土壤含水量產(chǎn)品降尺度方法的過程。

在該實例中，被動微波土壤含水量產(chǎn)品選擇9km空間分辨率SMAP L4數(shù)據(jù)，光學(xué)圖像選擇500m空間分辨率的MODIS 09A1數(shù)據(jù)。

獲取數(shù)據(jù)后，首先對MODIS數(shù)據(jù)利用ENVI軟件進行投影轉(zhuǎn)換、鑲嵌、裁剪、重采樣等預(yù)處理。

然后進行最優(yōu)光譜庫的選取，通過目視識別的方式分別建立土壤、植被類型光譜庫，利用COB(Count-based Endmember Selection)、EAR(Endmember Average RMES)、MASA(Minimum Average Spectral Angle)計算優(yōu)選最能代表土壤和植被的光譜曲線，建立最優(yōu)光譜庫，最終選擇出土壤最優(yōu)光譜13條，如附圖2所示，植被最優(yōu)光譜6條，如附圖3所示。

再將優(yōu)選的植被、土壤光譜集組合構(gòu)成3EM混合像元分析模型，對結(jié)果進行陰影歸一化處理，最終得到不同端元組合的混合像元分解結(jié)果。將500m空間分辨率的MODIS圖像進行重采樣，獲得9km、5km空間分辨的圖像。

對9km分辨率的MODIS圖像進行MESMA分析，如附圖4所示。計算得到9km分辨率土壤光譜的反射率、波段比值和曲率。隨機選取51個采樣點，如附圖5所示，其中，36個點作為校正集，15個點作為驗證集，分別提取51個采樣點的土壤光譜和SMAP土壤含水量，并將其作為參數(shù)輸入建立土壤光譜和含水量之間的GA-PLS模型。結(jié)果顯示，校正集決定系數(shù)為0.78，驗證集決定系數(shù)為0.72，如附圖6所示。

GA-PLS模型作為降尺度模型的基礎(chǔ)，建立逐級空間降尺度模型。逐級降尺度選擇空間間隔t＝4km，先將土壤含水量降至5km，需要對重采樣至5km分辨率MODIS數(shù)據(jù)進行MESMA處理分析，如附圖7所示，進一步計算得到5km尺度的土壤光譜的反射率、波段比值和曲率，利用上述公式(2)得到5km空間分辨率的土壤含水量，經(jīng)驗證5km分辨率含水量數(shù)據(jù)與實測數(shù)據(jù)相關(guān)性較9km并未降低，與站點得到的實測數(shù)據(jù)進行對比分析的精度R²為0.54，如附圖8所示。

重復(fù)上一步操作，在5km尺度下建立土壤光譜和含水量的GA-PLS模型，共選取了89個采樣點，如附圖5所示，其中62個校正集樣本，27個驗證集樣本，結(jié)果顯示，5km尺度下校正集決定系數(shù)為0.71，驗證集決定系數(shù)為0.67，如附圖9所示。將500m分辨率MODIS數(shù)據(jù)進行MESMA處理分析，如附圖10所示，并計算得到500m尺度的土壤光譜的反射率、波段比值和曲率，利用5km尺度建立的GA-PLS模型，再次運用公式(2)得到500m空間分辨率土壤含水量數(shù)據(jù)。利用實測土壤含水量與500m土壤含水量進行相關(guān)性分析，與站點得到的實測數(shù)據(jù)進行對比分析的精度R²為0.49，如附圖11所示。

本發(fā)明利用被動微波遙感數(shù)據(jù)和光學(xué)遙感數(shù)據(jù)二者在時空分辨率上的優(yōu)勢，將二者有效整合進行土壤含水量產(chǎn)品的降尺度研究，實現(xiàn)了流域尺度土壤含水量實時或準(zhǔn)實時的動態(tài)監(jiān)測。

以上所述，僅是本發(fā)明的較佳實施例而已，并非對本發(fā)明作任何形式上的限制，本領(lǐng)域技術(shù)人員利用上述揭示的技術(shù)內(nèi)容做出些許簡單修改、等同變化或修飾，均落在本發(fā)明的保護范圍內(nèi)。