• 中文核心期刊
  • 中国科技核心期刊
  • RCCSE中国核心学术期刊

开采沉陷动态预测时间函数参数变化规律研究——以Weibull为例

喻成林, 张宏贞, 范洪冬

喻成林,张宏贞,范洪冬. 开采沉陷动态预测时间函数参数变化规律研究——以Weibull为例[J]. 煤矿安全,2024,55(3):175−180. DOI: 10.13347/j.cnki.mkaq.20230032
引用本文: 喻成林,张宏贞,范洪冬. 开采沉陷动态预测时间函数参数变化规律研究——以Weibull为例[J]. 煤矿安全,2024,55(3):175−180. DOI: 10.13347/j.cnki.mkaq.20230032
YU Chenglin, ZHANG Hongzhen, FAN Hongdong. Study on change law of time function parameters for dynamic prediction of mining subsidence: Taking Weibull as an example[J]. Safety in Coal Mines, 2024, 55(3): 175−180. DOI: 10.13347/j.cnki.mkaq.20230032
Citation: YU Chenglin, ZHANG Hongzhen, FAN Hongdong. Study on change law of time function parameters for dynamic prediction of mining subsidence: Taking Weibull as an example[J]. Safety in Coal Mines, 2024, 55(3): 175−180. DOI: 10.13347/j.cnki.mkaq.20230032

开采沉陷动态预测时间函数参数变化规律研究——以Weibull为例

基金项目: 国家自然科学基金资助项目(42274054)
详细信息
    作者简介:

    喻成林(1977—),男,高级工程师,硕士,从事土地整治、采煤塌陷地综合治理技术及自然资源管理方面的工作。E-mail:ychenglin007@163.com

  • 中图分类号: TD325.4

