投稿问答最小化  关闭

万维书刊APP下载

永定河流域地表沉降监测及生态补水影响分析

2022/12/26 16:26:56  阅读:316 发布者: 来源:

永定河流域地表沉降监测及生态补水影响分析

田枭1,林敦灵2,刘洋1,韩森1,董梦2,李航昊1

1. 武汉大学 测绘学院,武汉 430079

2. 北京市勘察设计研究院有限公司,北京 100038

摘要:针对永定河流域(北京段)地表沉降及其与地下水位变化、生态补水之间的关系问题,利用小基线集合成孔径雷达干涉测量技术对覆盖研究区域的32景升轨雷达影像进行处理,提取了地表形变时空分布特征,分析了地表形变与地下水位变化的关系,并探讨了生态补水对区域地表沉降的影响。结果表明,流域上部具有一定的沉降区域,主要沉降速率为520 mm/a。地表形变与地下水位变化之间存在显著相关性:流域上部地下水位整体呈下降趋势,区域地表沉降相对严重;流域中下部地下水位整体呈上涨趋势,区域地表沉降得到遏制。生态补水对区域地表沉降有显著影响,但具体类型与区域地形地质、河道特征相关。

0 引言

永定河流域(北京段)位于华北平原北部,流域内城镇人口密集、工业发达、用水需求量巨大,在20世纪多次进入枯水期。该流域内地下水开采过量,地下水位长期下降,并出现了较为严重的地表沉降现象[1]。目前,已有学者对北京市地表沉降现象进行过研究,发现过量开采地下水导致地下水位长期下降是北京市地表沉降的主要因素[2-5]

21世纪以来,北京市采取系列政策限制地下水开采,并且实施了永定河河道的生态补水,由此,地表沉降现象较20世纪有了一定程度的好转[6-7]。生态补水是指补充生态系统被破坏所缺失的水,修复生态系统的结构、功能以及自我调节的能力,创造一个适宜人类生活的生态环境并使得生态环境和人类社会能良好健康发展,是一种恢复生态平衡、改善水环境的重要方法[8]。地表水和地下水之间存在密切的联系。生态补水不但能够显著影响地表水环境和生态环境,而且能够明显影响地下水状态与地表形变状态。

地表形变监测可以提供流域内地表形变的时空特征,为分析地表形变与地下水位变化之间的关系提供数据支持,从而为河道的补水工程提供参考信息[9]。永定河流域面积广阔、河道纵横,传统地表形变监测手段(如GPS、水准等)难以全面获得流域内地表形变的时空分布特征,且测量成本较高[10]20世纪90年代以来,合成孔径雷达干涉测量(interferometric synthetic aperture radarInSAR)技术飞速发展,该技术具有全天时、高空间分辨率、大范围等优点[11],逐渐成为大面积地表高程和地表形变监测的重要手段。小基线集合成孔径雷达干涉测量技术(small baseline subset interferometric synthetic aperture radarSBAS-InSAR)可以克服传统InSAR技术易受时空失相干、大气延迟等影响的难题,能够显著提高地表沉降监测的精度和可靠性[12],从而成为地表形变监测的主要手段之一。

目前,对北京市地表沉降的研究主要集中于北京市整体或者北京市局部区域,如朝阳区、通州区等[5,13-14],而有关永定河流域(北京段)地表沉降及其与地下水位变化、生态补水之间关系的研究较少。本文基于SBAS-InSAR技术,利用20181月—20216月期间的32景升轨Sentinel-1A合成孔径雷达(synthetic aperture radarSAR)影像提取并分析了永定河流域(北京段)地表形变的时空分布特征,结合收集到的永定河流域(北京段)10个地下水位监测点的监测数据,分析、探讨了地表形变与地下水位变化之间的关系以及生态补水对区域地表沉降的影响,能够为生态补水防治地表沉降等提供参考。

1 研究区域概况

北京市位于华北平原北部,西北部多为山地,东南部多为平原。永定河北端起始于阴山和太行山支脉恒山所包围的高原,出官厅水库后进入北京市,在北京市内自西北向东南流淌,流域覆盖北京市的门头沟区、石景山区、丰台区、房山区、大兴区等[15]。北京市城区人口密集且工业众多,大气降水难以满足用水需求。区域内长期过量开采地下水,地下水位长期下降,地表沉降严重[16]

