论华北地区的均衡状态(一)--方法和局部补偿

冯锐1, 王均2, 郑书真2, 黄桂芳2, 严惠芬2, 周海南1, 张若水3

冯锐1, 王均2, 郑书真2, 黄桂芳2, 严惠芬2, 周海南1, 张若水3. 1987: 论华北地区的均衡状态(一)--方法和局部补偿. 地震学报, 9(4): 406-416.
引用本文: 冯锐1, 王均2, 郑书真2, 黄桂芳2, 严惠芬2, 周海南1, 张若水3. 1987: 论华北地区的均衡状态(一)--方法和局部补偿. 地震学报, 9(4): 406-416.
FENG RUIup, WANG JUNup2, ZHENG SHUZHENup2, HUANG GUIFANGup2, YAN HUIFENup2, ZHOU HAINANup, ZHANG RUOSHUIup3dstyle. 1987: ISOSTATIC STATUS IN NORTH CHINA (1)--METHOD AND LOCAL COMPENSATION. Acta Seismologica Sinica, 9(4): 406-416.
Citation: FENG RUIup, WANG JUNup2, ZHENG SHUZHENup2, HUANG GUIFANGup2, YAN HUIFENup2, ZHOU HAINANup, ZHANG RUOSHUIup3dstyle. 1987: ISOSTATIC STATUS IN NORTH CHINA (1)--METHOD AND LOCAL COMPENSATION. Acta Seismologica Sinica, 9(4): 406-416.

论华北地区的均衡状态(一)--方法和局部补偿

