购买
下载掌阅APP,畅读海量书库
立即打开
畅读海量书库
扫码下载掌阅APP

1.2 研究现状

1.2.1 水体遥感理论基础

遥感(remote sensing),即遥远的感知,是指在不接触被观测目标的情况下,通过获取其反射、辐射或散射的电磁波(光、热、无线电等)、力场(重力、磁力等)、声波等信息,并进行分析应用的一门科学与技术 [36] 。不同的被观测目标具有不同的波谱特征,这是遥感的基本原理。

太阳光在水体中辐射传输过程的结果即对应水体波谱特征,即水体的表观光学特性(表观光学量,apparent optical properties,AOPs)和固有光学特性(固有光学量,inherent optical properties,IOPs)。水体的表观光学量会随光照条件及外部环境因素的变化而变化;而固有光学量仅与水体组分有关,主要包括水分子的吸收系数、散射系数和散射相函数,水中叶绿素的吸收系数、散射系数和前向后向散射系数,有色可溶性有机物的吸收系数,悬浮泥沙的吸收系数、散射系数和前向后向散射系数。总之,水体遥感反射信息包含了来自于底质、水体中内部光源、水体表面的镜面反射以及水体组分(悬浮泥沙、叶绿素、有色可溶性有机物)的吸收和散射等多方面特征。而各类传感器接收到的水体遥感辐射通量还包含了太阳光在大气传输过程中受到的大气分子、水汽和气溶胶等的散射、吸收和漫反射等的影响,如图1.1所示。

图1.1 太阳光在水体中的辐射传输过程 [37]

对太阳光在水体中的辐射传输过程而言,国内外专家学者以辐射传输理论为基础,通过大量的研究,建立了水体表观光学量与固有光学量之间定量的数学关系-辐射传输方程 [38,39] 。该方程严谨地描述了太阳光在水体中传播的多重散射和吸收过程,如式(1.1)所示。

式中: L ϕ φ , z )为离水辐射; 为散射相函数; c z )为衰减系数; S ϕ φ , z )为内部光源。

考虑到辐射传输方程所涉及的变量较多、形式复杂、不易实用的特点,许多研究基于一定的假设条件,对辐射传输方程做了许多简化和改进,得到了广泛的应用。

大量的水体悬浮泥沙现场试验和研究结果表明,在一定的悬浮泥沙含量范围内,水体在可见光波段的反射率会随着悬浮泥沙含量的增加而增加,超过一定的范围之后水体在可见光波段的反射率就会保持稳定或变化很小;而近红外及更长波波段则对高浓度悬浮泥沙水体敏感 [4,12,40] 。随着悬浮泥沙含量的增加,水体遥感反射率反射峰在向长波方向移动,反射峰宽度逐渐变窄。一般地,水体光谱反射率最大值小于7% [41,42] 。现场水体试验时,采用一定的测量方法可使水体的遥感信息中水体组分的贡献占绝大部分比例 [36] 。图1.2为现场水体实验实测结果,可以看出不同悬浮泥沙含量的水体光谱特征差异显著。

图1.2 现场水体光谱试验(不同悬浮泥沙含量的光谱反射率)

对水色水体遥感而言,常用的遥感信息集中在可见光到中远红外之间(见图1.2和表1.1)。既有的卫星遥感平台携带了不同的遥感传感器,这些遥感传感器在波段范围,时间、空间和光谱分辨率等方面各有异同,能够获取包括水量、深度、温度、水色及水体组分含量等多种水体信息。

(1)水量,即水体的空间分布(面积)信息,是相关研究工作的基础。简单来讲就是要准确地区分水陆边界。根据实际需求,在不同时空尺度上识别水量信息所利用的遥感数据也不同。全球或区域尺度的水体分布多基于MODIS(moderate resolution imaging spectroradiometer)、AVHRR (advanced very high resolution radiometer)等中低空间分辨率的遥感影像数据;而内陆水体(湖泊、河流、水库和人工水域等)或河口近岸的水量信息监测则多采用中高空间分辨率的遥感影像数据,如Landsat [43] ,SPOT [44] ,HJ星数据,GF数据等。

表1.1常用遥感电磁波谱

(2)深度,是水体一个重要的基础信息,与社会生产生活密切相关。水体深度信息常采用激光雷达或微波等手段获取,相关方法也多种多样。本书中涉及的水体深度信息相对较少,不作深入论述。

