利用MODIS探测华北平原浅层地下水水位变化的初步研究

刘 杰，何报寅，丁 超，张 文

（1.中国科学院 测量与地球物理研究所，武汉 430077；2.中国科学院 研究生院，北京 100049）

近年来，由于降水量减少；河道渗流补给减少；地下水开采量增大，华北平原形成了以石家庄－保定为代表的七个浅层地下水降落漏斗［1］.因此研究石家庄－保定地区地下水水位变化情况有着十分重要的意义.F.Alkhaier［2］等研究表明浅层地下水水位变化对地表温度有着一定的影响；K.Sreenivas［3］等研究了白天地表温度、夜间地表温度、地表温差等参数与浅层地下水水位的相关关系；付俊娥［4］等选用 NDVI，LST数据，应用支持向量机回归，建立了西北干旱地区地下水位遥感监测模型；塔西甫拉提·特依拜［5］等对不同深度地温与地下水埋深进行了相关性分析；徐亮等［6］利用Visual Modflow模型模拟了洪湖地区地下水剖面的变化情况.地下水对地表温度的影响主要是通过对浅层岩土含水量的影响来实现的，而浅层岩土含水量与地物的固有属性热惯量之间有着十分密切的联系.

热惯量的定义如下［7］：

其中，K 为热传导率（J·m－2·K－1·s－1）；ρ为物质密度（kg·m－3）；C为比热（J·m－2·K－1）.

热惯量是不能被直接测量的.同时如果要根据热传导率、比热、物质密度等因子来计算热惯量也是比较困难的.遥感技术作为获取热惯量值的一个重要手段，最普遍的应用是利用热惯量反演浅层岩土含水量.余涛等［8］、李星敏［9］等以卫星遥感数据为数据源，建立了热惯量与浅层岩土含水量之间的相关关系模型.因此可以通过研究热惯量来研究地下水水位的变化情况，即探寻热惯量与浅层地下水水位之间的关系.郭茜［10］等研究表明热惯量和浅层岩土含水量之间存在着密切的一元线性相关关系.因此影响浅层岩土含水量动态变化的因素也即影响热惯量的因素.影响浅层岩土含水量的因素有很多，在降水量丰富的地区主要受到降水、河湖径流、地下水等的影响.在干旱半干旱地区，浅层岩土含水量较小、降水量低、河湖较少，则浅层岩土含水量的年际变化主要受到地下水动态变化、灌溉、土地利用变化等的影响.因此在干旱半干旱地区需要尽量降低降雨以及植被覆盖对热惯量的影响.

为了最大限度的降低降雨和植被覆盖的影响，本文选取了相对于石家庄－保定地区而言，降雨量较低、植被覆盖较少的时间段，主要集中在冬季.在此情况下，热惯量年际变化的影响因素则是不断波动的浅层地下水埋深.

1 研究区与数据源

1.1 研究区概况

研究区（37°13′～39°10′N，114°22′～116°7′E）如图1所示位于华北平原北部、河北省的南部，包括石家庄、保定、衡水3个地级市以及下辖的32个县区，总面积2.07万km2，占河北省面积的10.6%.该区属于我国北方比较缺水的地区，人均水资源量只有450m3左右，是我国经济发展中受水资源制约最严重的地区之一［11］.

图1 研究区位置图Fig.1 Location of study area

图2 观测井位置图Fig.2 Location of observation well

该研究区是华北平原典型的山前平原区，冬季降雨量少、植被覆盖率低、地表水较少，以地下水为主要供水水源，且在研究区内，已探明的有石家庄和保定两个较大的地下水降落漏斗［1］.所以该区满足实验的基本条件，即：半干旱地区、以开采浅层地下水为主要的供水水源、冬季降雨量少、植被覆盖率低、研究区内有较大的地下水漏斗区.同时该研究区属于华北平原的代表区域，研究该区域的地下水变化情况有着十分重大的意义.

1.2 数据来源

MODIS传感器分别搭载在Terra（上午星）和Aqua（下午星）卫星上，有36个波段，其中1～20波段为反射波段、21～25和27～36为发射波段、26波段为反射发射并存波段［11］.MODIS数据1～2天可以获得一次全球观测数据（包括白天的可见光图像及白天／夜间的红外图像）.我国北方的观测时间一般为白天10∶30～12∶00，夜间21∶30～23∶00［13］.其优势是有较高的时间分辨率和适中的空间分辨率.

本文选取的MODIS数据有：MOD11A2（上午星地表温度八天合成产品）、MYD11A2（下午星地表温度八天合成产品）、MCD43B3（上下午星合成地表反照率参数Albedo八天合成产品）.对于以上3类产品采集了从2004年～2007年，每年11月1日～12月1日的数据，总计16个时相的数据.由于研究区在四景图像的交叉处，并且其投影是基于NASA的sinusoidal投影，所以需要对数据进行拼接、重投影以及裁切.其中拼接和重投影是利用MRT在DOS环境下进行批处理实现的.对于MOD11A2、MYD11A2、MCD43B3经过重投影后为科学数据集，还需进行辐射校正，将地表温度产品中的LST＿Day＿1km（1km白天地表温度）和LST＿Night＿1km（1km夜晚地表温度）转化为地表温度（单位：摄氏度），Albedo转换到0～1之间（无量纲）.