Study on change law of time function parameters for dynamic prediction of mining subsidence: Taking Weibull as an example

  • 摘要:

    针对开采沉陷动态预测模型中时间函数参数主要以最大下沉点的下沉时间序列为基础反演获取,在工作面上方空间位置变化规律研究不足的问题;依据工作面推进过程中覆岩受力、垮落情况,在平面方向上将采空区分为竖向整体受压区、中间区、动态完全区3个区间;以Weibull时间函数为例,分析了动态预测mk 2个参数的含义、空间分布规律及特点,并构建2个参数的分区数学模型及选取准则;分别以实测最终值、概率积分法预测值为基准,通过实例验证了其正确性,丰富了开采沉陷动态预测参数研究。

    Abstract:

    In view of the fact that the time function parameters in the dynamic prediction model of mining subsidence are mainly obtained by inversion based on the subsidence time series of the maximum subsidence point, its research on the spatial position change law above the working face is insufficient. According to the stress and collapse of the overburden during the advancing process of the working face, the goaf is divided into three sections in the plane direction: vertical overall compression zone, middle zone and dynamic full zone; taking Weibull time function as an example, the meaning, spatial distribution law and characteristics of the parameters m and k of dynamic prediction are analyzed, and the partition mathematical model and selection criteria of the two parameters are constructed; based on the measured final value and the predicted value of probability integral method, the correctness is verified through examples, which enriches the research on parameters of mining subsidence dynamic prediction.

  • 地下资源开采后,破坏其原有平衡状态,经过覆岩传递到地表,是1个复杂的动态移动过程,涉及地质、采矿、空间等多方面因素,许多学者在开采沉陷动态变形原理、模型、精度等方面进行了研究,并取得了大量成果。开采沉陷动态预测模型主要有动态力学模型和时间函数模型2类,目前主要采用时间函数构建开采沉陷的动态变形过程;常用的时间函数有Knothe时间函数模型、双曲线时间函数模型、Gompertz时间函数模型、logistic 曲线时间函数模型及 Weibull 曲线时间函数模型等。

    王悦汉等[1]将动态移动过程分为顶板初次垮落前、顶板垮落未充满采空区、顶板垮落已充满采空区、表土移动4个阶段,采动岩体分为连续介质、拟连续介质、非连续介质3种介质,分阶段建立了动态力学预测模型;吴侃等[2]提出了基于时间序列分析的动态预测模型;崔希民等[3]建立基于Knothe时间函数的概率积分预计方法,根据临界开采尺寸确定时间影响系数;张兵等[4-5]建立了1种新的优化分段Knothe时间函数,并提出了基于实测数据的反算时间函数对比求参法、基于概率积分参数的直接计算法的2种优化分段Knothe时间函数参数求取方法;刘玉成等[6-7]分析了下沉动态过程的Knothe时间函数模型的不足,提出了幂指数的Knothe时间函数,分析了常见地表沉陷时间函数模型的下沉、下沉速度、下沉加速度与时间关系,认为Weibull 曲线时间函数模型能完整地描述地表沉陷的动态过程;孙闯等[8]建立了考虑松散层下沉特点的双因素Knothe时间函数;王军保等[9]借鉴岩石流变力学中非定常流变模型,构建了变时间影响系数的Knothe的地表动态沉陷预测模型;许国胜等[10]回归分析得到地表动态移动变形参数与地质及开采技术参数之间的关系;李春意等[11-12]进行了基于Logistic时间函数、正态分布时间函数的动态沉陷预测;张永胜等[13]提出了以地表点最大下沉速度时刻为分界点,结合偏差改正、生长函数模型的1种分段Weibull时间函数;刘东海[14]根据常村观测资料,指出Weibull时间函数参数不固定,两参数变化规律为距边界变化类偏态分布;JAROSZ等[15]考虑开采产生的空间收敛,建立了双参数Sroka-Schober时间函数;KOWALSKI[16]在Knothe和Sroka-Schober双参数时间函数基础上,提出了广义时间函数;González Nicieza等[17]在Knothe时间函数基础上,引入正态分布函数,建立了正态分布时间函数。

    综上,总体认为地表动态下沉量与时间的关系曲线呈“S”形曲线,下沉速度为0→Vmax→0过程,下沉加速度为0→+amax→-amax→0过程。

    工作面上方的覆岩及地表在井下资源开采后,经历弯曲、垮落、再次压缩等过程,但在工作面不同位置的岩层完整厚度、位移、受力等经历过程不同,根据工作面回采、覆岩受力与垮落情况,将工作面回采结束后的地表分为3个部分,工作面动态分析分区图如图1

    图  1  工作面动态分析分区图
    Figure  1.  Partition diagram of dynamic analysis of working face

    1)竖向整体受压区(Ⅰ区)。处于工作面外侧煤柱区域,工作面开采边界至开采影响边界之间。此部分区域在整个回采期间,覆岩经历了弯曲,但没有垮落、再次压缩等过程,在竖向上是完整的。在整个开采期间竖向方向整体是受压状态。

    2)中间区(Ⅱ区)。处于工作面内测靠近煤柱区域,工作面开采边界至覆岩动态完全区边界(覆岩垮落角确定边界与工作面开采边界之间)。此部分区域在整个回采期间,覆岩经历了弯曲、垮落、非完全再次压缩过程等部分或全部过程,主要是由于岩梁的存在导致垮落非充分及垮落岩石受力非充分,在竖向上部分破裂,且离工作面边界(煤柱)越远,破裂高度越大,即完整覆岩厚度减小。在整个回采期间主要经历如下几个阶段:①受到影响到工作面推到下方前,在竖直方向上整体是受到压缩变形,竖向受力增大;②推过后至覆岩垮落前,在竖直方向上整体是受到拉伸变形、竖向受力减小;③垮落后至移动稳定,在竖直方向上整体是受到压缩变形,竖向受力增大至稳定。垮落岩石部分再次受到压缩变形。

    3)动态完全区(Ⅲ区)。处于工作面内测远离煤柱区域,覆岩动态完全区(覆岩垮落角确定边界至采空区中心)。此部分区域在整个回采期间,覆岩经历了弯曲、垮落、完整再次压缩过程,在竖向上部分破裂,破裂高度一致。在整个回采期间主要经历如下几个阶段:①受到影响到工作面推到下方(垮落前),在竖直方向上整体是受到压缩变形,竖向受力增大;②到再次受到压缩前,在竖直方向上整体是受到拉伸变形,竖向受力减小;③到移动稳定,在竖直方向上整体是受到压缩变形,竖向受力增大至稳定。

    根据开采沉陷岩梁相关理论,地表下沉可简化为各种类型岩梁弯曲,其弯曲下沉值与刚度(弹性模量与惯性矩之积)成反比,其中惯性矩与岩梁高度的三次方成正比。结合图1,在同种情况下,上覆岩层刚度、惯性矩变化规律为:Ⅰ区(整体一致,最大)>Ⅱ区(与离开切眼距离成反比)>Ⅲ区(整体一致,最小)。

    目前以时间函数为基础的各种地表动态沉陷预测模型参数主要包括:最大变形值、时间函数(时间函数不同,其参数不同),各类参数与覆岩、开采等因素有关,体现其差异,但目前时间函数都采用相同参数,即在工作面不同位置(煤柱内外)的时间函数参数都一样,这与处在工作面不同位置的受力过程、变形过程、岩层破坏情况及岩石力学性质等明显不符。

    Weibull时间函数f(t)数学表达式为:

    $$ f\left(t\right)=\left\{\begin{array}{cc}km{t}^{m-1}{{\mathrm{e}}}^{-k{t}^{m}}& t\geqslant 0\\ 0& t < 0\end{array}\right. $$ (1)

    式中:km为参数,正数;t为时间。

    据此建立开采沉陷过程中地表下沉时间函数、下沉速度时间函数、下沉加速度时间函数,开采沉陷动态下沉过程曲线如图2

    图  2  Weibull时间函数下沉、速度、加速度变化图
    Figure  2.  Weibull time function subsidence, velocity and acceleration variation diagram

    下沉时间函数W(t):

    $$ W\left(t\right)={W}_{0}\left(1-{{\mathrm{e}}}^{-k{t}^{m}}\right) $$ (2)

    式中:W0为该点下沉的最终值。

    下沉速度时间函数v(t):

    $$ v\left(t\right)={W}_{0}km{t}^{m-1}{{\mathrm{e}}}^{-k{t}^{m}} $$ (3)

    下沉加速度时间函数a(t):

    $$ a\left(t\right)={W}_{0}km\left(m-1-km{t}^{m}\right){t}^{m-2}{{\mathrm{e}}}^{-k{t}^{m}} $$ (4)

    根据理论分析及Weibull时间函数上述公式可知:① 当地面达到最大下沉速度时,下沉加速度为0,此时时间tvmax满足${t}_{vmax}=\sqrt[\leftroot{-1}\uproot{8}{{m}}]{{(m-1)}/{km}} $;② 由下沉加速度变化规律可知$m > 1 $;③ k反映了加速时间、加速度大小,k越小,加速度越小,到达最大下沉速度的时间越大;④ m越小,到达最大下沉速度的时间越小,最大下沉速度越小。

    根据前面分析可知:工作面不同位置在回采期间,其覆岩移动、破坏、受力、岩梁刚度等不同,则其动态时间函数参数应不同,mk决定了Weibull时间函数的形态变化过程,因此mk在工作面不同位置其存在差异,以垮落区边界为原点、指向煤柱方向为正,根据前述岩梁弯曲特征简述及相关岩石力学理论,建立其与地面点距离垮落区域边界距离S的函数,工作面不同位置mk变化情况图如图3

    图  3  工作面不同位置mk变化情况图
    Figure  3.  Changes of m and k at different positions of working face

    在煤柱侧区域(Ⅰ区)mk为:

    $$ m={a}_1 {{\mathrm{e}}}^{-b_1\frac{s}{H}}+c_1 $$ (5)
    $$ k={a}_2 \frac{{{S}}_{0}}{H}+b_2 $$ (6)

    在垮落区外侧至工作面边界(Ⅱ区)mk为:

    $$ {m}={a}_1 {{\mathrm{e}}}^{-b_1\frac{S}{H}}+c_1 $$ (7)
    $$ {k}={a}_2 \frac{{s}}{{H}}+{b}_2 $$ (8)

    在垮落区内侧(Ⅲ区)mk为:

    $$ {m}={a}_1+c_1 $$ (9)
    $$ {k}={b}_2 $$ (10)

    式中: a1b1c1a2b2为系数;S0为煤柱边界距离;H为工作面采高。

    mk计算公式及图3中可以看出:

    1)在垮落区外侧为负指数关系,以工作面边界为界,此区域分为2部分;Ⅰ区(煤柱范围)m变化较小,基本为常数;Ⅱ区(中间区)m变化较大,在垮落区内侧,m为一常数。

    2)以工作面边界、垮落区边界为界,此区域分为3部分;Ⅰ区(煤柱范围)k为常数;Ⅱ区(中间区)m为线性变化;在垮落区内侧,k为一常数。

    由于不同位置上覆岩层经历的受力、破坏、变形过程不同,根据开采沉陷相关研究,可采用如下方法确定垮落范围d。① 方法1:采用覆岩力学理论,根据岩梁的相关破裂理论确定垮落边界;② 方法2:采用概率积分法拐点偏移距;③ 方法3:采用顶板覆岩垮落角φ确定,$d={H}/{\mathrm{tan}\;\;\varphi } $。

    某工作面采用长壁一次性采全高采煤法回采。走向为774 m,倾斜宽为220.7 m。工作面开采深度为470~515 m。根据走向观测数据,采用上面的计算模型,获得的反演模型系数见表1

    表  1  反演模型系数
    Table  1.  Coefficients of inversion model
    点号 k m R2
    Z17 0.271 3 0.294 3 0.469 7
    Z18 0.163 1 0.386 6 0.629 3
    Z19 0.621 7 0.430 6 0.882 6
    Z20 0.131 8 0.401 2 0.945 1
    Z21 0.293 1 0.322 1 0.899 5
    Z22 0.271 7 0.365 3 0.880 4
    Z23 0.135 8 0.530 0 0.925 4
    Z24 0.009 9 1.133 0 0.924 2
    Z25 0.002 2 1.494 0 0.928 5
    Z26 0.001 4 1.545 0 0.952 1
    Z27 0.001 7 1.449 0 0.941 8
    Z28 0.002 5 1.301 0 0.927 0
    Z29 0.002 7 1.261 0 0.939 5
    Z30 0.001 0 1.414 0 0.941 2
    下载: 导出CSV 
    | 显示表格

    根据工作面开采情况,以工作面内侧90 m为垮落边界,以此确定各点距离边界的S值,根据表1反演动态时间函数$\varphi \left(t\right) $如下:

    $$ \varphi \left(t\right)=1-{{\mathrm{e}}}^{-k{t}^{m}} $$ (11)

    在采空区垮落内侧mk为:

    $$ m=1.419\;4, k=0.001\;908 $$ (12)

    在采空区垮落外侧mk为:

    $$ m\left(S\right)=1.099\cdot {{\mathrm{e}}}^{-12.97\times \frac{S}{495}}+0.320\;4 \;\;\;(R^{2}=0.921\;6) $$ (13)
    $$ k=1.517\times \frac{{S}}{495}+0.001\;908 \;\;\;(R^{2}=0.919\text{,}S<0.2H) $$
    $$ k=0.305 \qquad\qquad (S>0.2H)$$ (14)

    根据工作面开采情况及参数确定情况,分别以开采结束后实测结果作为最终值、概率积分法预测值作为最终值,进行开采第523 d的动态预测,实测结果、概率积分法预测值作为最终值的动态预测与实测对比如图4(图中距离以工作面边界为原点,采空区方向为正)。

    图  4  拟合与实测对比图(523 d)
    Figure  4.  Comparison between fitting and actual measurement (523 d)

    以实测结果作为最终值:动态预测中误差为110 mm,为最大下沉值的3.5%;以概率积分法预测结果作为最终值:动态预测中误差为130 mm,为最大下沉值的4.2%。

    1)分析了工作面上方不同位置的覆岩在开采过程中在竖向上的受力、变形、破坏等变化过程,根据煤柱、垮落情况将覆岩与地表划分为竖向整体受压区、中间区、动态完全区3个区域。

    2)以Weibull时间函数动态预测模型为例,分析了mk 2个参数的空间分布规律及特地,构建了2个参数的分区数学模型及选取准则。

    3)分别以实测结果、概率积分法预测结果作为工作面开采结束后的最终值,依据本文动态预测模型及参数确定方法,动态预测中误差分别为110、130 mm,为最大下沉值的3.5%、4.2%,证明本模型的可靠性。

  • 图  1   工作面动态分析分区图

    Figure  1.   Partition diagram of dynamic analysis of working face

    图  2   Weibull时间函数下沉、速度、加速度变化图

    Figure  2.   Weibull time function subsidence, velocity and acceleration variation diagram

    图  3   工作面不同位置mk变化情况图

    Figure  3.   Changes of m and k at different positions of working face

    图  4   拟合与实测对比图(523 d)

    Figure  4.   Comparison between fitting and actual measurement (523 d)

    表  1   反演模型系数

    Table  1   Coefficients of inversion model

    点号 k m R2
    Z17 0.271 3 0.294 3 0.469 7
    Z18 0.163 1 0.386 6 0.629 3
    Z19 0.621 7 0.430 6 0.882 6
    Z20 0.131 8 0.401 2 0.945 1
    Z21 0.293 1 0.322 1 0.899 5
    Z22 0.271 7 0.365 3 0.880 4
    Z23 0.135 8 0.530 0 0.925 4
    Z24 0.009 9 1.133 0 0.924 2
    Z25 0.002 2 1.494 0 0.928 5
    Z26 0.001 4 1.545 0 0.952 1
    Z27 0.001 7 1.449 0 0.941 8
    Z28 0.002 5 1.301 0 0.927 0
    Z29 0.002 7 1.261 0 0.939 5
    Z30 0.001 0 1.414 0 0.941 2
    下载: 导出CSV
  • [1] 王悦汉,邓喀中,吴侃,等. 采动岩体动态力学模型[J]. 岩土力学与工程学报,2003,22(3):352−357

    WANG Yuehan, DENG Kazhong, WU Kan, et al. On the dynamic mechanics model of mining subsidence[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(3): 352−357

    [2] 吴侃,靳建明,王卷乐. 开采沉陷时序预测技术[J]. 测绘工程,2001,10(4):20−23.

    WU Kan, JIN Jianming, WANG Juanle. The technique of time series prediction on mining subsidence[J]. Engineering of Surveying and Mapping, 2001, 10(4): 20−23

    [3] 崔希民,缪协兴,赵英利,等. 论地表移动过程的时间函数[J]. 煤炭学报,1999,24(5):453−455. doi: 10.3321/j.issn:0253-9993.1999.05.002

    CUI Ximin, MIAO Xiexing, ZHAO Yingli, et al. Discussion on the time function of time dependent surface movement[J]. Journal of China Coal Society, 1999, 24(5): 453−455 doi: 10.3321/j.issn:0253-9993.1999.05.002

    [4] 张兵,崔希民. 开采沉陷动态预计的分段Knothe时间函数模型优化[J]. 岩土力学,2017,38(2):541−548.

    ZHANG Bing, CUI Ximin. Optimization of segmented Knothe time function model for dynamic prediction of mining subsidence[J]. Rock and Soil Mechanics, 2017, 38(2): 541−548.

    [5] 张兵,崔希民,赵玉玲,等. 优化分段Knothe时间函数求参方法[J]. 煤炭学报,2018,43(12):3379−3386.

    ZHANG Bing, CUI Ximin, ZHAO Yuling, et al. Parameter calculation method for optimized segmented Knothe time function[J]. Journal of China Coal Society, 2018, 43(12): 3379−3386.

    [6] 刘玉成,庄艳华. 地下采矿引起的地表下沉的动态过程模型[J]. 岩土力学,2009,30(11):3406−3416. doi: 10.3969/j.issn.1000-7598.2009.11.031

    LIU Yucheng, ZHUANG Yanhua. Model for dynamic process of ground surface subsidence due to underground mining[J]. Rock and Soil Mechanics, 2009, 30(11): 3406−3416. doi: 10.3969/j.issn.1000-7598.2009.11.031

    [7] 刘玉成,曹树刚,刘延保. 可描述地表沉陷动态过程的时间函数模型探讨[J]. 岩土力学,2010,31(3):925−931. doi: 10.3969/j.issn.1000-7598.2010.03.044

    LIU Yucheng, CAO Shugang, LIU Yanbao. Discussion on some time functions for describing dynamic course of surface subsidence due to mining[J]. Rock and Soil Mechanics, 2010, 31(3): 925−931 doi: 10.3969/j.issn.1000-7598.2010.03.044

    [8] 孙闯,徐乃忠,刘义新,等. 基于双因素时间函数的松散地层条件下地表点动态沉降预计[J]. 岩土力学,2017,38(3):821−826.

    SUN Chuang, XU Naizhong, LIU Yixin, et al. Prediction of dynamic subsidence of ground point under loose stratum based on two-factor time function[J]. Rock and Soil Mechanics, 2017, 38(3): 821−826.

    [9] 王军保,刘新荣,刘小军. 开采沉陷动态预测模型[J]. 煤炭学报,2015,40(3):516−521

    WANG Junbao, LIU Xinrong, LIU Xiaojun. Dynamic prediction model for mining subsidence[J]. Journal of China Coal Society, 2015, 40(3): 516−521.

    [10] 许国胜,李德海,侯得峰,等. 厚松散层下开采地表动态移动变形规律实测及预测研究[J]. 岩土力学,2016,37(7):2056−2062.

    XU Guosheng, LI Dehai, HOU Defeng, et al. Measurement and prediction of the transient surface movement and deformation induced by mining under thick alluvium[J]. Rock and Soil Mechanics, 2016, 37(7): 2056−2062.

    [11] 李春意,赵亮,李铭,等. 基于Logistic时间函数地表动态沉陷预测及优化求参研究[J]. 安全与环境学报,2020,20(6):2202−2210.

    LI Chunyi, ZHAO Liang, LI Ming, et al. Prediction of surface progressive subsidence and optimization of predicting model parameters based on the Logistic time function[J]. Journal of Safety and Environment, 2020, 20(6): 2202−2210.

    [12] 李春意,高永格,崔希民. 基于正态分布时间函数地表动态沉陷预测研究[J]. 岩土力学,,2016,37(S1):108−116

    LI Chunyi, GAO Yongge, CUI Ximin. Progressive subsidence prediction of ground surface based on the normal distribution time function[J]. Rock and Soil Mechanics, 2016, 37(S1): 108−116.

    [13] 张永胜,胡海峰,蔡音飞,等. 分段Weibull时间函数开采沉陷动态预计研究[J]. 太原理工大学学报,2021(3):343−349.

    ZHANG Yongsheng, HU Haifeng, CAI Yinfei, et al. Study on dynamic prediction of mining subsidence with segmented Weibull time function[J]. Journal of Taiyuan University of Technology, 2021(3): 343−349.

    [14] 刘东海. 常村煤矿采区动态沉陷预测模型研究[D]. 西安:西安科技大学,2020.
    [15]

    JAROSZ A, KARMIS M, SROKA A. Subsidence development with time-experiences from longwall operations in the Appalachian coalfield[J]. International Journal of Mining and Geological Engineering, 1990, 8(3): 621−673

    [16]

    KOWALAKI A. Surface subsidence and rate of its increments based on measurements and theory[J]. Archives of Mining Science, 2001, 46(4): 391−406.

    [17]

    GONZALEZ-Nicieza C, ALVAREZ-Fernandez M I, MENENDEZ-Diaz A, et al. The influence of time on subsidence in the Central Austrian Coalfield[J]. Bulletin of Engineering Geology and the Environment, 2007, 66(3): 319−329 doi: 10.1007/s10064-007-0085-2

  • 期刊类型引用(5)

    1. 李鹏. 基于FLAC3D模拟的14103辅运顺槽支护参数设计. 江西煤炭科技. 2025(01): 90-92+96 . 百度学术
    2. 马清水,郭瑞,韩伟,宋永明,梁燕翔,刘耀,王佳明. 基于长短期记忆网络的煤巷支护设计参数预测. 能源与环保. 2025(02): 235-240 . 百度学术
    3. 李迎富,宋武,杨子仪,赵江梅,曾凡刚. 掘进巷道锚杆支护优化研究. 煤矿机械. 2024(05): 116-118 . 百度学术
    4. 陈万辉,郭瑞,韩伟,宋永明,梁燕翔,刘耀,王佳明,许娜,孟波. 煤矿巷道支护方案智能设计研究. 工矿自动化. 2024(08): 76-83+90 . 百度学术
    5. 彭启友,缪易辰,潘文,尚江淮. 吸能锚杆数值模拟研究现状与展望. 建筑技术. 2024(20): 2546-2551 . 百度学术

    其他类型引用(1)

图(4)  /  表(1)
计量
  • 文章访问数:  15
  • HTML全文浏览量:  1
  • PDF下载量:  3
  • 被引次数: 6
出版历程
  • 收稿日期:  2023-01-12
  • 修回日期:  2023-02-15
  • 刊出日期:  2024-03-19

目录

/

返回文章
返回