ISSN 1003-8035 CN 11-2852/P
  • 中国科技核心期刊
  • CSCD收录期刊
  • Caj-cd规范获奖期刊
  • Scopus 收录期刊
  • DOAJ 收录期刊
  • GeoRef收录期刊
欢迎扫码关注“i环境微平台”

基于GIS和GDAL的地质灾害汛期风险预警产品信息自动化生成算法的构建与实现

张像源

张像源. 基于GIS和GDAL的地质灾害汛期风险预警产品信息自动化生成算法的构建与实现[J]. 中国地质灾害与防治学报,2022,33(5): 76-82. DOI: 10.16031/j.cnki.issn.1003-8035.202205036
引用本文: 张像源. 基于GIS和GDAL的地质灾害汛期风险预警产品信息自动化生成算法的构建与实现[J]. 中国地质灾害与防治学报,2022,33(5): 76-82. DOI: 10.16031/j.cnki.issn.1003-8035.202205036
ZHANG Xiangyuan. Construction and implementation of an automatic algorithm for generating information of geological disaster floor risk warning products based on GIS and GDAL[J]. The Chinese Journal of Geological Hazard and Control, 2022, 33(5): 76-82. DOI: 10.16031/j.cnki.issn.1003-8035.202205036
Citation: ZHANG Xiangyuan. Construction and implementation of an automatic algorithm for generating information of geological disaster floor risk warning products based on GIS and GDAL[J]. The Chinese Journal of Geological Hazard and Control, 2022, 33(5): 76-82. DOI: 10.16031/j.cnki.issn.1003-8035.202205036

基于GIS和GDAL的地质灾害汛期风险预警产品信息自动化生成算法的构建与实现

详细信息
    作者简介:

    张像源(1970-),男,甘肃兰州人,硕士,教授级高级工程师。主要研究方向为地质环境信息化、系统研发和地质灾害预警研究。E-mail:14561996@qq.com

  • 中图分类号: P694