在此需要说明的是分别选取上午星的白天温度产品和下午星的夜间温度产品的原因有二：其一是因为8天合成温度产品精度较高，免去了计算地表温度的工作；其二是因为上午星和下午星均有白天温度、夜间温度产品，在对比这4类温度产品同一时间同一地点的温度值时发现，上午星的白天温度产品高于下午星的白天温度产品，下午星的夜间温度产品低于上午星的夜间温度产品，而计算过程中需要的是最高温和最低温，因此选择了温度值最高的上午星白天温度产品和温度值最低的下午星夜间温度产品.

华北平原地下水水位变化可划分为3个类型区，第一类型区：太行山、燕山山前平原地下水强烈开采区；第二类型区：黄河侧渗、引黄灌溉地下水补给区；第三类型区：中部、东部平原水质差、水量少区［14］.本文所研究的区域属于第一类型区，地下水位下降大，是地下水位降落漏斗的形成区，区内的河道几乎全部干涸.本文所选取的地下水数据是从中国地质调查局获取的2004年到2007年研究区内5个地下水监测点的潜水观测井的地下水水位数据.观测井的具体位置如图2所示.

其他数据：全国基础地理信息中心所获取的全国县市行政区划图、水体分布图、地级市分布图等数据.

2 热惯量的计算

以傅立叶级数表达的温度周日变化的热学模式为边界条件，解热传导方程［15］，得出如下关系：

C1为太阳赤纬和地理纬度的函数，表示为：

解上式得到热惯量的近似方程：

其中，S0为太阳常数（1367W·m2），V 为大气透过率（%），A为地表反照率（%），为地球自转角频率（rad／s），B为土壤发射率、比湿及气象参数的函数，Thigh、Tlow分别为昼夜地表最高、最低温度（K），二者之差用ΔT表示.

（3）式中第一项中所包含的A和ΔT可以通过遥感资料获取，它被定义为表观热惯量（ATI）即［15］：

B为太阳赤纬、地理维度、日照时数、日地距离的函数，表示为：

ATI是无量纲的相对值，反映了热惯量的相对大小.Price［16］等人在热惯量模型中，假定大气透过率V是常数（0.75），但是其实际情况因各地而异.S0VC1是到达地面的太阳总辐射.

地表反照率和地表温差是影响热惯量模型的两个关键因素，也是可以从遥感数据中得到的参量.如果反照率保持不变，随着温差的增大，热惯量将下降.如果地表温差大于10K，则热惯量随着温差的增大，热惯量变化趋于缓慢.如果温差恒定，热惯量将随着地表反照率的增加而减少.对于某一温差而言，地表反照率变化幅度相同的情况下，其对热惯量的影响不变［17］.

对于表观热惯量而言，最重要的两个因素是地表反照率和温差，为了便于计算，在只考虑地表反照率和温差，不考虑地理维度、太阳赤纬、日照时数、日地距离的情况下，可将表观热惯量简化为［10］：

为了尽量降低降雨对热惯量的影响，本文选取了一年中降雨较少的11月.由于利用的是MODIS八天合成产品，所以每年的11月均可以得到4个时间点的地表温度以及反照率数据，对每年这4个时间点的表观热惯量进行平均即得到了2004年～2007年每年11月的表观热惯量值.

需要说明的是在热惯量的提取过程中取的是观测井所在像元以及周围8像元总计9像元的平均得到的平均表观热惯量值，这是为了降低由于重投影以及几何校正误差的影响.

3 结果和讨论

如图3所示为根据表观热惯量公式计算得到的表观热惯量2004年～2007年分布图.从图中可以发现，虽然地表水已大部分干涸，其河道在图中仍然清晰可见.这是由于干涸的河道形成时间不长，河道的下垫面还没有完全被周围不同的土地覆盖所同化.

表1为表观热惯量值与地下水水位对比表.为了探寻表观热惯量与地下水水位之间的相关关系，计算表1中各个观测井表观热惯量与地下水的相关关系，得到表2.Person指数代表的是变量间的相关性，取值范围为－1～1，其绝对值越接近1说明相关性越好.从表3可以发现藁城、晋州、定州3个点表观热惯量与地下水水位有较好的负相关性，而正定县和涿州则表现出较低的相关性.分析原因可能是正定县距河北省省会石家庄较近，而石家庄地区的大气气溶胶的不规律变化容易影响到该区域附近的地表反照率，从而使得表观热惯量呈现出不规律性；涿州观测点离大清河较近，在大清河的影响范围内，容易受到地表水系的水量变化的影响，从而弱化了地下水水位变化对表 观热惯量的影响.

表1 观测井表观热惯量与地下水水位对比表Tab.1 Comparison between apparent thermal inertia and water table in the place of observation well

