利用MODIS探測(cè)華北平原淺層地下水水位變化的初步研究

劉 杰，何報(bào)寅，丁 超，張 文

（1.中國科學(xué)院 測(cè)量與地球物理研究所，武漢 430077；2.中國科學(xué)院 研究生院，北京 100049）

近年來，由于降水量減少；河道滲流補(bǔ)給減少；地下水開采量增大，華北平原形成了以石家莊－保定為代表的七個(gè)淺層地下水降落漏斗［1］.因此研究石家莊－保定地區(qū)地下水水位變化情況有著十分重要的意義.F.Alkhaier［2］等研究表明淺層地下水水位變化對(duì)地表溫度有著一定的影響；K.Sreenivas［3］等研究了白天地表溫度、夜間地表溫度、地表溫差等參數(shù)與淺層地下水水位的相關(guān)關(guān)系；付俊娥［4］等選用 NDVI，LST數(shù)據(jù)，應(yīng)用支持向量機(jī)回歸，建立了西北干旱地區(qū)地下水位遙感監(jiān)測(cè)模型；塔西甫拉提·特依拜［5］等對(duì)不同深度地溫與地下水埋深進(jìn)行了相關(guān)性分析；徐亮等［6］利用Visual Modflow模型模擬了洪湖地區(qū)地下水剖面的變化情況.地下水對(duì)地表溫度的影響主要是通過對(duì)淺層巖土含水量的影響來實(shí)現(xiàn)的，而淺層巖土含水量與地物的固有屬性熱慣量之間有著十分密切的聯(lián)系.

熱慣量的定義如下［7］：

其中，K 為熱傳導(dǎo)率（J·m－2·K－1·s－1）；ρ為物質(zhì)密度（kg·m－3）；C為比熱（J·m－2·K－1）.

熱慣量是不能被直接測(cè)量的.同時(shí)如果要根據(jù)熱傳導(dǎo)率、比熱、物質(zhì)密度等因子來計(jì)算熱慣量也是比較困難的.遙感技術(shù)作為獲取熱慣量值的一個(gè)重要手段，最普遍的應(yīng)用是利用熱慣量反演淺層巖土含水量.余濤等［8］、李星敏［9］等以衛(wèi)星遙感數(shù)據(jù)為數(shù)據(jù)源，建立了熱慣量與淺層巖土含水量之間的相關(guān)關(guān)系模型.因此可以通過研究熱慣量來研究地下水水位的變化情況，即探尋熱慣量與淺層地下水水位之間的關(guān)系.郭茜［10］等研究表明熱慣量和淺層巖土含水量之間存在著密切的一元線性相關(guān)關(guān)系.因此影響淺層巖土含水量動(dòng)態(tài)變化的因素也即影響熱慣量的因素.影響淺層巖土含水量的因素有很多，在降水量豐富的地區(qū)主要受到降水、河湖徑流、地下水等的影響.在干旱半干旱地區(qū)，淺層巖土含水量較小、降水量低、河湖較少，則淺層巖土含水量的年際變化主要受到地下水動(dòng)態(tài)變化、灌溉、土地利用變化等的影響.因此在干旱半干旱地區(qū)需要盡量降低降雨以及植被覆蓋對(duì)熱慣量的影響.

為了最大限度的降低降雨和植被覆蓋的影響，本文選取了相對(duì)于石家莊－保定地區(qū)而言，降雨量較低、植被覆蓋較少的時(shí)間段，主要集中在冬季.在此情況下，熱慣量年際變化的影響因素則是不斷波動(dòng)的淺層地下水埋深.

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

研究區(qū)（37°13′～39°10′N，114°22′～116°7′E）如圖1所示位于華北平原北部、河北省的南部，包括石家莊、保定、衡水3個(gè)地級(jí)市以及下轄的32個(gè)縣區(qū)，總面積2.07萬km2，占河北省面積的10.6%.該區(qū)屬于我國北方比較缺水的地區(qū)，人均水資源量只有450m3左右，是我國經(jīng)濟(jì)發(fā)展中受水資源制約最嚴(yán)重的地區(qū)之一［11］.

圖1 研究區(qū)位置圖Fig.1 Location of study area

圖2 觀測(cè)井位置圖Fig.2 Location of observation well