本文的研究区域主要包括永定河流经的门头沟区、石景山区、丰台区、房山区、大兴区(图1),经纬度范围为(116.07°~ 116.35°E39.71°~ 40.00°N),区域中部和东南部多为平原,西北部多为山地。区域地下水按赋存条件可分为平原区第四系水资源、基岩水资源,以大气降雨为主要补给来源,以蒸发、开采等方式排泄;按地貌特征可分为山区、平原区两个水文地质区[2]

1 研究区域

2 数据及其处理

2.1 实验数据

本文选取欧洲空间局Sentinel-1A卫星在20181月—20216月期间获取、覆盖永定河流域(北京段)的32景升轨SAR影像进行地表形变监测,其空间覆盖范围见图1。欧洲空间局Sentinel-1A卫星发射于201443日,雷达波为C波段,波长为5.63 cm,包含宽幅干涉模式(interferometric wide swathIW)、条带模式(strip map modeSM)、超幅宽模式(extra wide swathEW)以及波模式(wave modeWM4种工作模式。本文选取宽幅干涉模式,空间分辨率为5 m,平均入射角为39°,极化方式为VV极化。此外,选取Sentinel-1A卫星的精密轨道数据用以估计和校正轨道相位,采用30 m分辨率的航天飞机雷达地形测绘使命(shuttle radar topography missionSRTM)数字高程模型(digital elevation modelDEM)数据估计并去除地形相位。

2.2 SBAS-InSAR技术及数据处理

SABS-InSAR技术最早由Berardino提出[17],该技术设置空间基线阈值,将SAR影像组合成多个子集合。其中,子集合内SAR影像间的基线距小,子集合之间的基线距大,利用奇异值分解法(singular value decompositionSVD)将多个子集合进行联合求解[14,17]。本文设置空间基线小于300 m,时间基线小于365 d,对32幅雷达影像进行组合与干涉处理,其时空基线如图2所示,共得到228个干涉对。

2 干涉图的时空基线

假设采集的SAR影像为N+1幅,选择其中一幅作为主影像,并设置相应的时空基线阈值,再进行差分干涉处理,得到M幅差分干涉图。

在数据处理中,通常需要消除地形误差相位、大气相位及噪声相位等以提高形变监测精度和可靠性。根据地形误差相位与地形误差的线性关系,利用最小二乘法估计并消除地形误差相位[20];根据不同误差在时空域中的差异性,利用时间域高通滤波、空间域低通滤波减免或消除大气相位[21-22];利用滤波算法和多视处理减免或消除噪声相位[12]。上述非形变相位消除后,基于干涉图中相干性较高的点建立方程组,利用SVD法对方程组进行求解,进而得到研究区域的地表形变[23]

3为相应的SBAS-InSAR数据处理流程。基于上述处理,可获得永定河流域(北京段)的地表形变(图4),正值、负值分别表示地表点靠近、远离雷达卫星(即地表抬升、地表沉降)。其中,在干涉相位图生成中,采用距离向、方位向多视比4:1进行多视处理,采用30 m分辨率SRTM DEM数据去除地形相位,采用Goldstein滤波方法进行自适应滤波[24],采用最小费用流(minimum cost flowMCF)方法进行相位解缠[25]

3 基于SBAS-InSAR技术的数据处理流程 [26]

3 地表形变的时空特征

3.1 地表形变的空间分布特征

从图4可以看出,20181月—20216月期间,永定河流域(北京段)整体地表形变比较稳定,大部分区域地表形变速率小于5 mm/a。根据形变速率图,可以将研究区域分为沉降区域、稳定区域以及抬升区域。流域上部具有一定的沉降区域,地表沉降现象相对严重,主要的沉降速率大小为520 mm/a,最高可达34.9 mm/a。同时,流域上部存在局部稳定区域(图4中左上角的局部黄色区域),主要的形变速率大小约0 mm/a,其可能与局部地形地质与周围区域存在显著差异相关,该局部区域较周围区域平坦,在4.2节分析中认为地形平坦区域生态补水对地下水位的影响更大,因此该区域地表沉降得到遏制。流域中下部整体上呈轻微的地表抬升,主要的抬升速率大小为04 mm/a,最高可达9.9 mm/a。此外,区域内存在小范围的局部沉降现象(如图4A4点),其可能与区域范围内人类活动[27]或影像噪声等相关。

4 研究区域的形变速率

3.2 地表形变的时序特征

本文提取了10个地下水位监测点(图4P1-P10100 m范围内形变的均值时间序列(图5)。其中,P2位于沉降区域,P7P8P9位于抬升区域,其它监测点位于稳定区域。从图中可以看出,P2点的沉降趋势最为明显,最大沉降达到7.1 mm,与形变速率图(图4)中相应区域的沉降现象相一致;P7P8P9点处的时序形变均呈抬升趋势,P9点的最大抬升达6.6 mm,该3个点均位于流域中下部,与形变速率图(图4)中相应区域的抬升现象相一致。P1位于局部稳定区域,时序形变呈波动变化,波动幅度较小,与形变速率图(图4)相符合。

本文同时提取了5个典型点(图4A1-A5100 m范围内形变的均值时间序列(图6)。从图中可以看出,A2点位于流域中下部的丰台区,时序形变呈抬升趋势,最大抬升达8.7 mm;其余4个点均呈沉降趋势,其中,A1A5分别位于门头沟区北部、南部沉降较显著的区域,最大沉降达30.4 mmA3位于海淀区,最大沉降达23.1 mmA4位于房山区局部沉降区域,最大沉降达15.8 mm。从图5、图6可以看出,在沉降、抬升比较明显的区域形变趋势未出现显著波动。

5 地下水位监测点的形变时间序列

6典型点的形变时间序列

4 讨论与分析

4.1 地表形变与地下水位变化的关系

为分析地表形变与地下水位变化的关系,本文将10个地下水位监测点的水位变化时间序列值(监测时间与SAR数据获取时间相近)的前后各10个数据取平均、作差,进而给出了地下水位变化等值线(图7)。结合图4、图7可以看出,地下水位变化等值线和地表形变的空间分布特征相近:流域上部地下水位整体上呈下降趋势,最大值约0.9 m,区域地表沉降现象相对严重,主要的沉降速率大小为520 mm/a;流域中下部地下水位整体上呈上涨趋势,最大值约5.3 m,区域地表沉降现象得到遏制,整体上呈轻微的地表抬升,主要的抬升速率大小为04 mm/a。此外,流域上部P1点区域地下水位下降趋势较P2点区域要小,地表形变也较P2点区域要小,主要的形变速率大小约0 mm/a

本文同时选取典型点(P2P9)地下水位变化时间序列与监测点100 m范围内形变均值时间序列进行对比(图8),以进一步分析二者的相关性。其中,P2点是水位监测点中地表沉降最严重、地下水位下降最大的点位,P9点是地表抬升最显著、地下水位上涨最大的点位。可以看出,P2点地下水位下降幅度不大但长期位于低值,最低水位为-3.7 m,形变时间序列长期为下降趋势,形变时间序列与水位变化的整体趋势相一致;P9点地下水位上涨明显,尤其是2019年、2020年出现两个显著的上涨高峰,最高水位为7.5 m,形变时间序列长期为抬升趋势,形变时间序列与水位变化的整体趋势相一致。综上可知,区域地表形变与地下水位变化具有显著相关性,区域地表沉降可能与地下水位呈下降趋势相关[28-29]

7 2018-2021年间地下水位变化等值线

aP2

bP9

8 P2P9点的形变和地下水位时间序列比较

4.2 生态补水对区域地表沉降的影响

由于过量开采地下水,北京市地区地下水位长期下降。近年来,北京市大力推行限制地下水开采的政策,并且实施了永定河河道生态补水方案。相关措施在一定程度上遏制了该地区地下水位长期下降的趋势并减缓了地表沉降现象的加重[30-31]

生态补水使得部分水通过河床渗透来补给地下水,渗透总水量主要受河流接纳水量、河床渗透性等影响[32]。河流上部、中下部的地形地质通常不同,河道的水文地质构造也不同,以及水利工程的分布,导致生态补水对不同区段地下水位变化的影响也不同。由图9可以看出,永定河流域上部多为山地,地形落差较大,地质坚实,水流湍急,并且三家店拦水闸到卢沟桥拦水闸的河道存在人工底部防渗处理,渗水性较弱,河道渗水对地下水的影响较小,实施生态补水后,其地下水位整体上仍呈下降趋势(图7);中下部多为平原,地形落差较小,地质松软,水流平缓,湖泊较多,部分河道面积较大,河道多为采砂石坑[33],渗水性较强,河道渗水对地下水的影响较大,实施生态补水后,其地下水位整体上呈上涨趋势(图7)。根据4.1小节,地下水位变化对地表形变状态影响较大,流域上部整体上地表沉降相对严重,流域中下部地表沉降得到遏制。该分布特征与区域地形地质、生态补水特征相一致。综上可知,生态补水对区域地表沉降有显著影响,但具体影响类型与区域地形地质、河道特征相关。

9 研究区域地形

5 结束语

本文利用SBAS-InSAR技术处理覆盖永定河流域(北京段)、20181月—20216月期间的32景升轨Sentinel-1A SAR影像获取了区域地表形变时空分布特征,结合10个地下水位监测点的水位变化数据,分析了地表形变与地下水位变化的相关性,进而探讨了生态补水对区域地表沉降的影响。

1)流域上部具有一定的沉降区域,主要的沉降速率大小为520 mm/a;流域中下部整体上呈轻微的地表抬升,主要的抬升速率大小为04 mm/a;区域内存在小范围的局部形变现象。

2)区域地表形变与地下水位变化具有显著相关性。地下水位变化等值线与地表形变具有相近的空间分布特征,区域地表沉降可能与地下水位呈下降趋势相关。

3)生态补水对区域地表沉降有显著影响,但具体影响类型与区域地形地质、河道特征相关。永定河流域上部多为山地,局部河道存在底部防渗处理,渗水性较弱,生态补水对地下水影响较小,其地下水位整体上呈下降趋势,区域地表沉降相对严重;永定河流域中下部多为平原,河道多为采砂石坑,渗水性较强,生态补水对地下水影响较大,地下水位整体上呈上涨趋势,区域地表沉降得到遏制。

