基于Landsat8和Sentinel-1A數(shù)據(jù)的焦作廣利灌區(qū)夏玉米土壤墑情監(jiān)測(cè)方法研究

郭二旺，郭乙霏，羅蔚然，王文婷

(1.焦作市水利勘測(cè)設(shè)計(jì)院，河南 焦作 454003；2.東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院，哈爾濱 150003；3.鄭州大學(xué)水利與環(huán)境學(xué)院，鄭州 450001；4.焦作市抗旱防汛通訊站，河南 焦作 454150)

0 引 言

土壤墑情是水文學(xué)、氣候?qū)W、農(nóng)業(yè)學(xué)、生態(tài)環(huán)境學(xué)等各個(gè)研究領(lǐng)域中的重要指標(biāo)，在生態(tài)系統(tǒng)水循環(huán)中扮演著非常重要的角色，是農(nóng)作物旱情監(jiān)測(cè)、估產(chǎn)等研究的重要參數(shù)[1-4]。土壤墑情在時(shí)間、空間范圍上變化較大，傳統(tǒng)的土壤墑情監(jiān)測(cè)方法是通過設(shè)立監(jiān)測(cè)站點(diǎn)進(jìn)行人工監(jiān)測(cè)，不能滿足較大范圍上對(duì)土壤墑情進(jìn)行高精度、實(shí)時(shí)監(jiān)測(cè)的要求[5]。隨著微波遙感理論基礎(chǔ)經(jīng)過多年不斷的發(fā)展，利用微波數(shù)據(jù)對(duì)土壤水分監(jiān)測(cè)的研究已取得了突破性進(jìn)展，特別是合成孔徑雷達(dá)(Synthetic Aperture Radar, 簡(jiǎn)稱SAR)技術(shù)不斷的改進(jìn)，以其多頻率、多極化(或者全極化)、多角度、可變工作模式等能力，在土壤水分研究中變得日益重要[6]。根據(jù)微波遙感反演對(duì)象不同可分為裸露地表區(qū)和植被覆蓋區(qū)。國內(nèi)外研究現(xiàn)狀將針對(duì)國內(nèi)外裸露地表的反演和植被覆蓋區(qū)的反演進(jìn)行陳述。對(duì)于裸露地表土壤墑情反演，F(xiàn)ung[7]等構(gòu)建了IEM 模型(integral equation model)。IEM 模型是基于電磁波輻射傳輸方程的地表散射模型，能在一個(gè)很寬的地表粗糙度范圍內(nèi)再現(xiàn)真實(shí)地表后向散射情況，已經(jīng)被廣泛應(yīng)用于微波地表散射、輻射的模擬和分析。Chen 等[8]對(duì) IEM 模型存在的不足之處進(jìn)行分析，從而提出AIEM 模型。對(duì)于微波遙感反演植被覆蓋區(qū)的研究中，Ulaby的貢獻(xiàn)是毋庸置疑的，其主要突出貢獻(xiàn)為反演植被覆蓋區(qū)域奠定了扎實(shí)的基礎(chǔ)：一是水云模型[9](Water-Cloud model)，該模型是Ulaby和Attema等以植被覆蓋區(qū)作為研究對(duì)象，通過簡(jiǎn)化植被層的散射機(jī)制，得到的適用于低矮植被覆蓋地表的模型；一是密歇根微波植被散射模型(Michigan Microwave Canopy Scattering model, MIMICS)[10]，該模型是 Ulaby等人于 1990 年基于微波輻射傳輸式子的一階解建立的描述植被覆蓋地區(qū)散射特性的模型。

本研究采用Sentinel-1A作為雷達(dá)遙感監(jiān)測(cè)的數(shù)據(jù)源，結(jié)合Landsat8數(shù)據(jù)來進(jìn)行土壤墑情反演，基于改進(jìn)的水云模型研究植被的后向散射作用，嘗試消除植被效應(yīng)得到地表后向散射系數(shù)，根據(jù)地表反射率和地表參數(shù)關(guān)系實(shí)現(xiàn)研究區(qū)土壤墑情的反演。

1 研究區(qū)概況與數(shù)據(jù)選取

1.1 研究區(qū)簡(jiǎn)介