(3)温度,水温深刻地影响着发生在水体中各种生物物理化学过程,特别是进行光合作用的藻类、浮游植物等的生长状态 [45] 。相关研究中,常用各类卫星遥感传感器的热红外波段数据来反演水体温度。

(4)水色,即水体的颜色。水体组分,指水体关键参数含量,包括浮游植物、悬浮泥沙和有色溶解有机物(chromophoric dissolved organic matter,CDOM),盐度、富营养化参数(N、P、NH 3 -N)等。悬浮泥沙既是水体水色的主要影响因素之一,同时又是关键的水体组成成分。因此,悬浮泥沙是最为显著的水生环境影响要素。

1.2.2 河口海岸悬浮泥沙遥感研究进展

自遥感技术应用于河口海岸泥沙监测研究以来,国内外专家学者纷纷利用Landsat、SPOT、SeaWiFS、AVHRR、MERIS、MODIS、HICO、FY、HJ和GF等多源卫星遥感数据,结合现场实测数据在全球许多区域,进行了大量的研究工作,如美国的Apalachicola海湾 [5,46-48] ,法属Guyana海湾 [49] ,比利时的Scheldt河 [49] ,法国的Quiberon海湾 [50] 、Gironde河 [8,51,52] ,南美洲的la Plata河口海岸 [53] 、Mekong河 [54-56] 、Bassac河口 [56] 、Amazon流域 [57,58] 、Solimões河 [59] 、格棱兰岛近海 [60] 、Mississippi河及其支流 [61,62] ,加拿大的Minas河流域 [63] 、Saint John河 [64] 、 Macuse河口 [65] ,印度的Chilika湖 [66] 、Brahmani河 [67] ,以色列的Liverpool海湾 [68] ,伊朗的Chabahar海湾 [69] ,Burkina Faso的Bagre水库 [70] 、Danube河口海岸 [71] ,巴西的Tiete河流域 [72] ,我国的长江 [73-83] 、黄河 [84-86] 、珠江 [87-90] 以及沿海海域等 [4,13,91-96]

1.2.2.1 河口海岸悬浮泥沙遥感反演模型

根据悬浮泥沙模型建立的理论基础可将其分为三类:经验模型,辐射传输理论模型,半分析/半经验遥感反演模型 [30,94,97-99]

经验模型是在现场实测和采样数据基础上,基于水体辐射量(遥感反射比或离水辐射亮度等)与水体中悬浮泥沙含量之间的数学统计关系而建立起来的 [100] 。即通过对同步或准同步遥感数据与实测悬浮泥沙数据的相关性分析,确定两者间的相关系数及大致形式,进而建立经验模型;或者通过实验仪器测量多通道水体光谱曲线、水体表观光学量、大气光学特征以及水体中各组分的浓度而建立的地面光谱模型 [98,101] 。在悬浮泥沙研究的早期和当前,经验模型的应用都非常广泛。

基于现场实验结合同步(准同步)Landsat系列卫星影像的利用经验模型研究水体悬浮泥沙含量的有很多 [35,69,72,90,96,102] ,如Zhang等在黄河河口 [86] 、钟凯文等在珠江流域西江干流 [103] 、张毅博等在新安江水库 [104] 、Zheng等在洞庭湖 [105] 都反演分析了研究区悬浮泥沙含量;Vantrepotte等在法国的Guiana海湾 [106] 、Park等在亚马孙河流域借助MODIS数据分析了研究区悬浮泥沙的时空变化 [58] ;利用其他数据源进行悬浮泥沙反演的研究还有很多,如MERIS数据 [56,107] 、HJ星数据 [108] 、Hyperion影像 [109] 、GOCI影像 [73] 、航空高光谱数据 [61] 、Radarsat数据 [110,111] 和实测数据 [112,113] 等;Wu等还发展了一种在缺少星、地同步实验情况下的悬浮物遥感反演方法 [6]

经验模型易于操作实现,反演模型参数所需要数据源比较容易获取得到,可适用于大多数便于进行试验测量的区域水体,并且在特定的范围内可获得较高的反演精度,应用广泛。但是,经验模型的参数多具有区域性和季节性的特点,所以其普适性不足 [16,99] 。因经验模型的参数受研究对象的时空异质性影响较大,所以在推广应用到其他区域、时段的时候,需要进行重新检验验证 [97]

