多源数据的面向对象国际河流土地覆被分类研究

孔博, 邓伟, 李爱农, 杨勇

孔博, 邓伟, 李爱农, 杨勇. 多源数据的面向对象国际河流土地覆被分类研究[J]. 武汉大学学报 ( 信息科学版), 2015, 40(7): 943-949. DOI: 10.13203/j.whugis20130107
引用本文: 孔博, 邓伟, 李爱农, 杨勇. 多源数据的面向对象国际河流土地覆被分类研究[J]. 武汉大学学报 ( 信息科学版), 2015, 40(7): 943-949. DOI: 10.13203/j.whugis20130107
KONG Bo, DENG Wei, LI Ainong, YANG Yong. Object-oriented Landcover Classification of Multi-source Remote Sensing Data in International Trans-boundary River[J]. Geomatics and Information Science of Wuhan University, 2015, 40(7): 943-949. DOI: 10.13203/j.whugis20130107
Citation: KONG Bo, DENG Wei, LI Ainong, YANG Yong. Object-oriented Landcover Classification of Multi-source Remote Sensing Data in International Trans-boundary River[J]. Geomatics and Information Science of Wuhan University, 2015, 40(7): 943-949. DOI: 10.13203/j.whugis20130107

多源数据的面向对象国际河流土地覆被分类研究

基金项目: 中国科学院重点部署资助项目(KZZD EW 08 01);“一三五”方向性资助项目(sds 135 1205 03);中国科学院战略性先导科技专项(B类)资助项目(XDB03030507);国家自然科学基金资助项目(41301094);国土资源部地学空间信息技术重点实验室开放研究基金资助项目(KLGSTT2014 06)
详细信息
    作者简介:

    孔博,助理研究员,主要从事生态遥感研究。

    通讯作者:

    邓伟,研究员,博士。

  • 中图分类号: P237;P951

Object-oriented Landcover Classification of Multi-source Remote Sensing Data in International Trans-boundary River

