第37卷第3期农业工程学报V ol.37 No.3 132 2021年2月Transactions of the Chinese Society of Agricultural Engineering Feb. 2021
徐春萌1，田芷源2※，陈 威3,4，刘佳佳3,4，白 洁1
（1. 北京佳格天地科技有限公司，北京100190；2. 中国科学院南京土壤研究所，南京210008；
3. 中国农业科学院农业信息研究所，北京100081；
4. 农业农村部农业信息服务技术重点实验室，北京100081）
摘要：开展基于作物模型的大面积作物产量估测研究，可以为及时掌握全球重点地区农作物的生产情况提供数据支撑。 该研究以大豆为监测作物，选取中国吉林省和美国爱荷华州作为研究区域，基于DSSAT作物估产模型中的SOYGRO大豆模型，利用分辨率为0.5°×0.5°的生育期气象要素以及500 m×500 m绿叶绿素植被指数，进行遥感数据融合作物模型估测大豆单位产量的模拟与验证研究。结果显示，2008－2017年，美国爱荷华州大豆单位产量模拟值的平均误差为16.8%，均方根误差为762.8 kg/hm2，平均偏差为107.2 kg/h
m2；中国吉林省大豆单位产量估测的平均误差为36.3%，均方根误差为1 088.4 kg/hm2，平均偏差为−237.9 kg/hm2。在县域尺度下，大豆单位产量模拟值与调查值的拟合度较好，尤其在产量较低的年份，其中美国爱荷华州的产量相关系数最高可达0.78，中国吉林省的相关系数偏小，为0.59，表明对美国爱荷华州大豆单位产量的估测精度优于中国吉林省。研究所建立的大豆单位产量估测技术路线，可以为中美两国主产区作物单位产量的大面积有效估测提供参考。
关键词：遥感；作物；模型；大豆；产量；绿叶绿素植被指数
doi：10.11975/j.issn.1002-6819.2021.03.016
中图分类号：S127, TP79 文献标志码：A 文章编号：1002-6819(2021)-03-0132-08
徐春萌，田芷源，陈威，等. 基于DSSAT作物模型的中美大豆主产区单产模拟与验证[J]. 农业工程学报，2021，37(3)：132-139. doi：10.11975/j.issn.1002-6819.2021.03.016
Xu Chunmeng, Tian Zhiyuan, Chen Wei, et al. Simulations and validations of the soybean yields per unit area using DSSAT crop model in the major soybean producing areas of China and America[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(3): 132-139. (in Chinese with English abstract) doi：
10.11975/j.issn.1002-6819.2021.03.016
0 引 言
鉴于中国与国外的粮食贸易规模不断加大，对世界主要农作物产量估测的需求日益突出，对估产方法高精度与时效性的要求越来越高[1]，开展客观、及时、大面积的农作物产量预测研究，可以减少对人工统计数据的依赖性，对政府决策具有重要意义。农作物单位产量（单产）的预测研究方法，包括基于统计回归模型和相似年份聚类法的纯气象研究手段，建立植被指数与作物单产关系的纯遥感估产方法，以及逐步发展形成的作物模型同化遥感数据的集合预测方法[2-4]，尽管研究已提出多个作物模型[5-7]。但由于作物模型是基于单个采样点进行开发的，在大面积区域应用时仍面临模型区域化处理的问题。卫星遥感技术可进行大范围地表监测，并获取其中的植被指数[8-10]。遥感-作物模型同化是将遥感数据与基于单点尺度的作物模型进行互补，在考虑大豆生育期生长规律的基础上，借助遥感植被指数反演叶面积，融合
收稿日期：2020-05-25 修订日期：2020-12-02
基金项目：农业农村部农业信息服务技术重点实验室开放基金项目（CAAS-AII NYXXJSFW 2019-001）
作者简介：徐春萌，主要研究方向为农业气象应用与作物估产。
Email：******************
※通信作者：田芷源，博士，主要研究方向为遥感数据应用、土壤发生与侵蚀过程。Email：***********************生育期的气象时空分布数据，实现大面积范围内的大豆单产估测。
农业技术转让决策支持系统（the Decision Support System for Agro-technology Transfer，DSSAT）是由美国乔治亚大学组织开发的[11]，并于20世纪末引入中国，在中国不同地区进行了适用性验证[12-14]。DSSAT模型可以对众多作物种类进行生长模拟，包括玉米、大豆、棉花等重要粮食经济作物[15]。该作物模型采用绿叶绿素植被指数（Green Chlorophyll Vegetation Index，GCVI）作为模拟植被生长情况的指标。相比于常用的归一化植被指数（Normalized Difference Vegetation Index，NDVI），GCVI具有当叶面积很高时不会发生指标饱和的优势[16]。
在国内，科研工作者在中国以及国外地区（比如美国Corn Belt）开展了使用遥感信息及作物模型获取作物生物量及产量的探索研究[17-19]，积累了一定的技术方法。随着作物估产需求的提升，近年来越来越多的研究工作关注于同化遥感信息获取详细产量数据（如市县级、地块级别）的应用方法[20-21]。目前针对主粮作物（玉米和小麦）的试验较多，而对于其他重要粮食作物，比如大豆估产模型的研究相对缺乏。考虑到近年中国大豆进口需求的迅速增长，而美国是世界上重要的大豆生产国和中国大豆的进口来源[22]，本研究选取大豆作为作物估产研究对象，将中国吉林省与美国爱荷华州作为研究区域，在统一的估产技术路线基础上，对中美两国代表产区的
第3期徐春萌等：基于DSSAT作物模型的中美大豆主产区单产模拟与验证133
大豆单产进行研究，比较方法在不同地区的适用性。本研究通过使用大豆产量统计数据和DSSAT作物模型，并融合长时间序列的气象要素与MODIS中分辨率遥感影像，开展中美代表产区的大豆单产估测。通过在大面积范围内建立大豆单产估测模型，以期为全球重点区域的作物产量分析提供技术参考。
1 材料与方法
1.1 研究区概况
爱荷华州位于美国中北部黄金带，在北纬40°～44°，西经90°～97°之间，其耕地资源丰富，是美国大豆、玉米、小麦等经济及粮食作物的主产区之一。全州土地面积14.47万km2，大豆种植面积约4万km2。绝大部分地区地形平坦，从西北向南地势略微降低，水资源丰富[23]。爱荷华州属于温带大陆性气候，全年温差较大，降雨多集中在夏季且空间分布不均衡，大豆生育期间降雨总量为190～580 mm，日平均太阳辐射量为17～24 MJ/m2。
吉林省位于中国东北部黑土区，在北纬41°～47°，东经122°～131°之间，是中国重要的粮食生产地区之一。全省土地面积18.74万km2，大豆种植面积约2 000 km2。由东南向西北地势降低，东部为山地，
中西部为平原地貌，河流众多。吉林属于温带大陆性季风气候，四季分明，雨热同期，大豆生育期间降雨总量为150～900 mm，日平均太阳辐射量为14～23 MJ/m2。
1.2 数据来源
气象参数的模拟资料采用NASA世界能源预测项目（Prediction of Worldwide Energy Resources，简称POWER）所反演的降雨量、日最高温度以及太阳辐射数据，原始气象观测资料来源于POWER数据网的区域数据资料日值集（power.v/data-access-viewer/），空间分辨率为0.5°×0.5°。本研究使用ArcGIS10.3软件以格点数据为基础进行计算，分别获得爱荷华州和吉林省2008－2017年间每年6－8月的累积降雨量（mm），每年6－8月的日平均太阳辐射量（MJ/m2），以及每年8月的日最高温度平均值（℃），并采用双线性空间插值法将空间尺度细化至500 m×500 m。
爱荷华州的大豆分类数据来自美国农业部（United States Department of Agriculture，USDA）每年发布的农田数据层（Cropland Data Layers，CDL），其空间分辨率为30 m×30 m；使用基于2017年Sentinel-2遥感影像，辅助结合高分1号数据，提取中国吉林省的大豆种植分布数据，其空间分辨率为10 m×10 m，分类结果通过对全省各市县3 802个实地采样点进行验证，精度大于87%，以上数据由北京佳格天地科技有限公司与吉林省农委合作完成。研究中假设吉林省2008－2017年大豆种植区域不变，对其他年份进行统一的大豆分布图掩膜处理。以上大豆分类数据均重采样至500 m×500 m分辨率。
植被指数使用遥感影像MODIS MOD09A1数据集计算GCVI值，其空间分辨率为500 m×500 m。GCVI的计算式如下：
GCVI=NIR/G−1 （1）式中NIR和G指近红外光波段（光谱范围841～876 nm）、绿光波段（光谱范围545～565 nm）反射率。
土壤参数来自于国际土壤资料与信息中心（/）土壤数据集。根据DSSAT模型中大豆模型SOYGRO需要的基础土壤数据进行提取，包括土层深度、各层质地（黏粒、粉粒、砂粒含量）、pH值、有机质含量、全氮、容重以及持水性能等指标。
田间管理信息选取研究区域典型大豆品种作为研究对象，通过收集美国农业部发布的爱荷华州大豆种植数据（quickstats.v/）与吉林省6个气象站点（榆树、双阳、敦化、辽源、桦甸、延吉）记录的生育期数据，确定大豆各个生育期的日期。氮肥施用水平从50～300 kg/hm2不等。由于2个研究区域普遍不具备灌溉条件，因此在进行大范围估产的情况下将灌溉统一默认为无。
2008－2017年美国爱荷华州大豆单产数据来源于美国农业部调查的爱荷华州大豆产量历史数据（quickstats.v/），以县为统计单位；中国吉林省大豆产量数据来源于吉林省统计局公布的统计年鉴（v/tjsj/tjnj/），包含各市县的大豆播种面积和总产量，由此计算大豆单产数值。
1.3 单产模拟方法
windows server 20031.3.1 DSSAT模型调参
本研究选取DSSAT模型中的大豆模型SOYGRO作为模拟工具[24]。按照DSSAT模型的标准与格式，使用历史年份数据（2008－2017年）分别建立气象文件、土壤与田间管理文件、生育期及产量文件，并通过该模型中的最大似然估计（Generalized Likelihood Uncertainty Estimation，GLUE）模块对大豆遗传参数进行率定，分别对爱荷华州和吉林省构建了区域尺度的大豆模型参数，结果如表1所示。
利用历史年份（2008－2017年）数据，通过将大豆生长发育过程中的关键决定参数进行任意组合，包括种植日期、气象数据、氮肥水平播种密度、品种参数等，进而构建多场景的大豆生长情况，进行逐日模拟并绘制叶面积指数（Leaf Area Index, LAI）曲线（模拟起始日期为每年第100天，结束日期为第365天）。通过分析不同情境下LAI随时间的变化曲线，设置了2个28 d大豆生长窗口期：生育早期（每年第201天至第228天）及生育晚期（第229天至第256天）。利用大豆LAI与GCVI 的换算关系，如式（2）所示[25]，可将LAI随时间变化曲线转化为GCVI随时间变化的曲线。
GCVI=1.4LAI1.3+1.1 （2）结合历史年份（2008—2017年）的气象要素及遥感GCVI指数，获得研究区域内的大豆单产模拟向量系数表。向量系数表由14列（生育日期、气象要素、植被指数等相关参数），784行（28 d × 28 d）数据构成，可通过耦合大豆生长窗口期内的气象要素及遥感植被信息计
算得到大豆单位面积产量的估测值。为评价向量系数表
农业工程学报（）2021年134
的准确性，将2008—2017年气象站点的历史调查数据作
为单产真值与单产模拟值（通过向量系数表计算所得）进
行拟合，结果显示两者之间的线性回归决定系数R2大于
0.80。以上操作是将大豆模型SOYGRO链接到Python编
译器，通过Python调用SOYGRO的功能库完成。
卫慧作品表1 爱荷华州和吉林省的作物模型参数
Table 1 Parameters of crop model in Iowa State and Jilin
FL-SH、FL-SD、SD-PM、FL-LF分别为出芽到开花、开花到结荚、开花到鼓粒、鼓粒到成熟、开花到展叶的时间；LFMAX为最高光合速率；SLAVR 为比叶面积；SIZLF为三小叶期全叶最大尺寸；XFRT为最大种壳比；WTPSD 为每粒最高质量；SFDUR为种子灌浆时间；SDPUV为每荚平均种子数；PODUR为最终结荚天数；THRSH为最大脱净率（种子/（种子+种壳））；SDPRO为种子蛋白质含量；SDLIP为种子油质含量。
Note: CSDL, Critical Short Day Length below which reproductive development progresses with no daylength effect; PPSEN, slope of the relative response of development to photoperiod with time; EM-FL, FL-SH, FL-SD, SD-PM, FL-LF is time between plant emergence and flower appearance, first
flower and first pod, first flower and first seed, first seed and physiological maturity, time between first flower and end of leaf expansion, respectively; LFMAX is maximum leaf photosynthesis rate; SLAVR is specific leaf area of cultivar under standard growth conditions; SIZLF, maximum size of full leaf (three leaflets); XFRT, maximum fraction of daily growth that is partitioned to seed + shell; WTPSD, maximum weight per seed; SFDUR, seed filling duration for pod cohort; SDPUV, average seed per pod; PODUR, time required for cultivar to reach final pod load under optimal conditions; THRSH, threshing percentage, the maximum ratio of (seed/(seed+shell)) at maturity; SDPRO, fraction protein in seeds; SDLIP, fraction oil in seeds.
1.3.2 大豆单产模拟
将美国爱荷华州和中国吉林省作为研究区域，分别获取2008—2017年基于MODIS影像反演的逐年GCVI 植被指数，以及大豆生育期的气象要素，假设土壤性质和田间管理模式不变，利用研究区域所对应的大豆向量系数表，采用式（3）回归模型对大豆单产进行逐年模拟[26]。
Y=β0,d + β1,d·W + β2,d·RM d + β3,d·W·RM d （3）式中Y为大豆单产，kg/hm2；W为生育期内的天气特征（6－8月累积降雨量（mm），6－8月日平均太阳辐射量（MJ/m2），以及8月日最高温度（℃）向量，RM d为生育期内的遥感反演植被指数GCVI向量，系数β为向量系数表中大豆生育早期（第201天至第228天）与生育晚期（第229天至第256天）某一日期组合下的参数。
利用上述回归模型，以谷歌地球引擎（Google Earth Engine）作为遥感数据计算云平台，输入历史年份（2008－
2017年）的气象要素及遥感GCVI指数，分别进行2008－2017年美国爱荷华州和中国吉林省每年的大豆单产估测，输出单产模拟值的空间分辨率为500 m×500 m。
1.3.3 大豆单产精度验证
大豆单产精度验证所使用的数据集为2008－2017年爱荷华州各县（n=99）及吉林省各县（n=32～42）模拟大豆单产的平均值，验证集为中美两国政府发布的统计数据。对模拟大豆单产结果的误差分析，采用了3项计算指标进行评价：平均误差（Mean Percentage Error, MPE），均方根误差（Root Mean Square Error, RMSE），平均偏差（Mean Bias Error, BSE）。计算式如（4）～（6）所示。
1
[(||/)100%]
MPE
n
i i i
i
m d d
n
（4）
RMSE （5）
导航系统1
()
BSE
n
i i
i
m d
n
（6）式中n为统计县个数；m i为作物模型所得到的大豆单产，kg/hm2；d i为调查数据代表的大豆单产，kg/hm2。
同时，关于大豆产量的分县数据，在各年份作物模型的产量值与调查数据之间分别建立了线性回归模型，使用相关系数（r）验证两者之间的相关性。
2 结果与分析
2.1 美国爱荷华州大豆单产估测
美国爱荷华州2008－2017典型年份的大豆单产估测值空间分布如图1所示。从图中可以看出，爱荷华州大豆产量的空间变异性在不同年份存在差异：单产较高的大豆（>4 000 kg/hm2）在2008年分布于
东部和中部地区，在2011年主要位于东北地区，在2016年则分布于东部及西南地区；单产较低的大豆（<2 000 kg/hm2）主要分布于南部地区，但在2013年北部地区的大豆单产也较低。
美国爱荷华州2008－2017年大豆单产的年际变化如图2所示，模拟结果表明，大豆单产在2012年最低，仅为（2 139±193）kg/hm2，而在2010年大豆单产最高，为（4 766±970）kg/hm2，其余年份存在一定的产量波动。调查统计数据表明，爱荷华州大豆单产的最低值为2012年的（3 001±430）kg/hm2，而大豆单产最高值为2016年的（3 991±230）kg/hm2。从10 a间的单产波动来看，模拟值与调查值均呈现出大豆单产在2009年上升，于2011及2012年间下降，随后在2013至2016年上升，并在2017年再次下降的特征。美国爱荷华州大豆单产模拟值与调查数据的年际变化趋势较为一致，但作物模型估测的美国爱荷华州大豆单产的年际变化与统计调查数据相比，增减产的变化幅度较大，表现为单产预测值在高产年份比调查值更高，但在低产年份比调查值更低（图2）。
第3期 徐春萌等：基于DSSAT 作物模型的中美大豆主产区单产模拟与验证
东吴农村商业银行135
图1 美国爱荷华州典型年份大豆单产估测值空间分布
Fig.1 Spatial distribution of estimated soybean yields from typical years in Iowa State, USA
注：误差线是各县区单产的标准差。下同。
Note: Error bar represents the standard deviation among yield of counties. Same as below.
图2 美国爱荷华州2008－2017年大豆平均单产年际变化 Fig.2 Inter-annual variation of average soybean yields from 2008
to 2017 in Iowa State, USA
美国爱荷华州2008－2017年大豆单产估产的验证结果如表2所示。
表2 美国爱荷华州2008－2017年大豆单产模拟值与调查值误
差比较
Table 2 Errors of estimated soybean yields compared to investigated values from 2008 to 2017 in Iowa State, USA
年份Year n MPE /% RMSE /(kg·hm -2)
MBE /(kg·hm -2)
2008 99 22.0 441.6 274.7 2009 99 11.8 531.6 355.0 2010 99 40.7 1618.2 1 373.8 2011 99 8.7 369.5 −116.0 2012 99 27.8 920.9 −863.7 2013 99 15.5 539.1 −475.4 2014 99 9.1 383.4 102.9 2015 99 16.3 727.4 552.8 2016 99 12.5 627.1 233.2 2017
99
13.6
614.3
−371.2
平均值Average
16.8 762.8 107.2 注：n 为县的个数；MPE 为平均误差；RMSE 为均方根误差；MBE 为平均偏差。下同。
Note: n is the number of counties; MPE is mean bias error; RMSE is root mean square error; MPE is mean percentage error. Same as below.
综合来看，美国爱荷华州大豆模拟的单产区间略大于调查值，2008－2017大豆单产的模拟值为2 139～4 766
小电影qvod播放kg/hm 2
，而调查单产为3 002～3 991 kg/hm 2（图2），大豆单产估产的平均误差为8.7%～40.7%，估测值的均方根
误差为369.5～1 618.2 kg/hm 2
，平均偏差为−863.7～1373.9 kg/hm 2（表2）。其中，2010年模型估测的误差较大，大豆单产的估测值大于单产调查值1 373.8 kg/hm 2；而2011、2014年估测精确度较高，单产的平均偏差略高
于100 kg/hm 2。从10 a 平均来看，2008－2017年美国爱荷华州大豆单产估产的验证结果为：平均误差为16.8%，均方根误差为762.8 kg/hm 2，平均偏差为107.2 kg/hm 2，且估测值比调查值偏大（表2）。
对比2008－2017年美国爱荷华州大豆单产模拟数据与调查数据，在县级尺度上对估产结果的平均值进行了精度验证，典型年份线性回归模型的拟合结果如图3所示。遥感信息结合DSSAT 作物模型对美国爱荷华州2008－2017年大豆单产估测值与统计调查的产量值之间的相关系数（r ），从0.57到0.78不等。其中，相关程度较好的为2013年和2008年，最差的年份为2016年。整体来看，各县的大豆单产模拟值与调查值相关程度较好。当爱荷华州大豆单产较高时，回归模型的相关系数较小，当单产较低时，回归模型的相关系数较大，这可能与模型所适用的大豆产量区间有关。u活
注：r 为相关系数。
Note: r is correlation coefficient.
图3 美国爱荷华州典型年份大豆产量模拟值与调查值 Fig.3 Comparison of estimated with investigated soybean yields
from typical years in Iowa State, USA
2.2 中国吉林省大豆单产估测
中国吉林省2008－2017典型年份大豆单产估测值的空间分布如图4所示。从图中可以看出，吉林省大豆产量的空间变异性在不同年份略有差异。单产较高的大豆（＞3 000 kg/hm 2）在2012年分布在北部、中部和东部地区，在其它年份则主要分布于东部地区；而单产较低的大豆（＜2 000 kg/hm 2）除2012年以外主要分布于西部地区，部分东部地区单产也较低（图4）。
农业工程学报（ ） 2021年
136
图4 中国吉林省典型年份大豆产量估测值空间分布 Fig.4 Spatial distribution of estimated soybean yields from
typical years in Jilin Province, China
中国吉林省2008－2017年大豆单产的年际变化如图5所示。其中，吉林省大豆单产的估测数据显示，2009年的大豆单产最低，仅为（1 653±544）kg/hm 2，2017年大豆单产最高，为（2 766±933）kg/hm 2。从10 a 变化趋势来看，中国吉林省大豆单产模拟值在2009年显著下降，然后在2010年单产回升，直到2013年以后单产再次下降，并在2017年出现单产回升（图5）。吉林省的调查统计数据表明，大豆单产最低值为2009年的（1 997±667）kg/hm 2，最高值为2013年的（2 797±930）kg/hm 2。其中，大豆单产调查值在2009年下降，然后在2010年回升至稳定，直到2013年单产略微上升，并于2016年开始下降（图5）。综合来看，作物模型估测中国吉林省大豆单产的最低年份与调查值最低年份相吻合（2009年），且二者反映的高产年份相同（2008、2013、2015年）。除2017年模型认为产量上升以外，估测与调查所得大豆单产的年际变化波动趋势几乎一致。
图5 中国吉林省2008—2017年大豆平均单产年际变化 Fig.5 Inter-annual variation of average soybean yields from 2008
to 2017 in Jilin Province, China
中国吉林省2008－2017年运用遥感与作物模型计算所得的大豆单产估测值与调查值验证结果如表3所示。综合来看，2008－2017年中国吉林省大豆单产的模拟值为1 653～2 766 kg/hm 2，其产量区间非常接近于调查产量
1 997～
2 797 kg/hm 2
（图5），各年间大豆单产估产的平
均误差为24.5%～64.4%，估测单产的均方根误差为637.8～1 869.7 kg/hm 2，单产平均偏差为−655.2～534.1 kg/hm 2（表3）。其中，2015、2017年的大豆单产模型估测误差较大，而2008年的模型估测精确度较高，大豆单产的平均偏差仅为36.0 kg/hm 2。总结2008－2017年的10 a 平均值，大豆估产在中国吉林省的验证结果为：单产的平均误差为36.3%，均方根误差为1 088.4 kg/hm 2，平均偏差为−237.9 kg/hm 2，且估测值比调查值偏小（表3）。
表3 中国吉林省2008－2017年大豆产量模拟值与
调查值误差比较
Table 3 Errors of estimated soybean yields compared to investigated values from 2008 to 2017 in Jilin Province, China
年份Year n MPE /% RMSE /(kg·hm -2)
MBE /(kg·hm -2)
2008 32 24.8 837.4 36.0 2009 35 26.0 637.8 −293.2 2010 35 30.1 824.4 71.9 2011 36 31.4 1 132.7 −597.0 2012 38 24.5 735.3 −178.0 2013 42 31.6 1 056.2 −210.1 2014 38 40.5 1 460.7 −549.9 2015 36 51.5 1 869.7 −537.1 2016 37 38.6 1 262.2 −655.2 2017 37 64.4 1 067.7 534.1 平均值Average
36.3
1 088.4
−237.9
以县为统计单元，对比了中国吉林省2008－2017年间大豆单产的模拟值与调查数据，并对大豆估产结果在各县的平均值进行了精度验证，其中典型年份线性回归模型的拟合结果如图6所示。遥感信息融合DSSAT 作物模型对中国吉林省2008－2017年大豆单产估测值与统计调查值之间的相关系数（r ）在0.5～0.6，其中r 在2012年最高，且在2009、2016、2017年呈现出较好的相关性。整体来看，中国吉林省大豆单产模拟值与调查值的拟合程度较稳定，但其回归模型的相关系数相比于美国爱荷华州偏低。值得注意的是，吉林省模型的回归方程均位于1∶1线右下方，表现为大豆高产时模拟单产比调查单产小，低产时模拟单产比调查单产高（图6）。
图6 中国吉林省典型年份各县大豆单产模拟值与调查值 Fig.6 Comparison of estimated with investigated soybean yields
from typical years in Jilin Province, China
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议