参考文献

[1] 任永强.新时代北京地下水可持续发展工作探究[J].城市地质,202116(2)168-172(REN YongqiangDiscussion on sustainable groundwater development in new-era Beijing[J]Urban Geology202116(2)168-172)

[2] 贾三满,王海刚,赵守生,等.北京地面沉降机理研究初探[J].城市地质,20072(1)20-26(JIA SanmanWANG HaigangZHAO Shoushenget alA tentative study of the mechanism of land subsidence in Beijing[J]City Geology20072(1)20-26)

[3] 陶灿,曹成度,滕焕乐,等.基于长时序水准监测数据的北京东北部地区地面沉降趋势分析[J].铁道勘察,202046(4)17-23(TAO CanCAO ChengduTENG Huanleet alAnalysis and study of land subsidence trend in northeast Beijing based on long time series leveling monitoring data[J]Railway Investigation and Surveying202046(4)17-23)

[4] 冉培廉.基于时序InSAR技术的北京市地面沉降监测与影响因素分析[D].成都:成都理工大学,2020(RAN PeilianMonitoring and influence factors analysis of land subsidence in Beijing based on TS-InSAR technology[D]ChengduChengdu University of Technology2020)

[5] 周超凡,宫辉力,陈蓓蓓,等.北京地面沉降时空分布特征研究[J].地球信息科学学报,201719(2)205-215(ZHOU ChaofanGONG HuiliCHEN Beibeiet alStudy of temporal and spatial characteristics of land subsidence in Beijing[J]Journal of Geo-Information Science201719(2)205-215)

