长时间序列1984年—2020年密云水库水面信息遥感监测与分析

长时间序列1984年—2020年密云水库水面信息遥感监测与分析 Spatial-temporal water area monitoring of the Miyun Reservoir using remote sensing imagery from 1984 to 2020 1 2 3 first-author 刘畅，研究方向为水体遥感制图和应用。E-mail： lc_iecas@sina.com 刘畅，研究方向为水体遥感制图和应用。E-mail： lc_iecas@sina.com lc_iecas@sina.com 1 2 3 1 3 corresp 计璐艳，研究方向为陆表水体遥感制图。E-mail： jily@mail.ustc.edu.cn 计璐艳，研究方向为陆表水体遥感制图。E-mail： jily@mail.ustc.edu.cn jily@mail.ustc.edu.cn 1 2 3 中国科学院空天信息创新研究院, 北京 100094 中国科学院空天信息创新研究院, 北京 100094 Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China 中国科学院大学 电子电气与通信工程学院, 北京 100049 中国科学院大学 电子电气与通信工程学院, 北京 100049 School of Electronic, Electrical and Communication Engineering, University of Chinese Academy of Sciences, Beijing 100049, China School of Electronic, Electrical and Communication Engineering, University of Chinese Academy of Sciences, Beijing 100049, China 中国科学院空间信息处理与应用系统技术重点实验室, 北京 100190 中国科学院空间信息处理与应用系统技术重点实验室, 北京 100190 Key Laboratory of Technology in Geo-Spatial Information Processing and Application System, Beijing 100190, China Key Laboratory of Technology in Geo-Spatial Information Processing and Application System, Beijing 100190, China 密云水库在防洪、灌溉、发电、养殖、旅游、供应城市用水等方面产生了巨大效益，准确地完成水库制图并分析其特征时空变化对密云水库生态环境监测及管理具有重要意义。本文以密云水库为研究区域，收集了1984年—2020年所有的Landsat 5和Landsat 8遥感影像，利用水体的光谱、地形、时空特征，研究并解决了水库水体动态制图中的云、阴影、冰雪干扰，同物异谱及混合像元等难点问题，提出了一套自动化且高精度的水库水体动态制图算法，在无云无冰雪的情况下的总体精度高达98.2%，并引入了景观分裂度指数分析密云水库形态变化。基于制图结果发现水库水面面积变化可划分为5个时期：“增长期”（1984年—1993年）、“巅峰期”（1994年—1999年）、“衰落期”（2000年—2003年）、“保护期”（2004年—2014年）、“恢复期”（2015年—2020年），期间水库东、西两区经历了“三分三合”的过程，水库面积最大达到151.6 km²（1995年），最小仅为57.3 km²（2004年）；每年的最大面积出现在8—9月份，最小面积出现在5月份。此外，还讨论了密云水库水面面积变化与城市建设、政策、降水等多种驱动力的联系，可为水库的监测和管理提供参考。 Miyun Reservoir has produced huge benefits in flood control, agricultural irrigation, power generation, aquaculture, tourism, and urban water supply. Accurate water mapping is of great significance to the ecological environment monitoring of the Miyun Reservoir and the management of the South-to-North Water Diversion Project. The purpose of this research is to design a long-term dynamic water mapping method for the Miyun Reservoir by solving difficult problems such as cloud and snow interference, terrain shadows, and mixed pixels encountered in the mapping to realize the analysis and monitoring of water surface information changes in the Miyun Reservoir. The research successfully applied the tasselled cap transformation to the cloud detection of Landsat series data. DEM data are used to remove terrain shadows during processing. The water index WI, which is very suitable for long-term dynamic mapping, is introduced for water extraction, and the local unmixing method is used to make the extracted water contours and subpixel small targets more accurate. In addition, the research innovatively uses the landscape separation index to conduct a macroscopic analysis of the water surface morphology of the Miyun Reservoir. The algorithm in this paper has completed the dynamic water map of the Miyun Reservoir from 1984 to 2020, and the accuracy of the mapping result is high. The overall accuracy of direct verification is 98.2% under the condition of no clouds and snow, which is comparable to the existing water system map product (improved FROM-GLC). The overall consistency of cross-validation is as high as 99.4%. In addition, a long-term analysis of water surface information, such as the area, coverage, and morphological characteristics of the Miyun Reservoir, has been performed. (1) The area of Miyun Reservoir changed tremendously during the 37 years from 1984 to 2020, with the largest being 151.6 km² and the smallest being 57.3 km². The area changes are mainly concentrated in three areas, including the northern area, the area where the Chao and Bai Rivers enter the reservoir, and the island in the center of the reservoir. (2) Based on the changes in the area of Miyun Reservoir, it is divided into five periods: “growth” (1984—1993), “peak” (1994—1999), “decline” (2000—2003), “protection” (2004—2014), and “recovery” (2015—2020). (3) The changes in Miyun Reservoir during one year mainly occurred in 4 areas, including the northern area, the area where the Chao River enters the reservoir, the island in the center of the reservoir, and the West Stone Camel subdam area. Due to the release of water from the reservoir in May before the arrival of the rainy season each year, the area is the smallest. The area reached the largest in a year at the end of the rainy season in August or September because a large amount of rainwater was accumulated. (4) The Miyun Reservoir’s landscape separation changed greatly during the 37 years. When the water volume is small, it splits into two reservoir areas in the east and west. The east and west reservoirs underwent a process of three divisions and three consolidations from 1984 to 2020. The three split periods include 1984—1986, 2003—2005, and 2014—2015, and the other years are in the consolidation period. 光学遥感 动态水体制图 时间序列分析 密云水库 水面信息 景观指数 Landsat optical remote sensing dynamic water mapping time series analysis Miyun Reservoir water area landscape index Landsat 1　引言 1 北京是中国的首都，是政治、文化中心，但可利用的淡水资源却极其匮乏，水资源供给压力巨大。密云水库是北京最大的饮用水源供应地（ 李苗苗 等，2004 李苗苗 等，2004 ），上千万人的日常饮用水皆靠密云水库供给保障。同时，密云水库还承担着许多厂企的工业供水职责，包括燕化总厂、石电、首钢等。密云水库建成后，在防洪、灌溉、发电、养殖、旅游、供应城市用水等多方向产生了巨大效益（ 胡艳霞 等，2007 胡艳霞 等，2007 ），准确地完成密云水库制图并分析其特征时空变化，对监测及管理水库的生态环境具有重要意义。 中国水库生态环境保护工作起步晚、不成熟，缺乏科学的监测和治理方法，导致水库工程不能充分地发挥作用（ 张克锭，2019 张克锭，2019 ）。水库管理工作与信息技术的结合度较低，测绘、采样等监测管理工作依然主要由人工完成，这抑制了中国水库生态环境保护管理的实效（ 李秀歌，2019 李秀歌，2019 ）。得益于光学遥感技术的发展和高质量数据的积累，利用遥感数据进行长时间动态分析已成为水库监控的重要手段。 目前针对密云水库的遥感研究主要以水库水质、上游流域景观生态、库滨带植被分析等为主（ 胡晓静 等，2019 胡晓静 等，2019 ； 张微微 等，2020 张微微 等，2020 ； 马丰魁 等，2020 马丰魁 等，2020 ），对于水库水面信息监测的研究几乎是空白，仅有 曹荣龙等（2008） 曹荣龙等（2008） 采用水体指数法，利用7景Landsat影像实现了1984年—2005年密云水库每隔4年制图，但其制图结果陈旧且制图频次低，提取方法简单，精度有限。而随着2014年底南水北调中线工程通水，密云水库水面信息又呈现出新的特征，需要对水库进行更加科学的监测与管理（ 李秀歌，2019 李秀歌，2019 ）。 利用遥感数据进行水体提取，主要是利用了水体的光谱特征，即水体与植被、土壤等其他地物在可见光—短波红外波段的反射率差异。关于水库水体信息提取的现有方法，大致分为3类：单波段阈值法、多波段法（包括谱间关系法和水体指数法）、分类后提取法（ 刘桂林 等，2013 刘桂林 等，2013 ）。单波段阈值法利用遥感影像的短波红外波段阈值分割提取出水库水面，难以消除很多地物的影响，精度有限。分类后提取的方法包括监督分类、非监督分类等，大多着重于针对单景影像水体的提取，十分依赖人工解译。目前应用最广泛的是多波段法，其中谱间关系法通过构建不同波段的优势组合增强水体信息，从遥感影像中分离出水体，但谱间关系法的模型和阈值选取受数据源、地物谱间关系影响较大，且容易将云与水体混淆；水体指数法通过对波段进行算术运算，增强水体与其他地物的反差度，结合阈值分割进行水体制图，对于提取长时序水体信息，水体指数法是较为快捷可靠的手段（ 刘海 等，2016 刘海 等，2016 ）。在水体指数方面，目前常用的包括归一化差异水体指数（ McFeeters，1996 McFeeters，1996 ）、修正的归一化差异水体指数（ Xu，2006 Xu，2006 ）、归一化植被和水体综合指数（ Lu等，2011 Lu等，2011 ）、自动水体提取指数（ Feyisa等，2014 Feyisa等，2014 ）等。基于上述方法的现有研究主要针对单时相影像，利用长时间序列的遥感影像进行水库水面信息监测的研究相对较少。 已有的覆盖密云水库区域的水体制图产品主要有：（1）20世纪60年代—2008年中国湖泊产品（ 马荣华 等，2011 马荣华 等，2011 ），该产品每年只有1—2景，因此无法有效捕捉密云水库的动态特性，且产品依靠人工解译的方法，无法快速更新；（2）1984年—2015年30 m分辨率全球水体逐月制图产品（ Pekel等，2016 Pekel等，2016 ），该产品未处理云、冰雪以及混合像元；（3）2000年全球30 m分辨率内陆水体数据集（ Feng等，2016 Feng等，2016 ），该数据集只有2000年一期；（4）全球地表覆盖制图FROM-GLC（Finer resolution observation and monitoring of the global land cover）改进的水体制图层（Ji等，2015），原始FROM-GLC（ Gong等，2013 Gong等，2013 ）采用Landsat影像，利用监督分类方法（支持向量机）将地表分类为10个类别（其中第6类是水体），Ji等（2015）对其中的水体制图进行了改进，减少了云和山体阴影的影响，该制图产品只有2010年一期；（5）2000年—2016年全球水体逐日制图产品（ Ji等，2018 Ji等，2018 ），该产品基于500 m分辨率MODIS数据生产，相较于研究区域密云水库，空间分辨率较低，精度难以满足需求，且时间跨度较短。此外，上述水体产品均为大范围制图，在密云水库区域的提取精度难以充分保证。 因此，本文基于1984年—2020年Landsat反射率数据，研究解决了水库动态制图中的云和冰雪干扰、地形阴影及混合像元等难点问题，提出了一套自动化、高精度水库制图方法，实现了对密云水库的长时间、高频次、精确制图，并基于密云水库历史制图结果，分析了密云水库1984年—2020年的面积、覆盖率、景观形态等水面信息变化，讨论了密云水库面积变化的驱动力因素，能够为水库监测和管理提供一定的参考。 2　研究区域和数据 1 2.1　研究区域 2 密云水库（40°23′N—40°34′N，116°49′E—117°4′E）位于北京市密云县城北13 km处的燕山丘陵之中（ 图1 图1 ），在北京城区的东北方向，是北京最大的饮用水源供应地，是亚洲最大的人工湖，有“燕山明珠”之称（ 燕声，1981 燕声，1981 ）。水库建成于1960年9月，建成时面积约180 km 2 ，库容43.75亿m 3 ，平均水深30 m，分为东、西两个库区。密云水库乘纳潮白河两大支流水量，是重要的调蓄水库，具有强大蓄洪能力。白河起源于河北省沽源县，经赤城县、延庆县、怀柔区入库，潮河起源于河北省丰宁县，经滦平县、自古北口入库，流域内降水丰富且具有明显的季节性变化，7—8月降水占全年的70%，多为大雨和暴雨，冬季则降水稀少（ 李苗苗 等，2004 李苗苗 等，2004 ）。 （a） 北京市 （b） 密云水库 （a） Beijing （b） Miyun Reservoir 2.2　数据 2 2.2.1　多光谱反射率时间序列 3 Landsat系列卫星持续了长达40年的对地观测，对全球的水资源、食品、森林和其他自然资源的监测和管理起到了重要作用（ 姜高珍 等，2013 姜高珍 等，2013 ）。目前在服役的卫星为2013年发射的Landsat 8，而Landsat 5已于2013年6月退役，Landsat 7发射于1999年，在2003年发生设备故障。本次研究收集了1984年—2020年密云水库地区所有的Landsat 5和Landsat 8影像（PATH=123，ROW=032），包括485景Landsat 5影像和166景Landsat 8影像，共计651景。 图2 图2 为密云水库地区Landsat数据按年及按月统计的云覆盖率，可以看到1984年—2020年间，有10年数据的云覆盖度中值超过50%，从 图2 图2 （b）还可以发现，云的问题在夏季尤为严重。由于云覆盖度过高的影像上有效地面信息较少，本文筛除了云覆盖度超过80%的Landsat影像。此外还去除了有几何定标错误及条带状坏损的影像，最后实际使用数据474景。 （a） 年分布 （b） 月分布 （a） Annual （b） Monthly Landsat系列影像采用Worldwide Reference System-2（WRS-2）path/row参考系统，本文使用的全部波段空间分辨率为30 m，重访周期为16 d，波段范围为0.45—2.35 µm，Landsat 8和Landsat 5传感器对应波段设置如 表1 表1 所示，两者略有不同。 3 Landsat 8 3 Landsat 5 波段号 波段名 波长/µm 波段号 波段名 波长/µm 1 深蓝 0.43—0.45 2 蓝色 0.45—0.51 1 蓝色 0.45—0.52 3 绿色 0.53—0.59 2 绿色 0.52—0.60 4 红色 0.64—0.67 3 红色 0.63—0.69 5 近红外 0.85—0.88 4 近红外 0.76—0.90 6 短波红外1 1.57—1.65 5 短波红外1 1.55—1.75 7 短波红外2 2.11—2.29 7 短波红外2 2.08—2.35 2.2.2　数字高程模型 3 从 图1 图1 （b）可以看到，密云水库被群山包围，因此Landsat图像上存在大量的山体阴影。 图3 图3 给出了密云水库地区不同地物的光谱反射率曲线（基于秋、冬季节Landsat影像每类地物采集200个样本点绘制，曲线为平均值，误差棒为最小值—最大值），可以看到山体阴影和水体反射率均较低，且光谱形状有一定的相似性，因此山体阴影容易被误分类为水体。诸多研究表明，利用地形特征信息能够有效地去除山体阴影的影响（ 庞科臣等，2016 庞科臣等，2016 ），因而本文使用了“坡度阈值法”，有效地解决了这一问题。SRTM DEM（Shuttle Radar Topography Mission Digital Elevation Model）覆盖范围广（60°N—60°S），信噪比高，在遥感探测领域应用广泛（ 汪凌，2000 汪凌，2000 ），本研究利用该DEM数据计算密云水库区域坡度（ 图4 图4 ），用于密云水库附近阴影的消除。 2.2.3　验证数据 3 本文采取两种方式验证算法精度：（1）直接验证，即选择人工解译的样本作为测试样本；（2）交叉验证，即以公开发表的水体产品作为测试样本。针对每种验证方法，我们分别在密云水库及周边地区随机采样1000个样本点。具体地： （1）直接验证数据。为综合检验本文提出的算法在无云无冰雪、有云及有冰雪覆盖情况下的水体分类精度，分别选取了3景对应的影像：1）无云无冰雪影像（2013-10-31）；2）有云影像（2011-06-24）；3）有冰雪影像（2013-12-22）。然后，基于高分辨率Google Earth影像、Landsat 8多光谱数据、MODIS全球水体逐日制图产品（ Ji等，2018 Ji等，2018 ）等产品判断1000个像素点的类别，制作了3份测试样本（ 表2 表2 ），需要说明的是，这3份样本点的地理坐标都是相同的，随机样本点的位置分布如 图5 图5 所示。 图像类型 水 非水 无云无冰雪 394 606 有云 386 614 有冰雪 401 599 （2）交叉验证数据。本研究还选取了2010年度改进的FROM-GLC水体制图（Ji等，2015）对本算法结果进行交叉验证。类似直接验证，随机采样了1000个像素，制作了交叉验证样本，其中水体327个，非水体673个。 3　研究方法 1 本研究旨在研究密云水库长时间序列动态制图方法，解决制图中遇到的云和冰雪干扰、地形阴影、混合像元等难点问题，以实现对密云水库水面信息变化的分析和监控。本文成功地将缨帽变换应用到Landsat系列数据的云识别上，利用局部解混方法使提取的水体轮廓和亚像素级小目标更准确，并创新性地利用景观分离度指数对密云水库水面形态进行宏观分析。最终提出了密云水库长时间序列动态制图的整体处理流程（ 图6 图6 ），主要包括以下6个步骤： （1）针对所有Landsat影像进行云识别； （2）利用水体指数提取水库水体； （3）利用地形坡度数据，去除山体阴影； （4）基于水体提取结果，利用可见光波段阈值分割识别冰雪； （5）利用局部解混方法对水陆边界混合像元解混； （6）针对云及冰雪像素，通过时间序列插值处理，还原被覆盖的地物； 最后基于水库制图结果分析其水面信息变化，并引入了景观指数分析景观形态。 3.1　云识别 2 Landsat系列数据经常受到云层影响，致使地表反射率数据丢失，云的问题在夏季的影响最为严重，如 图2 图2 （b）所示。实现对地表的连续制图，需要准确地识别数据中的云并进行时间序列处理。 已有研究表明，缨帽变换TC（Tasseled Cap Transformation）第四分量可被用于Landsat影像的云层识别（ Joshi等，2019 Joshi等，2019 ）。缨帽变换也称K-T变换（ Kauth和Thomas，1976 Kauth和Thomas，1976 ），是将原6个光谱波段（蓝色，绿色，红色，近红外，短波红外1，短波红外2）线性变换为特定的6个复合波段的经验方法，并定义了6个复合波段的物理意义，其中前四个波段分别定义为亮度、绿度、湿度、噪声。缨帽变换系数只取决于传感器的波段设置及光谱响应特征，对同传感器的数据是固定的，不受图像的空间、时相影响，计算方法主要有施密特正交化和基于统计思想的主成分分析两种（ 王帅 等，2018 王帅 等，2018 ）。 Crist和Cicone（1984） Crist和Cicone（1984） 利用主成分变换和坐标轴旋转的方法提出了Landsat TM的缨帽变换系数，第四分量表达式为 <math display="block" id="M1"><mi>T</mi><mi>C</mi><msub><mrow><mn mathvariant="normal">4</mn></mrow><mrow><mi mathvariant="normal">T</mi><mi mathvariant="normal">M</mi></mrow></msub><mo>=</mo><mo>-</mo><mn mathvariant="normal">0.8242</mn><mo>×</mo><mi>B</mi><mo>+</mo><mn mathvariant="normal">0.0849</mn><mo>×</mo><mi>G</mi><mo>+</mo><mn mathvariant="normal">0.4392</mn><mo>×</mo><mi>R</mi><mo>-</mo><mn mathvariant="normal">0.058</mn><mo>×</mo><mi>N</mi><mi>I</mi><mi>R</mi><mo>+</mo><mn mathvariant="normal">0.2012</mn><mo>×</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">1</mn><mo>-</mo><mn mathvariant="normal">0.2768</mn><mo>×</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">2</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558959&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558950&type= 157.31066895 3.80999994 Baig等（2014） Baig等（2014） 给出了Landsat 8 OLI数据的缨帽变换系数，第四分量表达式为 <math display="block" id="M2"><mi>T</mi><mi>C</mi><msub><mrow><mn mathvariant="normal">4</mn></mrow><mrow><mi mathvariant="normal">T</mi><mi mathvariant="normal">M</mi></mrow></msub><mo>=</mo><mo>-</mo><mn mathvariant="normal">0.8239</mn><mo>×</mo><mi>B</mi><mo>+</mo><mn mathvariant="normal">0.0849</mn><mo>×</mo><mi>G</mi><mo>+</mo><mn mathvariant="normal">0.4396</mn><mo>×</mo><mi>R</mi><mo>-</mo><mn mathvariant="normal">0.058</mn><mo>×</mo><mi>N</mi><mi>I</mi><mi>R</mi><mo>+</mo><mn mathvariant="normal">0.2013</mn><mo>×</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">1</mn><mo>-</mo><mn mathvariant="normal">0.2773</mn><mo>×</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">2</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558955&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558953&type= 157.31066895 3.80999994 式中， <math id="M3"><mi>B</mi><mo>,</mo><mi>G</mi><mo>,</mo><mi>R</mi><mo>,</mo><mi>N</mi><mi>I</mi><mi>R</mi><mo>,</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">1</mn><mo>,</mo><mi>S</mi><mi>W</mi><mi>I</mi><mi>R</mi><mn mathvariant="normal">2</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558974&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558956&type= 51.39266586 3.13266683 分别表示蓝色，绿色，红色，近红外，短波红外1，短波红外2波段的反射率数值，对于Landsat 5来说，对应第1—5和7波段，对于Landsat 8来说，对应第2—7波段。由式（ 1 1 ）、（ 2 2 ）可以发现，Landsat TM和OLI传感器的系数极为接近，这主要是由于两种传感器的中心波长和光谱响应特征十分接近（ 表1 表1 ）。 <math display="block" id="M4"><mi mathvariant="normal">云</mi><mo>:</mo><mtext> </mtext><mi>T</mi><mi>C</mi><mn mathvariant="normal">4</mn><mo>≤</mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">C</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558978&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558976&type= 19.55800056 4.23333359 式中， <math id="M5"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">C</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558982&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558980&type= 3.47133350 3.72533321 表示云识别阈值，本文利用Otsu阈值分割方法（Otsu，2007）求得最佳阈值 <math id="M6"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">C</mi></mrow></msub><mo>=</mo><mo>-</mo><mn mathvariant="normal">0.046</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558985&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558983&type= 18.37266731 3.80999994 。 3.2　水体指数 2 由 图3 图3 给出的密云水库地区不同地物光谱反射率曲线可以看到：（1）水体在可见光波段（400—700 nm）的反射率比较高，但是由于水体在短波红外波段的高吸收率，导致水体这一范围的反射率几乎为零，因此水体光谱在可见光—短波红外波段表现出“前高后低”的特征；（2）水体在各波段的反射率数值存在一定的波动范围，即存在“同物异谱”现象，这主要是由于密云水库不同区域、不同时期水深及水质的差异造成的。针对密云水库水体的这两个特性，本文利用Ji等（2015）提出的水体指数WI（Water Index）进行水体提取，直接比较短波红外波段与可见光波段的最高反射率，表达式如下， <math display="block" id="M7"><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable><mtr><mtd><mn mathvariant="normal">0</mn><mo>,</mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">V</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">S</mi><mo>≤</mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mtd></mtr><mtr><mtd><mn mathvariant="normal">1</mn><mo>,</mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">V</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">S</mi><mo>></mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558988&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558987&type= 44.53466415 9.31333351 式中， <math id="M8"><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558963&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558961&type= 6.01133299 2.70933342 VIS和 <math id="M9"><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558965&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558990&type= 14.81666565 2.70933342 分别表示可见光波段和短波红外波段的最大反射率，对于Landsat 5，maxVIS=max｛band1, band2, band3｝, maxSWIR=max｛band 5, band 7｝, Landsat 8则为maxVIS=max｛band2, band3, and4｝, maxSWIR=max｛band6, band7｝。从上述表达式可以看出，WI指数仅利用水体光谱“前高后低”的形状特征，而避免使用光谱数值本身，因此可以有效避免水体“同物异谱”带来的影响。另外，WI为二值指数，因此其结果可直接作为水体分类结果。 对于无云、无冰雪或其他干扰的图像，WI的提取结果较好，鲁棒性较高，但与其他水体指数一样，WI容易将冰雪、山体阴影等误分为水体，因此还需要进行后续一系列的处理。 3.3　去除山体阴影 2 山体阴影是非常容易被误分类为水体的目标，通过 图3 图3 光谱曲线可以看出，山体阴影和水体的光谱曲线极为相似，且两者的整体反射率均较低。实验发现，有接近80%的地形阴影被WI误分类为水体，因此水库制图需要进一步消除山体阴影的影响。 中国湿地（包括水体）坡度小于3°和8°的区域分别占湿地总面积的93.85%和99.17%（ 牛振国等，2009 牛振国等，2009 ）， 图8 图8 给出了根据2017年度FROM-GLC密云水库制图结果（ http：//data.ess.tsinghua.edu.cn/fromglc10_2017v01.html http://data.ess.tsinghua.edu.cn/fromglc10_2017v01.html ［2020-11-06］）统计得到的坡度分布曲线，可以看到密云水库的水体几乎全部分布在坡度4°以下的平坦区域，因此本文利用坡度阈值法消除水体结果中的山体阴影，即 <math display="block" id="M10"><mi mathvariant="normal">阴影</mi><mo>:</mo><mi>s</mi><mi>l</mi><mi>o</mi><mi>p</mi><mi>e</mi><mo>></mo><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558971&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558993&type= 22.86000061 4.91066647 式中， <math id="M11"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559005&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559003&type= 3.04800010 3.72533321 表示坡度阈值， <math id="M12"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">s</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">4</mn><mo>°</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558994&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559007&type= 10.92199993 4.14866638 ，即把坡度大于 <math id="M13"><mn mathvariant="normal">4</mn><mo>°</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559009&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558996&type= 3.38666677 3.04800010 的像素从水体的结果中剔除。 3.4　冰雪识别 2 图3 图3 给出的密云水库地区不同地物光谱反射率曲线可以看到，冰雪光谱同样存在“前高后低”的光谱特征，因此WI指数也会将其误分为水体，但是同样可以看到，冰雪在可见光波段的反射率比水体要高很多，因此本研究利用可见光波段反射率阈值分割的方法，将水体结果中被误分的冰雪像素提取出来。具体地，如果某个像素满足以下条件： <math display="block" id="M14"><mi mathvariant="normal">冰</mi><mo>/</mo><mi mathvariant="normal">雪</mi><mo>:</mo><mi mathvariant="normal">m</mi><mi mathvariant="normal">a</mi><mi mathvariant="normal">x</mi><mi mathvariant="normal">V</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">S</mi><mo>≥</mo><msub><mrow><mi>T</mi></mrow><mrow><mrow><mrow><mi mathvariant="normal">i</mi></mrow><mo>/</mo><mrow><mi mathvariant="normal">s</mi></mrow></mrow></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559019&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559002&type= 29.71800041 4.06400013 那么将其标记为冰雪，其中阈值设置为 <math id="M15"><msub><mrow><mi>T</mi></mrow><mrow><mrow><mrow><mi mathvariant="normal">i</mi></mrow><mo>/</mo><mrow><mi mathvariant="normal">s</mi></mrow></mrow></mrow></msub><mo>=</mo><mn mathvariant="normal">0.15</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559022&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559021&type= 14.98600006 3.55599999 。 3.5　局部解混 2 受Landsat传感器分辨率影响，水陆边界处存在大量的混合像元，这些混合像元影响水库面积的准确计算，干扰了水陆分界线的精确提取，还影响对水库边界景观形态的研究。通常情况下，混合像元问题采用线性光谱解混方法来解决，针对水陆混合像元，可以用水陆双端元模型线性表出，端元系数满足非负及和为1的约束条件， <math display="block" id="M16"><mi mathvariant="bold-italic">r</mi><mo>=</mo><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><msub><mrow><mi mathvariant="bold-italic">e</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>+</mo><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub><msub><mrow><mi mathvariant="bold-italic">e</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559014&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559024&type= 22.26733208 3.72533321 <math display="block" id="M17"><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>+</mo><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559033&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559016&type= 16.42533302 3.72533321 <math display="block" id="M18"><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>,</mo><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub><mo>≥</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559036&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559035&type= 13.03866673 4.31799984 式中， <math id="M19"><mi mathvariant="bold-italic">r</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559038&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559037&type= 1.43933344 2.62466669 表示混合像元反射率向量， <math id="M20"><msub><mrow><mi mathvariant="bold-italic">e</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>,</mo><msub><mrow><mi mathvariant="bold-italic">e</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559029&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559040&type= 9.22866726 4.23333359 分别表示水体和陆地端元反射率向量， <math id="M21"><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>,</mo><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">L</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559042&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559032&type= 8.97466660 4.23333359 表示单个混合像元中水体和陆地端元的丰度系数。如果获知每个混合像元所对应的水体和陆地端元，就可以通过光谱解混的方式计算混合像元中的水体和陆地丰度。因此如何为每个混合像元确定水体和陆地端元成为关键。光谱解混分析通常利用多种地物的纯像素端元对混合像元进行线性光谱分解（ Nascimento和Dias，2005 Nascimento和Dias，2005 ； Geng等，2013 Geng等，2013 ；Ji等，2015），但实际情况卫星数据地物类型众多，且受传感器分辨率影响，纯端元难以获取。以密云水库为例，水陆边界两边的纯水端元可能是浅滩水像元或深水像元，陆地端元存在裸地、农田、草地、森林、不透水层等多种类型，而利用Landsat 30 m分辨率数据获取这些端元光谱是非常困难的。为了解决端元复杂和难以获取的问题，本文使用局部解混方法（Ji等，2015）。 利用局部解混法求取水陆边界水体丰度的具体步骤如下：（1）确定水陆混合区域。在这里，将与水陆边界像素相连的8邻域像素定义为混合区域（ 图9 图9 中的绿色像素）；（2）搜索水体、陆地端元。在以混合像元为中心的5×5窗口中搜索水陆端元。为确保水和陆地端元尽可能是纯像素，选择窗口中反射率最低和最高的像素分别作为水和陆地端元；（3）通过计算全约束条件下的最小二乘分解可获得混合像元的水体丰度（ Geng等，2016 Geng等，2016 ）；（4）阈值分割。即满足 <math id="M22"><msub><mrow><mi>c</mi></mrow><mrow><mi mathvariant="normal">w</mi></mrow></msub><mo>≥</mo><mn mathvariant="normal">0.5</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559046&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559044&type= 11.93799973 3.80999994 的混合像元分类为水体，否则为非水体。 3.6　时间序列处理 2 以往的制图研究中，大多对云和冰雪处理较粗糙，例如 Pekel等（2016） Pekel等（2016） 发布的全球逐月水体产品JRC-water将云和冰雪当作无遥感数据处理，导致其在多云季节和冬季的制图结果大量缺失。本文为排除云和冰雪造成的干扰，采用时间序列插值的方法进行处理。首先，对于制图初步结果中的每个空间点，建立时间序列向量， X =［ x 1 ， x 2 ， …， x N ］， <math id="M23"><mi>N</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559052&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559067&type= 2.62466669 2.62466669 为图像时相数目， <math id="M24"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559097&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559096&type= 2.45533323 3.72533321 代表像素点在 <math id="M25"><mi>i</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559072&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559070&type= 1.10066664 2.62466669 时相的分类结果，对于 <math id="M26"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559097&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559096&type= 2.45533323 3.72533321 存在以下4种情况： <math display="block" id="M27"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mfenced open="{" close="" separators="|"><mrow><mtable columnalign="left"><mtr><mtd><mn mathvariant="normal">1</mn><mo>,</mo><mtext> </mtext><mi mathvariant="normal">水体</mi></mtd></mtr><mtr><mtd><mn mathvariant="normal">0</mn><mo>,</mo><mtext> </mtext><mi mathvariant="normal">非水</mi><mi mathvariant="normal">体</mi></mtd></mtr><mtr><mtd><mi>C</mi><mo>,</mo><mtext> </mtext><mi mathvariant="normal">云</mi><mtext>    </mtext></mtd></mtr><mtr><mtd><mi>I</mi><mo>,</mo><mtext> </mtext><mi mathvariant="normal">冰雪</mi></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559079&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559074&type= 23.45266724 21.92866707 本文认为前两者为有效数据（1和0），而后两者（ <math id="M28"><mi>C</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559077&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559075&type= 2.28600001 2.62466669 和 <math id="M29"><mi>I</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559095&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559093&type= 1.43933344 2.62466669 ）为无效数据，对于任意时刻的 <math id="M30"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559097&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559096&type= 2.45533323 3.72533321 ，如果 <math id="M31"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mi>C</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559100&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559099&type= 9.22866726 3.80999994 或者 <math id="M32"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mi>I</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559104&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559102&type= 8.38199997 3.72533321 ，则搜索其前后时间邻域分类结果，用时间相隔最近的有效分类结果（1或0）对其进行插值，这样就完成了对无效数据的时间序列最近邻插值。 3.7　景观指数 2 景观指数可以被用来描述景观格局，进而建立景观结构与过程或现象的联系，更好地解释与理解景观功能。近年来，湿地景观格局研究已成为景观生态学和全球变化研究领域中的热点问题，通过景观指数描述景观格局，比较和分析不同尺度上的格局特点，长期以来一直倍受景观生态学家们关注（ Fu，1995 Fu，1995 ； 杨永兴，2002 杨永兴，2002 ； 宋敏敏 等，2018 宋敏敏 等，2018 ）。密云水库分为东、西两个库区（ 梁秀娟 等，2005 梁秀娟 等，2005 ； 李慧敏 等，2007 李慧敏 等，2007 ），多年来库区形态受水量影响，东、西库区不断地分裂和统一。本文利用景观分离度指数（Division）描述水库水景的破碎程度，尝试捕捉水库东、西库区之间的变化。景观分离度是描述景观中两个随机选择的像素不在同一斑块中的概率的指数（ Li等，2019 Li等，2019 ），本研究采用8连通算法计算景观分离度 <math id="M33"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559145&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559133&type= 3.80999994 3.72533321 ，公式如下： <math display="block" id="M34"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1</mn><mo>-</mo><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>k</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>M</mi></mrow></munderover></mstyle><mrow><mfrac><mrow><msubsup><mrow><mi>S</mi></mrow><mrow><mi>k</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msubsup></mrow><mrow><msup><mrow><mi>S</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></mfrac><mo>,</mo><mtext> </mtext><mi>M</mi><mo>≥</mo><mn mathvariant="normal">1</mn></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559083&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559110&type= 34.54399872 8.38199997 式中， <math id="M35"><mi>S</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559086&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559084&type= 1.86266661 2.62466669 为景观总面积， <math id="M36"><msub><mrow><mi>S</mi></mrow><mrow><mi>k</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559087&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559114&type= 2.96333337 3.72533321 为景观中互不连通的第 <math id="M37"><mi>k</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559118&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559116&type= 1.69333339 2.62466669 个斑块的面积， <math id="M38"><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559147&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559135&type= 3.04800010 2.62466669 为斑块总数， 图10 图10 为几种不同的景观形态对应的 <math id="M39"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559145&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559133&type= 3.80999994 3.72533321 示例，可以看到，当景观形态为一个整体时， <math id="M40"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559127&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559125&type= 10.15999985 3.80999994 ；当景观形态为 <math id="M41"><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559147&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559135&type= 3.04800010 2.62466669 等分时， <math id="M42"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1</mn><mo>-</mo><mn mathvariant="normal">1</mn><mo>/</mo><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559142&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559140&type= 20.48933220 3.80999994 ，因此 <math id="M43"><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559147&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559135&type= 3.04800010 2.62466669 越大， <math id="M44"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559145&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559133&type= 3.80999994 3.72533321 越大，即景观分裂成越多块，分离度指数值越大；当景观被非等分为 <math id="M45"><mi>M</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559147&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559135&type= 3.04800010 2.62466669 个斑块时，0 < P D < 1-1/ M （此时 <math id="M46"><mi>M</mi><mo>≥</mo><mn mathvariant="normal">2</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559137&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559150&type= 9.39799976 2.62466669 ）。 3.8　精度评估 2 为了定量化评估水库制图结果，本研究计算了总体精度OA（Overall Accuracy）、精确率（Precision）、召回率（Recall），计算如式（ 12 12 ）—（ 14 14 ）所示，三者都可以根据混淆矩阵（ 表3 表3 ）计算得出。 总体精度： <math display="block" id="M47"><mi>O</mi><mi>A</mi><mo>=</mo><mfrac><mrow><mi>T</mi><mi>P</mi><mo>+</mo><mi>T</mi><mi>N</mi></mrow><mrow><mi>T</mi><mi>P</mi><mo>+</mo><mi>T</mi><mi>N</mi><mo>+</mo><mi>F</mi><mi>P</mi><mo>+</mo><mi>F</mi><mi>N</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559174&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559157&type= 41.40200043 7.36600018 精确率： <math display="block" id="M48"><mi>p</mi><mi>r</mi><mi>e</mi><mi>c</mi><mi>i</mi><mi>s</mi><mi>i</mi><mi>o</mi><mi>n</mi><mo>=</mo><mfrac><mrow><mi>T</mi><mi>P</mi></mrow><mrow><mi>T</mi><mi>P</mi><mo>+</mo><mi>F</mi><mi>P</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559180&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559177&type= 31.49600029 7.28133297 召回率： <math display="block" id="M49"><mi>r</mi><mi>e</mi><mi>c</mi><mi>a</mi><mi>l</mi><mi>l</mi><mo>=</mo><mfrac><mrow><mi>T</mi><mi>P</mi></mrow><mrow><mi>T</mi><mi>P</mi><mo>+</mo><mi>F</mi><mi>N</mi></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559186&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559183&type= 26.83933067 7.36600018 真实值（1） 真实值（0） 预测值（1） True Positive （ <math id="M50"><mi>T</mi><mi>P</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559164&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559160&type= 3.30200005 2.03200006 ） False Positive （ <math id="M51"><mi>F</mi><mi>P</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559199&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559166&type= 3.47133350 2.03200006 ） 预测值（0） False Negative （ <math id="M52"><mi>F</mi><mi>N</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559190&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559201&type= 3.72533321 2.03200006 ） True Negative （ <math id="M53"><mi>T</mi><mi>N</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559204&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559202&type= 3.55599999 2.03200006 ） 4　结果 1 4.1　水库制图结果 2 基于本文算法得到的水库制图如 图11 图11 所示。对比 图11（a） 图11（a） 、（b）可以看到，水库水体被较好地提取了出来，由 图11（c） 图11（c） 、（d）可知水库的水体和湖心岛被很好地区分，水体边界细节被较好地保留，同时山体阴影也被消除。 图11（e） 图11（e） 、（f）为水库北部区域，可以看到，水库中的狭长围堰（箭头标出）被提取了出来，其实际宽度最大不超过25 m（基于Google Earth Pro高分影像目视解译测量），小于图像分辨率（30 m），实际上局部解混对于其提取具有一定的帮助。 1,2;3,4,5,6; 对于有云的图像，WI水体指数会将云错误分类，WI提取结果如 图12 图12 （b）红色矩形框所示。本文采用的缨帽变换及阈值分割能有效地识别云（ 图12 图12 （c）），再经过时间序列插值处理后，能恢复出云层下的水体（ 图12 图12 （d）），但由于云识别方法存在误差，未识别出来的云被当作地物进行分类（ 图12 图12 （d）中虚线方框区域所示，导致制图的精度下降，有云图像的提取精度详见4.2小节。 （a） False color image with cloud （b） WI result （c） Cloud extracted in this paper （d）时间序列插值后的水体制图 4.2　精度验证 2 4.2.1　直接精度验证 3 本文算法在3种情况下（图像无云无冰雪、有云、有冰雪覆盖）的精度如 表4 表4 所示，可以看到， 总体精度/% 精确率/% 召回率/% 正常影像 98.2 98.4 97.0 云影像 91.4 89.8 88.0 冰雪影像 95.8 92.8 97.0 （1）无云无冰雪。总体精度较高，为98.2%，特别是精确率达到了98.4%，表明在无云无冰雪条件下，本文算法能十分准确地提取密云水库的水体； （2）有云。在3种情况中，云对水体提取精度影响最大，精确率和召回率均低于90%，这主要是由于云识别方法存在误差，使得一些水体像素未被提取出来的（如 图12 图12 （d）），真阳性（ TP ）样本减少，导致各种精度指标降低。 （3）有冰雪。召回率较高，达到97%，但是精确率相对较低，为92.8%，有冰雪影像的分类精确率降低是因为本文方法在处理有冰雪的图像时，在冰雪与水体阈值分割过程中，存在将部分冰雪分类为水体的现象，而这部分冰雪可能实际覆盖在陆地上，导致假阳性（ FP ）样本增多，从而导致精确率降低。 4.2.2　FROM-GLC交叉验证 3 表5 表5 为利用改进的FROM-GLC水体图层进行交叉验证得到的混淆矩阵，其中总体提取结果有99.4%是相同的，结果表明本文水体提取结果与FROM-GLC的一致性非常高，这实际上说明了本文方法可以与基于人工选择训练样本及监督分类（支撑向量机）的水体提取方法精度相当，而本文方法不需要针对每景影像人工选择样本，因此本文的方法更适用于长时间序列的高精度水库动态制图。 总体精度/% 精确率/% 召回率/% FROM-GLC2010水体制图 99.4 98.2 100 5　密云水库水面信息的分析与讨论 1 5.1　密云水库年际变化 2 图13 图13 给出了1984年—2020年密云水库水体覆盖率图，可以直观地看出密云水库水域发生过较大的变化，其中变化最大的区域为北部区域，潮、白河入库口，湖心岛。在丰水期，水库向北扩张，潮白河入库口被水域覆盖，湖心岛被大部分淹没而只剩下较小的面积；而在缺水期水库面积缩减严重，湖心岛也裸露出来。 图14 图14 给出了密云水库面积（潮白河入库和出库区域不计算在内，下文同）从1984年—2020年的变化（2012年由于无卫星服役而数据缺失），图中箱体属性为上界限、上四分位、中位数、下四分位、下界限。密云水库水域37年间变化巨大，其最大面积为151.6 km 2 （1995年），最小为57.3 km 2 （2004年），2019年恢复到148.0 km 2 （2019.09.18），与北京园林绿化局（ http://yllhj.beijing.gov.cn/sdlh/jjjlylhxtfz/201909/t20190929_531708.shtml http://yllhj.beijing.gov.cn/sdlh/jjjlylhxtfz/201909/t20190929_531708.shtml ［2020-11-06］）公布的监测结果基本一致（最大161.1 km 2 ，最小56.4 km 2 ，于2019年恢复到148.7 km 2 （2019-08-15））。基于密云水库水面面积变化，可将其划分为五个时期：“增长期”（1984年—1993年）、“巅峰期”（1994年—1999年）、“衰落期”（2000年—2003年）、“保护期”（2004年—2014年）、“恢复期”（2015年—2020年）。具体地： 第一时期：“增长期”（1984年—1993年）：水库面积快速增长，这是由于国家为了治理和保护水资源被过度消耗的密云水库，颁布了一系列的政策。1985年，北京市公布“两库一渠条例”（ 朱祖希，1986 朱祖希，1986 ），禁止建设水利工程以外的任何项目，大量企业被关闭，在建项目被叫停。密云水库在政策的影响下生态环境不断改善，水域面积不断扩大。 第二时期：“巅峰期”（1994年—1999年）：水库面积维持在峰值水平。政府继续施行库区保护政策的同时，还采取了库区移民措施改善水库生态环境。到20世纪90年代初，水库周边地区的人口已严重超载，1995年首批移民迁出，前后历时6年，从库区总计移民12484人，密云水库周边地区人口超载的问题得到了大大缓解，在此期间，密云水库生态环境、水质不断改善，受利于国家的保护政策，水库面积在1995年前后达到峰值，为151.6 km 2 ，并于此时期维持在峰值水平。 第三时期：“衰落期”（2000年—2003年）：在此时期，水库面积呈现急速下降的趋势，主要由气候干旱和城区扩张用水影响。北京市此时正经历着连续多年的干旱， 图15 图15 为利用国家地球系统科学数据中心数据（ Peng等，2019 Peng等，2019 ）绘制的北京市2000年—2018年降水量曲线，由图可见在衰落期内北京市年降水量（424.48 mm）远低于2000年—2018的年平均降水量（488.42 mm），实际上据北京地区20个气象站的观测结果，1999年—2003年5年的平均降水量为1961年以来的最低值，仅为442.7 mm，比1950年以来年平均降水量少141.9 mm( 徐宗学 等，2006 徐宗学 等，2006 ； 李永坤和丁晓洁，2013 李永坤和丁晓洁，2013 )，全市平均年降水量仅有多年平均降水量的2/3，地表水可供量约占多年平均地表水可供量的1/3，平原区地下水补给量约占多年平均地下水补给量的1/2，年均供需水缺口约15亿m 3 ( 王金如，2004 王金如，2004 )，因此北京市水资源形势十分严峻，气候干旱伴随着北京市非常大的用水量缺口，密云水库陷入了危机，水库面积最低只有57.3 km 2 。 第四时期：“保护期”（2004年—2014年）：在“保护期”内，密云水库面积维持在低水平且较为稳定，主要是由于再次受到了国家保护政策的影响。2004年，王金如在第八届海峡两岸水利科技交流研讨会议上做了北京市水资源现状及展望的分析报告( 王金如，2004 王金如，2004 )，为了解决北京市以及2008年奥运会供水问题，国家再次对密云水库实施保护政策，随后密云水库面积持续下降的趋势得到控制，在长达十年的时间里，水库的面积维持在比较稳定的水平（ 贾东民 等，2012 贾东民 等，2012 ）。 第五时期“恢复期”（2015年—2020年）：水库面积逐渐恢复，近期已接近峰值水平。2014年12月12日下午14时32分，南水北调中线工程正式通水，受江水入库影响，水库面积不断攀升；此外由 图15 图15 可知，恢复期内北京市降水也有所增加（2016年降水量超过了700 mm，除2016年其他年份平均降水量（505.58 mm）略高于2000年—2018年平均水平（488.42 mm）），因此也有助于水库水量的恢复。所以密云水库在恢复期内，主要受南水北调影响，同时由于降雨量的增加、短时间暴雨以及减少出库水量的保护政策，水库面积已经恢复到接近1995年时的峰值。 综上，密云水库面积长时间发生的变化主要受政策、水利工程、气候、城市建设、工业用水等因素影响，在密云水库37年间所处的5个时期中，增长期、巅峰期、保护期均为主要受政策影响；衰落期产生的原因主要为气候连年干旱，并且城市建设消耗了过量水资源；而恢复期主要受南水北调水利工程影响，也同时受到保护政策和连续暴雨的影响。因此，政策、水利工程和气候是影响密云水库面积变化的主要驱动力因素。 5.2　密云水库年内变化 2 为了重点研究密云水库年内水体面积变化。根据 图14 图14 ，我们分别选择了1985（增长期）、1995（巅峰期，面积最大）、2002（衰落期）、2004（保护期，面积最小），2017（恢复期）五年制作了各年年内水体覆盖图，如 图16 图16 所示。可以看到，水库在年内水量充足时会向北侧高程较低的区域扩张，同时潮河入库区域的面积也会有明显的增加，其中1985年（ 图16（a） 图16（a） ）和2002年（ 图16（c） 图16（c） ）内湖心岛面积出现了较大变化，2017（ 图16（e） 图16（e） ）年水库南部区域也出现了较为明显的变化。具体地： 1,2,3,4,5; 1985年（ 图16（a） 图16（a） ）由于水库水量少，密云水库分裂为东、西两个区域，但在水量最多的9月份，水库东、西区仍呈现合二为一。实际上，东区和西区之间有一条人工建造的围堰（ 图11（e） 图11（e） 、（f）），水库水位较低时，东西区可通过围堰的缺口进行水交换（ 李慧敏 等，2007 李慧敏 等，2007 ）。 1995年（ 图16（b） 图16（b） ），水库水量充足，在潮河入库区域、北部区域有较为明显的水域变化。 2002年（ 图16（c） 图16（c） ），除北部、潮河入库口，湖心岛也有面积变化，由于此时处于衰落期，水库面积连年缩小，因此湖心岛面积不断扩大，从而引起湖心岛周边水体覆盖变化。2004年（ 图16（d） 图16（d） ），水库处在面积最小的时间点，东西区分区明显，年内除潮河入口区域有明显变化外，东西区之间还通过一条狭长的引水渠进行水交换，此时由于水库面积过小，水位过低，已不能达到围堰高度，只能通过引水渠连通。 2017年（ 图16（e） 图16（e） ），除北部、潮河入库口外，由于恢复期水库水量增加，南部的西石骆驼副坝区也出现了明显的水域扩张。 进一步，统计了密云水库每年达到最大面积和最小面积的月份，从 表6 表6 可以看到，密云水库的面积在8、9月份达到最大的次数远大于其他月份，而在5月份达到最小的次数远大于其他月份，因此可以认为密云水库面积在8—9月份达到最大，在5月份达到最小。这与 Ji等（2018） Ji等（2018） 关于全球水体面积变化分析的结论基本一致，该文献指出全球水体的最大面积出现在9月份（ Ji等，2018 Ji等，2018 ）。但该文献又表明水体面积最小出现在2月份，而本文结果是在5月份，差异较大，主要原因是 Ji等（2018） Ji等（2018） 把冬天的冰雪认为是非水，在统计水体面积时没有加入，而本文将冬天密云水库结冰区域时间序列处理后也统计在内，因此最小面积并没有出现在冬季，而出现在5月份，显然这一结果更符合密云水库。由于密云水库处于暖温带大陆性季风气候区，雨季为7—8月份，其降水占全年的70%（ 李苗苗 等，2004 李苗苗 等，2004 ），因此当8—9月份雨季结束时，密云水库由于蓄积了大量雨水，面积达到一年中的最大。在每年雨季结束之后，随着降雨减少，且上游河道常年干旱无水（ 李慧敏 等，2007 李慧敏 等，2007 ），水库下游持续错峰放水（ 张建涛 等，2010 张建涛 等，2010 ），密云水库面积持续下降，加之在来年的雨季到来之前，密云水库会为了发挥雨季防汛功能降低水位，因此水库面积也会在每年5月份达到最小。 2 参数 12 月份 1 2 3 4 5 6 7 8 9 10 11 12 最大/次数 1 0 0 0 0 0 1 13 19 1 0 1 最小/次数 1 1 0 3 28 1 0 2 0 0 0 0 5.3　密云水库景观形态变化 2 图17 图17 为密云水库景观分离度指数的年际变化趋势，可以看到，多年来景观分离度指数 <math id="M54"><msub><mrow><mi>P</mi></mrow><mrow><mi mathvariant="normal">D</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559383&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42559370&type= 3.80999994 3.72533321 基本分布在0或50%附近，但数值略高一些。这是由于在水库缺水的年份（如2004年），水域面积缩小，水库从湖心岛区域断裂开来，分成东、西两个库区（ 图16（d） 图16（d） ），根据景观分离度指数的定义可知，指数接近50%，表明水库东、西库区面积相当；而在水库水量较多年份（例如1995年），水库东西区连为一个整体，此时景观分离度较小，指数值接近0。由 图17 图17 可知，密云水库东、西库区37年间经历了“三分三合”的过程，结合 图14 图14 可以看到，水库在大部分时期内都是“合”的形态，而“分”主要出现在保护期早期（1984年、1985年）和衰落期（2004年）。相比于水库面积变化，景观分裂度指数的变化能够进一步地监控水库东、西库区之间的分合状态，并能提供东西库区面积比例等额外的形态信息，而景观形态往往又是水库水量、生态环境等因素的宏观体现，因此景观分离度指数可以作为密云水库资源信息监控和管理的重要指标之一。 景观指数能够真实地反映水库面积的宏观变化过程，景观指数的引入，给未来水资源监测和管理提供了一种新的科学角度，希望未来有更多种类的景观指数被实际应用到水资源科学监测中来。 6　结论 1 本文以Landsat数据为基础，以密云水库为研究区域，研究长时间序列动态水体制图中遇到的云和冰雪干扰、地形阴影、混合像元等难点，成功地应用缨帽变换对Landsat反射率数据进行云识别，利用局部解混方法有效地解决了大规模制图时的混合像元问题，并创新性地结合景观分离度指数对水面形态变化进行分析，最终提出了一套完整且高效的水库制图流程，包括：（1）云识别；（2）水体提取；（3）阴影去除；（4）冰雪识别；（5）混合像元解混；（6）时间序列插值。利用该算法完成了1984年—2020年的密云水库动态制图，制图结果精度较高，在无云无冰雪情况下直接验证总体精度为98.2%，与现有水体制图产品（改进的FROM-GLC）交叉验证的总体一致性高达99.4%。通过对密云水库的面积、覆盖率、形态特征等水面信息做了长时间序列分析，得到如下结论： （1）密云水库面积在1984年—2020年的37年间变化巨大，最大为151.6 km²，最小为57.3 km²，面积变化主要集中在3个区域，包括：北部区域，潮、白河入库区域，湖心岛； （2）基于密云水库面积变化，将其划分为五个时期：“增长期”（1984年—1993年）、“巅峰期”（1994年—1999年）、“衰落期”（2000年—2003年）、“保护期”（2004年—2014年）、“恢复期”（2015年—2020年）； （3）密云水库年内变化主要出现在4个区域，包括：北部区域、潮河入库区域、湖心岛、西石骆驼副坝区域。每年最大面积出现在8—9月份，最小出现在5月份，是由于在每年雨季到来之前的5月份水库泄水，导致面积最小，在8—9月份雨季结束时，蓄积了大量雨水，面积达到一年中的最大； （4）密云水库37年间景观分离度变化较大，水量较小时分裂为东、西两个库区，东、西库区从1984年—2020年经历了“三分三合”的过程，3个分裂期包括1984年—1986年、2003年—2005年、2014年—2015年，其余年份处于合并期。 后继研究中，将应用本文的水库制图方法持续地关注密云水库，并尝试推广到其他水库，更大规模地监测国水库的水面信息变化。 参考文献（References） Baig M H A , Zhang L F , Shuai T and Tong Q X . 2014 . Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance . Remote Sensing Letters , 5 ( 5 ): 423 - 431 [ DOI: 10.1080/2150704X.2014.915434 http://dx.doi.org/10.1080/2150704X.2014.915434 ] Derivation of a tasselled cap transformation based on Landsat 8 at-satellite reflectance Cao R L , Li C J , Liu L Y , Wang J H and Yan G J . 2008 . Miyun reservoir extracting miyun reservoir's water area and monitoring its change based on a revised normalized different water index . Science of Surveying and Mapping , 33 ( 2 ): 158 - 160 Miyun reservoir extracting miyun reservoir's water area and monitoring its change based on a revised normalized different water index 曹荣龙 , 李存军 , 刘良云 , 王纪华 , 阎广建 . 2008 . 基于水体指数的密云水库面积提取及变化监测 . 测绘科学 , 33 ( 2 ): 158 - 160 [ DOI: 10.3771/j.issn.1009-2307.2008.02.054 http://dx.doi.org/10.3771/j.issn.1009-2307.2008.02.054 ] 基于水体指数的密云水库面积提取及变化监测 Crist E P and Cicone R C . 1984 . A physically-based transformation of thematic mapper data---the TM tasseled cap . IEEE Transactions on Geoscience and Remote Sensing, GE-22 ( 3 ): 256 - 263 [ DOI: 10.1109/TGRS.1984.350619 http://dx.doi.org/10.1109/TGRS.1984.350619 ] A physically-based transformation of thematic mapper data---the TM tasseled cap Feng M , Sexton J O , Channan S and Townshend J R . 2016 . A global, high-resolution (30-m) inland water body dataset for 2000: first results of a topographic–spectral classification algorithm . International Journal of Digital Earth , 9 ( 2 ): 113 - 133 [ DOI: 10.1080/17538947.2015.1026420 http://dx.doi.org/10.1080/17538947.2015.1026420 ] A global, high-resolution (30-m) inland water body dataset for 2000: first results of a topographic–spectral classification algorithm Feyisa G L , Meilby H , Fensholt R and Proud S R . 2014 . Automated water extraction index: a new technique for surface water mapping using Landsat imagery . Remote Sensing of Environment , 140 : 23 - 35 [ DOI: 10.1016/j.rse.2013.08.029 http://dx.doi.org/10.1016/j.rse.2013.08.029 ] Automated water extraction index: a new technique for surface water mapping using Landsat imagery Fu B J . 1995 . Spatial pattern analysis of agricultural landscape in the loess area . Acta Ecologica Sinica , 15 ( 2 ): 113 - 120 Spatial pattern analysis of agricultural landscape in the loess area 傅伯杰 . 1995 . 黄土区农业景观空间格局分析 . 生态学报 , 15 ( 2 ): 113 - 120 [ DOI: 10.3321/j.issn:1000-0933.1995.02.001 http://dx.doi.org/10.3321/j.issn:1000-0933.1995.02.001 ] 黄土区农业景观空间格局分析 Geng X R , Ji L Y and Sun K . 2016 . Non-negative matrix factorization based unmixing for principal component transformed hyperspectral data . Frontiers of Information Technology and Electronic Engineering , 17 ( 5 ): 403 - 412 [ DOI: 10.1631/fitee.1600028 http://dx.doi.org/10.1631/fitee.1600028 ] Non-negative matrix factorization based unmixing for principal component transformed hyperspectral data Geng X R , Xiao Z Q , Ji L Y , Zhao Y C and Wang F X . 2013 . A Gaussian elimination based fast endmember extraction algorithm for hyperspectral imagery . ISPRS Journal of Photogrammetry and Remote Sensing , 79 : 211 - 218 [ DOI: 10.1016/j.isprsjprs.2013.02.020 http://dx.doi.org/10.1016/j.isprsjprs.2013.02.020 ] A Gaussian elimination based fast endmember extraction algorithm for hyperspectral imagery Gong P , Wang J , Yu L , Zhao Y C , Zhao Y Y , Liang L , Niu Z G , Huang X M , Fu H H , Liu S , Li C C , Li X Y , Fu W , Liu C X , Xu Y , Wang X Y , Cheng O , Hu L Y , Yao W B , Zhang H , Zhu P , Zhao Z Y , Zhang H Y , Zheng Y M , Ji L Y , Zhang Y W , Chen H , Yan A , Guo J H , Yu L , Wang L , Liu X J , Shi T T , Zhu M H , Chen Y L , Yang G W , Tang P , Xu B , Giri C , Clinton N , Zhu Z L , Chen J and Chen J . 2013 . Finer resolution observation and monitoring of global land cover: first mapping results with landsat TM and ETM+ data . International Journal of Remote Sensing , 34 ( 7 ): 2607 - 2654 [ DOI: 10.1080/01431161.2012.748992 http://dx.doi.org/10.1080/01431161.2012.748992 ] Finer resolution observation and monitoring of global land cover: first mapping results with landsat TM and ETM+ data Hu X J , Zhang Y F , Wu J D and Cheng J H . 2019 . Characteristics of vegetation change in the Miyun Reservoir before and after the South-to-North Water Transfer Project //Liu H L, ed. Beijing : China Water Resources and Hydropower Press: 140 - 147 Characteristics of vegetation change in the Miyun Reservoir before and after the South-to-North Water Transfer Project 胡晓静 , 张耀方 , 吴敬东 , 程金花 . 2019 . 南水北调来水前后密云水库库滨带植被变化特征分析 // 刘洪禄 . 北京水问题研究与实践( 2018 年). 北京 : 中国水利水电出版社: 140 - 147 南水北调来水前后密云水库库滨带植被变化特征分析 Hu Y X , Zhou L D , Yan M C , Dong X B , Wu Z Q and Wei C S . 2007 . The valuation and compensation policies suggestions of main capitals and services of Miyun reservoir . Journal of Natural Resources , 22 ( 4 ): 497 - 506 The valuation and compensation policies suggestions of main capitals and services of Miyun reservoir 胡艳霞 , 周连第 , 严茂超 , 董孝斌 , 吴志强 , 魏长山 . 2007 . 北京密云水库生态经济系统特征、资产基础及功能效益评估 . 自然资源学报 , 22 ( 4 ): 497 - 506 [ DOI: 10.3321/j.issn:1000-3037.2007.04.001 http://dx.doi.org/10.3321/j.issn:1000-3037.2007.04.001 ] 北京密云水库生态经济系统特征、资产基础及功能效益评估 Ji L Y , Geng X R , Sun K , Zhao Y C and Gong P . 2015a . Modified N-FINDR endmember extraction algorithm for remote-sensing imagery . International Journal of Remote Sensing , 36 ( 8 ): 2148 - 2162 [ DOI: 10.1080/01431161.2015.1034895 http://dx.doi.org/10.1080/01431161.2015.1034895 ] Modified N-FINDR endmember extraction algorithm for remote-sensing imagery Ji L Y , Gong P , Geng X R and Zhao Y C . 2015b . Improving the accuracy of the water surface cover type in the 30 m FROM-GLC product . Remote Sensing , 7 ( 10 ): 13507 - 13527 [ DOI: 10.3390/rs71013507 http://dx.doi.org/10.3390/rs71013507 ] Improving the accuracy of the water surface cover type in the 30 m FROM-GLC product Ji L Y , Gong P , Wang J , Shi J C and Zhu Z L . 2018 . Construction of the 500-m resolution daily global surface water change database (2001-2016) . Water Resources Research , 54 ( 12 ): 10270 - 10292 [ DOI: 10.1029/2018WR023060 http://dx.doi.org/10.1029/2018WR023060 ] Construction of the 500-m resolution daily global surface water change database (2001-2016) Jia D M , Gao X Y and Hao L J . 2012 . Discussion on water environment safety protection measures in the Miyun Reservoir basin . Beijing Water , ( 1 ): 7 - 10 Discussion on water environment safety protection measures in the Miyun Reservoir basin 贾东民 , 高训宇 , 郝丽娟 . 2012 . 密云水库流域水环境安全保护措施探讨 . 北京水务 , ( 1 ): 7 - 10 [ DOI: 10.3969/j.issn.1673-4637.2012.01.004 http://dx.doi.org/10.3969/j.issn.1673-4637.2012.01.004 ] 密云水库流域水环境安全保护措施探讨 Jiang G Z , Han B , Gao Y B and Yang C J . 2013 . Review of 40-year earth observation with Landsat series and prospects of LDCM . Journal of Remote Sensing , 17 ( 5 ): 1033 - 1048 Review of 40-year earth observation with Landsat series and prospects of LDCM 姜高珍 , 韩冰 , 高应波 , 杨崇俊 . 2013 . Landsat系列卫星对地观测40年回顾及LDCM前瞻（英文） . 遥感学报 , 17 ( 5 ): 1033 - 1048 Landsat系列卫星对地观测40年回顾及LDCM前瞻（英文） Joshi P P , Wynne R H and Thomas V A . 2019 . Cloud detection algorithm using SVM with SWIR2 and tasseled cap applied to Landsat 8 . International Journal of Applied Earth Observation and Geoinformation , 82 : 101898 [ DOI: 10.1016/j.jag.2019.101898 http://dx.doi.org/10.1016/j.jag.2019.101898 ] Cloud detection algorithm using SVM with SWIR2 and tasseled cap applied to Landsat 8 Kauth R J and Thomas G S . 1976 . The Tasselled Cap—a Graphic Description of the Spectral-Temporal Development of Agricultural Crops as Seen by LANDSAT The Tasselled Cap—a Graphic Description of the Spectral-Temporal Development of Agricultural Crops as Seen by LANDSAT Lehner B and Döll P . 2004 . Development and validation of a global database of lakes, reservoirs and wetlands . Journal of Hydrology , 296 ( 1/4 ): 1 - 22 [ DOI: 10.1016/j.jhydrol.2004.03.028 http://dx.doi.org/10.1016/j.jhydrol.2004.03.028 ] Development and validation of a global database of lakes, reservoirs and wetlands Li H M , Meng F Y , Du G S , Liu X D , Liu J and Song F . 2007 . Analysis on the phytoplankton and water quality in eastern and western Miyun Reservoir . Journal of Lake Sciences , ( 2 ): 146 - 150 Analysis on the phytoplankton and water quality in eastern and western Miyun Reservoir 李慧敏 , 孟凡艳 , 杜桂森 , 刘晓端 , 刘静 , 宋福 . 2007 . 密云水库东西库区的水质与浮游藻类分析 . 湖泊科学 , ( 2 ): 146 - 150 密云水库东西库区的水质与浮游藻类分析 Li M M , Wu B F , Yan C Z and Zhou W F . 2004 . Estimation of vegetation fraction in the upper Basin of Miyun reservoir by remote sensing . Resources Science , 26 ( 4 ): 153 - 159 Estimation of vegetation fraction in the upper Basin of Miyun reservoir by remote sensing 李苗苗 , 吴炳方 , 颜长珍 , 周为峰 . 2004 . 密云水库上游植被覆盖度的遥感估算 . 资源科学 , 26 ( 4 ): 153 - 159 [ DOI: 10.3321/j.issn:1007-7588.2004.04.022 http://dx.doi.org/10.3321/j.issn:1007-7588.2004.04.022 ] 密云水库上游植被覆盖度的遥感估算 Li X G . 2019 . Difficulties and countermeasures of ecological environment protection of reservoirs in my country . Agriculture of Jilin , ( 10 ): 39 Difficulties and countermeasures of ecological environment protection of reservoirs in my country 李秀歌 . 2019 . 我国水库生态环境保护面临的困境与对策 . 吉林农业 , 10 : 39 [ DOI: 10.14025/j.cnki.jlny.2019.10.016 http://dx.doi.org/10.14025/j.cnki.jlny.2019.10.016 ] 我国水库生态环境保护面临的困境与对策 Li Y K and Ding X J . 2013 . Analysis on the characteristics of precipitation changes in Beijing . Beijing Water , ( 2 ): 9 - 12 Analysis on the characteristics of precipitation changes in Beijing 李永坤 , 丁晓洁 . 2013 . 北京市降水量变化特征分析 . 北京水务 , ( 2 ): 9 - 12 [ DOI: 10.3969/j.issn.1673-4637.2013.02.003 http://dx.doi.org/10.3969/j.issn.1673-4637.2013.02.003 ] 北京市降水量变化特征分析 Li Z C , Feng Y J , Gurgel H , Xu L , Dessay N and Gong P . 2019 . Use of spatial autocorrelation and time series Landsat images for long-term monitoring of surface water shrinkage and expansion in Guanting Reservoir, China . Remote Sensing Letters , 10 ( 12 ): 1192 - 1200 [ DOI: 10.1080/2150704X.2019.1671634 http://dx.doi.org/10.1080/2150704X.2019.1671634 ] Use of spatial autocorrelation and time series Landsat images for long-term monitoring of surface water shrinkage and expansion in Guanting Reservoir, China Liang X J , Xiao C L , Yang T X , Wang J and Liu X D . 2005 . Research on the influence factors of nitrogen distribution and migration in Miyun Reservoir . Science in China Series D , 35 ( S1 ): 272 - 280 Research on the influence factors of nitrogen distribution and migration in Miyun Reservoir 梁秀娟 , 肖长来 , 杨天行 , 王静 , 刘晓端 . 2005 . 密云水库中氮分布及迁移影响因素研究 . 中国科学 D辑 地球科学 , 35 ( S1 ): 272 - 280 密云水库中氮分布及迁移影响因素研究 Liu G L , Zhang L C , Liu J and Li G Y . 2013 . Water body information extraction based on Landsat TM remote sensing imagery . Journal of University of Chinese Academy of Sciences , 30 ( 5 ): 644 - 650 Water body information extraction based on Landsat TM remote sensing imagery 刘桂林 , 张落成 , 刘剑 , 李广宇 . 2013 . 基于Landsat TM影像的水体信息提取 . 中国科学院大学学报 , 30 ( 5 ): 644 - 650 [ DOI: 10.7523/j.issn.2095-6134.2013.05.011 http://dx.doi.org/10.7523/j.issn.2095-6134.2013.05.011 ] 基于Landsat TM影像的水体信息提取 Liu H , Yin J , Chen J and Chen X L . 2016 . Dynamic change and casue of water area in Danjiangkou reservoir based on remote sensing image . Resources and Environment in the Yangtze Basin , 25 ( 11 ): 1759 - 1766 Dynamic change and casue of water area in Danjiangkou reservoir based on remote sensing image 刘海 , 殷杰 , 陈晶 , 陈晓玲 . 2016 . 基于遥感影像的丹江口水库水域面积动态变化与原因研究 . 长江流域资源与环境 , 25 ( 11 ): 1759 - 1766 [ DOI: 10.11870/cjlyzyyhj2016011015 http://dx.doi.org/10.11870/cjlyzyyhj2016011015 ] 基于遥感影像的丹江口水库水域面积动态变化与原因研究 Lu S L , Wu B F , Yan N N and Wang H . 2011 . Water body mapping method with HJ-1A/B satellite imagery . International Journal of Applied Earth Observation and Geoinformation , 13 ( 3 ): 428 - 434 [ DOI: 10.1016/j.jag.2010.09.006 http://dx.doi.org/10.1016/j.jag.2010.09.006 ] Water body mapping method with HJ-1A/B satellite imagery Ma F K , Jiang Q O , Xu L D , Liang Y , Wang R C and Su S . 2020 . Retrieval of water quality parameters based on bp neural network algorithm in Miyun reservoir . Ecology and Environmental Sciences , 29 ( 3 ): 569 - 579 Retrieval of water quality parameters based on bp neural network algorithm in Miyun reservoir 马丰魁 , 姜群鸥 , 徐藜丹 , 梁勇 , 王荣臣 , 苏帅 . 2020 . 基于BP神经网络算法的密云水库水质参数反演研究 . 生态环境学报 , 29 ( 3 ): 569 - 579 [ DOI: 10.16258/j.cnki.1674-5906.2020.03.017 http://dx.doi.org/10.16258/j.cnki.1674-5906.2020.03.017 ] 基于BP神经网络算法的密云水库水质参数反演研究 Ma R H , Yang G S , Duan H T , Jiang J H , Wang S M , Feng X Z , Li A N , Kong F X , Xue B , Wu J L and Li S J . 2011 . China’s lakes at present: Number, area and spatial distribution . Science China Earth Sciences , 54 ( 2 ): 283 - 289 China’s lakes at present: Number, area and spatial distribution 马荣华 , 杨桂山 , 段洪涛 , 姜加虎 , 王苏民 , 冯学智 , 李爱农 , 孔繁翔 , 薛滨 , 吴敬禄 , 李世杰 . 2011 . 中国湖泊的数量、面积与空间分布 . 中国科学 : 地球科学 , 41 ( 3 ): 394 - 401 [ DOI: 10.1007/s11430-010-4052-6 http://dx.doi.org/10.1007/s11430-010-4052-6 ] 中国湖泊的数量、面积与空间分布 McFeeters S K . 1996 . The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features . International Journal of Remote Sensing , 17 ( 7 ): 1425 - 1432 [ DOI: 10.1080/01431169608948714 http://dx.doi.org/10.1080/01431169608948714 ] The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features Nascimento J M P and Dias J M B . 2005 . Vertex component analysis: a fast algorithm to unmix hyperspectral data . IEEE Transactions on Geoscience and Remote Sensing , 43 ( 4 ): 898 - 910 [ DOI: 10.1109/TGRS.2005.844293 http://dx.doi.org/10.1109/TGRS.2005.844293 ] Vertex component analysis: a fast algorithm to unmix hyperspectral data Niu Z G , Gong P , Cheng X , Guo J H , Wang L , Huang H B , Shen S Q , Wu J Z , Wang X F , Wang X C , Wang Y Q , Liang L , Zhang L N , Wang L , Yao Q , Yang Z Z , Guo Z Q and Dai Y J . 2009 . Preliminary remote sensing mapping of wetlands in China and analysis of related geographic characteristics . Science in China : Series D , 39 ( 2 ): 188 - 203 Preliminary remote sensing mapping of wetlands in China and analysis of related geographic characteristics 牛振国 , 宫鹏 , 程晓 , 虢建宏 , 王琳 , 黄华兵 , 沈少青 , 吴昀昭 , 王晓风 , 王显成 , 王应清 , 梁璐 , 张丽娜 , 王雷 , 姚谦 , 杨镇钟 , 郭子祺 , 戴永久 . 2009 . 中国湿地初步遥感制图及相关地理特征分析 . 中国科学 D辑 地球科学 , 39 ( 2 ): 188 - 203 中国湿地初步遥感制图及相关地理特征分析 Otsu N . 1979 . A threshold selection method from gray-level histograms . IEEE Transactions on Systems, Man, and Cybernetics , 9 ( 1 ): 62 - 66 [ DOI: 10.1109/TSMC.1979.4310076 http://dx.doi.org/10.1109/TSMC.1979.4310076 ] A threshold selection method from gray-level histograms Pang K C , Chen L F and Wang S Y . 2016 . Water extraction method combining DEM for hillshade removal . Electronic Science and Technology , 29 ( 4 ): 76 - 78 Water extraction method combining DEM for hillshade removal 庞科臣 , 陈立福 , 王思雨 . 2016 . 融合DEM去除山体阴影的水体提取方法 . 电子科技 , 29 ( 4 ): 76 - 78 [ DOI: 10.16180/j.cnki.issn1007-7820.2016.04.020 http://dx.doi.org/10.16180/j.cnki.issn1007-7820.2016.04.020 ] 融合DEM去除山体阴影的水体提取方法 Pekel J F , Cottam A , Gorelick N and Belward A S . 2016 . High-resolution mapping of global surface water and its long-term changes . Nature , 540 ( 7633 ): 418 - 422 [ DOI: 10.1038/nature20584 http://dx.doi.org/10.1038/nature20584 ] High-resolution mapping of global surface water and its long-term changes Peng S Z , Ding Y X , Liu W Z and Li Z . 2019 . 1 km monthly temperature and precipitation dataset for China from 1901 to 2017 . Earth System Science Data , 11 ( 4 ): 1931 - 1946 [ DOI: 10.5194/essd-11-1931-2019 http://dx.doi.org/10.5194/essd-11-1931-2019 ] 1 km monthly temperature and precipitation dataset for China from 1901 to 2017 Song M M , Zhang Q F , Wu F Q , Wu B X and Wu B . 2018 . Landscape pattern changes and evaluation of ecological service values in a small watershed of the Loess Gully Region . Acta Ecologica Sinica , 38 ( 8 ): 2649 - 2659 Landscape pattern changes and evaluation of ecological service values in a small watershed of the Loess Gully Region 宋敏敏 , 张青峰 , 吴发启 , 吴秉校 , 吴驳 . 2018 . 黄土沟壑区小流域景观格局演变及生态服务价值响应 . 生态学报 , 38 ( 8 ): 2649 - 2659 [ DOI: 10.5846/stxb201705210935 http://dx.doi.org/10.5846/stxb201705210935 ] 黄土沟壑区小流域景观格局演变及生态服务价值响应 Wang J R . 2004 . Beijing water resources status and prospects // Cross-Strait Water Conservancy Science and Technology Exchange Seminar . Guangzhou : China Institute of Water Resources and Hydropower Research Beijing water resources status and prospects 王金如 . 2004 . 北京水資源現狀及展望 // 第八届海峡两岸水利科技交流研讨会 . 广州 : 中国水利水电科学研究院 北京水資源現狀及展望 Wang L . 2000 . A brief introduction to us shuttle radar topography mission . Bulletin of Surveying and Mapping , ( 12 ): 38 - 40 A brief introduction to us shuttle radar topography mission 汪凌 . 2000 . 美国航天飞机雷达地形测绘使命简介 . 测绘通报 , ( 12 ): 38 - 40 [ DOI: 10.3969/j.issn.0494-0911.2000.12.018 http://dx.doi.org/10.3969/j.issn.0494-0911.2000.12.018 ] 美国航天飞机雷达地形测绘使命简介 Wang S , Xu H Q and Shi T T . 2018 . Retrieval of tasseled cap transformation coefficients for GF-1 WFV2 sensor data . Advances in Earth Science , 33 ( 6 ): 641 - 652 Retrieval of tasseled cap transformation coefficients for GF-1 WFV2 sensor data 王帅 , 徐涵秋 , 施婷婷 . 2018 . GF-1 WFV2传感器数据的缨帽变换系数反演 . 地球科学进展 , 33 ( 6 ): 641 - 652 [ DOI: 10.11867/j.issn.1001-8166.2018.06.0641 http://dx.doi.org/10.11867/j.issn.1001-8166.2018.06.0641 ] GF-1 WFV2传感器数据的缨帽变换系数反演 Xu H Q . 2006 . Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery . International Journal of Remote Sensing , 27 ( 14 ): 3025 - 3033 [ DOI: 10.1080/01431160600589179 http://dx.doi.org/10.1080/01431160600589179 ] Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery Xu Z X , Zhang L and Ruan B Q . 2006 . Analysis on the spatiotemporal distribution of precipitation in the Beijing Region . Arid Land Geography , 29 ( 2 ): 186 - 192 Analysis on the spatiotemporal distribution of precipitation in the Beijing Region 徐宗学 , 张玲 , 阮本清 . 2006 . 北京地区降水量时空分布规律分析 . 干旱区地理 , 29 ( 2 ): 186 - 192 [ DOI: 10.3321/j.issn:1000-6060.2006.02.004 http://dx.doi.org/10.3321/j.issn:1000-6060.2006.02.004 ] 北京地区降水量时空分布规律分析 Yan S . 1981 . Protection of Miyun Reservoir, the Pearl of Yanshan Mountain . Environmental Protection , ( 5 ): 11 - 12 Protection of Miyun Reservoir, the Pearl of Yanshan Mountain 燕声 . 1981 . 保护燕山明珠密云水库 . 环境保护 , ( 5 ): 11 - 12 保护燕山明珠密云水库 Yang Y X . 2002 . New knowledge on the progress of international wetland science research and priority field and prospect of Chinese wetland science research . Advances in Earth Science , 17 ( 4 ): 508 - 514 New knowledge on the progress of international wetland science research and priority field and prospect of Chinese wetland science research 杨永兴 . 2002 . 国际湿地科学研究进展和中国湿地科学研究优先领域与展望 . 地球科学进展 , 17 ( 4 ): 508 - 514 [ DOI: 10.3321/j.issn:1001-8166.2002.04.008 http://dx.doi.org/10.3321/j.issn:1001-8166.2002.04.008 ] 国际湿地科学研究进展和中国湿地科学研究优先领域与展望 Zhang J T , Cong X H and Zhao Y F . 2010 . Study of super standard flood in the Chaobai river of Beijing . Beijing Water , ( 2 ): 23 - 25 Study of super standard flood in the Chaobai river of Beijing 张建涛 , 丛晓红 , 赵月芬 . 2010 . 北京市潮白河超标准洪水研究 . 北京水务 , ( 2 ): 23 - 25 [ DOI: 10.3969/j.issn.1673-4637.2010.02.009 http://dx.doi.org/10.3969/j.issn.1673-4637.2010.02.009 ] 北京市潮白河超标准洪水研究 Zhang K D . 2019 . Problems and countermeasures in reservoir management . Rural Economy and Science-Technology , 303 ( 14 ): 30 Problems and countermeasures in reservoir management 张克锭 . 2019 . 水库管理存在的问题及对策 . 农村经济与科技 , 303 ( 14 ): 30 水库管理存在的问题及对策 Zhang W W , Li X N , Wang C , Zhao C Q and Shi R S . 2020 . Water quality response to landscape pattern at different spatial scales in Baihe river in the upper reaches of the Miyun reservoir . Environmental Science , 41 ( 11 ): 4895 - 4904 Water quality response to landscape pattern at different spatial scales in Baihe river in the upper reaches of the Miyun reservoir 张微微 , 李晓娜 , 王超 , 赵春桥 , 史瑞双 . 2020 . 密云水库上游白河地表水质对不同空间尺度景观格局特征的响应 . 环境科学 , 41 ( 11 ): 4895 - 4904 [ DOI: 10.13227/j.hjkx.202003250 http://dx.doi.org/10.13227/j.hjkx.202003250 ] 密云水库上游白河地表水质对不同空间尺度景观格局特征的响应 Zhou Y , Xie G L , Wang S X , Wang F and Wang F T . 2014 . Information extraction of thin rivers around built-up lands with false NDWI . Journal of Geo-Information Science , 16 ( 1 ): 102 - 107 Information extraction of thin rivers around built-up lands with false NDWI 周艺 , 谢光磊 , 王世新 , 王峰 , 王福涛 . 2014 . 利用伪归一化差异水体指数提取城镇周边细小河流信息 . 地球信息科学学报 , 16 ( 1 ): 102 - 107 [ DOI: 10.3724/SP.J.1047.2014.000102 http://dx.doi.org/10.3724/SP.J.1047.2014.000102 ] 利用伪归一化差异水体指数提取城镇周边细小河流信息 Zhu Z X . 1986 . On the significance of the "interim measures for the protection and management of Water Sources in Beijing's Miyun Reservoir, Huairou Reservoir and Jingmi Diversion Canal ". Environmental Protection , ( 3 ): 11 - 14 On the significance of the "interim measures for the protection and management of Water Sources in Beijing's Miyun Reservoir, Huairou Reservoir and Jingmi Diversion Canal 朱祖希 . 1986 . 试论《北京市密云水库、怀柔水库和京密引水渠水源保护管理暂行办法》的重要意义 . 环境保护 , ( 3 ): 11 - 14 试论《北京市密云水库、怀柔水库和京密引水渠水源保护管理暂行办法》的重要意义 致谢:感谢中国科学院空天信息创新研究院的席彦新、朱亮亮、陈鑫帮助检查文章错误以及润色本文中的语言！感谢清华大学地球科学系统系提供遥感数据！ 10.11834/jrs.20220489 P2 1007-4619（2023）02-0335-16 刘畅，唐海蓉，计璐艳，赵永超.2023.长时间序列1984年—2020年密云水库水面信息遥感监测与分析.遥感学报，27（2）： 335-350 刘畅，唐海蓉，计璐艳，赵永超.2023.长时间序列1984年—2020年密云水库水面信息遥感监测与分析.遥感学报，27（2）： 335-350 Liu C，Tang H R，Ji L Y and Zhao Y C. 2023. Spatial-temporal water area monitoring of the Miyun Reservoir using remote sensing imagery from 1984 to 2020. National Remote Sensing Bulletin， 27（2）：335-350 Liu C，Tang H R，Ji L Y and Zhao Y C. 2023. Spatial-temporal water area monitoring of the Miyun Reservoir using remote sensing imagery from 1984 to 2020. National Remote Sensing Bulletin， 27（2）：335-350 第二次青藏高原综合科学考察研究（编号：2019QZKK0206）；国家自然科学基金（编号：61701477）；国家自然科学基金青年基金（编号：61805246） Supported by The Second Comprehensive Scientific Investigation of the Qinghai-Tibet Plateau (No. 2019QZKK0206); National Natural Science Foundation of China (No. 61701477); National Natural Science Foundation of China Youth Fund (No. 61805246) 2020-11-06 2023-02-07 2023-04-20 2021-12-27 遥感学报 2 27