該研究區(qū)是華北平原典型的山前平原區(qū)，冬季降雨量少、植被覆蓋率低、地表水較少，以地下水為主要供水水源，且在研究區(qū)內(nèi)，已探明的有石家莊和保定兩個(gè)較大的地下水降落漏斗［1］.所以該區(qū)滿足實(shí)驗(yàn)的基本條件，即：半干旱地區(qū)、以開采淺層地下水為主要的供水水源、冬季降雨量少、植被覆蓋率低、研究區(qū)內(nèi)有較大的地下水漏斗區(qū).同時(shí)該研究區(qū)屬于華北平原的代表區(qū)域，研究該區(qū)域的地下水變化情況有著十分重大的意義.

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

MODIS傳感器分別搭載在Terra（上午星）和Aqua（下午星）衛(wèi)星上，有36個(gè)波段，其中1～20波段為反射波段、21～25和27～36為發(fā)射波段、26波段為反射發(fā)射并存波段［11］.MODIS數(shù)據(jù)1～2天可以獲得一次全球觀測(cè)數(shù)據(jù)（包括白天的可見光圖像及白天／夜間的紅外圖像）.我國北方的觀測(cè)時(shí)間一般為白天10∶30～12∶00，夜間21∶30～23∶00［13］.其優(yōu)勢(shì)是有較高的時(shí)間分辨率和適中的空間分辨率.

本文選取的MODIS數(shù)據(jù)有：MOD11A2（上午星地表溫度八天合成產(chǎn)品）、MYD11A2（下午星地表溫度八天合成產(chǎn)品）、MCD43B3（上下午星合成地表反照率參數(shù)Albedo八天合成產(chǎn)品）.對(duì)于以上3類產(chǎn)品采集了從2004年～2007年，每年11月1日～12月1日的數(shù)據(jù)，總計(jì)16個(gè)時(shí)相的數(shù)據(jù).由于研究區(qū)在四景圖像的交叉處，并且其投影是基于NASA的sinusoidal投影，所以需要對(duì)數(shù)據(jù)進(jìn)行拼接、重投影以及裁切.其中拼接和重投影是利用MRT在DOS環(huán)境下進(jìn)行批處理實(shí)現(xiàn)的.對(duì)于MOD11A2、MYD11A2、MCD43B3經(jīng)過重投影后為科學(xué)數(shù)據(jù)集，還需進(jìn)行輻射校正，將地表溫度產(chǎn)品中的LST＿Day＿1km（1km白天地表溫度）和LST＿Night＿1km（1km夜晚地表溫度）轉(zhuǎn)化為地表溫度（單位：攝氏度），Albedo轉(zhuǎn)換到0～1之間（無量綱）.

在此需要說明的是分別選取上午星的白天溫度產(chǎn)品和下午星的夜間溫度產(chǎn)品的原因有二：其一是因?yàn)?天合成溫度產(chǎn)品精度較高，免去了計(jì)算地表溫度的工作；其二是因?yàn)樯衔缧呛拖挛缧蔷邪滋鞙囟取⒁归g溫度產(chǎn)品，在對(duì)比這4類溫度產(chǎn)品同一時(shí)間同一地點(diǎn)的溫度值時(shí)發(fā)現(xiàn)，上午星的白天溫度產(chǎn)品高于下午星的白天溫度產(chǎn)品，下午星的夜間溫度產(chǎn)品低于上午星的夜間溫度產(chǎn)品，而計(jì)算過程中需要的是最高溫和最低溫，因此選擇了溫度值最高的上午星白天溫度產(chǎn)品和溫度值最低的下午星夜間溫度產(chǎn)品.

華北平原地下水水位變化可劃分為3個(gè)類型區(qū)，第一類型區(qū)：太行山、燕山山前平原地下水強(qiáng)烈開采區(qū)；第二類型區(qū)：黃河側(cè)滲、引黃灌溉地下水補(bǔ)給區(qū)；第三類型區(qū)：中部、東部平原水質(zhì)差、水量少區(qū)［14］.本文所研究的區(qū)域?qū)儆诘谝活愋蛥^(qū)，地下水位下降大，是地下水位降落漏斗的形成區(qū)，區(qū)內(nèi)的河道幾乎全部干涸.本文所選取的地下水?dāng)?shù)據(jù)是從中國地質(zhì)調(diào)查局獲取的2004年到2007年研究區(qū)內(nèi)5個(gè)地下水監(jiān)測(cè)點(diǎn)的潛水觀測(cè)井的地下水水位數(shù)據(jù).觀測(cè)井的具體位置如圖2所示.