[6] 李海军,崔一娇,任永强,等.永定河2020年春季生态补水对北京地下水涵养效果分析[J].城市地质,202116(2)133-138(LI HaijunCUI YijiaoREN Yongqianget alAnalysis on the effect of ecological water replenishment of Yongding River in spring 2020 on groundwater conservation in Beijing[J]Urban Geology202116(2)133-138)

[7] .北京出台水污染防治条例[J].中国建设信息(水工业市场)2011(1)1(Beijing issuing water pollution prevention regulations [J]Information of China Construction (Water-Industry Market)2011(1)1)

[8] 王彪,姚旭初,石维新,等.永定河山峡段生态补水过程渗漏特征[J].南水北调与水利科技(中英文)202119(6)1208-1216(WANG BiaoYAO XuchuSHI Weixinet alEcological water supplement leakage characteristics in the canyon section of Yongding River[J]South-to-North Water Transfers and Water Science & Technology202119(6)1208-1216)

[9] 刘勇,李培英,丰爱平,等.黄河三角洲地下水动态变化及其与地面沉降的关系[J].地球科学,201439(11)1655-1665(LIU YongLI PeiyingFENG Aipinget alGroundwater dynamic evolutions and relationship between groundwater level and land subsidence in the Yellow River Delta[J]Earth Science201439(11)1655-1665)