ISOSTATIC STATUS IN NORTH CHINA (1)--METHOD AND LOCAL COMPENSATION

  • 摘要: 根据经典均衡的原理,本文重点分析了五种局部均衡补偿模式.在计算方法上利用频域三维重力场的理论公式,快速求得均衡校正值,据此计算了66个模型的均衡重力异常场,注意到对这许多不同参量的模型而言,它们的均衡异常场在形态与分布特征上基本一致,其中计入了地表和上下地壳密度差异分布的 Airy 模型具有最佳的补偿效果,它的均衡面在莫霍界面以下的上地幔中,标准深度50km.从整体上看,华北地台处于亚均衡状态,均衡异常的均值为1810-5m/s2.均衡重力异常的分布表现出明显的块体特征,正均衡异常区主要分布在东部胶辽地块和冀中平原北缘,在汾渭裂谷区存在负异常.模型对比表明,以莫霍界面作为均衡补偿面的模型是不可取的;Airy 模型比 Pratt 模型的补偿效果略好,这同地壳构造以层状为主而侧向变化有限的特征相符.有关复合补偿、均衡重力异常的基本特征和深部构造的关系等结果,将在文章的第二部分发表.
    Abstract: Follwing Airy and Pratt principles, five kinds of local compensation models are analysed with emphasis and a rapid 3-D gravity formula is utilized to calculate isostatic anomaly for sixty-six models with different parameters. It is noted that isostatic gravity maps appear nearly identical in their main patterns and features. The optimum compensation model in North China is one of modified Airy models in which the different density distribution in the surface, upper crust and lower crust are taken into account and the standard crustal thickness is 50 km. The position of complete compensation interface is located in the upper mantle. North China platform as a whole is under the sub-isostatic equilibrium status with isostatic anomaly of 1810-5 m/s2 on an average. The distribution of isostatic gravity anomaly shows an obvious blockwise pattern. Most positive anomaly areas occur over the eastern part, the Jiao-Liao block, Mt. Yan block and northern margin of Hebei-Shandong block, whereas negative area occurs in the Shanxi graben. Comparison of models indicates that the Moho discontinuity is not suitable to be taken as a compensation interface, and the compensation effects in Airy model are better than in Pratt model, which is consistent with the feature of dominantly layered structure with less lateral inhomogeneity in crust. Some results about composite compensation, the basic characteristics of isostatic gravity anomaly and deep structure will be published later in the second paper.
  • 近年来,随着空间对地观测技术的发展,对地震电磁电离层扰动的研究已成为热点,主要集中在对等离子体(电子浓度、离子浓度、电子温度、离子温度等)和低频电磁辐射参量的观测与研究. 然而,围绕地震电离层扰动的研究目前仍存在一些争议,部分原因是目前对强震电离层异常震例资料的积累仍不够充分,尤其是天基资料. 通过搭载在卫星上的测量仪器对电离层扰动区域进行测量是目前电离层变化监测最直接的方法. 许多发达国家已经发射了用于地震监测的电磁观测卫星,我国也计划于2016年底发射第一颗电磁监测试验卫星(申旭辉等,2007).

    电磁卫星载荷性能直接影响电磁观测数据的精准性及其在地震预测应用中的可靠性,但是复杂的电离层环境将会影响卫星载荷的性能. 目前,电磁卫星载荷性能分析多以实验室检测(Chao,Su,1999Klenzing et al,2008)和有限次野外测试为主(Sauvaud et al,2006). 这些分析方法由于受检测条件、人力物力等因素的制约而不能进行高频次的检测和分析,因此考虑对不同卫星搭载同一载荷或者同一卫星搭载同类型载荷的同步探测数据进行对比分析,进而间接地评价载荷性能(Friedel et al,2000王馨悦等,2008王春琴等,2013). 受限于卫星运行周期、在轨时间及覆盖范围,卫星载荷探测数据存在一定的时间差异,会影响探测数据的同步对比分析,因此,基于电磁探测数据的时序分析模型研究,可以解决多源探测数据间的时间匹配问题,有利于开展不同期在轨卫星探测数据的交叉检验分析(李传荣等,2015). 同时,通过时序分析建模建立电离层参量的背景场,可以克服传统的均值法和统计法所建立的背景场无法反映参量随时间变化规律的缺点,有效地分离地震等空间事件所引起的电离层扰动,从而开展基于空间事件同步响应分析的载荷探测数据交叉检验研究.

    本文首先将基于DEMETER卫星郎缪尔探针(instrument sonde de Langmuir,简写为ISL)探测的Ne数据,分析国际参考电离层(international reference ionosphere,简写为IRI)模型在电离层平静期预测Ne数据的相对误差随纬度的变化特征;然后,基于Ne数据明显的季节性变化特征,构建Ne数据的自回归移动平均(autoregressive integrated moving average,简写为ARIMA)模型;最后,对两种不同的时间序列分析预测模型模拟和预测电磁载荷探测数据的相对误差变化特征进行分析.

    本文采用DEMETER卫星ISL探测的Ne数据作为时序分析数据源. 法国DEMETER卫星于2004年6月29日发射,于2010年12月9日停止科学数据接收,可探测全球65°S—65°N范围内的电离层环境参量,为地震电磁耦合机理和地震前兆信息研究提供了大量宝贵的探测资料(徐方舟等,2012). 除了郎缪尔探针,DEMETER卫星上还搭载有感应式磁力仪、电场探测仪、等离子体分析仪和高能粒子探测仪等探测载荷(Parrot,2002).

    考虑到卫星在运行周期内的轨道变化特征(Parrot et al,2006),本文选取了2006年1月—2010年11月期间的Ne探测数据. 由于太阳辐照和空间磁暴影响电离层参量波动,为了更加真实地反映电离层参量的长时间变化特征,文中按“Kp≥5-Dst_min≤-50 nT”筛选标准排除磁暴影响,并选取DEMETER卫星升轨数据为研究对象(升轨数据对应地方时的夜晚). 同时,将全球划分成90×180个2°×2°经纬度网格,将落在网格内的探测数据月均值形成时间序列数据,以用于后续研究.

    IRI模型是在国际空间研究委员会和国际无线电委员会的联合资助下,IRI工作组根据大量的探测资料和多年积累的电离层研究成果编制开发的全球电离层经验模型(方涵先等,2012). IRI模型基于电离层垂测仪、非相干散射雷达、卫星和探空火箭等探测资料,融合了多个大气参数模型,引入太阳活动参数和地磁活动指数,给出了无极光情况下电离层在地磁平静条件下特定时间、特定地点上空50—2000 km范围内的电子浓度、电子温度、离子温度和电子含量等参数值(Bilitza, 19952001Bilitza,Reinisch,2008方涵先等,2012),因此可以利用IRI模型实现Ne时序数据的模拟与预测,进而进行不同卫星探测数据的时间匹配研究,从而满足非成像遥感数据的交叉检验需求.

    电离层Ne数据具有明显的季节变化和年变化特征(Zhang et al, 200820102012Liu et al,2010刘静等,2013),因此本文引入统计学中经典的时间序列分析法建立ARIMA模型. 该模型的建模思想是首先将预测对象随时间推移而形成的数据序列视为一个随机序列,然后以时间序列的自相关分析为基础,用一定的数学模型来近似描述该序列. 这个模型一旦被识别后,就可由时间序列的过去值和现在值来预测未来值(刘隽等,2011徐方舟等,2012). 其具体原理为:假设时间序列为xt,则xt的现在值可以用该序列滞后值的加权有限组合以及有限的随机扰动项加权组成. 该时间序列可表示为

    也可表示为

    其中,

    式中:φ1φ2,…,φp为自回归参数;θ1θ2,…,θq为移动平均参数;pq为模型阶数;ξt为白噪声过程.

    ARIMA模型是针对平稳时间序列建模的,针对非平稳Ne时间序列数据则需要经过季节差分和低阶差分转变为平稳序列进而进行建模,这个过程为时间序列平稳化. 季节差分算子定义为

    式中,s为时间序列周期. 对差分获得的平稳时间序列,建立周期为sP阶自回归Q阶移动平均季节时间序列模型(徐方舟等,2012):

    式中,ut为季节差分后的时间序列. 当ut非平稳且存在ARIMA成分时,ut可描述为

    将式(8)代入式(7),可得季节差分自回归移动平均模型,即

    式中:P,Qp,q分别表示季节和非季节的自回归移动平均算子的最大滞后阶数;Dd分别表示季节和非季节的差分次数(徐方舟等,2012);AP(Ls)和BQ(Ls)分别表示季节差分后时间序列的自回归成分和移动平均成分.

    本文利用最新版本的IRI-2012模型预测DEMETER卫星轨道高度处2006年1月—2010年11月的Ne月均值数据,并与DEMETER卫星ISL探测的Ne数据进行对比分析. 随机选取30°E上沿60°N—60°S间隔10°的不同纬度点为研究点,模拟该经纬度上的Ne数据,分析IRI模型预测的Ne数据与DEMETER卫星探测的Ne数据的相对误差. 图 1给出了(30°E,60°N)处IRI模型预测Ne数据与DEMETER卫星ISL探测Ne数据的对比图. 图 2给出了30°E和40°E上不同纬度处IRI模型预测Ne数据与DEMETER卫星探测Ne数据的平均相对误差分布. 可以看出,(30°E,60°N)处IRI模型预测Ne数据与DEMETER卫星ISL探测Ne数据的整体趋势较为符合,但也存在一定偏差; IRI模型预测Ne数据与DEMETER卫星ISL探测Ne数据的平均相对误差在高纬度地区较低,为20%左右,而在中低纬度地区则较高,表明IRI模型预测Ne数据相对于ISL探测Ne数据的偏差较大,这是由于中低纬度地区电离层受地磁环境和空间环境影响较大所致.

    图  1  (30°E,60°N)处IRI模型预测Ne数据与ISL探测Ne数据的对比
    Figure  1.  Comparison of Ne data predicted by IRI model to those detected by ISL at (30°E, 60°N)
    图  2  30°E(a)和40°E(b)上不同纬度处IRI模型预测Ne数据与ISL探测Ne数据的平均相对误差分布
    Figure  2.  Distribution of average relative error of Ne data predicted by IRI model relative to that detected by ISL in different latitude along longitude 30°E (a) and 40°E (b)

    本文首先以DEMETER卫星ISL在(30°E,60°N)处2006年1月—2009年12月探测的Ne月均值数据所建立的ARIMA模型来分析未来Ne数据的变化,然后利用2010年1—11月的Ne月均值数据验证已建模型的预测精度. 由于原始序列具有周期性,明显不平稳,故需对原始序列进行一次季节差分和一次非季节差分,差分前后的时间序列如图 3所示. 可以看出,差分后的时间序列在零均值上下波动.

    图  3  差分前(a)、后(b) ISL探测的Ne数据月均值时间序列
    Figure  3.  The time series of monthly average Ne data detected by ISL before (a) and after (b) difference

    图 4给出了差分后时间序列的自相关和偏自相关图. 可以看出,自相关图和偏自相关图均有明显的截尾特征,故差分后的时间序列为平稳序列. 根据图 4所确定的模型阶数进行建模,并运用最佳准则(AIC准则)函数定阶法对模型的阶数和相应参数同时给出一组最佳估计,最后选出的最佳模型阶数为 ARIMA(2,1,0)×(1,1,1)12. 根据所选定的模型阶数,所建ARIMA模型为

    图  4  差分后时间序列的自相关图(a)和偏自相关图(b)
    Figure  4.  Autocorrelation plot (a) and partial autocorrelation plot (b) of time series after difference

    图 5显示了ARIMA模型对原始序列的拟合程度及2010年1月—11月的预测值与序列原始值的比较. 可以看出,ARIMA模型对原始时间序列数据拟合得较好,并在短期预测中预测值相对真实值偏差较小,随着预测时间的增长,预测值相对真实值偏差较大,具体的定量分析将在下节进行.

    图  5  ARIMA模型预测Ne数据与ISL探测Ne数据的对比
    Figure  5.  Comparison of Ne data predicted by ARIMA model to those detected by ISL

    利用上节分析方法,在30°E和40°E经度上不同纬度点处构建ARIMA模型并分析该模型在空间上对原始Ne数据的拟合程度. 在部分测点处有个别月份数据缺失的时间序列采用相邻时间点取均值的方法,部分测点时间序列数据连续缺失较多的测点采用相邻网格点数据取均值的方法,对30°E和40°E经度上不同纬度点处ARIMA模型预测Ne数据与ISL探测Ne数据的平均相对误差进行了统计分析,结果如图 6所示. 可以看出,在不同纬度测点处所构建的ARIMA模型均能很好地拟合原始数据,最大误差为0.15左右,相对误差较小,同时在中低纬度地区的误差相对高纬度地区大,这是由于在中高纬度地区时间序列的季节性变化特征更明显,故所构建的ARIMA模型能更好地提取数据特征,其拟合精度更高.

    图  6  30°E(a)和40°E(b)上不同纬度处ARIMA模型预测Ne数据与ISL探测Ne数据的平均相对误差分布
    Figure  6.  Distribution of average relative error of Ne data predicted by ARIMA model relative to that detected by ISL in different latitude along longitude 30°E (a) and 40°E (b)

    为进一步分析IRI模型和ARIMA模型预测Ne数据的相对误差变化特征,将两种模型预测(30°E,60°N)处2010年1—11月的Ne数据与卫星实际探测的Ne数据进行对比,并分析其相对误差.

    为了更加直观地表现两种时间序列预测模型的预测精度随预测时长的变化特征,图 7给出了这两种时间序列预测模型预测Ne数据与DEMETER卫星ISL实际探测Ne数据的相对误差变化图. 可以看出,ARIMA模型在短期预测中的预测精度较高(预测4个月内的相对误差在10%以内),在11个月后预测的相对误差超过60%. 产生该现象的原因是ARIMA模型预测是依赖于过去时刻的序列值,而预测时刻的值与预测时刻前一段时刻的值有关,由于每次预测都会产生误差,所以预测时间越长,预测误差累积就越大. 而IRI模型的预测相对误差没有随着预测时间的增长而增大(基本保持在20%左右),这是由于IRI模型是集成了多种观测资料的经验模型,在同一地点预测平静时期电离层参量的预测精度不会有太大变化.

    图  7  2010年1—11月IRI模型和ARIMA模型预测Ne数据与ISL探测Ne数据的相对误差变化图
    Figure  7.  The variation of relative error of Ne data predicted by IRI model and ARIMA model relative to that detected by ISL in January to November, 2010
    表  1  IRI模型和ARIMA模型预测Ne数据与ISL实际探测Ne数据的对比及其相对误差
    Table  1.  Comparison of Ne data predicted by IRI model and ARIMA model relative to that detected by ISL as well as their relative errors
    月份ISL探测Ne数据IRI模型ARIMA模型
    预测Ne相对误差预测Ne相对误差
    16.296×1095.120×1090.1875.859×1090.069
    26.886×1097.998×1090.1626.869×1090.002
    31.203×10109.426×1090.2161.184×10100.016
    42.056×10101.592×10100.2261.883×10100.084
    52.621×10102.227×10100.1502.170×10100.172
    62.705×10102.371×10100.1232.304×10100.148
    72.515×10102.191×10100.1292.038×10100.190
    82.259×10101.875×10100.1701.637×10100.275
    91.650×10101.422×10100.1381.245×10100.245
    101.014×10101.190×10100.1746.461×1090.363
    116.072×1097.243×1090.1939.950×1090.639
    下载: 导出CSV 
    | 显示表格

    本文基于DEMETER卫星ISL载荷2006年1月—2010年11月探测的Ne月均值数据,首先验证了IRI模型在不同纬度上模拟Ne数据的精度,其次根据Ne数据的年变化特征和季节变化特征构建了ARIMA模型,最后比较分析了两种时间序列分析模型预测Ne数据的相对误差变化特征.

    对试验结果进行分析可知:IRI模型在高纬度地区模拟预测Ne数据的相对误差保持在20%左右,其在中低纬度地区的相对误差则较高;ARIMA模型对原始序列拟合得较好,短期内的预测相对误差在10%以内,证明了该模型能够反映电离层探测参量的时间变化特征,可以用于电离层参量背景场的构建.

    两种时间序列分析模型均有一定的局限性,如:IRI模型模拟预测精度较低,尤其是在中低纬度区域,只能反映平静状态下电离层参量的变化特征;ARIMA模型在短期内预测精度较高,但其预测误差会随着预测时间的增长而增大,同时ARIMA模型是根据时间序列的过去值特征所构建,提取时间序列的变化特征时用于构建模型的过去值越多,特征越明显,预测精度就越高,若进行分析研究的时间序列无明显变化特征或无足够的过去值则无法进行建模.

    在后期研究中我们将结合实际观测数据对模型进行改进,进一步提高预测精度,为研究不同期在轨卫星探测数据的时间匹配提供基础,进而开展交叉检验方法的研究.

  • [1] Heiskanen, W. A. and Vcning——Meinesz, F. A., The Earth and its Gravity Field, McGraw——Hill, New York,1958.

    [2] Stacey, F. D.,地球物理学,中国科学技术大学地球物理教研室译,傅承义校,地震出版社,1981,

    [3] Banks, R. J., R. L. Parker, and S. P. Huestis, Isostatic compensation on a continental scale, Local versus regional mechanisms, Geophys. J. R. astr. Soc., 51, 431——452, 1977.

    [4] Dorman, L. M. and B. T. R. Lewis, Experimental isostasy. I:Theory of the determination of the Earth's isostatic response to a concentrated load, J. Geophys. Rer., 75, 3357——3365, 1970.

    [5] Dehlinger, P.,海洋重力学,詹贤均、高玉芬等译,张少泉校,海洋出版社,1981,

    [6] 雷受显,重力广义地形改正值和均衡改正值的一种计算方法,海洋地质与第四纪地质,4, 101——111,1984,

    [7] 王悉基、程振炎,均衡异常与地壳结构,物化探研究报道,9,1——13, 1982,

    [8] Artyushkov, E. V., Rheological properties of the crust and upper mantle according to data on isostatic movements, J. Geophys. Rer., 76, 1376——1396, 1971.

    [9] 顾功叙、曾融生、张忠撒,中国境内208处重力加速度测定之海陆均衡变差(一),地球物理学报,1, 101——139, 1949,

    [10] 顾功叙、曾融生,中国境内208处重力加速度测点之海陆均衡变差(二),地球物理学报,2, 14——26, 1950,

    [11] 王憋基、程家印、程振炎,我国地壳深部构造的区域特征,物探与化探,5,193——204, 1981,

    [12] 魏梦华、王启鸣、史志宏、殷秀华、刘占坡、张玉梅,中国大陆自由空气重力场的初步研究,地震地质,3,3, 47——60, 1981.

    [13] 股秀华、史志宏、刘占坡、张玉梅,华北北部均衡重力异常的初步研究,地震地质,4,4, 27——34, 1982,

    [14] 袁宝印、孙建中,华北新生代沉积与断块构造,华北断块区的形成与发展,221——229,科学出版社,1980,

    [15] Parker, R. L., The rapid calculation of potential anomalies, Geophysics, 39, 447——455, 1973.

    [16] 冯锐,中国地壳厚度及上地慢密度分布(三维重力反演结果),地震学报,7, 143—— 157, 1985,

    [17] 冯锐,三维物性分布的位场计算,地球物理学报,29, 399——406, 1986,

    [18] 冯锐、严惠芬、张若水,三维位场的快速反演方法及其程序设计,地质学报,60, 390——403, 1986.

    [19] 冯锐、郑书真,对华北岩石圈构造的综合研究,科学通报,32,1987,

    [20] 李润江、伶淑娟、孟令顺、张琳、张凤华,东北地区重力均衡异常特征的初步研究,地震地质,8,3,61——68,1986,

    [21] 冯锐、周海南、姚政生、孙克忠,面波的频散、反演和层析成象,中国地震,3, 15——28, 1987,

    [22] Kahle, H. G. and D. Werner, A geophysical study of the Rhinegraben——II, Gravity anomalies and geothermal implications, Geophys. J. R. astr. Soc., 62, 631——647, 1985.

    [1] Heiskanen, W. A. and Vcning——Meinesz, F. A., The Earth and its Gravity Field, McGraw——Hill, New York,1958.

    [2] Stacey, F. D.,地球物理学,中国科学技术大学地球物理教研室译,傅承义校,地震出版社,1981,

    [3] Banks, R. J., R. L. Parker, and S. P. Huestis, Isostatic compensation on a continental scale, Local versus regional mechanisms, Geophys. J. R. astr. Soc., 51, 431——452, 1977.

    [4] Dorman, L. M. and B. T. R. Lewis, Experimental isostasy. I:Theory of the determination of the Earth's isostatic response to a concentrated load, J. Geophys. Rer., 75, 3357——3365, 1970.

    [5] Dehlinger, P.,海洋重力学,詹贤均、高玉芬等译,张少泉校,海洋出版社,1981,

    [6] 雷受显,重力广义地形改正值和均衡改正值的一种计算方法,海洋地质与第四纪地质,4, 101——111,1984,

    [7] 王悉基、程振炎,均衡异常与地壳结构,物化探研究报道,9,1——13, 1982,

    [8] Artyushkov, E. V., Rheological properties of the crust and upper mantle according to data on isostatic movements, J. Geophys. Rer., 76, 1376——1396, 1971.

    [9] 顾功叙、曾融生、张忠撒,中国境内208处重力加速度测定之海陆均衡变差(一),地球物理学报,1, 101——139, 1949,

    [10] 顾功叙、曾融生,中国境内208处重力加速度测点之海陆均衡变差(二),地球物理学报,2, 14——26, 1950,

    [11] 王憋基、程家印、程振炎,我国地壳深部构造的区域特征,物探与化探,5,193——204, 1981,

    [12] 魏梦华、王启鸣、史志宏、殷秀华、刘占坡、张玉梅,中国大陆自由空气重力场的初步研究,地震地质,3,3, 47——60, 1981.

    [13] 股秀华、史志宏、刘占坡、张玉梅,华北北部均衡重力异常的初步研究,地震地质,4,4, 27——34, 1982,

    [14] 袁宝印、孙建中,华北新生代沉积与断块构造,华北断块区的形成与发展,221——229,科学出版社,1980,

    [15] Parker, R. L., The rapid calculation of potential anomalies, Geophysics, 39, 447——455, 1973.

    [16] 冯锐,中国地壳厚度及上地慢密度分布(三维重力反演结果),地震学报,7, 143—— 157, 1985,

    [17] 冯锐,三维物性分布的位场计算,地球物理学报,29, 399——406, 1986,

    [18] 冯锐、严惠芬、张若水,三维位场的快速反演方法及其程序设计,地质学报,60, 390——403, 1986.

    [19] 冯锐、郑书真,对华北岩石圈构造的综合研究,科学通报,32,1987,

    [20] 李润江、伶淑娟、孟令顺、张琳、张凤华,东北地区重力均衡异常特征的初步研究,地震地质,8,3,61——68,1986,

    [21] 冯锐、周海南、姚政生、孙克忠,面波的频散、反演和层析成象,中国地震,3, 15——28, 1987,

    [22] Kahle, H. G. and D. Werner, A geophysical study of the Rhinegraben——II, Gravity anomalies and geothermal implications, Geophys. J. R. astr. Soc., 62, 631——647, 1985.

计量
  • 文章访问数:  1210
  • HTML全文浏览量:  21
  • PDF下载量:  93
  • 被引次数: 0
出版历程
  • 发布日期:  2011-09-01

目录

/

返回文章
返回