其他數(shù)據(jù)：全國基礎(chǔ)地理信息中心所獲取的全國縣市行政區(qū)劃圖、水體分布圖、地級(jí)市分布圖等數(shù)據(jù).

2 熱慣量的計(jì)算

以傅立葉級(jí)數(shù)表達(dá)的溫度周日變化的熱學(xué)模式為邊界條件，解熱傳導(dǎo)方程［15］，得出如下關(guān)系：

C1為太陽赤緯和地理緯度的函數(shù)，表示為：

解上式得到熱慣量的近似方程：

其中，S0為太陽常數(shù)（1367W·m2），V 為大氣透過率（%），A為地表反照率（%），為地球自轉(zhuǎn)角頻率（rad／s），B為土壤發(fā)射率、比濕及氣象參數(shù)的函數(shù)，Thigh、Tlow分別為晝夜地表最高、最低溫度（K），二者之差用ΔT表示.

（3）式中第一項(xiàng)中所包含的A和ΔT可以通過遙感資料獲取，它被定義為表觀熱慣量（ATI）即［15］：

B為太陽赤緯、地理維度、日照時(shí)數(shù)、日地距離的函數(shù)，表示為：

ATI是無量綱的相對(duì)值，反映了熱慣量的相對(duì)大小.Price［16］等人在熱慣量模型中，假定大氣透過率V是常數(shù)（0.75），但是其實(shí)際情況因各地而異.S0VC1是到達(dá)地面的太陽總輻射.

地表反照率和地表溫差是影響熱慣量模型的兩個(gè)關(guān)鍵因素，也是可以從遙感數(shù)據(jù)中得到的參量.如果反照率保持不變，隨著溫差的增大，熱慣量將下降.如果地表溫差大于10K，則熱慣量隨著溫差的增大，熱慣量變化趨于緩慢.如果溫差恒定，熱慣量將隨著地表反照率的增加而減少.對(duì)于某一溫差而言，地表反照率變化幅度相同的情況下，其對(duì)熱慣量的影響不變［17］.

對(duì)于表觀熱慣量而言，最重要的兩個(gè)因素是地表反照率和溫差，為了便于計(jì)算，在只考慮地表反照率和溫差，不考慮地理維度、太陽赤緯、日照時(shí)數(shù)、日地距離的情況下，可將表觀熱慣量簡(jiǎn)化為［10］：

為了盡量降低降雨對(duì)熱慣量的影響，本文選取了一年中降雨較少的11月.由于利用的是MODIS八天合成產(chǎn)品，所以每年的11月均可以得到4個(gè)時(shí)間點(diǎn)的地表溫度以及反照率數(shù)據(jù)，對(duì)每年這4個(gè)時(shí)間點(diǎn)的表觀熱慣量進(jìn)行平均即得到了2004年～2007年每年11月的表觀熱慣量值.

需要說明的是在熱慣量的提取過程中取的是觀測(cè)井所在像元以及周圍8像元總計(jì)9像元的平均得到的平均表觀熱慣量值，這是為了降低由于重投影以及幾何校正誤差的影響.

3 結(jié)果和討論

如圖3所示為根據(jù)表觀熱慣量公式計(jì)算得到的表觀熱慣量2004年～2007年分布圖.從圖中可以發(fā)現(xiàn)，雖然地表水已大部分干涸，其河道在圖中仍然清晰可見.這是由于干涸的河道形成時(shí)間不長，河道的下墊面還沒有完全被周圍不同的土地覆蓋所同化.

表1為表觀熱慣量值與地下水水位對(duì)比表.為了探尋表觀熱慣量與地下水水位之間的相關(guān)關(guān)系，計(jì)算表1中各個(gè)觀測(cè)井表觀熱慣量與地下水的相關(guān)關(guān)系，得到表2.Person指數(shù)代表的是變量間的相關(guān)性，取值范圍為－1～1，其絕對(duì)值越接近1說明相關(guān)性越好.從表3可以發(fā)現(xiàn)藁城、晉州、定州3個(gè)點(diǎn)表觀熱慣量與地下水水位有較好的負(fù)相關(guān)性，而正定縣和涿州則表現(xiàn)出較低的相關(guān)性.分析原因可能是正定縣距河北省省會(huì)石家莊較近，而石家莊地區(qū)的大氣氣溶膠的不規(guī)律變化容易影響到該區(qū)域附近的地表反照率，從而使得表觀熱慣量呈現(xiàn)出不規(guī)律性；涿州觀測(cè)點(diǎn)離大清河較近，在大清河的影響范圍內(nèi)，容易受到地表水系的水量變化的影響，從而弱化了地下水水位變化對(duì)表 觀熱慣量的影響.