(1)地理位置。廣利灌區(qū)位于河南省焦作市，地理位置為東經(jīng)112°37′～113°13′，北緯34°55′～35°11′，灌區(qū)范圍在濟(jì)源市、沁陽市沁南地區(qū)、溫縣城、武陟沁南幾個(gè)地區(qū)都有覆蓋。灌區(qū)的南面是黃河，北面挨著沁河，屬于黃河流域，引沁河水灌溉。灌區(qū)總面積3.4 萬hm2，如圖1所示。灌區(qū)內(nèi)主要種植小麥、玉米等，夏玉米的播種期為5月，5月中到7月為幼苗期，9月底成熟。

圖1 廣利灌區(qū)總體布置圖Fig.1 Distribution map of guangli irrigation area

(2)水文特征。規(guī)劃面積34 000 hm2，其中灌溉面積20 666 hm2、補(bǔ)源面積13 333 hm2，設(shè)計(jì)流量28.0 m3/s?，F(xiàn)有總干渠1條，長29.1 km；灌溉干支渠34條，長214.25 km；補(bǔ)水干支渠19條，長134.06 km。灌區(qū)多年平均降水量590.66 mm，降水量年內(nèi)分配極不均勻，多集中在汛期夏秋兩季之間，其中汛期6-9月份降水量占全年降水量的70%以上。由于受到降雨的影響，地下水位從4月中下旬開始由3 m以下逐漸上升，到9月中旬地下水位逐漸降至3 m以下。年內(nèi)蒸發(fā)量的變化基本與氣溫變化一致，變化規(guī)律呈正態(tài)分布，一年中以6月份最大，平均306.2 mm，11月份最小，平均86.7 mm。

1.2 遙感數(shù)據(jù)選取

本文采用的是2017年7月的Sentinel-1A數(shù)據(jù)和Landsat8數(shù)據(jù)，用SNAP軟件對(duì)Sentinel-1A影像的VV和VH影像進(jìn)行輻射定標(biāo)，將影像上的像元亮度值轉(zhuǎn)化為后向散射系數(shù)；對(duì)得到的結(jié)果去噪聲處理后，做投影轉(zhuǎn)換和坐標(biāo)系定義，將該影像設(shè)置為和Landsat8影像相同的投影坐標(biāo)系，即UTM投影、WGS-84坐標(biāo)系。研究中所用到的是Landsat8的OLI傳感器影像，對(duì)分辨率30 m的Landsat8數(shù)據(jù)利用廣利灌區(qū)的矢量邊界裁剪得到研究范圍的影像數(shù)據(jù)；利用ENVI5.2軟件對(duì)得到的研究區(qū)范圍的Landsat8影像進(jìn)行輻射定標(biāo)和Flaash大氣校正，并利用遙感影像計(jì)算NDVI，根據(jù)閾值分割法提取夏玉米分布，結(jié)果如圖2所示。

圖2 夏玉米空間分布圖Fig.2 Spatial distribution map of summer corn

1.3 實(shí)測(cè)野外樣品數(shù)據(jù)

選擇衛(wèi)星過境的同一時(shí)間進(jìn)行實(shí)地采樣，使用GPS記錄樣點(diǎn)坐標(biāo)。本研究選取2017 年7 月15 日的位于廣利灌區(qū)的24個(gè)實(shí)測(cè)點(diǎn)的土壤平均含水量數(shù)據(jù)、植被含水量數(shù)據(jù)及地表粗糙度等數(shù)據(jù)，圖3所示。①土壤墑情數(shù)據(jù) ：在焦作市廣利灌區(qū)進(jìn)行實(shí)測(cè)土壤墑情時(shí)，采用方法為土鉆取土，然后用烘箱進(jìn)行烘干。廣利灌區(qū)的土壤為粉質(zhì)壤土，土質(zhì)分布均勻，容重為1.40 g/m3，田間持水率為28%(重量含水率)。在野外采樣的過程中，每個(gè)采樣點(diǎn)取深度為0～10 cm處的土壤，用鋁盒收集樣品并及時(shí)記錄總重量，對(duì)樣本進(jìn)行烘干處理獲得烘干后的樣本鋁盒的重量并記錄，計(jì)算土壤含水量數(shù)據(jù)。圖3是研究區(qū)野外采樣點(diǎn)的分布情況。②植被含水量數(shù)據(jù)：測(cè)量采樣點(diǎn)植被的高度，葉片長度、寬度，分別記錄采樣點(diǎn)植株葉和根烘干前和烘干后的重量，計(jì)算鮮重含水量和干重含水量。