辐射传输理论严谨地描述了表观光学量、水体固有光学量(吸收、散射、光束衰减)在光场过程中的作用。辐射传输理论模型建立一般分为两个过程:以上行辐射与水体中光学活性物质特征吸收和后向散射特征之间的关系为基础,利用水体反射率反演水体中各组分的特征吸收系数和后向散射系数(表观光学量与固有光学量的过程);通过水体中各组分浓度与其特征吸收系数、后向散射系数相关联,反演水体中各组分的含量(固有光学量与水体中各组分浓度过程)。

基于辐射传输理论,已有许多学者 [114-116] 分析了水体表观光学量与固有光学量的关系,即发展推导了遥感反射比与水体总吸收系数、后向散射系数之间的函数关系,见式(1.2)。

式中: R 为遥感反射比; f Q 为太阳天顶角函数; t 为水气界面透射率; n 为水体折射系数, b b λ )为水体总的后向散射系数; a λ )为水体总吸收系数。

式(1.2)中, a λ )可以看作是纯水[ a w λ )]、CDOM [ a g λ )]、浮游植物[ a ph λ )]和悬浮物[ a x λ )]的吸收系数的线性组合;而 b b λ )主要包括悬浮物[ b bx λ )]和纯水[ b bw λ )]的后向散射,见式(1.3)和式(1.4)。

联合式(1.2)~式(1.4),经过波段代入求算水体固有光学量,进而反演水体各组分含量。

Giardino等在意大利Garda湖尝试了基于辐射传输理论模型的水质参数反演,重点关注了叶绿素和浑浊度 [117] ;Zhang等、孙德勇等在我国太湖 [118-119] ,戴永宁等在巢湖 [120] ,陈莉琼等在鄱阳湖 [121] 分别开展了辐射传输理论和水体固有光学特性的研究;李云梅等采用基于辐射传输理论的分析模型反演了太湖悬浮颗粒物浓度 [122] 。相关研究成果都为水质参数遥感监测的理论发展和实践提供了重要的参考和帮助。

辐射传输理论模型具有明确的物理意义和良好的普适性 [99] 。然而水体参数的反演是一个非常复杂的物理过程,辐射传输方程需要输入大量的固有光学参数以及地表参数。而目前这些参数的获取比较困难,辐射传输方程计算起来也比较复杂,限制了辐射传输方程的实际应用。鉴于辐射传输理论模型一些参数难以获取以及计算量大的问题,后人发展了简化的半分析模型。

半分析/半经验模型介于经验模型与辐射传输理论模型之间。半分析/半经验模型是在辐射传输理论模型的两个主要过程的某一阶段简化利用了经验统计分析方法的模型。一种形式是,首先建立固有光学量(吸收系数和后向散射系数)与悬浮泥沙的经验关系,而固有光学量则是基于辐射传输理论通过表观光学量(水体反射率)求解;另一种形式则相反,固有光学量通过表观光学量(水体反射率)通过近似简化或一些经验关系逐步逼近求解,而悬浮泥沙含量则通过水体中各组分浓度与固有光学量相关联建立的理论模型求解得到;半分析模型都是某一形式的两个过程相结合,最终建立表观光学量与悬浮泥沙的关系模式。因为半分析/半经验模型是简化了的辐射传输理论模型,有一定的理论基础和较为明确的物理意义,相对简便,应用也愈加广泛。

Chen等研究了我国渤海、长江入海口、浙江省瓯江入海口3个不同混浊度水体的泥沙含量 [123] ,建立了一个基于HJ-1A/CCD数据三波段组合的泥沙反演半分析模型,整体误差小于29%。Mao等为了改善常规模型不能准确测定我国东海高悬浮泥沙含量的问题 [13] ,以复杂的光谱指数代替生物光学模型输入中单一的波段。结果显示,光谱指数与泥沙含量的相关系数达0.912,高于任何单一波段的结果,模型反演结果的平均相对误差为23%。Nechad等发展了一个适用于多种海洋卫星遥感数据悬浮泥沙反演的生物光学模型 [124] 。仅以红光到近红外区间的单一波段作为模型输入,平均相对误差可优于30%。陈燕等基于渤海湾的站点观测数据,结合该区域水体光谱特性,首先利用多波段半分析算法(quasi-analytical algorithm,QAA)计算了研究区水体的固有光学特性,建立了渤海湾悬浮泥沙遥感反演的半经验模型 [17] 。结果表明,该模型的精度更高、普适性更好,相对误差18.28%。

