Attribution analysis of precipitous decrease of sediment loads in the Hekou-Longmen section of Yellow River since 2000
GAO Haidong1,2, LIU Han2, JIA Lianlian3, PANG Guowei4, WANG Jie1收稿日期:2018-12-5修回日期:2019-07-28网络出版日期:2019-09-25
基金资助: |
Received:2018-12-5Revised:2019-07-28Online:2019-09-25
Fund supported: |
作者简介 About authors
x高海东(1983-),男,内蒙古乌审旗人,博士,讲师,中国地理学会会员(S110010718M),研究方向为土壤侵蚀与水土保持E-mail:hdgao@msn.cn。

摘要
关键词:
Abstract
Keywords:
PDF (4237KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
高海东, 刘晗, 贾莲莲, 庞国伟, 王杰. 2000-2017年河龙区间输沙量锐减归因分析. 地理学报[J], 2019, 74(9): 1745-1757 doi:10.11821/dlxb201909004
GAO Haidong.
1 引言
2000年以来,在气候变化和人类活动共同影响下,黄河输沙量呈显著减少趋势[1,2]。黄河水沙变化关系到整个流域的经济安全、能源安全、粮食安全和生态安全[3]。因此,气候变化以及黄土高原高强度生态建设对黄河水沙的影响机理研究,受到国内外****的高度关注。研究表明,大规模生态建设活动是黄河输沙量减少的主因[4]。认识黄河水沙变化特征,科学评估生态建设减沙贡献率,对厘清黄河水沙变化原因,预测未来黄河水沙变化趋势,具有重要意义。对于黄河输沙量变化的阶段性与基准期确定问题上,可采用双累积曲线法(DMC)、累积距平法(CAM)、有序聚类方法(OCM)等方法进行识别[5]。黄河水沙变化特征及其影响因子识别,国内外****开展了大量研究。研究方法有:气候弹性系数法[6,7]、水保法[8,9,10]、水文模拟法[11,12]、累积斜率变化率法[13,14]以及神经网络模型等[15]。如傅伯杰等[1]研究认为坝库、梯田等工程措施是20世纪70-90年代黄土高原产沙减少的主要原因,占54%,2000年以来,随着退耕还林(草)工程的实施,植被措施成为了土壤保持的主要贡献者,占57%。刘晓燕[16]研究发现,2007-2014年,潼关以上地区水库及淤地坝的实际减沙量为3.2亿t/年,林草和梯田等坡沟要素合计减沙12.54亿~14.11亿t/年。
对于气候变化和人类活动减沙贡献率的定量识别,目前已有多种方法。但是对于植被、淤地坝以及梯田等各项人类活动的减沙贡献率,还缺少有效的识别手段,有待发展新的计算方法。本文以河龙区间为研究对象,构建了包含气象数据、水沙数据、植被数据以及淤地坝数据在内的河龙区间地理数据库,根据建立的各项人类活动减沙贡献率计算公式分析2000-2017年间河龙区间输沙量变化特征,并揭示了2000年以来河龙区间输沙量锐减原因,从而提出河龙区间治理对策。研究结果可为黄土高原生态文明建设提供借鉴。
2 材料与方法
2.1 研究区概况
河口—龙门区间位于108°E~113°E、35°N~41°N,是黄土高原的主要组成部分,总面积11万km2(图1a)。河龙区间黄河干流长725 km。流域面积在1000 km2以上的较大支流有21条,主要有皇甫川、窟野河、秃尾河、无定河、延河等(图1b)。研究区属温带干旱半干旱地区,地貌以黄土丘陵地貌为主,主导植被类型为草原。土壤类型主要是黄绵土,土质疏松,易侵蚀。年降水量变化于300~800 mm,集中于7-8月,次降雨过程历时短、强度大,导致土壤侵蚀严重。为了遏制严重的土壤侵蚀,在此区域进行了大规模的生态建设活动,截至2011年底共修建大型淤地坝3786座(图1c),截至2015年底修建梯田6002 km2,经过退耕还林和封禁等措施,2017年植被覆盖度达到60%。经过治理,这一区域土壤侵蚀强度呈下降趋势,区间输沙量锐减。图1

图1研究区范围及大型淤地坝空间分布图
Fig. 1Location of the study area and distribution of large-sized dams
2.2 数据来源与处理
(1)降水量数据。本文共收集了河龙区间及周边46个气象站的逐日降水量数据,其中来自中国气象科学数据共享服务网(http://data.cma.cn/)的站点有21个,来自黄河水利委员会的站点有25个。采用章文波等[17]提出的方法计算降雨侵蚀力(R)值,公式为:式中:Ri为第i个半月的降雨侵蚀力值(MJ·mm/(hm2·h));k为该半月的天数(d);Dj为半月时段内第j天的侵蚀性日降雨量(mm)(研究区侵蚀性日降雨量标准为≥ 12 mm);Pd12为侵蚀性日降雨量平均值(mm);Py12表示日降雨量≥ 12 mm的年平均降雨量。
计算得到的降雨侵蚀力(R)值,在ArcGIS支持下,采用反距离权重插值(IDW)获得整个河龙区间的降雨侵蚀力(R)值。
(2)植被、梯田以及淤地坝数据。植被数据分别采用了GIMMS NDVI数据集(时间序列1982-2011年,空间分辨率8 km)和MODIS NDVI数据集(时间序列2000-2015年,空间分辨率250 m)。GIMMS NDVI用于确定权重系数(ωi),其中1950-1981年NDVI值以1982-1985年平均值近似。MODIS NDVI用于计算土壤侵蚀模数。根据MODIS NDVI数据集统计,河龙区间NDVI值在2000-2015年间呈增加趋势(图2a)。
图2

图2研究区NDVI、梯田面积及大型淤地坝数量变化
Fig. 2Change of NDVI, terrace area, and number of large-sized dams in the study area
使用像元二分模型计算植被覆盖度(FV):
式中:NDVImax和NDVImin分别为河龙区间内最大和最小NDVI值。NDVImax在NDVI累积概率95%确定,NDVImin在5%处确定。
作物覆盖与管理因子(C)采用江忠善等[18]提出的公式计算:
梯田数据来源于河龙区间各区县统计年鉴,2000-2015年梯田面积逐步增加(图2b)。大型淤地坝数据来源于第一次全国水利普查,调查了每座大型淤地坝的修建年份、控制面积、总库容以及截至2011年底的淤积量等数据,精度较高,是目前淤地坝拦沙效果评估的主要数据源。截至2011年底,河龙区间共有大型淤地坝3786座(图2c)。中小型坝数据来源于2008年淤地坝安全大检查。
(3)输沙量与径流量数据。输沙量与径流量数据来源于中华人民共和国水文年鉴和黄河水资源公报,时间尺度为年。河龙区间输沙量(径流量)定义为龙门水文站输沙量(径流量)与头道拐水文站输沙量(径流量)的差值。年引水量数据来源于黄河水资源公报。
(4)数字高程模型(DEM)和土壤数据。数字高程模型(DEM)来源于地理空间数据云平台(http://www.gscloud.cn),土壤类型数据来源于“黑河计划数据管理中心”(http://westdc.westgis.ac.cn)。
2.3 研究方法
2.3.1 水沙变化趋势分析 使用Mann-Kendall(M-K)趋势检验判断输沙量变化趋势[19],使用双累积曲线结合Pettitt检验来确定序列跃变时间[5, 20]。2.3.2 植被、梯田以及淤地坝减沙贡献率 根据生态措施的分布特征及其对侵蚀输沙影响分析可知,植被和梯田影响坡面的土壤侵蚀过程,而淤地坝影响流域泥沙输送。基于此,本文首先使用修正通用土壤流失方程(RUSLE)计算无梯田影响的坡面土壤侵蚀量,其次确定梯田的减沙贡献率,最后计算淤地坝的拦沙量。
(1)植被对土壤侵蚀影响。
RUSLE是被广泛使用的坡面土壤侵蚀模数计算工具[21],表达式为:
式中:A是年平均土壤侵蚀模数(t·/(hm2·a));R是降雨侵蚀力因子(MJ·mm/(hm2·h·a)),计算方法见公式(1);K是土壤可蚀性因子(t·hm2·h/(hm2·MJ·mm));S和L分别是坡度和坡长因子;C是作物覆盖—管理因子,计算方法见公式(5);P是水土保持措施因子。
(2)梯田减沙贡献率计算。
梯田减沙贡献率采用刘晓燕[16]提出的公式计算:
式中:ΔWs是减沙幅度(%);Ti是梯田比,指研究区梯田面积占轻度及以上土壤侵蚀强度面积的比例(%),轻度及以上土壤侵蚀强度面积根据公式(6)确定。
(3)淤地坝对产沙影响。
首先对大型淤地坝的泥沙总淤积量逐年还原;其次确定中型和小型淤地坝的淤积量。由于淤地坝的淤积速率取决于淤地坝控制范围内的土壤侵蚀强度,对于同一地点来讲,土壤可蚀性因子(K)、坡度坡长因子(SL)以及水土保持措施因子(P)较稳定,因此,年际土壤侵蚀强度的差异主要受降雨侵蚀力因子(R)及作物覆盖—管理因子(C)影响。因此,采用公式(8)确定权重系数(ωi),使用该权重系数,根据公式(9)对大型淤地坝的总淤积量进行逐年还原:
式中:BSKi为第i年大型淤地坝拦沙量(万t);Vt为大型淤地坝总拦沙量(万t);i为年份。
黄土高原淤地坝建设以坝系为基本单元,配置模式以大型坝为主体,内部布设中型坝和小型坝。因此,根据定义的淤地坝系结构系数(fs),确定坝系内的中型坝和小型坝拦沙量:
式中:AK为研究区大型坝的平均控制面积,根据调查,其值取5.00 km2;AM为中型坝控制面积,取1.65 km2;AS为小型坝控制面积,取0.70 km2;nK(i)、nM(i)、nS(i)分别为具有拦沙能力的大型坝、中型坝及小型坝数量。
有拦沙能力的淤地坝是指未淤满的淤地坝,淤地坝淤满时间使用公式(11)确定:
式中:ts为淤地坝淤满时间(a);Le为可拦沙库容占总库容的比例,根据文献[22]研究成果,大型淤地坝取0.77,中型坝取0.88;V为淤地坝总库容(m3);1.35为泥沙容重(t/m3);A为淤地坝控制面积(km2);SEM为坝控流域土壤侵蚀模数(t/(km2·a));c为淤地坝修建年代。根据调查,黄土高原小型淤地坝平均总库容为5万m3,平均控制面积0.70 km2。土壤侵蚀模数按照7500 t/(km2·a)计算,小型坝拦沙寿命为10年,根据统计的各个年代小型坝数量,按照平均10年拦沙运行期,估算各年份的有效小型坝数量。
根据淤地坝系结构系数,淤地坝的总拦沙量(BSi)计算公式为:
(4)其他措施减沙量。
河龙区间减沙因素还有水库淤积和引水减沙。根据刘晓燕[16]计算结果,研究区水库年均淤积量为0.88 亿t。引水减沙量采用引水量和年均含沙量的乘积确定,数据来自黄河水资源公报。
3 结果与分析
3.1 河龙区间水沙变化趋势
根据M-K趋势检验,河龙区间年降水量无显著变化趋势(p > 0.05),河龙区间年径流量(Z = -7.42,p < 0.001)和年输沙量(Z = -6.87,p < 0.001)均呈极显著减少趋势(图3)。图3

图3河龙区间降水量、径流量及输沙量变化趋势
Fig. 3Variation trends of precipitation, runoff, and sediment load in the Hekou-Longmen (HL) section
使用双累积曲线及Pettitt检验,河龙区间输沙量突变年份为1979年(p < 0.001)和1999年(p < 0.001)(图4)。因此,结合河龙区间生态建设实际,河龙区间输沙量变化可分为3个典型时期:基准期(P1,1952-1979年),此时段区间输沙量较大,治理强度低;工程措施影响期:(P2,1980-1999年),输沙量变化主要受梯田及淤地坝等工程措施影响;综合影响期:(P3,2000-2017年),输沙量变化受工程措施和植被措施综合影响。
图4

图4河龙区间降水量—输沙量双累积曲线图
Fig. 4Double cumulative curve of precipitation and sediment load in the HL section
相比于基准期(P1)时段,工程措施影响期(P2)和综合影响期(P3)时段的径流量分别减少了40.84%和67.47%,输沙量分别减少了54.84%和88.92%,输沙量减幅高于径流量减幅。P3时期,河龙区间平均年输沙量为1.03亿t(表1)。
Tab. 1
表1
表1各时段降水量、径流量及输沙量统计表
Tab. 1
时段 | 年均降水量(mm) | 年均径流量(亿m3) | 年均输沙量(亿t) |
---|---|---|---|
P1(1952-1979年) | 477 | 66.50 | 9.30 |
P2(1980-1999年) | 437 | 39.34 | 4.20 |
P3(2000-2017年) | 448 | 21.63 | 1.03 |
新窗口打开|下载CSV
3.2 大规模生态建设对输沙量影响
根据2.3.2节中的减沙量计算公式,首先计算降雨和植被影响下的坡面土壤模数,其次确定梯田减沙量,最后计算淤地坝拦沙量。坡面土壤侵蚀量减去梯田减沙量、淤地坝拦沙量、水库拦沙量和引水减沙量,得到计算输沙量,拟合显示,实测输沙量和计算输沙量之间的决定系数R2为0.7997,纳什效率系数(NSE)为0.80(图5)。图5

图5研究区输沙量计算值与输沙量实测值拟合曲线
Fig. 5Fitting curve of calculated annual sediment load and measured annual sediment load
3.2.1 坡面措施对土壤侵蚀影响 不考虑梯田影响,2000-2015年河龙区间坡面土壤侵蚀量变化介于2.61亿~5.83亿t/a,呈减少趋势,年均减少0.23亿t。土壤侵蚀模数较高的地区主要是研究区的中部和南部地区(图6)。
图6

图62000-2015年植被影响下的河龙区间土壤侵蚀模数变化
Fig. 6Change of soil erosion modulus in the HL section under the effect of vegetation in 2000-2015
随着坡面土壤侵蚀模数的降低和梯田面积增加,梯田比逐年增加,梯田减沙幅度从14.02%(2000年)增加到25.99%(2015年)。同时,受梯田和植被共同影响,研究区坡面土壤侵蚀量变化介于1.90亿~5.13亿t/a,呈降低趋势(表2)。
Tab. 2
表2
表22000-2015年河龙区间坡面土壤侵蚀量变化
Tab. 2
年份 | 植被影响下的坡面年土壤侵蚀量(亿t) | 梯田比(%) | 梯田减沙贡献率(%) | 植被和梯田共同作用下的坡面年土壤侵蚀量(亿t) |
---|---|---|---|---|
2000 | 3.97 | 6.92 | 14.02 | 3.42 |
2001 | 5.83 | 6.42 | 11.96 | 5.13 |
2002 | 5.82 | 6.61 | 12.74 | 5.08 |
2003 | 5.38 | 6.80 | 13.51 | 4.65 |
2004 | 3.92 | 7.25 | 15.42 | 3.32 |
2005 | 3.39 | 7.94 | 18.52 | 2.76 |
2006 | 3.80 | 8.05 | 18.98 | 3.08 |
2007 | 4.06 | 7.43 | 16.22 | 3.40 |
2008 | 3.13 | 8.13 | 19.38 | 2.53 |
2009 | 3.69 | 8.31 | 20.22 | 2.94 |
2010 | 3.17 | 8.32 | 20.27 | 2.52 |
2011 | 3.09 | 9.02 | 23.54 | 2.36 |
2012 | 3.59 | 8.13 | 19.39 | 2.90 |
2013 | 4.09 | 8.16 | 19.50 | 3.29 |
2014 | 2.61 | 9.78 | 27.21 | 1.90 |
2015 | 2.90 | 9.53 | 25.99 | 2.15 |
新窗口打开|下载CSV
3.2.2 沟道治理工程拦沙贡献率 1952-2011 年,河龙区间淤地坝逐年拦沙量变化介于0~3.91亿t,累积拦沙量为73.17亿t,1952-2011年平均拦沙量1.22亿t(图7)。在P1时期,年均拦沙量0.93亿t,P2时期年均拦沙量为1.53亿t,P3时期年均拦沙量1.38亿t。根据计算,P2时期淤地坝减沙贡献率为30%,P3时期降至17%。
图7

图7河龙区间淤地坝拦沙量计算结果
Fig. 7The amount of sediment reduction by check dams in the HL section
综上,P1时期区间输沙量为9.30亿t,P3时期,降至1.03亿t,降幅8.27亿t。进一步可以确定各措施贡献率。结果显示,植被贡献率为54%。梯田和淤地坝各贡献了17%,二者合计贡献了34%。水库拦沙和引水取沙贡献了12%。植被恢复是河龙区间输沙量减沙的主要原因。
4 讨论
4.1 植被恢复对水沙过程影响
关于黄河泥沙减少原因,目前取得的研究共识主要有两个:一是人类活动是黄河泥沙减少的主要原因[23,24],二是植被恢复是主要原因[1, 16]。这与本文研究结果一致。另外,时鹏等2019年在无定河流域的研究表明,淤地坝减沙贡献率为11.7%,淤地坝减沙贡献比例与本研究接近[25]。进一步分析植被恢复对水沙过程影响,计算NDVI和径流量(Q)、输沙量(S)以及含沙量(Cs)之间的相关系数(表3)。发现NDVI与径流量之间相关性较差,NDVI与输沙量之间的相关关系好于其与径流量的相关关系,且为负相关,而NDVI与含沙量的相关性最高,即植被恢复主要导致径流含沙量降低。河龙区间以黄土丘陵区沟壑区地貌为主,从分水岭至沟底,由峁边线、沟坡线将流域划分为3个部分:梁峁坡、沟谷坡以及沟谷底。梁峁坡地形条件较好,植被恢复速度和覆盖度一般好于沟谷坡,根据图2a显示,研究区植被呈增加趋势,随着坡面植被恢复加快,坡面径流含沙量减少,下沟径流挟沙能力增大,可能会加大沟谷坡侵蚀。
Tab. 3
表3
表3河龙区间主要支流NDVI与径流量、输沙量以及含沙量相关系数
Tab. 3
流域(水文站) | NDVI与Q | NDVI与S | NDVI与Cs |
---|---|---|---|
皇甫川(皇甫) | 0.116 | -0.215 | -0.657* |
孤山川(高石崖) | -0.015 | -0.626* | -0.897** |
窟野河(温家川) | 0.390 | -0.645* | -0.783** |
秃尾河(高家川) | -0.537* | -0.454 | -0.439 |
佳芦河(申家湾) | 0.402 | 0.161 | -0.166 |
无定河(白家川) | 0.556* | -0.404 | -0.479 |
清涧河(延川) | -0.045 | -0.280 | -0.448 |
延河(甘谷驿) | 0.112 | -0.332 | -0.504 |
新窗口打开|下载CSV
4.2 泥沙输移比变化
淤地坝是黄土高原沟壑整治的主要措施,坡面措施主要对土壤侵蚀过程产生影响,而沟道措施主要影响流域产沙过程,即对泥沙输移比(SDR)造成影响。一般认为,没有坝库工程,黄土丘陵区泥沙输移比将接近于1[26,27,28]。选择河龙区间皇甫川、孤山川、佳芦河以及大理河4条支流为分析对象,分析淤地坝建设对泥沙输移比影响。这4条流域,主河道均没有水库工程,除佳芦河外,其他流域内没有中型及以上水库,已建水库控制面积小于10%,可以最大程度减少水库对泥沙输移比的影响。根据计算结果(表4),4条流域泥沙输移比均减少,如皇甫川流域,从1965-1979年的0.95降至2000年后的0.47。进入2000年以来,4条流域的泥沙输移比均不超过0.5。
Tab. 4
表4
表4河龙区间典型流域泥沙输移比变化
Tab. 4
流域 | 项目 | 1965-1979年 | 1980-1999年 | 2000-2011年 |
---|---|---|---|---|
皇甫川 | 输沙量 | 0.61 | 0.34 | 0.08 |
拦沙量 | 0.03 | 0.05 | 0.09 | |
泥沙输移比 | 0.95 | 0.87 | 0.47 | |
孤山川 | 输沙量 | 0.29 | 0.12 | 0.02 |
拦沙量 | 0.02 | 0.03 | 0.02 | |
泥沙输移比 | 0.94 | 0.80 | 0.50 | |
佳芦河 | 输沙量 | 0.25 | 0.06 | 0.02 |
拦沙量 | 0.11 | 0.11 | 0.05 | |
泥沙输移比 | 0.69 | 0.38 | 0.29 | |
大理河 | 输沙量 | 0.51 | 0.22 | 0.19 |
拦沙量 | 0.15 | 0.25 | 0.19 | |
泥沙输移比 | 0.77 | 0.47 | 0.50 |
新窗口打开|下载CSV
4.3 河龙区间未来治理对策
根据2017年NDVI反演计算,河龙区间平均植被覆盖度为60%,整个河龙区间植被覆盖度可达68%[29]。河龙区间40%的区域年降水量大于450 mm,整个研究区植被还有一定的提升空间,根据公式(5),植被覆盖度每提高1%,流域土壤侵蚀模数约降4%。可以预见,未来随着植被恢复,流域土壤侵蚀模数会进一步降低。退耕还林(草)工程是河龙区间植被恢复的主要外动力,目前,河龙区间仍有一定数量的坡耕地,根据在黄土高原调查,坡耕地单产一般为1000 kg/hm2左右,黄土高原460万hm2坡耕地,粮食总产量为460万t,2013年,黄土高原粮食总产量为4712万t,坡耕地以40%的比例贡献了大约10%的产量。因此,仍应持续推进退耕还林(草)工程,力争做到河龙区间无坡耕地的理想状态。根据地形条件和筑坝材料限制,整个河龙区间适宜修建淤地坝的面积为5.47万km2。按照大型淤地坝控制面积不超过8 km2估算,整个河龙区间大型淤地坝布设潜力为6838座。将流域现存大型淤地坝数量和大型淤地坝布设潜力的比值称为淤地坝布设强度。按照25%、50%以及75%划分标准,可以看出河龙区间中部淤地坝布设强度高(图8)。对于淤地坝布设强度较高的地区,淤地坝工作的重点为除险加固,而布设强度较低地区,应该合理规划,适当推进淤地坝建设。
图8

图8河龙区间淤地坝布设强度与限制区分布
Fig. 8Construction strength of check dams and restricting areas in the HL section
5 结论
(1)河龙区间年输沙量呈极显著下降趋势,其突变年份为1979年和1999年。1952-1979年区间年均输沙量为9.30亿t,至2000-2017年,年均输沙量大幅减少至1.03,减幅89%。(2)基于各生态建设措施的减沙量计算方法,本文评估了各措施的减沙贡献率,植被贡献达54%,梯田和淤地坝合计贡献了34%。
(3)NDVI与径流量之间相关性较差,而与含沙量的相关性较好,植被恢复主要导致径流含沙量降低。淤地坝建设降低了流域泥沙输移比,进入2000年以来,河龙区间典型支流泥沙输移比均不超过0.5。
(4)河龙区间未来植被还有一定的提升空间,仍应持续推进退耕还林(草)工程,力争做到无坡耕地的理想状态。随着植被恢复,流域土壤侵蚀模数会进一步降低。对于淤地坝工程,应因地制宜,布设强度较高地区,应加强除险加固;布设强度较低地区,应合理规划,适当推进淤地坝建设。
参考文献 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子
[本文引用: 3]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 2]
[本文引用: 2]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 4]
[本文引用: 4]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
DOI:10.2307/2346729URL [本文引用: 1]

Non-parametric techniques are introduced for the change-point problem. Exact and approximate results are obtained for testing the null hypothesis of no change. The methods are illustrated by the analysis of three sets of data illustrating the techniques for zero-one observations, Binomial observations and continuous observations. Some comparisons are made with methods based on CUSUMS. [ABSTRACT FROM AUTHOR]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]
[本文引用: 1]