高级检索

近距离煤层群协调开采支架工作阻力计算方法与系统

李杨, 任玉琦, 李铁峥, 杨坤鹏, 金向阳

李 杨,任玉琦,李铁峥,等. 近距离煤层群协调开采支架工作阻力计算方法与系统[J]. 煤炭科学技术,2023,51(7):268−277

. DOI: 10.13199/j.cnki.cst.2022-0768
引用本文:

李 杨,任玉琦,李铁峥,等. 近距离煤层群协调开采支架工作阻力计算方法与系统[J]. 煤炭科学技术,2023,51(7):268−277

. DOI: 10.13199/j.cnki.cst.2022-0768
LI Yang,REN Yuqi,LI Tiezheng,et al. Calculation and system of support resistance of shield for contugous-multiple coal seams with coordinated mining[J]. Coal Science and Technology,2023,51(7):268−277. DOI: 10.13199/j.cnki.cst.2022-0768
Citation: LI Yang,REN Yuqi,LI Tiezheng,et al. Calculation and system of support resistance of shield for contugous-multiple coal seams with coordinated mining[J]. Coal Science and Technology,2023,51(7):268−277. DOI: 10.13199/j.cnki.cst.2022-0768

近距离煤层群协调开采支架工作阻力计算方法与系统

基金项目: 

国家自然科学基金资助项目(52074293);河北省自然科学基金资助项目(E2020402041);中央高校基本科研业务费研究生科研创新能力提升资助项目(2022YJSNY07)

详细信息
    作者简介:

    李杨: (1982—),男,河北唐山人,教授,博士生导师。E-mail:liyangcumtb@163.com

    通讯作者:

    任玉琦: (1993—),男,山西大同人,博士研究生。E-mail:ryq2019cumtb@126.com

  • 中图分类号: TD82

Calculation and system of support resistance of shield for contugous-multiple coal seams with coordinated mining

Funds: 