Zhang等基于太湖现场测量光谱和悬浮物浓度数据,利用近红外两个波段开发了适用于较深、悬浮物含量高的水体的生物光学模型 [16] ,验证结果平均相对误差仅为13%。Sun等通过对我国太湖、巢湖、滇池、三峡库区等内陆水体的水体特性实地观测,发展了反演中国内陆水体悬浮物组分的半分析模型 [125] ,并结合新型的光学卫星数据HJ-1A/HIS影像分析了太湖水体的悬浮颗粒组成及比例。Giardino等以航空高光谱数据为基础,结合实测数据与生物光学模型 [126] ,研究了意大利中部一个浅而浑浊的Trasimeno湖水体悬浮物浓度和水生植物的空间分布,发现水生植物的存在是悬浮物保持低浓度积极因素。金鑫等 [127] 根据实测数据,在确定了水体固有光学参数的基础之上,构建了巢湖的悬浮物浓度生物光学遥感反演模型。结果表明,随着悬浮物浓度的增加,模型精度越高,平均相对误差为17.25%。张红等在巢湖的研究也有类似发现 [128] 。徐京萍等分析了吉林省长春市石头门水库400~1200nm范围内的水体光谱特性,发现高浓度悬浮物含量对水体总的反射率贡献较大 [97] ,并分别以 808nm、873nm和1067nm处的光谱反射率建立了反演悬浮物浓度的生物光学模型。施坤等根据生物光学原理,基于对太湖、巢湖、三峡库区实测数据的分析 [99] ,建立了浑浊湖泊水体总悬浮物浓度的单波段估算模型,平均相对误差低于24%,均方根误差低于18mg/L。

半分析/半经验模型结合了经验模型与辐射传输理论模型的优点,普适性较好且具有一定的物理意义。半分析模型不需要大量的实地采样数据作为支撑,也不同于辐射传输理论模型需要大量的固有光学参数、地面参数输入,根据已有的离水辐射高度或遥感反射比加上水体中各组分的吸收和散射信息就可以直接较为准确地估算出各组分的含量。

然而,半分析/半经验模型也存在一定的不足之处 [30] 。一般而言,河口海岸水域深度较深,水动力较强,水体更新周期短,水生环境复杂易变,且受河流上游、沿岸生产活动、洋流潮汐等诸多因素的影响;内陆湖泊、水库则不同,相对而言其水深较浅,水动力较弱,水体更新周期长,水生环境简单稳定,对其影响因素相对较少。建立半分析模型需要已知水体各组分的吸收系数、后向散射系数、体散射系数等相关参数,而这些参数需要较多的一致度、逼真度高的数据模拟得出。河口海岸和内陆湖泊、水库自身的特点决定了其可被模拟性、逼近程度、参数本身的复杂度、可被认知程度及实验获取难易等问题。显然,半分析/半经验模型在内陆湖泊、水库应用更多,也更为成功。对此,最大浑浊带作为河口海岸的一个“较为稳定”的动态变化区域,几乎是所有相关因素在河口海岸的综合作用“信息承载体”,拓展半分析/半经验模型的应用可以考虑优先开展河口海岸最大浑浊带的理论分析、现场实验、数据模拟、参数测定等相关研究。加强对最大浑浊带的研究,不仅可以直接准确把握、分析河口海岸悬浮泥沙的动态变化,而且能够促进对其他水质参数的遥感反演研究。另外,半分析/半经验方法在模型建立求解的一些近似或简化过程中,相关前提假设条件缺乏相应的验证,可能会带入未知的误差。

悬浮泥沙遥感反演依据不同的理论,过程各有异同,数据源也不尽相同,包含波段的各种变换组合(如差值、比值、对数等)形式也千差万别,但最终仍可以归结为以下六种。

(1)线性形式。线性形式是最常见的泥沙含量估算形式之一,包括一元和多元的,适合于简单水生环境的情况。其一般形式如式(1.5)所示。

式中:TSS (total suspended solid)为泥沙含量; R 1 R 2 为光谱信息; a b c 为待定常量参数。