[10] 杨国创,段春磊.联合InSARGNSS监测沿海地区地面沉降[J].测绘通报,2021(6)76-82(YANG GuochuangDUAN ChunleiJoint InSAR and GNSS to monitor land subsidence in coastal areas[J]Bulletin of Surveying and Mapping2021(6)76-82)

[11] 李瑞峰,常乐,秦海.InSAR监测技术与水准测量技术对比研究[J].工程质量,202139(3)72-76(LI RuifengCHANG LeQIN HaiComparative study on InSAR monitoring technology and leveling technology[J]Construction Quality202139(3)72-76)

[12] 蒲川豪,许强,赵宽耀,等.利用小基线集InSAR技术的延安新区地面抬升监测与分析[J].武汉大学学报·信息科学版,202146(7)983-993(PU ChuanhaoXU QiangZHAO Kuanyaoet alLand uplift monitoring and analysis in Yan'an new district based on SBAS-InSAR technology[J]Geomatics and Information Science of Wuhan University202146(7)983-993)

[13] 张双成,许强,罗勇,等.时序InSAR解译20172020年北京地面沉降时空变化[J].大地测量与地球动力学,202242(1)48-53(ZHANG ShuangchengXU QiangLUO Yonget alTemporal and spatial variation of land subsidence in Beijing from 2017 to 2020 interpreted by time series InSAR[J]Journal of Geodesy and Geodynamics202242(1)48-53)

[14] 李永生,张景发,李振洪,等.利用短基线集干涉测量时序分析方法监测北京市地面沉降[J].武汉大学学报(信息科学版)201338(11)1374-1377(LI YongshengZHANG JingfaLI Zhenhonget alLand subsidence in Beijing City from InSAR time series analysis with small baseline subset[J]Geomatics and Information Science of Wuhan University201338(11)1374-1377)

[15] 杨柠.永定河引黄生态补水长效机制初步探索[J].水利发展研究,202020(2)13-16(YANG NingPreliminary exploration on the long-term mechanism of ecological water replenishment from the Yellow River in Yongding River [J]Water Resources Development Research202020(2)13-16)

[16] 叶超,田芳,罗勇,等.北京地面沉降控制区划及防控措施[J].南京大学学报(自然科学)201955(3)440-448(YE ChaoTIAN FangLUO Yonget alControl zoning and prevention measures on land subsidence in Beijing[J]Journal of Nanjing University (Natural Science)201955(3)440-448)

[17] BERARDINO PFORNARO GLANARI Ret alA new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]IEEE Transactions on Geoscience and Remote Sensing200240(11)2375-2383

[18] 李珊珊,李志伟,胡俊,等.SBAS-InSAR技术监测青藏高原季节性冻土形变[J].地球物理学报,201356(5)1476-1486(LI ShanshanLI ZhiweiHU Junet alInvestigation of the seasonal oscillation of the permafrost over Qinghai-Tibet plateau with SBAS-InSAR algorithm[J]Chinese Journal of Geophysics201356(5)1476-1486)

[19] CASU FMANZO MLANARI RA quantitative assessment of the SBAS algorithm performance for surface deformation retrieval from DInSAR data[J]Remote Sensing of Environment2006102(3/4)195-210

[20] 张子文,杨帆,吴文豪,等.地下水开采与地面沉降关系的短基线集分析[J].测绘科学,201641(6)64-69(ZHANG ZiwenYANG FanWU Wenhaoet alStudy of land subsidence and groundwater based on SBAS time-series analysis[J]Science of Surveying and Mapping201641(6)64-69)

[21] MORA OMALLORQUI J JBROQUETAS ALinear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J]IEEE Transactions on Geoscience and Remote Sensing200341(10)2243-2253

[22] HANSSEN RRadar interferometryData interpretation and error analysis[Z].[S.l.]:Kluwer Academic Publishers,2001.

