单春芝,来庆广,刘焱春,石锦坤
(1.青岛地质工程勘察院(青岛地质勘查开发局),山东 青岛 370200;2.自然资源部 滨海城市地下空间地质安全重点实验室,山东 青岛 370200)
受黄河泥沙输送以及海洋动力作用影响,黄河三角洲海岸带的形态一直处于动态变化的状态[1]。笔者在多年持续的业务工作中发现,近年来黄河三角洲形态变化明显,原因未知。因此,分析黄河三角洲近年来形态变化特征,把握其变化规律,不论是对黄河的调水调沙工作,还是对研究近年来黄河三角洲变化加剧的原因都具有重要意义。
对于宏观区域的演化分析,遥感手段是一种有效、精准的方法。已有研究主要从确定海岸线的角度分析滨海地区形态变化,主要有一般高潮线法、低潮线法等[2~3]。薛允传等[4]采用平均高潮线法解译黄河三角洲25 a的岸线信息,分析其岸线变化特征及时空发育规律。杨江平等[1]以低潮线为基线,研究了1996—2011年黄河口河道、岸线演变特征及近岸淤蚀状况。上述方法主要是利用多期遥感影像,通过遥感影像自动分类、人机交互识别[5-9]的方法确定岸线位置,分析历史变化规律。此外,黄李冰等[10]从黄河行河河道的洪枯季变化的角度,研究黄河现行入海流路的摆动规律及其趋势。程慧等[11]采用遥感技术及数理统计法对岸线及面积变化进行监测计算,分析了黄河三角洲4个不同阶段的冲淤演变及影响因素。
滨海区域演变过程分析中的海岸线,是基于遥感影像人工判读、加工出的“一条线”,而从遥感影像获取的最直观、最精确的瞬时水边界数据并未得到直接利用,也没有体现高潮、低潮间海岸带的变化。本研究收集2003年以来黄河三角洲地区遥感影像,根据临近海域潮汐的变化规律,利用农历日潮汐的变化规律推算影像成像时间对应的潮位,将影像根据潮位分为高潮、低潮两组。通过人机交互的方式解译两组瞬时水边界数据,分析黄河三角洲整体及局部变化过程、新生沙洲以及入海口位置、方向的变化规律,用于表征黄河三角洲海岸带形态变化特征。
研究区位于黄河三角洲区域1号与2号点之间,1号点坐标为东经119°05′48.31″、北纬37°51′09.23″,2号点坐标为东经119°01′19.87″、北纬37°41′28.67″,如图1所示,影像为2020年10月24日Landsat 8低潮时遥感影像。1976—1996年黄河向东南方向入海,1996—2007年改向东流,2007年再次改道向北入海,2012年河道继续向北迁移,并且开始分叉。至2020年,沿主河道向北淤积严重,已经与西侧孤东大堤形成U形海湾,主河道入海口处形成两处沙洲。
图1 研究区范围
收集29景研究区美国陆地资源系列卫星影像(http://ids.ceode.ac.cn/获取),空间分辨率30 m,补充1景国产高分系列卫星影像(中国资源卫星应用中心http://www.cresda.com/CN/获取),空间分辨率16 m,时间跨度为2003—2020年。该影像既能满足海岸带变化分析的精度要求,又能满足覆盖范围的要求,适用性较好。
根据罗丹等[12]的研究成果,黄河三角洲区域与莱州湾潮汐均为不规则半日潮。由于无法获得黄河三角洲区域详尽的历史潮汐数据,因此利用莱州湾潍坊港的潮汐数据代表黄河口区域的潮汐数据。由于只能获得潍坊港2020年的数据,其他历史潮汐数据需要推算获得,因此利用影像获取时间对应的农历日时间,对照2020年相同农历日所对应的日期,从海事服务网(ht⁃tps://www.cnss.com.cn/tide/)查询潍坊港2020年潮位,推算影像获取时的潮位信息。根据海事服务网统计数据,潍坊港高潮为2.2 m左右、低潮为0.7 m左右,获得完全与高、低潮时间相匹配的影像较为困难,利用潮位推算结果,筛选潮位在1.5 m以上时对应的影像作为高潮影像,潮位在1 m以下时的影像作为低潮影像,见表1、表2。
表1 高潮影像信息
表2 低潮影像信息
2.2.1 数据处理
本研究所获取的卫星遥感影像为L1级数据,经过了系统辐射校正和地面控制点的几何校正。按照研究区范围裁剪影像,为便于后续假彩色合成,保留蓝光、绿光、红光以及近红外波段数据。为保证对比分析的可靠性,以2003年影像为基准影像,配准其他影像,配准误差在1个像元以内。长度、面积、角度计算在高斯-克吕格投影(3°分带,120°中央经线)坐标系进行。
本研究通过人机交互方式目视解译各期遥感影像瞬时水边界数据,为突出水陆分界,采用绿光、红光、近红外波段标准假彩色合成影像[13],解译过程可通过调节影像直方图、对比度、亮度等方式增强图像,提高解译效果。解译过程原则统一,在滩涂和水域明显区分的地方勾画水边界,潮滩上细小的潮沟和黄河入海口处,拉直勾画;对于处于零散状态、边界不够明显的区域,以呈片状分布的滩涂区域外边界勾画;沙洲单独勾画。最后形成高潮、低潮水边界矢量数据集,获得9条高潮瞬时水边界和10条低潮瞬时水边界,见图2。
图2 瞬时水边界解译结果及断面布设
2.2.2 断面计算
在海岸形态变化分析方面,采用较为成熟的美国地质调查局(USGS)研发的数字海岸线分析系统(DSAS,Digital Shoreline Analysis System)进行岸线形态的计算、分析。该系统主要是在研究区均匀设置一定密度的断面,通过计算断面点的变化来反映区域内海岸形态的变化。比如刘鹏等[14]通过获取黄河三角洲多年一般大潮高潮线,采用数字岸线分析系统和岸线分形分析法,定量分析了黄河三角洲岸线长度、形态及变化过程。
根据瞬时水边界解译结果,黄河入海口处岸线形态变化角度、距离较大,无法直接利用DSAS系统计算断面变化。分析可知黄河三角洲主要有8个典型变化区域,本研究参照DSAS方法,根据这些变化区域分别设计8条监测断面(高潮GC01~GC08、低潮DC01~DC08),断面尽量设计在变化波动最大的区域,鉴于高、低潮时变化位置有差别,部分断面位置在高、低潮时有调整,如图2所示。以最早一期的瞬时水边界作为基线,断面与基线交点为基点,分别计算高、低潮时其他各个时期水边界与断面交点(简称“断面交点”)至基点的距离,从而计算出断面上各个时期淤积或侵蚀的距离以及平均速率。
2.2.3 入海口位置及方向计算
黄河入海口的迁移变化是整个黄河三角洲海岸带形态变化最典型、最直观的反映,已有的研究重点则更多集中在海岸线的变化上,如杨江平等[1]在研究中定性给出了黄河入海的位置及方向。
为便于分析研究,本研究将河口处水边界的中点作为入海口位置,河口水边界上游1 km处河道中心点至河口水边界中心点所形成矢量的地理方位角作为黄河的入海方向,如图3所示。根据瞬时水边界的解译结果,进一步解译并定量计算黄河入海口的迁移过程。存在沙洲时,为了突出河道的变化,将沙洲与陆地相夹形成的河道作为入海口计算其入海方向。本文所说的入海口位置是指入海口处遥感解译水边界的中点,并不具有严格的地学意义。
图3 入海口位置及方向解译示意
3.1.1 总体变化情况
根据基期和末期的瞬时水边界矢量信息,通过线-面转换、矢量叠加分析,得到黄河三角洲海岸带区域淤积、侵蚀空间分布状况,如图4(a)所示。高潮时,淤积面积66.39 km2,年平均淤积速率3.91 km2/a,主要集中在三角洲东北部区域;侵蚀面积64.46 km2,年平均侵蚀速率3.79 km2/a,主要集中在三角洲西北、东南及南部区域。低潮时,淤积面积90.36 km2,年平均淤积速率5.32 km2/a,主要集中在三角洲北部、东部和西南部区域;侵蚀面积33.50 km2,年平均侵蚀速率1.97 km2/a,主要集中在老黄河口以及西北部区域。
为了比较海岸带的变化情况,根据2004年、2020年高潮、低潮瞬时水边界形成两期海岸带区域,如图4(b)所示,2004年海岸带面积82.99 km2,2020年海岸带面积83.73 km2,总面积并未发生明显变化,但在空间分布上发生了较为明显的变化,北部区域海岸带变窄,现行河口区域海岸带整体向海迁移。东部区域海岸带也整体向外迁移,但高潮、低潮线还能明显区分。老黄河区域高潮、低潮线整体向陆地迁移,海岸带宽度未发生明显变化。南部及西南部区域海岸带宽度明显变大。海岸带变化定量表征见3.1.3节。将本研究成果与黄河口区域已有的一般高潮线法[1]和低潮线法[2]成果进行比较,鉴于研究时间和范围的不同,本文仅对部分断面进行定性比较,结果见表3。可以看出,本研究方法与一般高潮线法和低潮线法在整体上具有相同的结论。
图4 淤积、侵蚀空间分布
表3 与已有研究成果的对比
3.1.2 区域变化情况
计算高潮、低潮时各断面交点至基点的距离d1、d2,保留整数米,见表4、表5,正值代表向海淤积,负值代表向陆侵蚀。结合高潮、低潮时8条断面交点至基点距离变化图,可以清晰看出黄河口区域淤积-侵蚀变化过程,如图5所示。
图5 断面交点至基点距离变化
表4 高潮断面交点至基点距离d1 m
表5 低潮断面交点至基点距离d2 m
第2、3、6条断面所在区域一直处于淤积状态,且第3条断面淤积速度最快。高潮时,其最大淤积距离为10.30 km,平均淤积速度0.61 km/a;低潮时,其最大淤积距离为9.62 km,平均淤积速度0.60 km/a。高潮、低潮时淤积趋势、速度相似。
第1、7条断面所在区域均处于逐渐侵蚀的状态,且第7条断面侵蚀速度最快。高潮时,其最大侵蚀距离为4.51 km,平均侵蚀速度0.27 km/a;低潮时,其最大侵蚀距离为3.91 km,平均侵蚀速度0.24 km/a。高潮、低潮时侵蚀趋势、速度相似。
第4条断面在高潮时呈现逐渐淤积的状态,低潮时先侵蚀后淤积。第5条断面在高潮、低潮时均为先淤积后侵蚀的过程,但高潮时整体处于淤积状态。第8条断面在高潮时处于波动侵蚀的状态,而在低潮时呈现逐渐淤积的状况,这导致该区域潮间带的宽度被拉大、滩涂面积增大。具体见图4。
3.1.3 海岸带变化表征指标
基于上述分析,提出表6所列指标,用于表征区域海岸带的变化趋势。
根据表6中的指标定义,计算黄河三角洲2003—2020年海岸带变化指标,结果见表7。总体指标中海岸带高潮时向陆迁移总面积为64.46 km2,向海迁移总面积为66.39 km2;低潮时海岸带向陆迁移总面积33.50 km2,向海迁移总面积90.36 km2。有3条断面高潮线向陆侵蚀,最大侵蚀距离4.51 km,5条断面高潮线向海淤积,最大淤积距离10.30 km;3条断面低潮线向陆侵蚀,最大侵蚀距离3.91 km,5条断面低潮线向海淤积,最大淤积距离9.62 km。海岸带宽度在断面8处变宽达8.13 km。海岸带中心点在断面3处向海迁移达9.96 km,在断面7处向陆迁移4.21 km。空间分布变化状况见图4。
表6 海岸带变化表征指标
表7 黄河三角洲海岸带形态变化结果
根据韩香举等[15]关于黄河三角洲冲淤时空演变的研究成果,在2009—2014年、2017—2018年现行河口口门处淤积量大幅增加。根据本研究影像解译结果,2013年11月高潮时在入海口处出现面积为0.31 km2的沙洲,河道分叉,随后向东北侧发育,如图6所示,2020年已发展为面积6.95 km2的稳定滩涂,滩面逐渐被互花米草稳固。2020年起,在北侧区域发育形成第二片沙洲,面积5.12 km2,再次分叉河道。
图6 沙洲面积变化
低潮时出现了类似的情况,但沙洲出现的时间较早,面积0.28 km2,在空间上较高潮时的沙洲偏北2.5 km,低潮时沙洲的演变方向为东南向,与高潮时不同。发展到2020年,面积达到11.53 km2,北部沙洲在低潮时面积也达到8.04 km2。
关于黄河入海口位置、方向变化定量精确计算的研究成果较少,为了便于分析入海口迁移变化过程,根据高低潮瞬时水边界的解译结果,按照2.2.3节的方法计算高潮、低潮时入海口位置及方向,结果见表8。东向、北向坐标为高斯投影(3°带,中央经线120°)下的平面坐标。高潮、低潮日期后的(1)(2)代表分叉河道的入海口位置,高潮时有3个时期出现河道分叉,低潮时有5个时期出现河道分叉。
表8 黄河入海口位置、入海方向
根据图7所示黄河入海口时间、空间迁移变化情况可知,高潮时,2007年之前入海方向在58°~99°之间,入海口向海推进约5 km。2007年9月解译结果显示,黄河河道向北分叉,分别向349°、33°方向入海,此后,在7°~26°形成主河道。2013年河口出现小规模沙洲,2015年沙洲规模变大,使得河道明显向5°和74°两个方向分叉。2020年北向河道再次向72°和331°方向分叉,至此黄河入海口处形成3处入海河道,比2003年河道分别向东北方向迁移约9.3、7.0 km,向西北方向迁移约11.4 km。
图7 入海口时间、空间迁移变化
低潮时,2007年入海口位置与2004年相比变化不大,2008年改向8°方向流出,空间位置上比2004年向西北迁移约6.2 km。2012年出现小规模沙洲,2014年沙洲面积扩大,形成较为明显的分叉河道,分别向88°和357°方向流出。此后,河道继续沿这两个方向向海推进,2020年北向河道再次发生分叉,分别向90°和336°方向流出,至此低潮时也形成3处入海河道,比2004年河道分别向东北方向迁移约8.9、5.4 km,向西北方向迁移约12.3 km。
本研究基于黄河三角洲历史影像以及潮位数据,解译三角洲区域高潮、低潮瞬时水边界,通过计算典型断面变化量,分析海岸带、新生沙洲以及入海口位置、方向变化趋势,提出表征海岸带变化的指标。研究结果表明,利用高潮、低潮两方面数据分析能够全面表征2003—2020年黄河三角洲海岸带形态演变趋势。利用海岸带变化指标能够定量计算区域内海岸带的具体变化量,证明该指标在海岸带演变的研究工作中能够发挥作用,笔者将在今后工作中继续完善、推广。