本书以下其他公式中所列参数含意与式(1.5)中相同。 R i i =1,2,3,…)可以是遥感影像某波段数据,也可以是高光谱或多光谱仪器测量值,或者是由其变换而成的其他组合。

(2)高次形式。高次形式是线性形式的延伸[见式(1.6)],提高了其应用范围。

(3)对数形式。对数形式包括对光谱数据和泥沙数据分别或单一作对数变换,在低泥沙含量的水体中反演精度更高。对数变换常见的有以10为底数的变换和以自然常数e为底数的变换[见式(1.7)],并且可以代入多组光谱数据组合计算。

(4)Gordon形式。Gordon形式由准单散射理论近似得到,比较著名、多被提及,见式(1.8)。但是,该形式假设水体光学性质稳定均一,因此反演精度不高、应用较少。学者对Gordon方法的适用性也存在不同见解,何青等的研究结果表明Gordon方法适合于低含量泥沙水体或者高含量泥沙水体 [129] ,而陈勇等则认为该方法对中等含量悬浮泥沙水体反演效果较好 [130] 。这样的差异可能是由于学者对悬浮泥沙高、低浓度的定性描述不同而引起的。

(5)指数形式。指数形式考虑了水体组分对光学性质的影响,较真实地反映了光谱反射率与泥沙的相关关系,应用范围广泛,在各类水体都有较好的适用性 [12] 。常见的指数形式多以自然常数e作为底数,见式(1.9)。

(6)综合形式。综合形式是对上述形式的组合,相关研究较少 [131] 。当然,同样存在以上各反演模式的反函数形式。

需要指出的是,以上6种悬浮泥沙反演模型同时又可以进一步归结为单波段模型和多波段组合模型两类 [30]

悬浮泥沙单波段遥感反演模型是指仅利用单一波段数据来反演悬浮泥沙含量。此类模型的遥感数据来源较多,因此在很多区域都有应用。然而,由于不同波段对不同悬浮泥沙含量水体的敏感性有很大的差别,此类模型难以适用于悬浮泥沙含量较大范围变化的水体 [35] 。例如,许多研究已经表明红光波段遥感反射率会随着水体悬浮泥沙含量的增加而升高,但超过一定的浓度范围(在高悬浮泥沙含量的水体),红光波段遥感反射率则不再升高,出现饱和现象 [12] 。与红光波段相比,近红外波段对高悬浮泥沙含量水体更为敏感 [4,12,40,46] 。因此,基于单一波段的悬浮泥沙遥感反演模型在悬浮泥沙含量大范围变化的水体应用时有较大的局限性。

悬浮泥沙多波段组合遥感反演模型是指利用多个波段数据来反演悬浮泥沙含量,可以较好地避免遥感反射饱和现象 [12,132-134] 。常见的多个波段的组合形式有波段比率 [132,135-137] 、泥沙指数和其他复杂的形式 [76,104,133,134,138] 。此类模型所对应的模型形式如上面提到的线性函数、指数函数、对数函数以及二次函数等。需要注意的是此类多波段组合模型的形式多数为单调函数。应用单调函数性质的悬浮泥沙遥感反演模型时,存在两个潜在的问题 [30,35] :一是,由于单调函数的性质,使得悬浮泥沙的反演结果随着遥感反射率的变化呈现固定的变化(如线性函数),这可能与实际情况不符;二是,对于指数函数或对数函数而言,在一定的区间,遥感反射率的微小改变可能会造成悬浮泥沙结果的巨大变化,产生过拟合现象。基于以上分析可知,非单调性质的多波段组合模型可能是悬浮泥沙遥感反演的最恰模型形式 [30,35] 。尽管基于遥感方法的悬浮泥沙反演研究已取得丰富的成果,但是,目前仍然缺乏一种可应用于多个河口海岸悬浮泥沙含量显著分异、精度高、鲁棒性好的定量遥感反演模型 [34,35]

1.2.2.2 河口海岸悬浮泥沙主要研究内容

根据既有研究的具体内容可将相关工作分为最大浑浊带、特殊气象条件效应和时空变化特征三个方面的研究。