[23] 许军强,马涛,卢意恺,等.基于SBAS-InSAR技术的豫北平原地面沉降监测[J].吉林大学学报(地球科学版)201949(4)1182-1191(XU JunqiangMA TaoLU Yikaiet alLand subsidence monitoring in north Henan plain based on SBAS-InSAR technology[J]Journal of Jilin University (Earth Science Edition)201949(4)1182-1191)

[24] GOLDSTEIN R MWERNER C LRadar interferogram filtering for geophysical applications[J]Geophysical Research Letters199825(21)4035-4038

[25] COSTANTINI MA novel phase unwrapping method based on network programming[J]IEEE Transactions on Geoscience and Remote Sensing199836(3)813-821

[26] 杨帆,张磊,张子文,等.利用短基线集InSAR技术监测抚顺市地面沉降[J].测绘通报,2018(3)84-88(YANG FanZHANG LeiZHANG Ziwenet alLand subsidence in Fushun City from InSAR with small baseline subset[J]Bulletin of Surveying and Mapping2018(3)84-88)

[27] 唐益群,宋寿鹏,陈斌,等.不同建筑容积率下密集建筑群区地面沉降规律研究[J].岩石力学与工程学报,201029(S1)3425-3431(TANG YiqunSONG ShoupengCHEN Binet alStudy of land subsidence rule of dense buildings under different floor area ratios[J]Chinese Journal of Rock Mechanics and Engineering201029(S1)3425-3431)

[28] 张雯,宫辉力,陈蓓蓓,等.北京典型区地面沉降演化特征与成因分析[J].地球信息科学学报,201517(8)909-916(ZHANG WenGONG HuiliCHEN Beibeiet alEvolution and genetic analysis of land subsidence in Beijing typical area[J]Journal of Geo-Information Science201517(8)909-916)

[29] 郭海朋,白晋斌,张有全,等.华北平原典型地段地面沉降演化特征与机理研究[J].中国地质,201744(6)1115-1127(GUO HaipengBAI JinbinZHANG Youquanet alThe evolution characteristics and mechanism of the land subsidence in typical areas of the North China Plain[J]Geology in China201744(6)1115-1127)

[30] 田芳,罗勇,周毅,等.北京地面沉降与地下水开采时空演变对比[J].南水北调与水利科技,201715(2)163-169(TIAN FangLUO YongZHOU Yiet alContrastive analysis of spatial-temporal evolution between land subsidence and groundwater exploitation in Beijing[J]South-to-North Water Transfers and Water Science & Technology201715(2)163-169)

[31] 孙冉,潘兴瑶,王俊文,等.永定河(北京段)河道生态补水效益分析与方案评估[J].中国农村水利水电,2021(6)19-24(SUN RanPAN XingyaoWANG Junwenet alAn analysis and evaluation of ecological water replenishment benefit of Yongding River(Beijing section)[J]China Rural Water and Hydropower2021(6)19-24)

[32] 翟远征,姜亚,夏雪莲,等.永定河生态补水下渗对地下水质量的影响[J].中国环境科学,202242(4)1861-1868(ZHAI YuanzhengJIANG YaXIA Xuelianet alImpact of infiltration of ecology water replenishment of the Yongding River on groundwater quality and the mechanism[J]China Environmental Science202242(4)1861-1868)

[33] 马尧,杨勇,胡国金,等.永定河永定河(北京段)生态补水对地下水的补给分析[J].北京水务,2020(4)22-27(MA YaoYANG YongHU Guojinet alAnalysis of groundwater replenishment caused by water supplement in Yongding River(Beijing section)[J]Beijing Water2020(4)22-27)

作者简介:田枭(1999—),男,湖南怀化人,硕士研究生,主要研究方向为InSAR数据处理及应用。

E-mail2752997936@qq.com

基金项目:国家自然科学基金项目(41874011);国家重点研发计划项目(2018YFC1503603

通信作者:刘洋 副教授 E-mail: Yang.Liu@sgg.whu.edu.cn

转自:“测绘学术资讯”微信公众号

如有侵权,请联系本站删除!


本文评论

暂无相应记录!

首页<<1>>尾页共0页共0条记录
  • 万维QQ投稿交流群    招募志愿者

    版权所有 Copyright@2009-2015豫ICP备2021036211号

     纯自助论文投稿平台    E-mail:1121090112@qq.com;eshukan@163.com