圖3 研究區(qū)樣點(diǎn)分布Fig.3 Sample point distribution in research area

(3)地表粗糙度:本研究中的粗糙度數(shù)據(jù)是采用剖面板測(cè)量法測(cè)量采樣點(diǎn)的數(shù)據(jù)，隨后利用離散數(shù)據(jù)的均方差高度以及表面相關(guān)長度計(jì)算公式來計(jì)算S(均方根高度)和CL(相關(guān)長度)，對(duì)于一維離散數(shù)據(jù)，表面高度的標(biāo)準(zhǔn)偏差(均方根高度)s為：

在針對(duì)城鄉(xiāng)進(jìn)行規(guī)劃和具體設(shè)計(jì)的時(shí)候，為了保證城鄉(xiāng)規(guī)劃的科學(xué)性和合理性，首先要做的一點(diǎn)就是要樹立具有生態(tài)文明特征的規(guī)劃理念，并且在實(shí)踐中可以將該理念落實(shí)到實(shí)處。也就是說在實(shí)踐中，要采取有針對(duì)性的措施，盡可能保證遵循人與自然環(huán)境相互協(xié)調(diào)的基礎(chǔ)發(fā)展原則。首先，要采取有針對(duì)性的措施，盡可能降低CO2的排放力度，同時(shí)還要使用一些具有清潔性功能特征的能源。這樣不僅能夠針對(duì)城鄉(xiāng)規(guī)劃相關(guān)措施具體落實(shí)過程中的一些污染問題進(jìn)行有針對(duì)性的處理，而且還能夠保證生態(tài)文明城市理念在實(shí)踐中的有效落實(shí)。

式中：N為采樣數(shù)目。

對(duì)于離散數(shù)據(jù)，相距x′=(j-1)x(j為自然數(shù))的歸一化自相關(guān)函數(shù)由下式給出：

當(dāng)相關(guān)函數(shù)ρ(x′)=1/e時(shí),間隔x′稱為表面相關(guān)長度。

為了在一定程度上避免采樣誤差，在進(jìn)行樣本采集時(shí)，對(duì)于每一個(gè)樣點(diǎn)附近分3個(gè)子樣點(diǎn)，將子樣本點(diǎn)的均值作為該點(diǎn)土壤含水量，同時(shí)考慮遙感影像的與實(shí)測(cè)數(shù)據(jù)之間的差異性,將測(cè)點(diǎn)對(duì)應(yīng)的影像點(diǎn)的鄰近像元值進(jìn)行均值化處理，使得采樣點(diǎn)與影像點(diǎn)能夠更好的匹配。

2 植被覆蓋地表土壤墑情反演模型

微波傳感器接收的回波信號(hào)中包含了土壤水分和植被水分的信息，地表植被對(duì)微波遙感反演土壤水分的精度有較大的影響。因此在進(jìn)行土壤墑情反演時(shí)，消除植被的影響是研究的重點(diǎn)。

2.1 植被散射模型

Attema和Ulaby等在1987年提出了以農(nóng)作物為研究對(duì)象的水云模型。在水云模型中，研究區(qū)植被層被理想化為云層一樣均勻分布在地表上方，只考慮植被散射和地表散射而忽略了植被和土壤表面的多次散射?？捎孟率奖硎荆?/p>

γ2(θ)=exp[-2Bmv/cos(θ)]

式中：γ2(θ)為植被雙層衰減因子；θ指入射角；A和B的值取決植被類型及入射電磁波頻率,A=0.001 2，B=0.091；mv表示植被含水量。

本文將采取歸一化水指數(shù)(NDWI)來估算植被含水量，Jackson(Jacson et al.，2004)[11]等人的研究證明了NDWI估算植被指數(shù)比NDVI有優(yōu)勢(shì)。 公式如下：

從實(shí)測(cè)的植被含水量數(shù)據(jù)中選取一部分來建立NDWI和植被含水量的關(guān)系，如圖4所示：

y=-49.027x+4.413 1

(5)

利用公式(5)經(jīng)波段運(yùn)算得到的研究區(qū)植被含水量分布。根據(jù)實(shí)測(cè)植被含水量數(shù)據(jù)檢測(cè)NDWI估算的VWC值的準(zhǔn)確度，結(jié)果表明模擬植被含水量和實(shí)測(cè)值R2=0.535 7，可以用來估算植被含水量。