Funds: TheImportDeployingProject,ChineseAcademyofSciences,No.KZZD EW 08 01;the“One Three Five”Directiv ityInnovationProjectofIMHE,No.sds 135 1205 03;theLeadStrategicProject,CAS,No.XDB03030507;theNationalNaturalSci
More Information
    Author Bio:

    KONG Bo: KONGBo,assistantprofessor,specializesinecologicalremotesensing.

    Corresponding author:

    DENG Wei: DENGWei,PhD,professor.

  • 摘要: 作为横跨3个国家(尼泊尔、印度、中国)的国际跨界河流———柯西河流域,地形高差巨大,土地覆被结构组成复杂,进行土地覆被的自动分类研究具有典型意义。基于面向对象方法多源遥感数据、训练规则、丰富的细节信息为复杂土地覆被自动分类研究提供了可能。选择合适的影像分割特征和最优分割尺度,按照数据挖掘中的规则顺序逐步进行各个土地覆被的提取。总体精度说明分类结果与野外点相一致的概率能达到90.05%,说明国际跨界河流土地覆被分类方法是可行的,分类结果是准确、可信的。
    Abstract: Landcoverclassificationisnoteasyforitsinnercomplexityresultingfromhugeterrainele vationinKosiRiverasinternationaltrans boundaryriver,ithaspossibletoclassificationdifficultiesflowingthroughthreecountries.Theemergencyofmulti sourceremotesensingimagesandtrainingalgorithmmakeitpossibletoclassifythelandcoverforitsampledetailsbasedonobject orientedmethod.Thepaperisaboutlandcoverclassificationmethodselectingfeatureandoptimalscaleofseg mentation,theinnovationofthismethodliesinselectionofproperscaleparameterresultingfromproperimagedataandcertainclassificationorder.Thetotalaccuracyhashighly90.05%comparedwithobjectorientedclassificationresultandactualsamplingpoint,whichisafeasiblemethod,andtheclassificationresultsaremoreaccurate.
  • 随着中国经济的飞速发展,工业、交通也迅猛发展,城市扩张,人口集中,各种环境问题也层出不穷,其中大气污染,尤其是颗粒物(气溶胶)污染具有广泛的空间分布,对人体健康产生较大危害[1]。空气动力学当量直径小于10 μm的颗粒物被称为PM10,直径小于2.5 μm的颗粒物被称为PM2.5,二者均具有可吸入性,前者通常可以沉积到上呼吸道,而后者可以进入肺泡,对身体危害极大;同时细颗粒物在大气中停留时间长,也影响了大气能见度和大气环境质量[1-3]

    在颗粒物污染监测中,近地表的监测备受人们关注,各个城市中普遍设立了空气质量监测站,同时监测颗粒物浓度(PM10、PM2.5)、硫氧化物、氮氧化物等参数,能够测量站点附近大气污染的详细参数,测量精度高,但受制于成本,不能广泛布站,覆盖区域的测量很难实现,无法全面反映气溶胶的分布。激光雷达作为一种遥感手段,能够进行长距离、高频次的探测,探测结果可以有效反映激光传输路径内的气溶胶分布情况,是当前遥感探测颗粒物污染的有效手段,被广泛应用于城市大气污染探测中[4-11]。国内外多家单位开展了激光雷达用于大气探测的研究,其中中科院安光所[4-7]、武汉大学[8]、西安理工大学[9]等多家国内单位,以及法国索邦大学[10]、日本国立环境研究所[11-12]、意大利[13]、美国NOAA[14]等国外研究机构均利用激光雷达进行了颗粒物探测。

    当前国内外的研究中,都是单独进行垂直探测[1, 5, 8-10, 12-17]或者水平探测[4, 6-7],鲜有结合三维探测的激光雷达数据进行相关研究。垂直探测可以确定探测地点上空的气团垂直分布和运动,但是不能获得区域性的颗粒物分布;水平探测可用于识别该区域的颗粒物水平分布,但牺牲了时间分辨率,无法获得颗粒物的传输规律;而同时结合垂直和水平探测组成的三维数据,即可实现两种探测方式的优势互补,从而实现对颗粒物更全面的监测和研究。本文即利用两台微脉冲激光雷达对天津市武清区中心城区进行了全天时的遥感探测,其中一台激光雷达进行垂直测量,探测颗粒物的传输、沉降;另外一台激光雷达架设于高楼顶,进行水平扫描,探测颗粒物的水平分布。通过对两台激光雷达的三维数据的分析,揭示区域性的颗粒物污染分布和传输情况。

    激光雷达是一种主动探测设备,它向大气中发射激光脉冲,激光束与大气气溶胶粒子和大气分子等发生相互作用,产生回波信号被接收望远镜接收,然后对回波信号进行采集、反演就可以得到大气气溶胶的光学特征信息。激光雷达探测大气颗粒物具有很高的空间和时间分辨率,是目前有效获得大气颗粒物垂直分布廓线的重要手段。

    本文研究中为了获取大气颗粒物的三维数据,同时采用了两台激光雷达分别进行大气垂直和水平扫描探测。垂直探测和水平扫描激光雷达的技术指标如表 1所示。在大气激光雷达探测中一般采用的是Nd:YAG激光器的基频光,即1 064 nm,因此本文采用的激光雷达的激光波长均为1 064 nm。同时对激光雷达进行了一体化设计,并加入窄带滤光片,使得其具有全天时、全天候工作能力,从而实现24 h连续监测。

    表  1  激光雷达技术参数
    Table  1.  Parameters of the Two LiDARs
    参数 垂直LiDAR 水平LiDAR
    波长/nm 1 064 1 064
    能量/μJ >100 >20
    脉冲频率/kHz 1 1
    望远镜口径/ mm 100 100
    距离分辨率/m 30 15
    滤光片带宽/nm 0.5 0.5
    下载: 导出CSV 
    | 显示表格

    激光雷达的布设区域选择在天津市武清区城区。武清城区处于京津冀的中心地带,其大气颗粒物污染受外地源影响较为明显,希望通过激光雷达遥感探测的手段找出污染的传输方向,并初步判定污染源。垂直探测激光雷达和水平探测激光雷达分别部署于同一小区中不同建筑楼顶,为保证水平探测路径内无遮挡,将水平探测激光雷达部署于城区内的一栋33层居民楼楼顶,具体地理位置在水平扫描结果中可见。

    大气探测激光雷达信号反演中一般采用Fernald算法[18],该算法主要是区分了大气气溶胶和大气分子,并引入激光雷达比的概念,进行大气气溶胶消光系数和后向散射系数的计算。具体如下:

    $$ {P(r) = EC{r^{ - 2}}\beta (r)\exp \left[ { - 2\int_0^r \sigma (r){\rm{d}}r} \right]} $$ (1)
    $$ {{S_1} = {\sigma _1}(r)/{\beta _1}(r)} $$ (2)
    $$ {{S_2} = {\sigma _2}(r)/{\beta _2}(r)} $$ (3)
    $$ \begin{array}{l} {\sigma _1}(I - 1) + \frac{{{S_1}}}{{{S_2}}}{\sigma _2}(I - 1) = X(I - 1) \times \\ \exp [ + A(I - 1, I)]/\left\{ {\frac{{X(I)}}{{{\sigma _1}(I) + \left( {{S_1}/{S_2}} \right){\sigma _2}(I)}} + } \right.\\ X(I)\Delta r + X(I - 1)\exp [ + A(I - 1, I)]\Delta r\} \end{array} $$ (4)

    式(1)为典型的大气激光雷达方程,是大气颗粒物探测的基本方程,所有利用激光雷达进行大气探测均是基于该方程进行的,其中,P为激光雷达回波信号强度功率;r代表探测距离;E为激光雷达固定参数,包括发射激光功率等;C为矫正常数;β为大气后向散射系数;σ为大气消光系数。式(2)和式(3)分别是大气气溶胶和大气分子的激光雷达比计算方程,是Fernald算法的核心,其中S1S2分别是大气气溶胶和大气分子的激光雷达比;大气分子的激光雷达比为固定值8π/3,下标1表示大气颗粒物/气溶胶,下标2表示大气分子。式(4)为利用Fernald算法进行消光系数的后向积分方程,其中I代表高度相关的计数值。

    通过确定激光雷达比、标定高度和标定高度处的大气气溶胶消光系数3个参数后,即可以利用Fernald算法进行气溶胶消光系数的计算。其中激光雷达比一般采用定值简化计算,本文设定为50;标定高度选择在相对洁净、大气气溶胶少的区域,本文选择高度10 km为标定高度rc,该处的气溶胶消光系数计算采用经验公式a1(rc)=(1.01-1)a2(rc)计算[18],其中分子的消光系数由大气模式给出。

    水平扫描激光雷达信号处理算法是在垂直探测算法的基础上完成的。目前一般采用的是斜率法或者Fernald法,但前者的问题在于其假定的大气分布均匀状态在实际情况下不成立,而后者的弊端在于无法获得标定高度(距离)及在此处的气溶胶消光系数[18-21]

    为减小反演误差,本文结合斜率法和Fernald算法,提出了一种新的水平扫描反演算法。斜率法不存在激光雷达比的问题,但水平方向上某小段距离内可以认为大气状态均匀,在本次实验中,先根据原始激光雷达信号确定大气分布均匀的一段距离,并以斜率法获取该处的气溶胶平均消光系数,将结果作为Fernald算法的标定高度(距离)和此处的消光系数标定值。

    由于获取的是中间距离的消光系数,因此利用Fernald算法时,分别进行前向和后向积分,从而得到较为准确的大气气溶胶消光系数反演结果。

    本文选取2018年7月16日水平和垂直激光雷达数据进行反演,揭示该地区的颗粒物三维分布特征,并同地面空气质量国控站点的PM数据进行校正和对比分析。

    对7月16日连续24 h的垂直探测结果,利用Fernald算法进行反演,标定高度选为10 km,图 1(a)为2018年7月16日凌晨00:30时刻的单条垂直激光雷达距离校正信号,可以看出激光传输距离较远,表明高层大气较为洁净,在高度0.5 km附近,信号强度先增加后减小是由于overlap效应造成的,即激光没有完全进入望远镜视场;在高度2.5 km处,激光雷达回波信号出现突增现象,表明此处有明显的云层,但厚度较薄,穿过云层后激光的衰减较小,继续传输至高空中。

    图  1  垂直激光雷达探测结果
    Figure  1.  Observation Results of Vertical LiDAR

    图 1(b)是7月16日连续24 h垂直探测得到的大气气溶胶消光系数,通过消光系数可以直接反映出大气气溶胶的浓度情况,其中高度2.5 km处的消光系数出现极大值,是由于云层造成的,不在本文的研究范围内。

    一般情况下,气溶胶主要集中在边界层内,由于探测时间为夏天,边界层高度较高,利于污染物扩散。根据查询到的气象数据,激光雷达探测从0时至早上8时左右,边界层以下的颗粒物浓度较低,且呈逐渐消散(扩散)的趋势;但从早上8时左右开始,上午低空的气溶胶消光系数持续增加,可以看出,颗粒物污染也持续增加,表明当地的人为颗粒物排放较为严重,同时颗粒物由低空向高空传输的特征明显,因此可以判定当天的大气气溶胶污染主要是本地源,即从地面向高空扩散;白天污染持续增加,到夜晚24时达到顶峰,这一现象主要是人为活动造成的。武清区是天津和北京的交界处,该处白天和夜晚的交通活动较强,同时夜晚气溶胶浓度持续增加的主要原因不排除当地企业夜间生产偷排的情况。

    通过布设在居民楼顶的水平扫描激光雷达(117.04°E,39.39°N)探测气溶胶在水平方向上的分布。为了保证信噪比,将单次扫描角度时间设置为20 s,扫描步进角度为2°,水平方向顺时针扫描一圈时间为1 h。利用本文提出的水平反演算法进行了消光系数反演,图 2(a)为7月16日凌晨00:30激光雷达单次扫描得到的距离校正信号,同图 1(a)相比,水平探测中,由于地表附近较高空气溶胶浓度高,因此回波信号随距离增加而衰减的程度较大;同时由于水平气溶胶分布的不均匀性,水平探测的回波信号波动也较大,采用斜率法进行反演会造成一定的反演误差。

    图  2  水平扫描激光雷达探测结果
    Figure  2.  Observation Results of Horizontal Scanning LiDAR

    图 2(b)为水平方向上1 h(16日凌晨00:25-01:25)探测得到的大气气溶胶消光系数分布图,其中A为激光雷达站点位置,图中消光系数较大的区域气溶胶浓度较高,此段时间内污染集中在以探测点为中心的主城区内,以及南部高铁站附近。由于是凌晨,气溶胶污染主要集中在城区附近,人口密度高,与生产、生活、交通的污染排放有关,而郊区的人类活动和生产较少,污染较低。

    图 2(c)为一天24 h连续探测后的数据融合结果,反映24 h内区域中的污染较大区域分布。采取的方法是提取高污染区域的颗粒物数据,剔除低污染部分,并融合到一张图中。可以看出,当天24 h内共有8个污染较重的区域,标号3~8的污染区域主要集中在城区中,表明城区人为活动是造成颗粒物污染的主要原因;1号、2号区域位于远城区,人类活动较少,但实际污染程度也较大,根据实地调查,两处地区有本地工厂,工业活动是造成该区域气溶胶浓度上升的主要原因。

    根据垂直激光雷达探测的结果,如图 1(a)所示,可知当天气溶胶污染源主要是本地源,结合图 2(c),充分表明在本地污染源起主导作用的情况下,大气气溶胶污染主要集中在城区部分;同时结合图 2(b)可知,城区地区发达的交通和密集的人口加重了武清区的夜间污染,且密集的建筑物使得气溶胶扩散不易。

    需要说明的是,激光雷达探测反演得到的颗粒物消光系数是相对值,可以在一定程度上反映出观测区域中的颗粒物浓度情况,一般情况下,这种探测结果结合地面国控站数据即可实现区域内的颗粒物密集监测,但无法直接将颗粒物的消光系数转换到绝对的颗粒物浓度值。研究表明,颗粒物浓度同激光雷达探测得到的消光系数间存在一定的线性关系,尤其是颗粒物消光系数和PM10之间的线性关系更为明显。因此在完成垂直和水平激光雷达探测实验后,为初步验证激光雷达探测结果的准确性,将激光雷达水平扫描结果同地面站点数据进行对比,地面站点数据采用的是扫描区域内的国控站点PM10数据。

    需要指出的是,由于激光雷达是位于高点处进行扫描的,其探测结果同地面国控站点间必定存在一定的偏差,因此在进行相关性计算时,需要考虑到当天的风速和风向,若风速较大,需要以下风处的气溶胶消光系数平均值作为对比值,为了降低对比误差,本文此次选取了风速较低的情况下的数据。

    图 3是激光雷达探测到的气溶胶消光系数同PM10间的相关性对比,由于水平激光雷达数据获取周期为1 h,因此对国控站的PM10数据进行小时平均,采用223组对比数据,即近10 d的数据,拟合得到的颗粒物消光系数反演PM10浓度的经验公式为:y=767.82x+10.08,R2为0.85,表明二者间的线性关系较为明显,即颗粒物的消光系数可以转换为精度较高的颗粒物浓度数据。若通过长期数据积累,可以得到更为准确的经验公式。

    图  3  气溶胶消光系数和PM10的相关性对比分析图
    Figure  3.  Comparison and Analysis of Correlation Between Aerosol Extinction Coefficient and PM10

    通过上述实际数据对比分析,可以认为气溶胶消光系数和颗粒物浓度(尤其是PM10)间的线性程度较高,因此,在获取一定时间内的经验公式后,可以利用激光雷达探测到的气溶胶消光系数的水平分布图,直接反映扫描区域内颗粒物浓度的分布情况。

    本文利用两台激光雷达分别进行垂直方向和水平扫描探测,获得了区域性的气溶胶三维分布结果。其中利用垂直探测得到气溶胶的垂直分布,并说明一天24 h内边界层高度变化下低层大气气溶胶的扩散情况,初步揭示当天的气溶胶污染源是本地或是外地;同时利用布设在高楼顶的水平扫描激光雷达获取气溶胶的水平分布结果,探明扫描区域内的气溶胶污染分布情况。三维数据的激光雷达相比于传统的被动站点布设,能以更少的成本反映区域的气溶胶分布情况以及变化趋势,为城市颗粒物监测和治理提供了有效的测量手段。

    需要指出的是,单独的激光雷达只能对区域气溶胶进行实时监测,若加上偏振通道,则可以区分本地和外源沙尘性气溶胶,但仍然无法实现对污染物的溯源、预测。而卫星遥感手段能够对大尺度的PM2.5、PM10进行监测,结合气象数据,构建后向轨迹模型,追溯气团运动轨迹,可以实现对颗粒物污染的追踪和预测,这也是下一步需要完善的工作。

  • [1] LiZhifei.Trans boundaryRiverIssuesandtheChi nesePeripherySecurity[J].犃犮犪犱犲犿犻犮犈狓狆犾狅狉犪 狋犻狅狀,2011(1):113(李志斐.跨国界河流问题与中国周边关系[J].学术探索,2011(1):113)[2] ChenYunhao,FengTong,ShiPeijun,etal.Clas sificationofRemotSensingImageBasedonObjectOrientedandClassRules[J].犌犲狅犿犪狋犻犮狊犪狀犱犐狀犳狅狉 犿犪狋犻狅狀犛犮犻犲狀犮犲狅犳犠狌犺犪狀犝狀犻狏犲狉狊犻狋狔,2006,31(4):316320(陈云浩,冯通,史培军,等.基于面向对象和规则的?杏跋穹掷嘌芯浚郏剩荩浜捍笱аПāば畔⒖蒲О?2006,31(4):316 320)[3] ChenQihao.ResearchandRealizationofMulti SourcesRemoteSensingDataonObjectOriented[D].Wuhan:ChinaUniversityofGeoseiences,2007(陈启浩.面向对象的多源遥感数据分类技术研究与实现[D].武汉:中国地质大学,2007)[4] SuWei,LiJing,ChenYunhao,etal.Object orien tedUrbanLand coverClassificationofMulti scaleImageSegmentationMethod———ACaseStudyinKualaLumpurCityCenter,Malaysia[J].犑狅狌狉狀犪犾狅犳狉犲犿犪狋犲狊犲狀狊犻狀犵,2007,11(4):521529(苏伟,李京,陈云浩,等.基于多尺度影像分割的面向对象城市土地覆被分类研究———以马来西亚吉隆坡市城市中心区为例[J].?醒Пǎ玻埃埃罚保保ǎ矗海担玻豹?29)[5] SunDanfeng,YangJihong,LiuShunxi.Applica tionofHigh SpatialIKNOSRemoteSensingImagesinLandUseClassficationandChangeMonitoring [J].犜狉犪狀狊犪犮狋犻狅狀狊狅犳狋犺犲犆犛犃犈,2002,18(2):160 164(孙丹峰,杨冀红,刘顺喜.高分辨率遥感卫星影像在土地利用分类及其变化监测的应用研究[J].农业工程学报,2002,18(2):160 164)[6] BlaschkET,LangS,LorupE,etal.Object orien tedImageProcessinginanIntegratedGIS/RemoteSensingEnvironmentandPerspectivesforEnviron mentalApplications[C].EnvironmentalInformationforPanning,PoliticsandthePublicMetropolisVer lag,Marburg,2000[7] MauroC,EufemiaT.AccuracyAssessmentofPer fieldClassificationIntegratingveryFineSpatialRes olutionSatelliteImagerywithTopographicData[J].犑狅狌狉狀犪犾狅犳犌犲狅狊狆犪狋犻犪犾犈狀犵犻狀犲犲狉犻狀犵,2001,3(2):127 134[8] BenzUC,HofmannP,IllhauckWG,etal.Multi resolution,Object orientedFuzzyAnalysisIsofRe moteSensingDataforGISReadyInformation[J].犐犛犘犚犛犑狅狌狉狀犪犾狅犳犘犺狅狋狅犵狉犪犿犿犲狋狉狔犪狀犱犚犲犿狅狋犲犛犲狀狊犻狀犵,2004,58:239258[9] LoboA,ChicO,CasteradA.ClassificationofMediterraneanCropswithMulti sensorData:Per pixelVersusPer objectStatisticsandImageSeg mentation[J].犐狀狋犲狉狀犪狋犻狅狀犪犾犑狅狌狉狀犪犾狅犳犚犲犿狅狋犲犛犲狀狊犻狀犵,1996,17:23582400[10] BaatzM,SchapeA.Object OrientedandMulti ScaleImageAnalysisinSemanticNetworks[C].The2ndInternationalSymposiumonOperationali zationofRemoteSensing,Ensehede,1999[11] BaatzM,SchapeA.Multi resolutionSegmenta tion:AnOptimizationApproachforHighQualityMulti scaleImageSegmentation[C].ProceedingsofAGIT’00,Heidelberg,Germany,2000[12] BhaskaranS,ParamanandaS,RamnarayanM.Per pixelandObject OrientedClassificationMethodsforMappingUrbanFeaturesUsingIKONOSSatelliteData[J].犃狆狆犾犻犲犱犌犲狅犵狉犪狆犺狔,2010,30:650665[13] YuHuan,ZhangShuqing,KongBo,etal.OptimalSegmentationScaleSelectionforObject orientedRe moteSensingImageClassification[J].犑狅狌狉狀犪犾狅犳犐犿犪犵犲犪狀犱犌狉犪狆犺犻犮狊,2010,15(2):352360(于欢,张树清,孔博,等.面向对象?杏跋穹掷嗟淖钣欧指畛叨妊≡裱芯浚郏剩荩泄枷裢夹窝Пǎ玻埃保埃?5(2):352 360)[14] LiHui,HuXiaomei.AnalysisandComparisonBe tweenID3AlgorithmandC4.5AlgorithminDeci sionTree[J].犠犪狋犲狉犚犲狊狅狌狉犮犲狊犪狀犱犘狅狑犲狉,2008,26(2):129132(李会,胡笑梅.决策树中ID3算法与C4.5算法分析与比较[J].水电能源科学,2008,26(2):129 132)
  • 期刊类型引用(3)

    1. 陶叶青,束明聪,陈浩. 基于分类估计的异质数据融合M_(split)模型. 淮阴师范学院学报(自然科学版). 2025(01): 45-51 . 百度学术
    2. 张双成,赵颖,张成龙,张菊清,樊茜佑,司锦钊,张雅斐,朱武,李振洪. 联合光学遥感和SAR影像青海玛多Ms 7.4地震同震形变场分析. 武汉大学学报(信息科学版). 2025(03): 469-482 . 百度学术
    3. 丁明涛,陈浩杰,李振洪,刘振江. 利用光学遥感影像光流场模型进行地表形变分析. 武汉大学学报(信息科学版). 2024(08): 1314-1329 . 百度学术

    其他类型引用(2)

计量
  • 文章访问数: 
  • HTML全文浏览量: 
  • PDF下载量: 
  • 被引次数: 5
出版历程
  • 收稿日期:  2013-05-01
  • 修回日期:  2015-07-04
  • 发布日期:  2015-07-04

目录

/

返回文章
返回