表1 觀測(cè)井表觀熱慣量與地下水水位對(duì)比表Tab.1 Comparison between apparent thermal inertia and water table in the place of observation well

圖3 石家莊保定地區(qū)2004年～2007年表觀熱慣量分布圖Fig.3 Apparent thermal inertia from 2004to 2007in Shijiazhuang－Baoding region

表2 表觀熱慣量與地下水水位相關(guān)性表Tab.2 Correlation coefficients between apparent thermal inertia and water table

為了得到表觀熱慣量與地下水水位變化的關(guān)系，選取相關(guān)性最好的定州觀測(cè)井進(jìn)行統(tǒng)計(jì)擬合，得到表觀熱慣量與地下水水位的簡(jiǎn)單統(tǒng)計(jì)關(guān)系模型，如下所示：

其中，得到的地下水水位是根據(jù)表觀熱慣量計(jì)算得到的相對(duì)地下水水位，此相對(duì)地下水水位并不能代表水位的絕對(duì)值，計(jì)算此相對(duì)水位的目的是為了得到地下水水位的年際變化情況，如圖4.

圖4為地下水水位變化圖，大于0表示地下水下降，小于0表示地下水上升.從圖4中可以發(fā)現(xiàn)石家莊東北部地區(qū)和石家莊的南部地區(qū)的地下水水位下降較快，這部分區(qū)域在這這四年內(nèi)下降了4～6m，速度十分驚人，這部分已探明為地下水漏斗區(qū)域，這說明石家莊地下水漏斗區(qū)域的地下水水位仍然在以較快的速率下降；同時(shí)從圖中可以看到地表水干涸河道影響區(qū)域地下水水位下降較慢，明顯將地下水下降較快區(qū)域隔離開；保定的南面地下水有略微的下降，下降幅度在0～1m范圍內(nèi)；保定東面地下水水位下降較慢，部分地區(qū)有輕微的上升；衡水的西南面地下水水位下降較慢，部分地區(qū)有輕微的上升；衡水東面地下水下降較快.

熱慣量方法監(jiān)測(cè)土壤含水量適用于裸土或有少量植被覆蓋的情況［9］，因此利用熱慣量監(jiān)測(cè)地下水水位變化應(yīng)當(dāng)降低植被覆蓋的影響.鄧輝［18］在利用遙感方法監(jiān)測(cè)土壤水人時(shí)認(rèn)為NDVI大于0.3的區(qū)域不適合用熱慣量法進(jìn)行監(jiān)測(cè)，而應(yīng)當(dāng)利用植被供水指數(shù)法.本文嘗試用NDVI＜0.3對(duì)研究區(qū)進(jìn)行掩膜處理，但是結(jié)果發(fā)現(xiàn)大部分區(qū)域均小于0.3，嘗試將閾值降低到0.2或者更低但是缺乏研究依據(jù).所以未對(duì)研究區(qū)進(jìn)行植被掩膜處理.在后續(xù)的工作中，為了更進(jìn)一步提高地下水監(jiān)測(cè)的精度將植被覆蓋的影響也考慮進(jìn)來，可以引入其他的遙感監(jiān)測(cè)方法，如前所提到的植被供水指數(shù)法、溫度植被指數(shù)法等進(jìn)行進(jìn)一步的研究.

圖4 2007年～2004年地下水水位變化分布圖（＞0表示地下水位下降，＜0表示上升）Fig.4 Differential chart of water table from 2004to 2007（＞0represent decline of water table down and＜0represent rising）

4 結(jié)論

利用表觀熱慣量研究地下水水位變化受到地表覆蓋類型的影響較大，適合于在裸地進(jìn)行研究，高植被覆蓋區(qū)地下水水位的研究應(yīng)當(dāng)引入其他的方法如植被供水指數(shù)法、溫度植被指數(shù)法等.

對(duì)表觀熱慣量而言，地表反照率與地表溫差是兩個(gè)最主要的因素.對(duì)于城市地區(qū)而言，由于氣溶膠受到人為干擾較大從而使得地表反照率的變化波動(dòng)性很強(qiáng)，因此對(duì)于城市地區(qū)而言很難通過表觀熱慣量法來研究地下水水位的變化.