(1)最大浑浊带研究。最大浑浊带是河口近岸一个动态变化的区域,与来水输沙、潮汐、沉积作用密切相关 [139] ,其悬浮物浓度稳定高于河口近岸其他区域,在河口近岸的物质循环过程中起到了非常重要的作用 [75,140] 。实测资料表明,最大浑浊带的范围与河口盐水楔锋面的进退区域及滞流点的上、下位移范围相对应 [141] 。然而,国内外相关研究较少,仅在Yenisei河 [140] 、长江 [75,142-44] 相关研究中提及,专门研究的有黄河 [145] 、鸭绿江 [146] 、珠江 [147] 。但是,专门研究也多为较早时期的结果,难以满足对包括位置、范围、含沙浓度、垂直分布等在内的最大浑浊带的特征进行全方面、全方位的研究,而这些特征受到河口海岸水动力与水化学环境,如潮流、径流作用,盐淡水混合作用等时间、空间变化的密切影响。

(2)特殊气象条件效应研究。在暴雨、台风等特殊气象条件下,研究区多被云层覆盖,卫星传感器很难有效获取到地表相关信息,造成此时段内研究区遥感数据源缺失或存在质量等问题,致使基于遥感方法研究悬浮泥沙在特殊气象条件效应下变化特征的机会非常难得。因此,国内外相关报道较少。其中,Chen等、Huang等分别在不同地区研究了飓风对海湾水质参数的影响 [48,148] ;Chen等 [47,149] 、Michael等 [150] 分析了水体泥沙含量随暴雨过程的变化;另外,还有研究监测了海水入侵对河口混浊度的影响 [151] 。我国自然灾害频发,特别是东南沿海一带,受极端天气影响尤为严重,尽管受限于客观条件等因素,当前相关研究比较缺乏 [152] ,所以更需要重点关注。既有的情景分析和数值模拟等方法值得借鉴 [22,153-155]

(3)时空变化特征研究。基于情景分析(又称物理模型、物理实验)和数值模拟(又称计算模型、数学模型)对泥沙运动输移、趋势变化的研究由来已久 [156-163] 。而当前悬浮泥沙时空演变遥感研究相关成果略显不足,亟待加强。已有研究多为在短时间或中等时间跨度上利用不同方法结合不同数据对不同区域的研究 [12,20,56,70,86,89,106,164-166] ,对把握泥沙时空演变规律有重要的借鉴意义和参考价值 [53,59,71,74,93,167-170] 。但是,因其时间跨度相对较短,难以发现悬浮泥沙含量的长期变化规律和趋势,并且其相关研究结论的准确性、可扩展性还有待讨论和进一步验证。虽然长江河口近岸已有近40年时间序列的研究结果 [130] ,但因卫星遥感影像获取或影像质量等问题造成其数据密度极低,可能会造成一些重要信息的缺失。

1.2.2.3 河口海岸悬浮泥沙遥感研究总结

目前,基于站点监测和遥感数据反演水体悬浮泥沙含量的相关研究已有很多。尽管传统的水文站点和行船作业调查具有较高的精度和一定的连续性,但耗时费力且代表范围有限,实效性低,当前多用于精确的对比验证。整体来看,基于遥感方法反演水体悬浮泥沙的相关研究虽然是水色遥感中方法最成熟、成果最丰硕的,包括了现场光谱数据采集、水体样品实验室分析、遥感影像预处理、泥沙反演模型的建立、泥沙时空变化相关因素关联关系分析、水体泥沙反演结果的分析及应用等方方面面。但是,既有的悬浮泥沙定量遥感模型的普适性不足,并且关于悬浮泥沙时空演变规律(最大浑浊带、特殊气象效应、时空变化特征)的研究和讨论还远不充分,亟待加强。

因此,建立一个可适用于多个河口海岸区域、水体悬浮泥沙含量大范围波动的定量遥感反演模型,并结合长时间序列遥感影像、大量现场水体实测数据及其他相关数据,准确反演区域内水体悬浮泥沙的空间分布及时序变化特征,分析其影响因素及响应机制,对理解和把握河口海岸悬浮泥沙时空演变规律具有重要意义!不仅可以为河口海岸区域的生态环境保护、资源科学开发利用提供重要的科技支撑和决策依据,同时也将有助于我们提高对全球物质循环、全球变化和气候变化等科学问题的认识和理解 [30,35,90,96] CAPXlbIuz3hMpcuQ6nkdAUgtOXWJMVyxqhXCzfNDvqvJJrNOzA3xqgn1/WKFZu1n

点击中间区域
呼出菜单
上一章
目录
下一章
×