图3 石家庄保定地区2004年～2007年表观热惯量分布图Fig.3 Apparent thermal inertia from 2004to 2007in Shijiazhuang－Baoding region

表2 表观热惯量与地下水水位相关性表Tab.2 Correlation coefficients between apparent thermal inertia and water table

为了得到表观热惯量与地下水水位变化的关系，选取相关性最好的定州观测井进行统计拟合，得到表观热惯量与地下水水位的简单统计关系模型，如下所示：

其中，得到的地下水水位是根据表观热惯量计算得到的相对地下水水位，此相对地下水水位并不能代表水位的绝对值，计算此相对水位的目的是为了得到地下水水位的年际变化情况，如图4.

图4为地下水水位变化图，大于0表示地下水下降，小于0表示地下水上升.从图4中可以发现石家庄东北部地区和石家庄的南部地区的地下水水位下降较快，这部分区域在这这四年内下降了4～6m，速度十分惊人，这部分已探明为地下水漏斗区域，这说明石家庄地下水漏斗区域的地下水水位仍然在以较快的速率下降；同时从图中可以看到地表水干涸河道影响区域地下水水位下降较慢，明显将地下水下降较快区域隔离开；保定的南面地下水有略微的下降，下降幅度在0～1m范围内；保定东面地下水水位下降较慢，部分地区有轻微的上升；衡水的西南面地下水水位下降较慢，部分地区有轻微的上升；衡水东面地下水下降较快.

热惯量方法监测土壤含水量适用于裸土或有少量植被覆盖的情况［9］，因此利用热惯量监测地下水水位变化应当降低植被覆盖的影响.邓辉［18］在利用遥感方法监测土壤水人时认为NDVI大于0.3的区域不适合用热惯量法进行监测，而应当利用植被供水指数法.本文尝试用NDVI＜0.3对研究区进行掩膜处理，但是结果发现大部分区域均小于0.3，尝试将阈值降低到0.2或者更低但是缺乏研究依据.所以未对研究区进行植被掩膜处理.在后续的工作中，为了更进一步提高地下水监测的精度将植被覆盖的影响也考虑进来，可以引入其他的遥感监测方法，如前所提到的植被供水指数法、温度植被指数法等进行进一步的研究.

图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 结论

利用表观热惯量研究地下水水位变化受到地表覆盖类型的影响较大，适合于在裸地进行研究，高植被覆盖区地下水水位的研究应当引入其他的方法如植被供水指数法、温度植被指数法等.

对表观热惯量而言，地表反照率与地表温差是两个最主要的因素.对于城市地区而言，由于气溶胶受到人为干扰较大从而使得地表反照率的变化波动性很强，因此对于城市地区而言很难通过表观热惯量法来研究地下水水位的变化.

不同的地下水水位在地表反映出不同的热学特性，表观热惯量作为地物的固有属性可以作为一个参量进行研究，同时还有许多其他的热学参量来反映这种变化特征，如地表最高温、地表最低温、地表温差等.这些不同的参量与地下水水位变化的关系将在后续的研究中进行.

［1］王仕琴，宋献方，王勤学.华北平原浅层地下水水位动态变化［J］.地理学报，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］.地球信息科学学报，2010，12（4）：466－472.

［5］塔西甫拉提·特依拜，崔建永，丁建丽.热红外遥感技术探测干旱区绿洲－荒漠交错带地下水分布的方法研究［J］.干旱区地理，2005，28（2）：252－257.

［6］徐 亮，陈世俭，向先富.洪湖地区浅层承压水动态模拟研究［J］.华中师范大学学报：自然科学版，2008，42（2）：299－303.

［7］张霄羽，毕于运，李召良.遥感估算热惯量研究回顾与进展［J］.地理科学进展，2008，27（3）：166－172.

［8］余 涛，田国良.热惯量法在监测表层浅层岩土含水量变化中的研究［J］.遥感学报，1（1）：24－31.

［9］李星敏，刘安麟，张树誉.热惯量法在干旱遥感监测中的应用研究［J］.干旱地区农业研究，2005，23（1）：54－59.

［10］郭 茜，李国春.用表观热惯量法计算浅层岩土含水量探讨［J］.中国农业气象，2005，26（4）：215－219.

［11］贾金生，刘昌明.华北平原地下水动态及其对不同开采量响应的计算：以河北省栾城县为例［J］.地理学报，2002，57（2）：201－209.

［12］刘荣高，刘 洋，刘纪远.MODI科学数据处理研究进展［J］.自然科学进展，2009，19（2）：147－147.

［13］赵英时.遥感应用分析原理与方法［M］.北京：科学出版社，2003.

［14］张兆吉，费宇仁，陈宗宇，等.华北平原地下水可持续利用调查评价［M］.北京：地质出版社，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数据的地表温度、热惯量反演研究及其在浅层岩土含水量、地气间热交换方面的应用［D］.北京：中国科学院研究生院，2006.

［18］邓 辉.基于MODIS数据的大区域土壤水分遥感监测研究 ［D］.北京：中国农业科学院研究生院，2004.