圖4 植被含水量和NDWI之間的關(guān)系圖Fig4.Map of the relationship between vegetation water content and NDWI

根據(jù)公式(1)～(3)運(yùn)用ENVI5.2中的波段運(yùn)算結(jié)合Sentinel-1A影像和Landsat8影像可以計(jì)算去除植被影響的后向散射。利用AIEM模型，研究在VV極化下地表土壤墑情對(duì)后向散射的影響。設(shè)置AIEM模型中的參數(shù)：雷達(dá)入射頻率采用Sentinel1A固有頻率f=5.4 GHz，VV極化，相關(guān)長度cl取值范圍5～20 cm，均方根高度s取值范圍0.3～3.3 cm，由于研究區(qū)地勢(shì)比較平坦，雷達(dá)入射角波動(dòng)范圍較小，整體入射角在θ=40°左右變化，因此在參數(shù)模擬過程中，入射角固定為40°。模擬結(jié)果如圖5所示。

圖5 不同地表粗糙度下土壤含水量對(duì)后向散射系數(shù)的影響Fig.5 Effect of soil water content on backscattering coefficient under different surface roughness

2.2 聯(lián)合地表參數(shù)反演土壤墑情

利用AIEM模型模擬的過程中需要輸入均方根高度、相關(guān)長度、入射角、介電常數(shù)等參數(shù)，采用Zribi[12]等提出的組合參數(shù)(Zs=s2/cl)來表征地表粗糙度參數(shù)。其中，用土壤墑情Mv和地表粗糙度Zs的函數(shù)來表示同極化雷達(dá)后向散射系數(shù)[13]：

聯(lián)合以上公式，組合粗糙度的反演公式為：

結(jié)合土壤水分和土壤粗糙度參數(shù)的式子：

得到結(jié)果：

圖6 實(shí)測(cè)值和模擬值關(guān)系曲線Fig.6 Curve of relationship between measured and simulated values

3 反演結(jié)果與精度分析

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

利用回歸公式對(duì)交叉極化組合反演的結(jié)果進(jìn)行計(jì)算，得到VV-VH極化下的土壤墑情空間分布圖，如圖7所示。在圖中可以看到土壤水分大多集中在20%～30%之間，并且土壤水分在空間上的分布也比較均勻。

圖7 交叉極化VV-VH的土壤水分反演圖Fig.7 Soil moisture inversion map of cross-polarized VV-VH

3.2 精度分析

對(duì)模擬得到的土壤水分?jǐn)?shù)據(jù)選取9個(gè)野外實(shí)測(cè)數(shù)據(jù)進(jìn)行精度驗(yàn)證，圖8是交叉極化VV-VH組合反演的土壤水分值與野外實(shí)測(cè)值的擬合關(guān)系，相關(guān)系數(shù)為0.7。

圖8 模擬數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的相關(guān)關(guān)系Fig.8 The correlation between simulated data and measured data

通過模擬數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)的相關(guān)關(guān)系可以得知本研究中使用的基于水云模型、AIEM模型和組合粗糙度參數(shù)建立的反演模型在廣利灌區(qū)有一定的適用性。但由于實(shí)測(cè)數(shù)據(jù)數(shù)量的限制和模型本身存在的問題使得精度不是很高，得到的土壤墑情分布圖有一定的參考作用。

4 結(jié) 語

本文研究區(qū)是河南省焦作市的廣利灌區(qū)，該研究區(qū)的植被類型主要是玉米、大豆、小麥等低矮植被。使用的數(shù)據(jù)是Sentinel-1A和Landsat8遙感數(shù)據(jù)，結(jié)合光學(xué)遙感和微波遙感的特性，考慮到灌區(qū)植被的生長和分布狀況，選用水云模型去除地表植被對(duì)后向散射的影響。研究利用水云模型得到地表后向散射系數(shù)之后，通過VV、VH交叉極化組合方式得到交叉極化差，結(jié)合組合粗糙度參數(shù)建立研究區(qū)的土壤水分反演模型，結(jié)果較為滿意。本文研究了植被體散射之間單次散射關(guān)系，并未考慮土壤--植被層之間的多次交互散射的影響。為進(jìn)一步提高土壤水分反演精度，應(yīng)考慮實(shí)際植被枝葉形態(tài)、朝向、體積以及介電常數(shù)等因素，進(jìn)而來獲取更高的反演精度。