留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测

李梦华 张路 董杰 蔡杰华 廖明生

李梦华, 张路, 董杰, 蔡杰华, 廖明生. 四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测[J]. 武汉大学学报 ● 信息科学版, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
引用本文: 李梦华, 张路, 董杰, 蔡杰华, 廖明生. 四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测[J]. 武汉大学学报 ● 信息科学版, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
LI Menghua, ZHANG Lu, DONG Jie, CAI Jiehua, LIAO Mingsheng. Detection and Monitoring of Potential Landslides Along Minjiang River Valley in Maoxian County, Sichuan Using Radar Remote Sensing[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
Citation: LI Menghua, ZHANG Lu, DONG Jie, CAI Jiehua, LIAO Mingsheng. Detection and Monitoring of Potential Landslides Along Minjiang River Valley in Maoxian County, Sichuan Using Radar Remote Sensing[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367

四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测

doi: 10.13203/j.whugis20210367
基金项目: 

国家自然科学基金 41774006

国家自然科学基金 41904001

四川省科技厅重点研发计划 2019YFS0074

测绘遥感信息工程国家重点实验室开放基金 21R02

详细信息

Detection and Monitoring of Potential Landslides Along Minjiang River Valley in Maoxian County, Sichuan Using Radar Remote Sensing

Funds: 

The National Natural Science Foundation of China 41774006

The National Natural Science Foundation of China 41904001

the Provincial Key Research and Development Program of Sichuan Ministry of Science and Technology 2019YFS0074

the Open Research Fund of State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing 21R02

More Information
    Author Bio:

    LI Menghua, PhD, specializes in geohazard monitoring using radar remote sensing. E-mail: menghuali@kust.edu.cn

    Corresponding author: ZHANG Lu, PhD, professor. E-mail: luzhang@whu.edu.cn
图(14) / 表(1)
计量
  • 文章访问数:  382
  • HTML全文浏览量:  169
  • PDF下载量:  108
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-07-01
  • 刊出日期:  2021-10-05

四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测

doi: 10.13203/j.whugis20210367
    基金项目:

    国家自然科学基金 41774006

    国家自然科学基金 41904001

    四川省科技厅重点研发计划 2019YFS0074

    测绘遥感信息工程国家重点实验室开放基金 21R02

    作者简介:

    李梦华, 博士, 主要研究方向为雷达遥感地质灾害监测。menghuali@kust.edu.cn

    通讯作者: 张路,博士,教授。luzhang@whu.edu.cn
  • 中图分类号: P237

摘要: 长期以来,中国四川省茂县地区受地质、地形条件和构造活动的影响,滑坡等地质灾害频发,给人民的生命财产和公路等基础设施安全带来了巨大的威胁,因此需要对滑坡隐患区域进行有效识别和监测。以时序哨兵1号A、B卫星(Sentinel-1A/1B)影像为数据源,利用时间序列合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术对茂县岷江河谷区段的潜在滑坡隐患开展识别监测,对重点区段进行了分析,同时分析探讨了InSAR滑坡监测中不同轨道数据的视线方向形变测量灵敏度差异。从实验结果中成功探测识别出了茂县岷江河谷沟口乡至石大关乡段的20余处滑坡隐患,现场实地考察验证了识别结果的准确性,证明了时序InSAR方法在高山河谷区滑坡隐患识别监测中的有效性。

English Abstract

李梦华, 张路, 董杰, 蔡杰华, 廖明生. 四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测[J]. 武汉大学学报 ● 信息科学版, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
引用本文: 李梦华, 张路, 董杰, 蔡杰华, 廖明生. 四川茂县岷江河谷区段滑坡隐患雷达遥感识别与形变监测[J]. 武汉大学学报 ● 信息科学版, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
LI Menghua, ZHANG Lu, DONG Jie, CAI Jiehua, LIAO Mingsheng. Detection and Monitoring of Potential Landslides Along Minjiang River Valley in Maoxian County, Sichuan Using Radar Remote Sensing[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
Citation: LI Menghua, ZHANG Lu, DONG Jie, CAI Jiehua, LIAO Mingsheng. Detection and Monitoring of Potential Landslides Along Minjiang River Valley in Maoxian County, Sichuan Using Radar Remote Sensing[J]. Geomatics and Information Science of Wuhan University, 2021, 46(10): 1529-1537. doi: 10.13203/j.whugis20210367
  • 2017-06-24发生在中国四川省阿坝州茂县叠溪镇的特大高位山体滑坡瞬间摧毁了坡脚的新磨村,掩埋了64户农房和近2 km的松坪沟河道,造成10人死亡和73人失踪,引起了国内外对该地区地质灾害问题的广泛关注[1-3]。茂县地处青藏高原东缘,构造运动强烈,地震频发,1933年叠溪地震、1976年松潘-平武地震和2008年汶川地震使得山体震裂松动,岩体破碎,为滑坡灾害形成提供了大量的物质基础。岷江自北向南穿过该地区形成高山峡谷地貌,河流的侵蚀加上陡峭的山坡为滑坡的发育提供了有利的地形条件。此外,短时强降雨也会进一步诱发滑坡体失稳成灾[1]。受这些因素的共同影响,长期以来茂县地区滑坡灾害频发,对当地居民人身安全以及沿岷江修建的213国道的运行安全构成了严重威胁。

    利用现代测绘遥感技术手段对潜在滑坡隐患进行早期识别、圈定高风险区域以便有针对性地制定科学防范措施,可有效减少滑坡灾害造成的人员伤亡和财产损失[4]。滑坡体表面形变是判断滑坡稳定性最重要、最直观的表观特征[5],传统的裂缝计、全球导航卫星系统(global navigation satellite system,GNSS)等地面测量手段需要耗费大量的人力、物力,在复杂山区很难开展大范围监测。星载合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)技术具有全天时、全天候、大范围地表形变监测能力[6],为滑坡等地质灾害隐患的广域识别与长时间变形监测提供了一种新的技术手段,在滑坡研究中得到了广泛应用[7-12]。时间序列InSAR技术经过十几年的发展,基本算法已较为成熟[13-16]。但传统InSAR技术所能测量到的形变是真实地表三维形变在雷达视线方向(line of sight,LOS)上的投影。在坡度坡向各异的复杂高山峡谷地区,InSAR方法对于不同坡体的形变测量敏感程度差异很大。因此,基于InSAR处理结果进行区域滑坡灾害隐患识别时,不能仅根据单一LOS形变速率来简单地判断评定不同坡体变形的潜在危险程度。

    本文以茂县岷江河谷区段为例,收集了46景升轨、42景降轨哨兵1号A、B卫星(Sentinel-1A/1B)数据,采用时序InSAR技术反演得到了岷江上游茂县沟口乡至石大关乡一线的区域地表形变场,同时基于卫星成像几何参数,结合地形信息估算了不同轨道InSAR测量的灵敏度,综合考虑升降轨InSAR测量灵敏度和视线向形变速率,对该区域内潜在的滑坡隐患进行了识别解译。

    • 时序InSAR技术通过提取观测时间段内具有稳定散射特征的高相干点集合,并依据不同相位信号分量的时空分布特性来实现误差分离,从而获得经过地形改正以及轨道、大气误差去除的形变速率场,代表性技术主要有永久散射体(persistent scatterer interferometry,PSI)技术[13]和小基线集(small baseline subsets,SBAS)技术[14],具体细节可参见文献[13-1416-17]。

    • 本文采用的时序InSAR处理流程[18]综合了PSI和SBAS两种方法,保留了PSI技术处理的全分辨率和高精度的特点,以及SBAS通过构成短基线集来获得高相干干涉对的处理方法,同时不预设形变模型来进行参数估计,这样能获得更多相干点,从而更可靠地提取地表形变量。该处理流程如图 1所示,从缠绕的差分干涉图出发,无需多视处理,利用距离向频谱滤波找出在短时间干涉像对中失相干很小的点,即滤波后相位微失相干(slowly-decorrelating filtered phase,SDFP)点。利用低通滤波器从干涉相位观测值中分离出空间相关的相位分量,包括地表形变相位、大气湍流混合相位、部分轨道误差以及空间相关的地形误差;由视角引起的空间不相关的相位,如空间不相关的地形误差和相位中心误差,可以用残余相位与垂直基线的线性关系拟合估计去除。去除掉这两部分相位之后得到像素的去相关噪声,根据噪声水平的大小设定阈值筛选出SDFP点。对SDFP上的缠绕相位进行时空三维解缠,恢复相位绝对值。解缠后,利用高通滤波器滤除大气误差,再减去空间相关和空间不相关的误差相位,就可以解算出最终的形变量。

      图  1  滑坡隐患识别监测时序InSAR处理流程

      Figure 1.  Workflow of Identification and Monitoring Potential Landslide Using Time-Series InSAR

    • InSAR测量得到的形变量是地表三维形变矢量(包含南北方向分量dnorth、东西方向分量deast和垂直方向分量dup)在雷达LOS上的投影[19]。以右视成像为例,视线方向形变dLOS可以表示为:

      dLOS=dnorthsinθsinαh-deastsinθcosαh+dupcosθ ]]>

      式中,αhθ分别代表飞行指向角和名义入射角。

      以卫星视角30°、飞行方向角345°的情况为例,如图 2(a)所示,假设InSAR测得的LOS形变为30 mm(缓慢变形滑坡体常见的年平均形变量),使用3 mm作为形变测量的标准差,理论上存在着无数种不同方向的三维形变矢量LOS投影量等于30 mm。3 mm的LOS测量精度对应于垂直方向的精度约为3.5 mm,对应于零多普勒面内的水平形变测量精度为6 mm。零多普勒面方向的水平形变可以看做是其他方向形变在该方向上的投影。零多普勒方向的测量精度最高,为6 mm,换算到其他各个方向的精度如图 2(b)所示。同时可以看出在此成像几何情况下,3 mm的LOS形变测量精度转换到南北方向精度降低约8倍,可以达到23.2 mm。传统InSAR技术不具有沿卫星飞行方向形变的测量能力,所以在卫星飞行方向σ=。值得注意的是,3 mm的LOS形变测量精度在城市地区通常可以达到,但是在复杂山区滑坡监测中,受相干性、大气延迟等因素影响,很难达到。滑坡研究中一般比较合理的精度要求在10 mm左右[20],此时对应的南北方向水平形变测量精度接近80 mm。

      图  2  InSAR LOS观测对滑坡三维运动灵敏度分析示意图

      Figure 2.  Diagram of the Sensitivity of InSAR Line-of-Sight Observation to the Real Motion of Landslides

      在山区滑坡监测中,不同坡体的形变方向差异很大。假设滑坡体沿着最大地形梯度方向(坡度φ、方向角αasp)发生位移,考虑到InSAR测量结果是滑坡体真实三维形变LOS上的投影,结合某一轨道的卫星成像几何参数(入射角θ和飞行方向角α)和区域地形,可以得到该轨道下InSAR形变测量结果dLOS与该区域真实形变dslope之间的定量化关系:

      dLOS=cosθ×dvertical-sinθcosα×deast+sinθsinα×dnorth ]]>
      dvertical=dslope×sinφ ]]>
      deast=-dslope×cosφsinαasp ]]>
      dnorth=-dslope×cosφcosαasp ]]>

      这样就可以估算得到不同坡体上InSAR测量的灵敏度,直观地判断不同坡体形变的真实量级和危险程度,更可靠地捕捉区域内的形变信号,实现对灾害隐患区域的准确探测。

    • 本文选择四川茂县岷江河谷典型区域作为研究区域,从沟口乡开始,沿岷江向上延伸至石大关乡,区域范围以及地形如图 3所示。茂县位于四川省西北部的阿坝藏族羌族自治州,地处青藏高原向四川盆地过渡地区,地貌以高山峡谷为主,地势从西北向东南倾斜。发源于岷山南麓的岷江由北向南贯穿四川省全境,途经松潘、茂县等地,从都江堰流入成都平原,最终在宜宾汇入长江,是长江的重要支流之一。岷江河谷两岸地势陡峭,大多呈V形。研究区域气候类型属于高原季风气候,气象条件受海拔影响显著,暴雨、泥石流等灾害集中多发于春夏季节。

      图  3  研究区域地形图

      Figure 3.  Topographic Map of the Study Area

      茂县地区地质构造复杂,分布有岷江断裂带、松坪沟断层及石大关断层,复杂的地质构造使得区域内岩石高度破碎,地层出露以砂岩、板岩、千枚岩等软质岩体为主。研究区处于南北地震带中段的龙门山地震带和松潘地震带上,地震活动频繁。本文使用了两组Sentinel-1A/1B时间序列数据集,分别为升轨和降轨观测获取,其基本参数如表 1所示,各自的时空基线分布如图 4所示。图 4中每个圆圈代表一景Sentinel-1A/1B影像,每条黑线连接的两景SAR影像组成一个干涉像对进行差分干涉处理。

      表 1  实验采用的SAR数据的基本参数

      Table 1.  Basic Parameters of SAR Datasets

      卫星 时间 方向 入射角/(°) 飞行方向角/(°) 数量/景
      Sentinel-1A/1B 2014-10-14—2017-07-06 升轨 37.35 347.23 46
      2014-10-09—2017-06-07 降轨 35.91 192.73 42

      图  4  Sentinel-1数据干涉对组合图

      Figure 4.  Interferometric Combinations for Sentinel-1 Data

    • 采用时序InSAR中的StaMPS方法对覆盖茂县岷江河谷区段的升轨和降轨Sentinel-1数据进行了处理,提取了2014—2017年该地区LOS的地表年平均形变速率场,叠加在谷歌地球

      (Google Earth)影像上(见图 5(a)图 6(a))。图 5图 6中红色表示沿着LOS远离卫星运动,蓝色表示地物朝向卫星运动。探测到的相干点主要集中在坡体中下部,靠近坡顶的区域受陡峭地形和植被覆盖的影响相干点较为稀疏。从形变的空间分布来看,升轨数据集探测到的形变区域较多,大多集中在岷江河谷的西岸山坡上,最大LOS形变速率超过50 mm/a。降轨数据集中只有擦耳岩村的区域形变较为明显。对比同一区域不同轨道观测获取的形变场,不难看出不同观测几何的数据在同一山区存在明显的灵敏度差异。采用前述的滑坡形变监测灵敏度估算方法,结合参考地形数据,计算得到该区域InSAR LOS形变测量灵敏度分布(见图 5(b)图 6(b))。升轨数据处理结果中箩兜寨到两河口一线密集分布的形变区在降轨数据集中几乎探测不到,这说明降轨观测数据对于该区域的形变测量灵敏度几乎为零。结合探测灵敏度分布图可以更准确地解译InSAR结果,判断高风险形变区域。

      图  5  Sentinel-1升轨数据集时序InSAR形变测量结果

      Figure 5.  Time-Series InSAR Results of Ascending Sentinel-1 Data

      图  6  降轨数据时序InSAR形变探测结果

      Figure 6.  Time-Series InSAR Results of Descending Data

    • 在区域形变探测结果和灵敏度分析的基础上,将该区域划分为沟口乡-飞虹乡、水草坪村-两河口、两河口3大块区域分别开展分析。

    • 茂县沟口乡位于两条峡谷的交汇处,图 7给出了从Sentinel-1升轨数据集中提取的沟口乡-飞虹乡区域年平均形变速率图及其对应区域的谷歌地球光学影像。飞虹乡擦耳岩村所在的坡体、鸡公寨所在山体以及沟口乡政府对面坡体这3块区域可以探测到明显的形变信号。其中鸡公寨滑坡形变量级最大,LOS最大形变速率接近50 mm/a,而沟口乡对面坡体的最大LOS形变速率仅为25 mm/a左右。结合图 5的形变测量灵敏度图可以发现,鸡公寨滑坡在升轨数据集中的形变测量灵敏度在0.8以上,升轨数据集的观测几何非常有利。而茂县沟口乡政府对面的坡体朝向偏南北向,在图 5(b)中该区域的灵敏度仅为0.4左右,实际沿坡面方向的形变速率远比InSAR结果的LOS形变速率大。

      图  7  沟口乡-飞虹乡区域Sentinel-1降轨数据形变情况

      Figure 7.  Deformation Map from Ascending Sentinel-1 Data in the Region from Goukou Township to Feihong Township

      飞虹乡政府对面擦耳岩村所在的斜坡也是呈近南北朝向,图 7(a)的升轨数据集处理结果中LOS形变速率为10 mm/a左右,但是其对应的形变测量灵敏度仅为0.2左右。图 8给出了擦耳岩滑坡区域降轨数据集的形变结果,在降轨数据集中的灵敏度约为0.4~0.5。降轨数据集LOS最大形变速率超过50 mm/a,对应实际沿坡面运动速率超过100 mm/a。图 8中同时给出了该斜坡上一个点对应的降轨时间序列形变量,可以看出该区域时间序列形变呈现周期性的变化,2015年和2017年都出现过形变放缓的趋势。

      图  8  飞虹乡擦耳岩滑坡Sentinel-1降轨数据形变图

      Figure 8.  Deformation Map from Descending Sentinel-1Data of Ca'eryan Landslide in Feihong Township

      图 9(a)是高分辨率光学卫星影像中擦耳岩滑坡的大致形态,图 9(b)是从茂县飞虹乡政府所在位置拍摄的擦耳岩滑坡体的现场实景照片。从图 9中可以看出整面山坡呈现圈椅状的典型滑坡地貌特征,斜坡上部坡度较陡,下部较缓,坡体上存在大量小型崩塌、落石留下的痕迹。滑坡体后缘可见滑动后形成的陡壁,前缘受岷江河水长期冲刷,坡度较陡,形成临空面,存在小型滑塌堆积挤压河道的现象。

      图  9  飞虹乡擦耳岩村斜坡

      Figure 9.  The Slope in Ca'eryan Village, Feihong Township

    • 茂县水草坪村到两河口之间的河谷西侧山体存在连续几个大型的不稳定斜坡。这一段岷江河道呈西北-东南走向,切割山体形成的坡体朝向与升轨Sentinel-1数据的雷达视线方向近似平行,成像几何有利于观测,探测灵敏度均在0.8以上,形变速率图如图 10所示。其中位于水草坪村和箩兜寨的两个滑坡在雷达LOS形变量最大,均超过50 mm/a,并且滑坡边界在形变速率图中较为清晰,该区域岷江东岸的斜坡无明显的整体变形信号。

      图  10  水草坪村-两河口区域哨兵1号升轨数据集形变速率图

      Figure 10.  Deformation Velocity from Descending Sentinel-1 Data in the Region from Shuicaoping to Lianghekou

      图 11是岷江西岸坡体形变探测结果的放大图,其中ABCD 4个区域分别取自4个不同的滑坡体且4块不稳定的区域是同一片山坡被沟壑分割而成的不同坡体。从2014—2017年Sentinel-1升轨数据集观测到的形变结果可以看出,B所在的滑坡体面积和LOS形变速率最大,滑坡后缘位于距离山顶约1/3处,前缘临江,长期受岷江河水冲刷,最大形变速率超过50 mm/a,出现在靠近坡脚的位置。A区域的滑坡体在面积和形变速率方面都比B区域滑坡小,前缘不临江。水草坪村和箩兜寨分别位于AB所在的滑坡体上,受到滑坡灾害的直接威胁。

      图  11  水草坪村和箩兜寨滑坡形变探测结果放大图

      Figure 11.  Enlarged View of the Deformation Map inShuicaoping Village and Luodou Village

      CD所在的两个滑坡体位于箩兜寨滑坡的上游,靠近两河口位置,其中C区域滑坡形变速率较大,超过50 mm/a,并且滑坡前缘临江。D区域的滑坡体相比于C区域形变范围和形变速率都更小。CD两个滑坡体区域附近均没有居民点,唯一的威胁对象是沿着岷江河谷伸展的213国道。需要注意的是,CD区域在升轨数据集中的形变探测灵敏度比AB区域低,在升轨LOS相同大小的形变速率情况下,CD区域的真实沿坡面形变速率更大。

      图 12给出了ABCD对应位置测量点集合的平均时间序列形变量。这4个区域2014-10—2017-07的累积形变量均超过100 mm,其中B区域累积形变超过150 mm,在2016-01前后形变放缓,随后加速,在2017-05之后又再次出现加速;AC区域大体上呈现线性形变特征;D区域则是在2015-03之后出现短暂的形变减缓。

      图  12  水草坪村和箩兜寨ABCD 4个区域时间序列形变

      Figure 12.  Time-Series Deformation of Regions A, B, C and D in Shuicaoping and Luodou Villages

      2019年11月初,对B区域所在的箩兜寨滑坡进行了实地考察,图 13中两张照片分别是在图 11白色三角形标注的P1P2点处拍摄的。P1点靠近箩兜寨滑坡的顶部,可以清晰地观察到滑坡后缘因坡体长期下滑形成的陡坎。在P2点处,山坡上通向箩兜寨山顶的水泥道路有多处受坡体局部下滑推挤影响出现明显的破损痕迹。

      图  13  箩兜寨滑坡现场考察照片

      Figure 13.  In-situ Photos of the Landslide in Luodou Village

    • 两河口区域位于多条峡谷交汇处,通往黑水县的302省道在黑水河汇入岷江处与213国道汇合。该区域的不稳定斜坡主要从升轨数据形变探测结果中识别得到,图 14给出了两河口区域的升轨数据LOS形变速率图。

      图  14  两河口区域Sentinel-1升轨数据集形变速率图

      Figure 14.  Deformation Velocity from Ascending Sentinel-1 Data in the Region of Lianghekou

      除了黑水河北岸白布村所在的南北向坡体外,岷江两岸山体的坡脚处分布着一些小型滑坡体。形变的范围和形变速率都相对较小,Sentinel-1升轨数据LOS形变速率大约在30 mm/a左右。白布村所在的山坡呈南北向,同时坡体表面沟壑多,局部地形十分复杂,形变探测的灵敏度差异也较大,所以在图 14的形变测量结果中表现为一部分区域朝向卫星运动,一部分区域远离卫星运动。此外,在时序InSAR结果中,白布村滑坡体的大片中心区域无法提取出有效的测量点。根据现场考察的情况看,该坡体上植被覆盖很少,推测这部分地区形变梯度很大,导致了严重的失相干现象。

    • 四川茂县岷江河谷地区受地质条件、地形地貌、地震活动、强降雨等多种因素的影响,发育分布有众多的滑坡灾害隐患,对当地居民的生命财产安全和国道省道等基础设施的安全运营构成了严重威胁。本文收集了46景升轨和42景降轨的Sentinel-1数据,利用时间序列InSAR技术在岷江河谷沟口乡到石大关乡一线区域开展了滑坡隐患识别,共识别出岷江及其支流黑水河沿岸大小滑坡隐患20余处。同时对不同区段的形变测量结果进行了详细分析,其中针对飞虹乡到箩兜寨一线的多个滑坡分别提取分析了时间序列形变,并实地开展了现场调查核实,验证了InSAR探测识别结果的有效性。

      本文研究结果表明,利用时序InSAR技术在坡度和坡向各异的复杂高山峡谷地区开展滑坡隐患识别和形变监测时,需要综合考虑SAR数据的成像几何条件与区域地形特征的共同影响分析不同轨道的InSAR测量灵敏度,避免造成对于滑坡体真实形变量级的误判。

参考文献 (20)

目录

    /

    返回文章
    返回