National Natural Science Foundation of China (52074293); Natural Science Foundation of Hebei Province (E2020402041); Basic Research Funds for Central Universities Graduate Research Innovation Ability Enhancement Project (2022YJSNY07)

  • 摘要:

    近距离煤层群协调开采模式下,由于煤层间的开采扰动影响,各煤层开采顶板结构特征不尽相同,因此各煤层开采支架工作阻力的计算方法也存在差异。为针对性地为各煤层开采支架额定工作阻力确定提供思路,综合应用理论分析、系统开发及现场实测方法,对各煤层开采支架工作阻力的计算方法进行了研究。主要成果有:①建立了煤层开采顶板砌体梁式、散体给定式及散体给定−砌体梁式平衡结构模型。砌体梁式平衡结构适用于未受开采扰动影响或受到的开采扰动影响较小的煤层开采;散体给定载荷平衡结构适用于煤层顶板岩层为单一岩层,且受到上部煤层开采影响的煤层开采;散体给定载荷−砌体梁式平衡结构适用于顶板岩层为多层岩层,且层间存在厚度较大,岩性较硬的岩层,同时受到上煤层开采影响岩层仍能保持连续性与完整性的煤层开采;②开发了适用于开滦集团“近距离煤层群开采覆岩破断及支架载荷评价系统”,提出了各煤层支架额定工作阻力的建议选择结果。通过各煤层开采支架压力实测,采用支架额定工作阻力经验选择结果后,各煤层开采支架载荷利用率普遍偏低,支架载荷富余量偏大。而采用各煤层支架额定工作阻力的建议选择结果后,各煤层开采支架载荷利用率明显提高。

    Abstract:

    Under the coordinated mining mode of close-multiple coal seams, due to the mining influence between coal seams, the roof structural characteristics are different after each coal seam extraction, so the calculation methods of support resistance shield of each coal seam are also different. In order to provide ideas for the setting load of shield determination in each coal seam, the calculation methods of the support capacity of shield in each coal seam is given by comprehensive use of theoretical analysis, system development and field measurement. The results show that: ①The voussoir beam, given load of loose body and voussoir beam with given load of loose body balance roof structure models after each coal seam extraction are established. Voussoir beam balance roof structure model is applicable to the coal seams extractions that are not affected by mining or are less affected by mining. Given load of loose body balance roof structure model is applicable to the coal seams extractions with a single roof stratum and is affected by the upper coal seam extraction. Voussoir beam with given load of loose body balance roof structure model is applicable to the coal seams extractions with multi-rock strata and within has a thick and hard lithology. At the same time, affected by the extraction of the upper coal seam, the rock stratum can still maintain continuity and integrity. ②The “overburden breaking and load evaluation system for close-multiple coal seams extraction” suitable for Kailuan Group is developed, and the recommended selection results of setting load of shield in each coal seam are put forward. Through the field measurement of support capacity of shield, the load utilization rate of shield in each coal seam is generally low and the load margin of shield is too large after using the empirical selection results of the setting load of shield. After adopting the recommended selection results of the setting load of shield in each coal seam, the load utilization rate of shield in each coal seam is significantly improved and the load margin of shield is significantly reduced.

  • 深部煤层气资源是中国煤层气产业规模性发展的重要基础,埋深小于2 000 m的煤层气地质资源量为29.82×1012 m3,其中埋深大于1 000 m的深部煤层气资源量为18.71×1012 m3,占比63%[1-2] ,资源潜力巨大。而当前我国煤层气开发的深度大多位于1 000 m 以浅[3-4],深部煤层气勘探开发研究较薄弱。随着浅部已探明可动用储量的减少,深部煤层气开采技术将成为非常规天然气勘探开发的一个新领域,具有重要战略意义[5]。目前对深部煤层气储层地质及关键影响因素分析[6-11]、资源特征及开发潜力[12-13]、富集成藏效应[14-15]、开发工艺及新理论的探索[2,16-17]、生产特征[12,18]已多有研究,同时,在鄂尔多斯、沁水等盆地也启动了试采工程并取得突破[19-20]。一系列的探索与成果,昭示着中国深部煤层气商业性开发具有可观前景,尤其是随着吉深6-7平01井达到日产气量10.1×104 m3 高产工业气流后,更是颠覆了对深部煤层气气勘探开发的许多传统认识[21],对其他区块深部煤层气的开发具有指导性意义。

    横岭区块前期煤层气资源评估显示区内煤层厚度大,含气量高,煤层气开发潜力大。笔者利用地震相控反演开展沁水盆地横岭区块储层空间分布、含气性、煤层厚度、煤层气有利区的精细预测,综合分析物性参数及围岩特征等,在研究区划分出4个煤层气有利区,为后续的煤层气勘探开发提供有力依据。

    沁水盆地目前沉积盖层主要构造线呈北北东向展布,在南、北边缘转折端部位,受边界构造的影响,构造线方向偏转为北东东向或近东西向。按断块学说,和顺横岭区块所在的大地构造Ⅱ级区划为华北断块,Ⅲ级区划为吕梁太行断块,Ⅳ级区划为沁水块坳。位于沁水盆地复向斜东翼的横岭区块,地层走向沿北北东方向,总体由一个背斜和一个向斜构成,在此基础上伴有次一级宽缓褶曲,一般倾角为3°~10°,局部因构造影响可达16°左右。区块因多期构造应力叠加作用,形成挤压变形与伸展拉张变形共存的格局。区内褶皱、断层发育,共发育有向背斜7个,断层7条,以及陷落柱2个;褶皱多呈NNE向,向斜与背斜紧邻交替出现,褶皱延伸从几千米到几十千米,具体如图1所示。

    图  1  沁水盆地横岭区块构造纲要
    Figure  1.  Structure outline map of Hengling block in Qinshui Basin

    横岭区块开发层系为石炭系上统太原组15号煤层,埋深在1 352~1 921 m,垂深多大于1 500 m,属于典型深部煤层气开采;基于岩芯观察,15号煤以光亮型−半亮型煤为主,半暗煤为辅;煤岩成分主要为亮煤,其次为暗煤,煤岩具玻璃光泽,未见裂隙,棱角状断口。试验分析资料表明煤岩演化程度较高,镜质组反射率Ro介于2.0%~3.0%,处于贫煤−无烟煤阶段,多为无烟煤;含气量15~33 m3/t,平均23 m3/t,整体具备良好的生气条件。

    地震波阻抗反演技术是储层反演和预测的重要技术之一。在实际生产中应用的地震反演技术多种多样[22],笔者采用了地震相控非线性随机反演技术,在充分吸取宽带约束反演与模型法反演优点的同时,将标准化或重构之后的测井资料与地震信息有机结合,采用非线性最优化理论、随机模拟算法等,保证了反演结果具有明确的地质意义又有较高的纵向分辨率和良好的预测性。

    利用本区块钻井资料、测井资料编制多口井层序划分与地层界面解释对比图。在对应的连井地震剖面图上解释出层序界面,建立层序或相控模型,进而可以在平面上和三维空间勾画出目的层等不同层序间的匹配关系,为地震相控约束反演奠定约束条件。地震剖面层序解释示意如图2所示。

    图  2  地震剖面层序解释示意
    Figure  2.  Sequence interpretation diagram of seismic section

    考虑地下地质的随机性,相控外推计算中采用多项式相位时间拟合方法建立道间外推关系。在相界面控制的时窗范围内从井出发,将测井资料得到的先验模型参数向量或井旁道反演出的模型参数向量,沿多项式拟合出的相位变化方向进行外推,参与下一地震道的约束反演[23]

    N为给定的正整数,给定数值${f}\left({-N}\right) {f}\left({-N+1}\right){,…,f}\left({N}\right)$,2N多项式拟合数据${f}{(}{x}{)}$,有:

    $$ f\left(x\right)={c}_{0}{p}_{0}\left(x\right)+{c}_{1}{p}_{1}\left(x\right)+ … {+c}_{{\rm{n}}}{p}_{{\rm{n}}}\left(x\right) $$ (1)

    这里每个${p}_{{{j}}}\left(x\right)({{j}}=\mathrm{0,1},2,\cdots ,n)$为xj次多项式,且满足:

    $$ \left\{\begin{array}{l}{p}_{0}\left(x\right)=1\\ \displaystyle\sum {p}_{{\rm{k}}}\left(x\right){p}_{{\rm{m}}}\left(x\right)=0\end{array}\right. $$ (2)

    ${p}_{{{k}}}\left(x\right)$与${p}_{{{m}}}\left(x\right)(k\ne m)$相互正交。由$ {p}_{0}\left(x\right)=1 $可以递推出全部的${p}_{{{j}}}\left(x\right)(j>0)$,用3次多项式拟合得到:

    $$ {c}_{{{k}}}=\sum _{k=-N}^{N}{p}_{{{k}}}\left(x\right)f\left(x\right)/\sum _{k=-N}^{N}{p}_{{{k}}}^{2}\left(x\right)({{k}}=\mathrm{0,1},2,\cdots ,n) $$ (3)

    利用变差函数来描述空间数据场中数据之间的相互关系,建立空间储层参数点之间的统计相关函数。区域化变量Z(x)在x和(x+h)两点处的增量的平方累加起来再除以2倍的t,得到的以两点间距h为变量的函数值为

    $$ G\left(h\right)=\frac{1}{2t}\sum _{i=1}^{t}{\left[{Z}_{{{i}}}\left(x\right)-{Z}_{{{i}}}(x+h)\right]}^{2} $$ (4)

    基于地震道非线性最优化反演的思想,将地震道与波阻抗关系的目标函数定义为式(5),即求解目标函数在最小二乘意义下的极小值,若假设岩石密度为常数,则波阻抗反演变换为速度反演。

    $$ g\left(v\right)=\sum _{i=0}^{n-1}{\left({S}_{{{i}}}^\Delta -D_{{i}} \right)}^{2}\to \mathrm{m}\mathrm{i}\mathrm{n} $$ (5)

    式中:v为速度,m/s;${S}_{{{i}}}^\Delta$为模型响应,即速度预测结果对应的合成地震记录,由地震子波与反射系数褶积得到;${D}_{{{i}}}$为实际地震记录,i为地震记录的采样点序号。

    根据Cook的广义线性反演思想,用Taylor公式将$\left({S}_{{{i}}}^\Delta -{D}_{{{i}}}\right)$在初始模型响应${S}_{{{i}}}$处展开:

    $$ {S}_{{{i}}}^\Delta -{{D}}_{{i}}={{S}}_{{i}}-{{D}}_{{i}}+\sum _{{{{k}}=0}}^{{n-}{1}}\Delta {v}_{{{{k}}}}\frac{\partial{{S}}_{{i}}}{\partial {{V}}_{{k}}}+\frac{{1}}{{2}}\Delta {v}_{{{{k}}}}^{{2}}\left(\sum _{{{{k}}=0}}^{{n-}{1}}\frac{{\partial}^{{2}}{{S}}_{{i}}}{\partial{v}_{{k}}\partial{v}_{{j}}}\right)+ \cdots $$ (6)

    式中:$ {{S}}_{{i}} $为速度初始模型对应的合成地震记录,$\Delta v$为模型参数摄动量。保留二次项,将以上的高次项略掉,即:

    $$ {{S}}_{{i}}^\Delta -{{D}}_{{i}}={{S}}_{{i}}-{{D}}_{{i}}+\sum _{{k=0}}^{{n-1}}\Delta {v}_{{k}}\frac{\partial {{S}}_{{i}}}{\partial {v}_{{k}}}+\frac{{1}}{{2}}\Delta {v}_{{k}}^{{2}}\left(\sum _{{k=}{0}}^{{n-}{1}}\frac{{\partial }^{{2}}{{S}}_{{i}}}{\partial {v}_{{k}}\partial {v}_{{j}}}\right) $$ (7)

    将式(7)对Δv求一阶导数,可得:

    $$ \begin{aligned} &\qquad \quad \frac{\partial {S}_{{{i}}}^{\Delta }}{\partial \Delta {v}_{{{j}}}}=\frac{\partial {S}_{{{i}}}}{\partial {v}_{{{j}}}}+\sum _{{{k}}=0}^{n-1}\Delta {v}_{{{k}}}\frac{{\partial }^{2}{S}_{{{i}}}}{\partial {v}_{{{j}}}\partial {v}_{{{k}}}} \\&\qquad ({{i}}=\mathrm{0,1},\cdots ,n-1;{{j}}=\mathrm{0,1},\cdots ,n-1) \end{aligned}$$ (8)

    将式(8)右端对$\Delta v$求一阶导数,并令该导数为0,可得:

    $$ \sum _{{{i}}=0}^{n-1}\frac{{\partial ({S}_{{{i}}}^{\Delta }-{D}_{{{i}}})}^{2}}{\partial \Delta {v}_{{{j}}}}=2\sum _{{{i}}=0}^{n-1}\left({S}_{{{i}}}^{\Delta }-{D}_{{{i}}}\right)\frac{\partial {S}_{{{i}}}^{\Delta }}{\partial \Delta {v}_{{{j}}}}=0 $$ (9)

    将式(7)和式(8)代入式(9),则有:

    $$ \begin{array}{c} 2 \displaystyle\sum\limits_{{{i}} = 0}^{n - 1} {\left\{ {\left[ {{S_{{i}}} - {D_{{i}}} + \displaystyle\sum\limits_{{{k}} = 0}^{n - 1} {\Delta {v_{{k}}}} \dfrac{{\partial {S_{{i}}}}}{{\partial {v_{{j}}}}} + \dfrac{1}{2}\Delta v_{{k}}^2\left( {\displaystyle\sum\limits_{{{k}} = 0}^{n - 1} {\dfrac{{{\partial ^2}{S_{{i}}}}}{{\partial {v_{{k}}}\partial {v_{{j}}}}}} } \right)} \right]} \right.} \\ \left. {\left( {\dfrac{{\partial {S_{{i}}}}}{{\partial {v_{{j}}}}} + \displaystyle\sum\limits_{{{k}} = 0}^{n - 1} \Delta {v_{{k}}}\dfrac{{{\partial ^2}{S_{{i}}}}}{{\partial {v_{{j}}}\partial {v_{{k}}}}}} \right)} \right\} = 0 \\[-20pt] \end{array}$$ (10)

    将式(10)两边除以2,再将左端展开后,省略部分高阶极小量简化为:

    $$\begin{array}{c} \displaystyle\sum\limits_{{{i}} = 0}^{n - 1} {\left[ {\left( {{S_{{i}}} - {D_{{i}}}} \right)\dfrac{{\partial {S_{{i}}}}}{{\partial {v_{{j}}}}} + \dfrac{{\partial {S_{{i}}}}}{{\partial {v_{{j}}}}}\displaystyle\sum\limits_{{{k}} = 0}^{n - 1} \Delta {v_{{k}}}\dfrac{{\partial {S_{{i}}}}}{{\partial {v_{{k}}}}} + } \right.} \\ \left. {\left( {{S_{{i}}} - {D_{{i}}}} \right)\displaystyle\sum\limits_{{{k}} = 0}^{n - 1} \Delta {v_{{k}}}\dfrac{{{\partial ^2}{S_{{i}}}}}{{\partial {v_{{j}}}\partial {v_{{k}}}}}} \right] = 0 \end{array} $$ (11)

    为便于理解,将式(11)用简单形式表示为:

    $$ A\Delta v+B\Delta v+C=0 $$ (12)

    其中:

    $$ \begin{aligned} A= & \sum _{i=0}^{n-1}\left[\frac{\partial {S}_{{{i}}}}{\partial {v}_{{{j}}}}\sum _{{{k}}=0}^{n-1}\frac{\partial {S}_{{{i}}}}{\partial {v}_{{{k}}}}\right]; B = \sum _{{{i}}=0}^{n-1}\left[\left({S}_{{{i}}}-{D}_{{{i}}}\right)\sum _{{{k}}=0}^{n-1}\frac{{\partial }^{2}{S}_{{{i}}}}{\partial {v}_{{{j}}}\partial {v}_{{{k}}}}\right] ; \\ C = & \sum _{{{i}}=0}^{n-1}\left[\left({S}_{{{i}}}-{D}_{{{i}}}\right)\frac{\partial {S}_{{{i}}}}{\partial {v}_{{{k}}}}\right] \\[-21pt] \end{aligned}$$ (13)

    利用式(12)求取模型摄动量$ \Delta v $时,一般多采用矩阵求逆的方法,这样很容易因矩阵奇异而无解。为此,笔者将矩阵求逆蜕变为一元一次方程求解来减少反演的多解性,增强其稳定性。

    由式(12)求出$ \Delta v $,通过式(14)迭代得到最终的反演速度v

    $$ {v}^{a+1}={v}^{a}+\Delta {v}^{a} $$ (14)

    式中:a为迭代次数。

    利用HL-T-01井参与叠后反演,利用HL-T-02、HL-T-03、HL-T-04、HL-T-06、HL-T-09、ZK0101和ZK0102井进行验证。另外,由于本区块用了以前的测线(OL2和OD1)的部分结果,为了保证不同工区的闭合性,OD1与SN02的交点、OL2与EW02的交点、OL2与SN03的重合部分的结果也用来约束本区块的反演。

    从反演成果中储层特征是否明显、反演结果与钻井结果是否吻合、反演结果是否稳定、其预测性是否强以及反演结果是否充分反映区块目的层地质、沉积特征等方面来衡量反演效果,本次采用了先进的相控非线性随机反演方法,因此区块目的层储层反演取得了较好的效果。

    1)反演分辨率高,反演剖面与测井数据吻合。

    图3是HL-T-01钻孔所在SN03测线与测井资料对比分析图,从图3可以看出,剖面储层特征明显,低速煤层(深蓝色)、高速石灰岩(黄褐色)、中高速度的砂(黄色)和泥页岩(浅蓝色)都能分辨出来,反演剖面对岩性识别效果非常好。同时,反演剖面在HL-T-01钻孔上的岩性与钻井资料对应较好。

    图  3  HL-T-01井所在SN03测线反演剖面与测井对比分析图
    Figure  3.  Comparative analysis of inversion profile and logging of sn03 logging line in HL-T-01 well

    2)反演剖面与地震波形特征一致。

    图4为区块EW06测线的反演速度剖面与地震剖面的叠合图,从4图中可看出反演与地震波形保持一致,储层岩性的横向细微变化与地震剖面上同相轴及其振幅的变化是同步的,地震特征变了,岩性随之发生变化,这表明反演结果中岩性、储层的横向变化主要与地震信息(属性)的变化有关,随地震属性(如振幅等)的变化而变化,说明反演结果忠实于地震资料,反演方法稳定。

    图  4  EW06测线反演剖面与地震剖面叠合图
    Figure  4.  Superposition of EW06 line inversion section and seismic section

    3)沉积特征明显。

    由于本次反演融入了相控反演的思想,使得反演结果沉积特征清楚。图5为EW02测线的反演剖面图,剖面上的地层跨过断层出现缺失,符合地质特征。剖面上的砂体岩性超复、尖灭、分叉以及相变特征非常明显,符合地质特征。

    图  5  EW02测线叠后反演剖面
    Figure  5.  Post stack inversion profile of EW02 survey line

    在反演速度剖面上,以地层对比为基础,对应区块内现有钻孔的储层数据进行标定,确定储层追踪解释目标和层数。在区块内共确定了3个煤储层小层:① 经过K10灰岩和K8灰岩以及第一段泥岩后是3号煤储层;② 过3号煤底-K4石灰岩顶的第二段泥页岩后是8号煤;③ 8号煤底经过K4等灰岩层后是15号煤层。钻孔小层位置示意如图6所示。

    图  6  HL-T-01钻孔小层位置示意
    Figure  6.  Location of borehole sublayer of HL-T-01 well

    根据测井分层资料确定钻孔上储层位置,继而从钻孔出发向两侧追踪解释。如图7所示,反演剖面上储层位置与测井一致。

    图  7  小层位置在剖面上的展示图
    Figure  7.  Display of sublayer position on section

    追踪完小层后,将各储层的顶和底相减得到的时间厚度值,进行剖面间插值,并将其与层速度相乘,得到深度域的厚度平面等值线图。

    利用叠后地震相控非线性随机反演技术预测了煤储层(3号、8号、15号煤层)的厚度分布。结果显示:8号和3号煤储层厚度分布较不均匀,厚度范围为0.5~2.5 m,沉积的稳定性较差;15号煤层整体分布为NNE方向,全区厚度分布稳定,厚度变化范围在4~7 m(图8)。煤层最厚处位于西北部,最薄处位于中西部。

    图  8  15号煤层厚度预测平面图
    Figure  8.  Coal seam thickness plan of 15# coal

    图9为EW05测线的叠后反演剖面与叠后含气性预测剖面的对比分析,上图为EW05测线的叠后反演剖面,下图为该测线的叠后含气性预测剖面。上图中有煤层发育的地方在下图中不一定能看到含气的存在,可见叠后反演可以有效识别储层分布,但无法反映储层中是否含气,上图中速度特低区在下图有明显的含气性异常。可见将叠后含气性检测与叠后速度反演交互,能够有效识别含气储层。除此之外,基于流体活动性属性分析技术从追踪的小层中提取储层的含气量,采用含气厚度的概念宏观表征储层含气性及含气量的相对大小,即厚度=速度×时间,具体是:首先根据钻孔测井曲线得到煤层的速度、拉梅系数乘密度这2个约束条件,优选含气煤层,再根据含气煤层采样点的速度乘以时间间隔,累计起来就是含气厚度。

    图  9  测线EW05叠后含气性检测剖面与叠后反演剖面对比
    Figure  9.  Comparison between post stack gas bearing detection section and post stack inversion section of survey line EW05

    叠后含气性预测与叠后反演剖面结合得到的15号煤层含气性预测平面分布图如图10所示,由图10可以看出,15号煤层在区块中部和北部含气特征明显,厚度最高可达4.5 m。区块内整体含气性较好。

    图  10  15号煤层叠后含气性预测平面分布
    Figure  10.  Plane distribution of post stack gas bearing prediction of No.15 coal seam

    在圈定目标区的含气有利区的工作中,对于目标区的精细勘探是必不可少的,因此对目标区的储层物性的研究显得尤为重要。研究过程中发现,区块目的层的孔隙度较低,在低孔隙的背景下寻找相对高孔隙部位对分析含气成藏机理寻找有利成藏区域以及后续的井位设计是很重要的一个环节。孔隙度反演的目的就是寻找到储层发育的高孔隙区分布。

    由于区块缺少预测孔隙度的资料,邻区榆社-武乡区块仅有十几个泥页岩岩心的试验结果,所以只能用其做拟合来反演整个区域的孔隙度。

    通过分析泥页岩孔隙度$\varphi$与速度${V}_{{\rm{p}}}$的关系,得到交会图(图11),从拟合关系可以得到如下拟合公式:

    图  11  泥页岩孔隙度与速度交会图
    Figure  11.  Cross plot of shale porosity and velocity
    $$\varphi =-0.009\,234 {V}_{{\rm{p}}}+37.88 $$ (15)

    利用式(15)对泥页岩孔隙度进行反演,因勘探区无煤层孔隙度样品数据,因此调研了勘探区附近沁水煤田其他区块的孔隙度数据。因煤层孔隙度数据比泥页岩孔隙度大,所以在采用泥页岩公式的基础上对其进行改动,得到煤层孔隙度反演公式为:

    $$ \varphi =-0.009\,234 {V}_{{\rm{p}}}+39.88 $$ (16)

    Gardner在1974年给出密度与速度的经验公式[24]

    $$ \rho =0.31{V}_{{\rm{p}}}^{0.25} $$ (17)

    进一步获得:

    $$ \varphi =-0.009\,234 \sqrt[0.25]{\rho /0.31} $$ (18)

    本次预测了3号、8号、15号煤储层的孔隙度分布。图12为15号煤层孔隙度预测平面分布图,通过15号煤层孔隙度与煤层厚度、煤层含气性对比可知:15号煤孔隙度分布整体较均匀,范围在4%~6%。煤层厚的区域孔隙度较大,含气高部位对应的孔隙度发育较好。

    图  12  15号煤层孔隙度分布
    Figure  12.  Porosity distribution of No.15 coal seam

    采用多尺度边缘检测技术实现对裂缝特征的研究,即利用图像边缘检测理论,对复杂地质体引起地震同相轴横向错断及地震振幅的突变能够有效突显,识别地震信号振幅纵横向强弱变化的边界,实现对小尺度裂缝的刻画识别。与地震相干体相比,多尺度边缘检测能够利用多尺度函数检测算子,不仅对大断裂等大尺度地质异常信号能很好地识别,而且对很多小尺度断裂及裂缝引起细微变化信息进行有效监测,极大的提高了地震裂缝的分辨率。

    1)二维地震小波变换

    通过将三维地震数据的水平切片表示,对$ G(x,y) $进行二维小波变换,首先设二维小波${\varPsi }_{{\rm{s}}}^{\left(1\right)}$(x,y)${\varPsi }_{{\rm{s}}}^{\left(2\right)}$(x,y)在尺度${{ s}}$下的伸缩为:

    $$ \left\{\begin{array}{c}{\varPsi }_{{\rm{s}}}^{\left(1\right)}(x,y)=\dfrac{1}{{\rm{s}}}{\varPsi }^{\left(1\right)}\left(\dfrac{x}{{\rm{s}}},\dfrac{y}{{\rm{s}}}\right) \\ {\varPsi }_{{\rm{s}}}^{\left(2\right)}(x,y)=\dfrac{1}{{\rm{s}}}{\varPsi }^{\left(2\right)}\left(\dfrac{x}{{\rm{s}}},\dfrac{y}{{\rm{s}}}\right)\end{array}\right. $$ (19)

    那么水平切片$ G(x,y) $的二维小波变换由式(20)定义:

    $$ \left\{\begin{array}{c}{\varPsi }_{\rm{s}}^{\left(1\right)}(x,y)=G(x,y){\varPsi }_{\rm{s}}^{\left(1\right)}(x,y)\\ {\varPsi }_{\rm{s}}^{\left(2\right)}(x,y)=G(x,y){\varPsi }_{\rm{s}}^{\left(2\right)}(x,y)\end{array}\right. $$ (20)

    二维小波变换思路在一维小波变化基础上进行,具体形式如式(20)。再进行二进离散化,就可以获得如下形式:

    $$ \left\{\begin{array}{c}{W}_{{2}'}^{\left(1\right)}G(x,y)=G(x,y){\psi }_{{2}'}^{\left(1\right)}(x,y)\\ {W}_{{2}'}^{\left(2\right)}G(x,y)=G(x,y){\psi }_{{2}'}^{\left(2\right)}(x,y)\end{array}\right. $$ (21)

    式(21)表明,$ G(x,y) $的边缘点就是由${W}_{{\rm{s}}}^{x}G(x,y)$与${W}_{{\rm{s}}}^{y}G(x,y)$的模同时取极大的点确定,而边缘方向实际上是给定点$ (x,y) $所在曲面$G(x,y)*{\theta }_{{\rm{s}}}(x,y)$的梯度方向。梯度向量的幅值为:

    $$ {M}_{{\rm{s}}}f(x,y)=\sqrt{{\left|{W}_{{\rm{s}}}^{x}G(x,y)\right|}^{2}+{\left|{W}_{{\rm{s}}}^{y}G(x,y)\right|}^{2}} $$ (22)

    幅角为:

    $$ \alpha ={A}_{{\rm{s}}}G(x,y)=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\frac{\left|{W}_{{\rm{s}}}^{x}G(x,y)\right|}{\left|{W}_{{\rm{s}}}^{y}G(x,y)\right|} $$ (23)

    式中:${W}_{{\rm{s}}}^{x}G(x,y)$与${W}_{{\rm{s}}}^{y}G(x,y)$分别为小波变换后沿x方向的梯度和y方向的梯度;${\theta }_{{\rm{s}}}(x,y)$为小波函数;${A}_{{\rm{s}}}$为一个整体,代表梯度方向。

    2)边缘检测正演成像试验

    地震裂缝识别与图像裂缝识别的原理一样,图像(图13)的每一个像素点都对应一个数值,黑色代表小值,白色代表大值。裂缝检测就是查找某个点与周围点的差异。如果某个点的值与周围点的值差异大于某个阈值,而且相邻点也具有类似特征,且这些点呈线性分布,就认为是一个裂缝。地震检测也是这样,首先把地震数据提取跟断层敏感的属性,比如相干、曲率、方差等属性。空间上每个点代表一个属性值,就可以按照图片中裂缝识别的原理进行地震裂缝预测。

    图  13  多尺度裂缝检测正演成像
    Figure  13.  Multiscale fracture detection forward imaging

    通过对裂缝图片进行裂缝检测正演试验(图13),多尺度边缘检测纵向检测、横向检测,分别能够很好地检测出照片中纵向、横向展布的裂缝;同时利用全方位检测,能更加全面地将纵向和横向的裂缝都检测出来,是一种行之有效的裂缝检测方法。

    横岭区块15号煤层裂隙主要分布在区块西部,断层发育的破碎带裂隙发育好。裂隙密度随着深度增加逐渐减少,煤层裂隙密度平面分布图如图14所示。

    图  14  15号煤层裂隙密度平面分布
    Figure  14.  Plane distribution of fracture density of No.15 coal seam

    前人研究结果显示,煤层埋藏深度、煤层顶板5 m内砂岩厚度和顶板岩性是影响勘探区煤层含气量的3个关键因素[25]。因此,本次研究利用地震反演获得的速度数据体,统计煤层上、下各4 m范围内的速度,用以区分煤层直接顶、底板的岩性。岩石物理分析显示,砂泥岩、砂质泥岩速度部分重叠,但是平均速度和趋势有明显差异,本次采用的是层平均速度来区分的。

    从顶底板岩性分布图(图15图16)可以看出,15号煤层直接顶板主要是砂质泥岩,尤其是区块北部,主要为泥岩,说明顶板对煤层气的封闭能力较强;直接底板主要是砂质泥岩封盖,局部地区为泥岩和砂岩,砂岩分布连片性较差,对煤层气逸散影响不大。总体来说,顶底板对煤层气封闭能力较强。

    图  15  15号煤层直接顶板岩性
    Figure  15.  Lithology of direct roof of No.15 coal seam
    图  16  15号煤层直接底板岩性
    Figure  16.  Lithology of direct floor of No.15 coal seam

    煤层气一般富集在含煤区构造高点、上斜坡区以及背斜向斜相间的褶皱区等构造有利部位,适量的构造裂隙可增加煤层中的裂隙,但对于构造非常复杂的地区,特别是张性断层的存在,将成为煤层气逸散的通道,不利于煤层气的保存,在选择煤层气有利区时,一般应避开构造特别复杂的部位。本次有利区的划分首先是避开构造特别复杂的部位,选择背斜向斜相间的褶皱区等构造有利部位,在此基础上选择储层厚度大、含气显示良好、孔隙度较高的区域作为有利区。

    根据15号煤层含气性分布及构造图,结合15号煤层煤层厚度和孔隙度图,15号煤层共圈定了含气有利区4块(图17),总面积51.41 km2,有利区内气层厚度大,孔隙度发育好。

    图  17  15号煤层有利区分布示意
    Figure  17.  Distribution diagram of favorable areas of No.15 coal seam

    有利区Ⅰ主要分布于区块的北部,SN01线以西,位于DS1向斜的西翼,面积约8.15 km2;有利区Ⅱ主要分布于区块的中北部,EW05线以北,SN01线以东,位于背斜和向斜相间的褶皱区,面积约10.78 km2;有利区Ⅲ分布于区块的西部,EW02线和EW04线之间,位于背斜和向斜相间的褶皱区,面积约15.56 km2;有利区Ⅳ分布于区块的东南部,EW01线和EW03线之间,位于背斜和向斜相间的褶皱区,面积约16.92 km2。除此之外,结合确定的8号、3号煤层气有利区分布情况,以储层厚度分布在6~10 m,气层厚度分布在3~8 m,孔隙度4%~10%之间为标准,圈定了煤层气综合有利区3处。有利区内煤层气有利区Ⅰ主要分布于区块西北部,EW06线以南、SN01线以西,位于DS1向斜的西翼,面积约3.34 km2;煤层气有利区Ⅱ主要分布于区块中西部边界附近,EW02线和EW03线之间,位于DS4背斜的两翼,面积约6.40 km2;煤层气有利区Ⅲ主要分布于区块东南部,位于DS5向斜、DS6背斜、DS7向斜相间的褶皱区,面积约17.27 km2。总面积27.01 km2图18)。

    图  18  煤层气综合有利区分布示意
    Figure  18.  Distribution diagram of favorable areas for coalbed methane

    1)相控非线性随机反演的方法在目的储层反演中取得效果显著,反演分辨率高,与测井吻合,反演剖面对岩性识别效果明显,且与地震波形特征一致,反演结果沉积特征清楚,符合地质特征。在反演速度剖面上,以地层对比为基础确定了3个煤储层小层。

    2)15号煤层整体分布为NNE方向,全区分布稳定,厚度范围为4~7 m;15号煤层在区块中部和北部含气特征明显,厚度最高可达4.5 m;顶底板岩性主要为砂质泥岩、泥岩,煤层气封闭能力较强。

    3)结合叠后反演和叠前反演,对煤层孔隙度进行了预测,15号煤储层的孔隙度分布在4%~6%。通过厚度,含气性和孔隙度图的对比,孔隙度高的地方煤层的厚度及含气性相对较高。

    4)在构造有利部位,综合考虑储层厚度、含气性、孔隙度等因素,圈定了15号煤层气有利区4块,总面积51.41 km2。根据储层厚度分布在6~10 m,气层厚度大,分布在3~8 m,孔隙度4%~8%,确定了煤储层综合有利区4处,总面积41.99 km2

  • 图  1   煤层群赋存示意

    Figure  1.   Occurrence of multiple coal seams

    图  2   煤系地层岩性与厚度

    Figure  2.   Thickness and lithology of coal-bearing strata

    图  3   砌体梁式平衡结构

    Figure  3.   Voussoir beam balanced structure

    图  4   散体给定载荷平衡结构

    Figure  4.   Balance structure with given load of loose body

    图  5   散体拱给定载荷平衡结构简化示意

    Figure  5.   Simplified diagram of given load of loose body.

    图  6   散体给定载荷−砌体梁式平衡结构

    Figure  6.   Voussoir beam balanced structure with given load of loose body

    图  7   近距离煤层群开采覆岩破断及支架载荷评价系统

    Figure  7.   Roof fracture and shield load evaluation system for close-multiple coal seams extraction

    图  8   各煤层支架最大工作阻力计算结果

    Figure  8.   Maximum calculation results of each coal seams extraction

    图  9   9号煤层支架工作阻力影响因素分析

    Figure  9.   Analysis on support capacity of shield in No. 9 coal seam with different influencing factors

    图  10   敏感性分析结果

    Figure  10.   Sensitivity analysis results

    图  11   各煤层支架压力实测结果

    Figure  11.   Field measurement of support capacity of shield in each coal seams

    图  12   支架载荷利用率

    Figure  12.   Load utilization of shield

    表  1   各煤层煤厚与间距

    Table  1   Thickness and distance of each coal seam

    煤层厚度/ m 间距/m
    最小最大平均最小最大平均
    505.401.4
    709.804.117.8863.6926.6
    805.121.8017.395.3
    90.195.131.90.6317.395.6
    12−10.308.163.410.2560.9329.5
    下载: 导出CSV

    表  2   顶板结构模型及其适用煤层

    Table  2   Roof structure model and its application in coal seam

    顶板结构适用煤层开采顺序原因
    砌体梁式

    平衡结构
    7第1顶板岩层完整、未受到开采扰动影响
    12−1第3与8号煤层间距较大,顶板岩层开采扰动影响程度小,顶板岩层仍能保持连续性与完整性
    5第5位于煤系地层顶部,与7号煤层间距大,顶板岩层受开采扰动小
    散体给定载荷
    平衡结构
    8第2与7号煤层间距小,层间岩层
    (顶板岩层)只有一层细砂岩,
    受到开采扰动影响
    散体给定载荷砌体
    梁式平衡结构
    9第4与8号煤层间距小,顶板岩层受到8号煤层开采扰动影响,层间具有一层基本顶硬岩,但仍能保持连续性与完整性
    下载: 导出CSV

    表  3   系统开发应用的主要控件

    Table  3   Main controls of system development and application

    序号控件
    英文
    控件
    中文
    介绍
    1Pushbutton按钮鼠标单击实现某种行为、
    调用相应的回调子函数
    2Text文本框控制用户编辑或修改
    字符串的文本域
    3Edit文本信息其他控件的标签
    4Axes坐标轴显示图片
    下载: 导出CSV

    表  4   各煤层支架工作阻力结果

    Table  4   Results of support capacity of shield in each coal seam

    煤层7号煤层8号煤层12−1号煤层9号煤层5号煤层
    最大工作阻力
    计算结果/kN
    3 948.55 4 018.32 4 101.63 3 560.03 4 015.3
    额定工作阻力建议
    选择结果/kN
    4 500 4 300 4 300 4 000 4 300
    额定工作阻力经验
    选择结果/kN
    6 400 4 800 4 800 4 800 5 000
    下载: 导出CSV

    表  5   各煤层实测结果分析

    Table  5   Analysis results of field measurement in each coal seam

    煤层支架型号实测支架压力
    /MPa
    实测支架工作阻力
    /kN
    额定工作阻力
    /kN
    支架载荷利用率/%
    最小值最大值最小值最大值经验
    结果
    建议
    结果
    经验结果建议结果
    7ZY6400-21/4523.126.9 3 714.574 325.63 6 4004 500 58.04~67.5982.55~96.13
    8ZY4800-13/3224.231.52 978.463 876.724 8004 30062.05~80.7769.27~90.16
    12−1ZY4800-19/4020.628.72 535.383 531.314 8004 30052.82~73.5958.96~82.15
    9ZY4800-13/3224.029.52 953.853 630.774 8004 00061.54~75.6473.85~90.77
    5ZY5000-09/20D16.428.42 019.703 497.545 0004 30040.39~69.9546.97~81.34
    下载: 导出CSV
  • [1] 彭赐灯. 矿山压力与岩层控制研究热点最新进展评述[J]. 中国矿业大学学报,2015,44(1):1−8.

    PENG Syd S. Topical areas of research needs in ground control: a state of the art review on coal mine ground control[J]. Journal of China University of Mining & Technology,2015,44(1):1−8.

    [2]

    PENG Syd S. Coal mine ground control (Third edition) [M]. Morgamtown: College of Engineering and Mineral Resources, 2008.

    [3] 王国法. 大采高技术与大采高液压支架的开发研究[J]. 煤矿开采,2009,14(1):1−4.

    WANG Guofa. Research on mining technology with high mining height and development of powered support for high mining height[J]. Coal Mining Technology,2009,14(1):1−4.

    [4] 张可斌,钱鸣高,郑朋强,等. 采场支架围岩关系研究及支架合理额定工作阻力确定[J]. 采矿与安全工程学报,2020,37(2):215−223.

    ZHANG Kebin,QIAN Minggao,ZHENG Pengqiang,et al. Relationship between support and surrounding rocks and determination of reasonable rated working resistance against support[J]. Journal of Mining & Safety Engineering,2020,37(2):215−223.

    [5] 陈炎光, 钱鸣高. 中国煤矿采场围岩控制[M]. 徐州: 中国矿业大学出版社, 1994.
    [6] 史元伟. 采煤工作面围岩控制原理和技术[M]. 徐州: 中国矿业大学出版社, 2003.
    [7] 郝宪杰,许家林. 综采支架工作阻力确定方法综述[J]. 神华科技,2009,7(4):12−16.

    HAO Xianjie,XU Jialin. Summary on determination of the reasonable working resistance of hydraulic support in fully-mechanized face[J]. Energy Science and Technology,2009,7(4):12−16.

    [8] 钱鸣高,缪协兴,何富连,等. 采场支架与围岩耦合作用机理研究[J]. 煤炭学报,1996,21(1):40−44.

    QIAN Minggao,MIAO Xiexing,HE Fulian,et al. Mechanism of coupling effect between supports in the workings and the rocks[J]. Journal of China Coal Society,1996,21(1):40−44.

    [9] 曹胜根,钱鸣高,刘长友,等. 采场支架-围岩关系新研究[J]. 煤炭学报,1998,23(6):17−21.

    CAO Shenggen,QIAN Minggao,LIU Changyou,et al. New research about support and surrounding rock relationship in working face[J]. Journal of China Coal Society,1998,23(6):17−21.

    [10] 王家臣,王 蕾,郭 尧. 基于顶板与煤壁控制的支架阻力的确定[J]. 煤炭学报,2014,39(8):1619−1624.

    WANG Jiachen,WANG Lei,GUO Yao. Determining the support capacity based on roof and coal wall control[J]. Journal of China Coal Society,2014,39(8):1619−1624.

    [11]

    WANG Jiachen,YANG Shengli,LI Yang,et al. A dynamic method to determine the supports capacity in longwall coal mining[J]. International Journal of Mining Reclamation and Environment,2015,29(4):277−288. doi: 10.1080/17480930.2014.891694

    [12] 王国法,庞义辉,李明忠,等. 超大采高工作面液压支架与围岩耦合作用关系[J]. 煤炭学报,2017,42(2):518−526.

    WANG Guofa,PANG Yihui,LI Mingzhong,et al. Hydraulic support and coal wall coupling relationship in ultra large height mining face[J]. Journal of China Coal Society,2017,42(2):518−526.

    [13] 于 雷,闫少宏,刘全明. 特厚煤层综放开采支架工作阻力的确定[J]. 煤炭学报,2012,37(5):737−742.

    YU Lei,YAN Shaohong,LIU Quanming. Determination of support working resistance of top coal caving in extra thick coal seam[J]. Journal of China Coal Society,2012,37(5):737−742.

    [14] 闫少宏,尹希文,许红杰,等. 大采高综采顶板短悬臂梁-铰接岩梁结构与支架工作阻力的确定[J]. 煤炭学报,2011,36(11):1816−1820.

    YAN Shaohong,YIN Xiwen,XU Hongjie,et al. Roof structure of short cantilever-articulated rock beam and calculation of support resistance in full-mechanized face with large mining height[J]. Journal of China Coal Society,2011,36(11):1816−1820.

    [15] 冯国瑞,杨文博,白锦文,等. 非等宽复合柱采区中部遗煤开采可行性分析[J]. 采矿与安全工程学报,2021,38(4):643−654.

    FENG Guorui,YANG Wenbo,BAI Jinwen,et al. Feasibility analysis of residual coal mining in the middle of unequal width composite pillar mining area[J]. Journal of Mining & Safety Engineering,2021,38(4):643−654.

    [16] 王 寅,付兴玉,孔令海,等. 近距离煤层群上行式开采悬空结构稳定性研究[J]. 煤炭科学技术,2020,48(12):95−100.

    WANG Yin,FU Xingyu,KONG Linghai,et al. Study on stability of dangling structure in ascending mining contiguous coal seams[J]. Coal Science and Technology,2020,48(12):95−100.

    [17] 王龙飞,常泽超,杨战标,等. 深井近距离煤层群采空区下回采巷道联合支护技术[J]. 采矿与安全工程学报,2018,35(4):686−692.

    WANG Longfei,CHANG Zechao,YANG Zhanbiao,et al. Combined support technology of roadway under mined gob of ultra-distance seams in deep mine[J]. Journal of Mining & Safety Engineering,2018,35(4):686−692.

    [18] 胡少轩,许兴亮,田素川,等. 近距离煤层协同机理对下层煤巷道位置的优化[J]. 采矿与安全工程学报,2016,33(6):1008−1013.

    HU Shaoxuan,XU Xingliang,TIAN Suchuan,et al. Optimization of roadway location in lower coal seam from synergy mechanism of contiguous seam mining[J]. Journal of Mining & Safety Engineering,2016,33(6):1008−1013.

    [19] 钱鸣高,缪协兴,何富连. 采场“砌体梁”结构的关键块分析[J]. 煤炭学报,1994,19(6):557−563.

    QIAN Minggao,MIAO Xiexing,HE Fulian. Analysis of key block in the structure of voussoir beam in longwall mining[J]. Journal of China Coal Society,1994,19(6):557−563.

    [20] 钱鸣高, 石平五, 许家林. 矿山压力与岩层控制[M]. 徐州: 中国矿业大学出版社, 2003.
    [21] 杨路林,李亚春,吴士良. 近距离煤层采空区下综采支架合理工作阻力研究[J]. 煤炭科学技术,2020,48(9):189−194.

    YANG Lulin,LI Yachun,WU Shiliang. Study on reasonable fully-mechanized mining support working resistance under gob of contiguous coal seam[J]. Coal Science and Technology,2020,48(9):189−194.

    [22] 李 杨,雷明星,郑庆学,等. 近距离“薄-中-厚”交错分布煤层群上行协调开采定量判别研究[J]. 煤炭学报,2019,44(S2):410−418.

    LI Yang,LEI Mingxing,ZHENG Qingxue,et al. Quantitative criterion on coordinated ascending mining in close multiple “thin- medium- thick” coal seams[J]. Journal of China Coal Society,2019,44(S2):410−418.

    [23]

    LI Yang,REN Yuqi,PENG Syd S,et al. Measurement of overburden failure zones in close-multiple coal seams mining[J]. International Journal of Mining Science and Technology,2021,31(8):43−50.

    [24]

    LI Yang,WANG Jiachen,CHEN Yiding,et al. Overlying strata movement with ground penetrating radar detection in close-multiple coal seams mining[J]. International Journal of Distributed Sensor Networks,2019,15(8):1550147719869852.

    [25] 李 杨,任玉琦,王 楠,等. 采空区垮落顶板形态及其演化特征[J]. 煤炭学报,2021,46(12):3771−3780.

    LI Yang,REN Yuqi,WANG Nan,et al. structure form and evolution characteristics of collapsed roof in goaf[J]. Journal of China Coal Society,2021,46(12):3771−3780.

    [26] 薛广哲,任博玲,冯宇峰. 近距离煤层下层煤顶板结构力学分析[J]. 煤矿安全,2011,42(2):143−144,149. doi: 10.13347/j.cnki.mkaq.2011.02.013
    [27] 李红涛,刘长友,汪理全. 上位直接顶“散体拱”结构的形成及失稳演化[J]. 煤炭学报,2008,33(4):378−381.

    LI Hongtao,LIU Changyou,WANG Liquan. Generating and destabilization evolutionary of granular arch structure of upper immediate roof[J]. Journal of China Coal Society,2008,33(4):378−381.

    [28] 朱 涛,张百胜,冯国瑞,等. 极近距离煤层下层煤采场顶板结构与控制[J]. 煤炭学报,2010,35(2):190−193. doi: 10.13225/j.cnki.jccs.2010.02.012

    ZHU Tao,ZHANG Baisheng,FENG Guorui,et al. Roof structure and control in the lower seam mining field in the ultra-close multiple seams[J]. Journal of China Coal Society,2010,35(2):190−193. doi: 10.13225/j.cnki.jccs.2010.02.012

    [29] 中国矿业大学(北京). 近距离煤层群开采覆岩破断及支架载荷评价系统[Z]. 中国: ZL32021SR0857193, 2021.
    [30] 天工在线. 中文版MATLAB2018从入门到精通(实战案例版)[M]. 北京: 中国水利水电出版社, 2018.
  • 期刊类型引用(1)

    1. 杨欢欢,王占辉,李依霖,李彩云,吴艳玲,董玉. 2020—2021年承德市生活饮用水中氟化物检测结果分析. 医学动物防制. 2024(05): 517-521 . 百度学术

    其他类型引用(0)

图(12)  /  表(5)
计量
  • 文章访问数:  157
  • HTML全文浏览量:  34
  • PDF下载量:  91
  • 被引次数: 1
出版历程
  • 收稿日期:  2022-09-10
  • 网络出版日期:  2023-06-16
  • 刊出日期:  2023-07-24

目录

/

返回文章
返回