Construction and implementation of an automatic algorithm for generating information of geological disaster floor risk warning products based on GIS and GDAL

  • 摘要: 我国大陆从2003年开始开展了汛期区域地质灾害气象预警工作,并形成预警产品面向公众的发布,在地质灾害防治中发挥着重要的作用。而地质灾害预警描述信息是地质灾害预警产品的重要组成,包括了预警范围、等级、时段和文字说明等内容,为了快速形成权威、科学、精准的这一信息,文章研究了一种可高效自动分析形成初步的预警产品描述信息的技术算法,并对其进行了验证分析。算法以研究区行政区划图和预警区划图为数据源,以设计的三级预警等级为主键的预警等级集合概念模型为基础,采用GIS技术,通过空间信息叠加与判别、预警等级分布面积占比求算、分布范围描述和单元分布方位求算等为关键方法,将图面内容快速地转换为直观的文本描述,可实现预警产品中描述信息的自动化获取。在验证环节,设计了11个虚拟乡镇为预警范围,预警区划图层以3 km×3 km为预警等级单元,预警单元总数为5194个,文件为Shapefile格式;验证过程基于GDAL库(Geospatial Data Abstraction Library)和C++语言,按照算法设计了软件工具,利用实验数据进行了验证与分析,结果表明,通过该算法形成的预警描述信息和预警区划图层的空间展现结果完全吻合,且计算过程高效快捷,可显著提高了预警精细化程度和工作效率,完全适合在实际工作推广使用。
    Abstract: Since 2003, mainland China has carried out regional geological disaster early warning work during the flood season and formed early warning products to be released to the public, which plays an important role in geological disaster prevention and control. The geological disaster warning description information is an important component of the geological disaster warning product, including the warning range, grade, time period, and text description. In order to quickly form authoritative, scientific, and accurate information, this paper studies a technical algorithm that can efficiently and automatically analyze and form preliminary early warning product description information and conducts verification analysis on the algorithm. The algorithm uses the administrative division map and early warning division map of the study area as data sources.Based on the designed three-level warning level as the main key of the set concept model of the warning level, Using GIS technology, through the use of spatial information superposition and discrimination, the calculation of the distribution area ratio of the warning level, the description of the distribution range of the warning level, and the calculation of the unit distribution orientation of the warning level as key methods, the information on the map can be quickly Converted to an intuitive text description, the automatic acquisition of early warning product description information can be realized.In the verification process, 11 virtual townships are designed as the warning range. The warning zoning layer takes 3×3KM as the warning level unit. The total number of warning units is 5194, and the file is in Shapefile format. The verification process is based on the GDAL library (Geospatial Data Abstraction Library)and the C++ language. A software tool is designed for the algorithm, and the experimental data is used for verification and analysis. The results show that the early warning description information formed by the algorithm is completely consistent with the spatial presentation results of the early warning zoning layer, and the calculation process is efficient and fast, which can significantly improve the precision of early warning. It is fully suitable for promotion and use in practical work.
  • 我国地质灾害种类多、分布广、危害大。为进一步减轻地质灾害风险,最大限度保障人民群众生命财产安全,从2003年起,我国大陆开始开展汛期区域地质灾害气象预警工作,并形成预警产品面向公众的发布,在地质灾害防治中发挥了 “消息树”和“发令枪”的作用,取得了明显的社会经济效益[1],极大地提升了公众社会对防范地质灾害的认知。随着该项工作的不断推进,预警产品的内涵从粗到细逐步走向规范化,包括了预警范围、等级、时段和文字说明等内容[2],且临灾发布工作要求快速、高效和精准。然而这些信息若靠人工获取并不是件易事,比如要经技术人员通过读取预警区划图并研判后,才能粗略得出空间分布范围等定性信息,这一过程耗时耗力,且难以达到定量描述的准确效果。因此自动计算预警结果并快速生成预警产品是预警系统的一个重要功能[2],围绕聚焦解决好地质灾害可能发生的地点、成灾范围等预警预报问题[3-4],考虑如何利用信息技术自动获取预警产品描述信息的实现势在必行。然而经检索发现国内直接进行相关研究还较少,可供参考的文献不多。为了快速形成权威、科学、符合实际的这一产品,本文提出一种可高效自动分析形成初步的预警产品描述信息的技术算法,为有关信息平台的功能模块的研发、支撑专家做出更详尽的预警产品研判提供理论支撑。

    以研究区行政区划图和预警区划图[2]为数据来源,应用GIS技术将空间位置信息和属性信息无缝结合,结合数学统计和地理知识等,精确获取预警等级在行政区划所处方位、区域占比和防治措施等精细描述的预警信息,将图面内容转换为直观的文本描述,为最终预警产品的快速生成奠定基础(图1)。

    图  1  技术路线
    Figure  1.  Technical route

    为了实现目标,在开始研究前,要准备必要的空间图层作为计算的基础,即将具有不同预警等级属性的矢量栅格地质灾害预警区划图层和反映行政区名称及空间范围的研究区行政区划图层作为数据源,二者要求具有相同的地图投影参数。其中,前者一般通过网格剖分后通过多因子要素叠加进行综合评价得到,这一过程在当前地质灾害空间评价预警研究中是普遍采用的方法[5-6],但因不是本文的重点,故不再赘述。

    很显然,前述空间图层叠加结果中,不同预警级别对应的地区预警信息处于离散状态,且一般情况下,这些信息数量较为可观,为了后续分析数据高效便捷,有规律可循,需要提前分析预警等级、地区、方位等数据间的关联性,通过聚类分析建立分类簇集合信息,这一预警等级集合概念模型设计见图2

    图  2  预警等级集合概念模型示意图
    Figure  2.  Schematic diagram of the conceptual model of early warning level aggregation

    图2中可以看出,预警等级集合以三级预警等级为唯一的主键,包括了这一等级下的地区列表预警信息子集合,而该子集合以地区名称为主键,包括了一对多的分布方位和面积占比列表,该列表以分布方位为唯一索引,在实际计算时,相同分布方位需要进行去重处理,所占的面积要进行求和运算。最终通过综合计算,形成该预警等级下的总体分布区域、分布方位等综合预警描述信息。

    空间关系描述是GIS系统的基本功能之一,GIS的技术支持的地质灾害风险区划的最终目的是划分不同灾害等级的区域,可为地质灾害预警提供依据[7-8]。通过将预警区划图层和行政区划分区图层进行空间叠加,遍历判断每个预警等级矢量栅格单元和行政区划单元的空间拓扑关系,确定该行政单元是否包括某预警级别(图3)。如果二者拓扑关系为不相离,说明该地区具备该等级,反之如果是包含、相交、穿越等非相离关系,则还需进行两两拓扑求交运算,并重新采用交集中的预警等级区参与面积计算更具科学性和精确性。

    图  3  空间判别示意图
    Figure  3.  Schematic diagram of spatial discrimination

    一个由N个拐点(xi, yi)确定的封闭多边形的面积如式(1)计算:

    A=12i=0N1(xiyi+1xi+1yi) (1)

    式中:i——拐点序号;

    N——拐点个数;

    xi——第i个拐点x坐标;

    yi——第i个拐点y坐标;

    A——封闭多边形的面积。

    相同方位的预警等级所占行政区划单元的面积占比(R)是对预警等级广泛程度的描述,算式如下:

    R=i=0N1AiSi (2)

    式中:i——拐点序号;

    N——拐点个数;

    R——某预警等级的面积占比;

    Ai——某预警等级单元的面积;

    Si——某行政区划单元面积。

    对分布范围则用绝大部分、大部分、局部、个别四级进行空间范围的广泛程度描述,判别指标为集合中单元格之和的占比(R),定义见表1

    表  1  预警等级占比描述表
    Table  1.  Description of the proportion of early warning levels
    某预警等级在行政单元中的占比程度描述
    R≥0.9基本全域
    0.7≤R<0.9绝大部分
    0.5≤R<0.7大部分
    0.2≤R<0.5局部
    R≤0.2个别
    下载: 导出CSV 
    | 显示表格

    空间方向的定性描述是用若干主方向粗略地描述空间方向。而定量描述则是用方位角来量测空间目标之间的方向关系,因此方位角是空间方向描述的一个重要手段[9]。要获取一个预警等级单元在行政区划中的分布方向,实际上是通过计算该等级相对于所处区域的方位角得到(图4),结合地理知识,根据实际情况共划分出了8个方位角区间和对应的分布关系(图5)。

    图  4  方位角示意图
    Figure  4.  Schematic diagram of the azimuth angle
    图  5  方位角与分布方向的映射关系图
    Figure  5.  The mapping relationship between azimuth and distribution direction

    1、方位角计算

    (1)形心求算

    本次方位角计算要获取行政区划单元和预警等级单元的几何形心。一个由N个拐点(xi, yi)确定的封闭多边形的中心如式(3)、式(4)计算:

    Cx=16Ai=0N1(xi+xi+1)(xiyi+1xi+1yi) (3)
    Cy=16Ai=0N1(yi+yi+1)(xiyi+1xi+1yi) (4)

    式中:i——拐点序号;

    N——拐点个数;

    A——多边形的面积,由式(1)得出;

    xi——第i个拐点x坐标;

    yi——第i个拐点y坐标;

    Cx——几何形心x坐标;

    Cy——几何形心y坐标。

    (2)方位角求算

    当在平面上2个点的坐标已知时,给出方位角(十进制度数)公式如式(5):

    α=180πarctan2(X2X1,Y2Y1) (5)

    式中:X2——平面上终点x坐标;

    X1——平面上起点x坐标;

    Y2——平面上终点y坐标;

    Y1——平面上起点y坐标;

    α——方位角。

    (3)分布方向描述

    当方位角得出后,便可以根据方位角与分布方向的映射关系图(图5)获取具体的分布方向。

    经过前述步骤后,将得到一个以预警等级为索引的数据集合,其中包括各个地区的预警信息子集合。将地区预警信息子集合按照预警方位的个数进行升序排列,为了简洁,可以选取前若干个方位作为主要的方位,其余的则以“等地区”代替,而该级别的防治措施则从表2中对照获取[2]

    表  2  预警等级防治措施描述表
    Table  2.  Early warning level control measures description table
    预警等级防治措施
    1请严密防范
    2请加强防范
    3请注意防范
    4请监测分析
    下载: 导出CSV 
    | 显示表格

    通过遍历预警等级集合,形成该地区综合的预警产品描述信息。最终某地区预警信息描述格式形成通用模板举例如下:风险大(Ⅱ级)主要分布于××地区的西北部、北部、东北等局部地区,请加强防范;风险较大(Ⅲ级)主要分布于××地区的西南、南部等个别地区,请注意防范。

    考虑到预警信息的范围划分取决于行政区划图,不同的划分将得到不同的预警产品描述信息,故可将行政区划图看作一个相对的变量,为了验证工作的灵活性和可扩展性,本次通过GIS软件形成了包括11个虚拟乡镇的地区作为预警范围实验数据进行工作(图6),该地区某日的矢量栅格地质灾害预警区划图层业已形成(图7),采用3 km×3 km为预警等级单元,预警单元总数为5194个。上述文件均为Shapefile格式,其投影参数为CGCS2000。

    图  6  某地区地质灾害风险预警范围
    Figure  6.  Geological hazard risk warning range in a region
    图  7  某日地质灾害预警区划图层
    Figure  7.  Geological disaster early warning zoning map on a certain day

    GDAL库(Geospatial Data Abstraction Library)是一个开源的用于栅格和矢量地理空间数据格式的C++转换器库,由开源地理空间基金会在MIT(麻省理工学院)风格的开源许可下发布,目前几乎所有的GIS和RS软件底层都使用GDAL来读写空间数据[10]。由于GDAL库能够很好地支持包括Shp格式在内的很多数据格式,因此本次工作基于GDAL和C++语言,在Visual Studio 2019平台下实现了整个研发验证过程。研发过程中,以本算法为理论基础,结合UML理论建立了空间分析类、因子获取类、信息处理类,并形成了类的方法、属性和事件(图8),最后编译形成了预警信息生成软件工具。

    图  8  模块构成图
    Figure  8.  Module composition diagram

    利用已有的实际数据图层和软件工具,耗时不足1分钟形成了预警产品信息(表3)。从表中可以看出,该日预警级别包括了Ⅱ级、Ⅲ级,主要是以风险较大(Ⅲ级)为主,分布于秦姜董镇、王家湾镇、白堆子镇等地,并且有方位和分布范围程度描述。依据表中信息,通过组合可形成规范格式的综合预警产品描述信息,见图9

    表  3  某地区预警产品描述信息表
    Table  3.  Table of product description for a specific region
    级别分布地区方位占比/%防治措施
    风险大(Ⅱ级)王家湾镇西北等个别地区8.40加强防范
    五里埔乡东北、东部等个别地区4.40
    秦姜董镇南部、东南、东北等个别地区1.70
    风险较大(Ⅲ级)秦姜董镇基本全域96.60注意防范
    孙各庄乡东北、西北、东南等31.88
    王家湾镇基本全域91.00
    别山镇东北、西南、西部等局部地区33.00
    白堆子镇基本全域96.60
    龚家庄镇西北、东南、东北等局部地区24.60
    雁儿湾镇绝大部分74.60
    邱家庄乡西南、东南、西部等局部地区34.30
    西岗镇大部分地区58.20
    五里埔乡东南、东部、东北等局部地区32.00
    上营镇东北、西南、西北等30.40
    下载: 导出CSV 
    | 显示表格
    图  9  预警产品描述信息
    Figure  9.  The description information of the warning product

    为了进一步验证算法的正确性和优越性,将风险预警范围图(图6)和预警区划图层(图7)进行叠加,通过逐一人工判断每个乡镇区域预警区划情况,可看出所有的乡镇的结果和预警区划图层(图7)的空间展现结果完全吻合,如王家湾镇预警产品信息为风险大(Ⅱ级)、分布于西北,占比不大(图10),这一目估结果与采用本算法生成的预警产品描述信息一致,其次应用计算过程耗时短暂,可见算法能完全满足汛期地质灾害风险预警之快速高效的需要,与传统人工定性分析相比,具有规范、快速、准确的特点,没有遗漏,其优越性显著。

    图  10  预警产品描述信息验证
    Figure  10.  Early warning product description information verification

    文章提出了基于GIS技术支撑下的地质灾害风险预警产品描述信息的自动化生成技术方法,并以随机的行政区划图和已有的地质灾害预警区划图为数据源,采用GDAL库,通过研发预警信息获取模块进行了实例验证。结果表明,采用本算法形成预警产品描述信息快速高效,结果科学规范、完整全面,进一步提高了预警信息精度,不仅大大降低了预警产品信息获取的繁琐程度,而且节省了时间成本,又兼顾了产品的准确性,可为相关信息系统功能研发及专家研判提供基础理论支撑,将显著提高预警精细化程度和工作效率,适合在实际工作推广使用。

  • 图  1   技术路线

    Figure  1.   Technical route

    图  2   预警等级集合概念模型示意图

    Figure  2.   Schematic diagram of the conceptual model of early warning level aggregation

    图  3   空间判别示意图

    Figure  3.   Schematic diagram of spatial discrimination

    图  4   方位角示意图

    Figure  4.   Schematic diagram of the azimuth angle

    图  5   方位角与分布方向的映射关系图

    Figure  5.   The mapping relationship between azimuth and distribution direction

    图  6   某地区地质灾害风险预警范围

    Figure  6.   Geological hazard risk warning range in a region

    图  7   某日地质灾害预警区划图层

    Figure  7.   Geological disaster early warning zoning map on a certain day

    图  8   模块构成图

    Figure  8.   Module composition diagram

    图  9   预警产品描述信息

    Figure  9.   The description information of the warning product

    图  10   预警产品描述信息验证

    Figure  10.   Early warning product description information verification

    表  1   预警等级占比描述表

    Table  1   Description of the proportion of early warning levels

    某预警等级在行政单元中的占比程度描述
    R≥0.9基本全域
    0.7≤R<0.9绝大部分
    0.5≤R<0.7大部分
    0.2≤R<0.5局部
    R≤0.2个别
    下载: 导出CSV

    表  2   预警等级防治措施描述表

    Table  2   Early warning level control measures description table

    预警等级防治措施
    1请严密防范
    2请加强防范
    3请注意防范
    4请监测分析
    下载: 导出CSV

    表  3   某地区预警产品描述信息表

    Table  3   Table of product description for a specific region

    级别分布地区方位占比/%防治措施
    风险大(Ⅱ级)王家湾镇西北等个别地区8.40加强防范
    五里埔乡东北、东部等个别地区4.40
    秦姜董镇南部、东南、东北等个别地区1.70
    风险较大(Ⅲ级)秦姜董镇基本全域96.60注意防范
    孙各庄乡东北、西北、东南等31.88
    王家湾镇基本全域91.00
    别山镇东北、西南、西部等局部地区33.00
    白堆子镇基本全域96.60
    龚家庄镇西北、东南、东北等局部地区24.60
    雁儿湾镇绝大部分74.60
    邱家庄乡西南、东南、西部等局部地区34.30
    西岗镇大部分地区58.20
    五里埔乡东南、东部、东北等局部地区32.00
    上营镇东北、西南、西北等30.40
    下载: 导出CSV
  • [1] 刘传正,刘艳辉,温铭生,等. 中国地质灾害气象预警实践:2003—2012[J]. 中国地质灾害与防治学报,2015,26(1):1 − 8. [LIU Chuanzheng,LIU Yanhui,WEN Mingsheng,et al. Early warning for regional geo-hazards during 2003-2012,China[J]. The Chinese Journal of Geological Hazard and Control,2015,26(1):1 − 8. (in Chinese with English abstract)
    [2] 中国地质灾害防治工程行业协会. 地质灾害区域气象风险预警标准(试行): T/CAGHP 039-2018[S]. 武汉: 中国地质大学出版社, 2018

    China Geological Disaster Prevention Engineering Association. Regional meteorological risk early warning standard for geological disaster (trial): T/CAGHP 039-2018[S]. Wuhan: ChinaUniversity of Geosciences Press, 2018. (in Chinese)

    [3] 汪民. 关于地质灾害防治需要关注的几个问题[J]. 中国地质灾害与防治学报,2022,33(1):1 − 5. [WANG Min. Several problems needing attention in the prevention and control of geological disasters[J]. The Chinese Journal of Geological Hazard and Control,2022,33(1):1 − 5. (in Chinese) DOI: 10.16031/j.cnki.issn.1003-8035.2022.01-01
    [4] 唐军峰, 唐雪梅, 周基, 等. 滑坡堆积体变形失稳机制—以贵州剑河县东岭信滑坡为例[J]. 吉林大学学报(地球科学版),2022,52(2):503 − 516. [TANG Junfeng, TANG Xuemei, ZHOU Ji, et al. Deformation and instability mechanism of landslide accumulation:A case study of donglingxin landslide accumulation in Jianhe County, Guizhou Province[J]. Journal of Jilin University (Earth Science Edition),2022,52(2):503 − 516. (in Chinese with English abstract)
    [5] 金亚兵, 傅旭东. 滑坡泥石流灾害风险评价研究及其应用[M]. 北京: 地质出版社, 2018.

    JIN Yabing, FU Xudong. Landslide and debris flow disaster risk assessment research and applications[M]. Beijing: Geological Publishing House, 2018. (in Chinese)

    [6] 曹洪洋, 袁颖, 贾磊. 区域降雨型滑坡灾害预警预报[M]. 北京: 地质出版社, 2017

    CAO Hongyang, YUAN Ying, JIA Lei. Early warning and forecast of regional rainfall landslide disaster[M]. Beijing: Geological Publishing House, 2017. (in Chinese)

    [7] 张以晨, 张继权, 张峰. 极端降雨诱发地质灾害风险评价、预警及管理对策研究以吉林省东南部山区为例[M]. 北京: 科学出版社, 2017

    ZHANG Yichen, ZHANG Jiquan, ZHANG Feng. Study on risk assessment, early warning and management of geological disasters caused by extreme rainfall[M]. Beijing: Science Press, 2017. (in Chinese)

    [8]

    TAN Q L,HUANG Y,HU J,et al. Application of artificial neural network model based on GIS in geological hazard zoning[J]. Neural Computing and Applications,2021,33(2):591 − 602. DOI: 10.1007/s00521-020-04987-4

    [9] 王晓芳,侯恩科,张志华. 坐标方位角算法在GIS中的应用[J]. 测绘科学,2011,36(2):239 − 240. [WANG Xiaofang,HOU Enke,ZHANG Zhihua. Application of azimuth angle algorithm in GIS[J]. Science of Surveying and Mapping,2011,36(2):239 − 240. (in Chinese with English abstract) DOI: 10.16251/j.cnki.1009-2307.2011.02.037
    [10] 李民录. GDAL源码剖析与开发指南[M]. 北京: 人民邮电出版社, 2014

    LI Minlu. GDAL source code analysis and development Guide[M]. Beijing: Posts & Telecom Press, 2014. (in Chinese)

  • 期刊类型引用(9)

    1. 董颐. 一种面向云原生GIS平台的算力网络计算任务设计. 自动化与仪表. 2025(05): 144-148 . 百度学术
    2. 周诗凯,刘正华,余丰华,朱浩濛,黄丽,佘恬钰. 浙江省地质灾害气象风险预警一体化建设的探索与实践. 中国地质灾害与防治学报. 2024(02): 21-29 . 本站查看
    3. 杨铭珂,张永富,邓朋朋,陈忠超,李娟,袁方芳,张波. 贵州省建设用地利用情况起底大调查系统的建设与应用. 测绘与空间地理信息. 2024(07): 105-107+114 . 百度学术
    4. 凌成. 基于三维GIS空间信息地质灾害预警模型的研究. 矿产勘查. 2024(07): 1254-1261 . 百度学术
    5. 张东宇. 基于三维GIS的省域气象防灾减灾系统设计与开发. 现代信息科技. 2024(17): 129-133 . 百度学术
    6. 曹满,王洋,陈曦,魏运芳,王传涛. 基于时空分布特征的地质灾害气象预警预报模型研究. 计算机测量与控制. 2024(11): 278-286 . 百度学术
    7. 隋嘉,孙皓,张丽华,辛倩男,孙永旺,王栋. 青海省地质灾害监测预警信息化平台建设与实现. 中国地质灾害与防治学报. 2023(02): 92-101 . 本站查看
    8. 唐书亮. 基于数字孪生的炼化企业风险智能预警系统设计. 信息与电脑(理论版). 2023(06): 143-145 . 百度学术
    9. 马娟,张鸣之,齐干,叶思卿,黄喆,丁帆. 地质灾害监测复杂场景下压电式雨量计精度标定及适宜性分析. 中国地质灾害与防治学报. 2023(05): 91-96 . 本站查看

    其他类型引用(1)

图(10)  /  表(3)
计量
  • 文章访问数:  334
  • HTML全文浏览量:  217
  • PDF下载量:  303
  • 被引次数: 10
出版历程
  • 收稿日期:  2022-05-23
  • 修回日期:  2022-07-03
  • 网络出版日期:  2022-09-26
  • 刊出日期:  2022-10-19

目录

/

返回文章
返回