不同的地下水水位在地表反映出不同的熱學(xué)特性，表觀熱慣量作為地物的固有屬性可以作為一個(gè)參量進(jìn)行研究，同時(shí)還有許多其他的熱學(xué)參量來反映這種變化特征，如地表最高溫、地表最低溫、地表溫差等.這些不同的參量與地下水水位變化的關(guān)系將在后續(xù)的研究中進(jìn)行.

［1］王仕琴，宋獻(xiàn)方，王勤學(xué).華北平原淺層地下水水位動(dòng)態(tài)變化［J］.地理學(xué)報(bào)，2008，63（5）：462－472.

［2］Alkhaier F，Schotting R J，Su Z.A qualitative description of shallow groundwater effect on surface temperature of bare soil［J］.Hydrol Earth Syst Sci，2009，12：1749－1756.

［3］Screenivas K，Dwivedi R S，Singh A N，et al.Detection of Sub－Surface Waterlogging using Terra－1MODIS Data［J］.J Indian Soc Remote Sens，2010，38∶119－132.

［4］付俊娥，蘇慶珣，潘世兵，等.基于支持向量機(jī)理論的地下水動(dòng)態(tài)遙感監(jiān)測(cè)模型與應(yīng)用［J］.地球信息科學(xué)學(xué)報(bào)，2010，12（4）：466－472.

［5］塔西甫拉提·特依拜，崔建永，丁建麗.熱紅外遙感技術(shù)探測(cè)干旱區(qū)綠洲－荒漠交錯(cuò)帶地下水分布的方法研究［J］.干旱區(qū)地理，2005，28（2）：252－257.

［6］徐 亮，陳世儉，向先富.洪湖地區(qū)淺層承壓水動(dòng)態(tài)模擬研究［J］.華中師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2008，42（2）：299－303.

［7］張霄羽，畢于運(yùn)，李召良.遙感估算熱慣量研究回顧與進(jìn)展［J］.地理科學(xué)進(jìn)展，2008，27（3）：166－172.

［8］余 濤，田國良.熱慣量法在監(jiān)測(cè)表層淺層巖土含水量變化中的研究［J］.遙感學(xué)報(bào)，1（1）：24－31.

［9］李星敏，劉安麟，張樹譽(yù).熱慣量法在干旱遙感監(jiān)測(cè)中的應(yīng)用研究［J］.干旱地區(qū)農(nóng)業(yè)研究，2005，23（1）：54－59.

［10］郭 茜，李國春.用表觀熱慣量法計(jì)算淺層巖土含水量探討［J］.中國農(nóng)業(yè)氣象，2005，26（4）：215－219.

［11］賈金生，劉昌明.華北平原地下水動(dòng)態(tài)及其對(duì)不同開采量響應(yīng)的計(jì)算：以河北省欒城縣為例［J］.地理學(xué)報(bào)，2002，57（2）：201－209.

［12］劉榮高，劉 洋，劉紀(jì)遠(yuǎn).MODI科學(xué)數(shù)據(jù)處理研究進(jìn)展［J］.自然科學(xué)進(jìn)展，2009，19（2）：147－147.

［13］趙英時(shí).遙感應(yīng)用分析原理與方法［M］.北京：科學(xué)出版社，2003.

［14］張兆吉，費(fèi)宇仁，陳宗宇，等.華北平原地下水可持續(xù)利用調(diào)查評(píng)價(jià)［M］.北京：地質(zhì)出版社，2009.

［15］Price J C.On the analysis of thermal infrared imagery：The limited utility of apparent thermal inertia［J］.Remote Sensing of Environment，1985，18：59－73.

［16］Price J C.Thermal Inertia Mapping：A New View of the Earth［J］.J Geophys Res，1977，82：2582－2590.

［17］蔡國印.基于MODIS數(shù)據(jù)的地表溫度、熱慣量反演研究及其在淺層巖土含水量、地氣間熱交換方面的應(yīng)用［D］.北京：中國科學(xué)院研究生院，2006.

［18］鄧 輝.基于MODIS數(shù)據(jù)的大區(qū)域土壤水分遙感監(jiān)測(cè)研究 ［D］.北京：中國農(nóng)業(yè)科學(xué)院研究生院，2004.