首页> 中国专利> 一种基于云模型元胞自动机的城市扩展预测方法

一种基于云模型元胞自动机的城市扩展预测方法

摘要

本发明公开了一种基于云模型元胞自动机的城市扩展预测方法,包括以下步骤:1)栅格化城市土地利用图,每一栅格对应一个元胞;2)对各元胞转换概率变量及非城市元胞距各空间变量的空间距离采用概念集“近、中、远”表达,并采用云模型定量表达概念集;3)构造规则发生器;4)通过规则发生器得到空间变量影响下的元胞转换概率;5)制定非城市元胞的行为转变规则,并根据非城市元胞的行为转变规则计算每个元胞的最终元胞转换概率。本发明方法构建了不确定性与模拟结果精度之间的定量模型,并将上述模型应用于城市扩展预测中,有助于更真实、准确地模拟城市扩展。

著录项

  • 公开/公告号CN102254105A

    专利类型发明专利

  • 公开/公告日2011-11-23

    原文格式PDF

  • 申请/专利权人 武汉大学;

    申请/专利号CN201110227421.4

  • 发明设计人 王海军;贺三维;张文婷;何青青;

    申请日2011-08-10

  • 分类号G06F19/00(20060101);G06N3/00(20060101);

  • 代理机构武汉科皓知识产权代理事务所(特殊普通合伙);

  • 代理人张火春

  • 地址 430072 湖北省武汉市武昌区珞珈山武汉大学

  • 入库时间 2023-12-18 03:43:07

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2020-07-24

    未缴年费专利权终止 IPC(主分类):G06Q10/04 授权公告日:20130605 终止日期:20190810 申请日:20110810

    专利权的终止

  • 2013-06-05

    授权

    授权

  • 2012-01-04

    实质审查的生效 IPC(主分类):G06F19/00 申请日:20110810

    实质审查的生效

  • 2011-11-23

    公开

    公开

说明书

技术领域

本发明涉及城市扩展预测方法,尤其涉及一种基于云模型元胞自动机的城市扩展预测方法。

背景技术

城市扩展是元胞自动机(cellular automata,简称CA)应用的一个热点领域。Tobler院士在20世纪70年代首次将元胞自动机方法引入城市研究,并采用元胞自动机模拟美国五大湖边底特律地区城市的迅速扩展。但是,这时元胞自动机在城市扩展的应用还处于一般性讨论与方法“地理化”阶段,标准元胞自动机模型应用在城市扩展领域的性能和实用性有待进一步完善。许多元胞自动机扩展模型相继被提出并得到了更为广泛的发展。首先,多准则判断[1]、Logistic回归[2]等数理统计方法被应用于元胞自动机转换规则构建中;其次,近年来数据挖掘[3]、神经网络[4]、多智能体[5]、核学习机[6]、支持向量机[7]、人工免疫系统[8]等方法被用于转换规则的智能获取。

以黎夏、叶嘉安和周成虎为代表的学者们在元胞自动机扩展模型上做了很多有意义的理论尝试,但是这些元胞自动机扩展模型过分依赖高效的模型算法,不能从根本上解决城市扩展元胞自动机的真实性问题,很少考虑城市扩展中的不确定问题。一般元胞自动机采用蒙特卡洛方法[9]、模糊逻辑[10-11]来反映城市扩展的不确定性,但上述方法仅反映了城市扩展不确定性的随机性或模糊性,不能完整表达不确定性。Zhang等[12]采用位置和数量Kappa系数来分析空间明确模型中的预测不确定性和模拟的精度问题,该研究的实质是对模拟结果进行敏感性分析,从模型构建上并未体现不确定性定量传播问题。

文中涉及的参考文献如下:

[1] Wu F, Webster CJ.  Simulation of land development through the integration of cellular automata and multicriteria evaluation. Environment and Planning B: Planning and Design. 1998, 25(1):103-126.

[2] Wu F. Calibration of stochastic cellular automata: the application to rural-urban land conversions. International Journal of Geographical Information Science. 2002. 16(8):795-818.

[3] Li X, Yeh AGO. Data mining of cellular automata's transition rules. International Journal of Geographical Information Science. 2004. 18(8):723-744.

[4] Li X, Yeh A.G.O. Neural-network-based cellular automata for simulating multiple land use changes using GIS. International Journal Of Geographical Information Science: 2002. 16(4): 323-343.

[5] Chen YM, Li X, Liu XP, et al. An agent-based model for optimal land allocation (AgentLA) with a contiguity constraint. International Journal of Geographical Information Science. 2010.  24(8): 1269-1288.

[6] Liu X, Li X, Shi X, et al. Simulating complex urban development using kernel-based non-linear cellular automata. Ecological Modelling. 2008. 211(1-2):169-181.

[7] Yang QS, Li X, Shi X. Cellular automata for simulating land use changes based on support vector machines. Computers & Geosciences. 2008. 34:592–602.

[8] Liu XP, Li X, Shi X, et al. Simulating land-use dynamics under planning policies by integrating artificial immune systems with cellular automata. International Journal of Geographical Information Science. 2010.24(5): 783-802.

[9] Li D, Cheunq D, Shi XM. Uncertainty  Reasoning  Based  on  Cloud  Models  in  Controllers. Computers and Mathematics with Applications. 1998. 35(3):99-123.

[10] Al-Kheder S, Wang J, Shan J. Fuzzy inference guided cellular automata urban-growth modelling using multi-temporal satellite images. International Journal Of Geographical Information Science. 2008. 22(11-12):1271-1293.

[11] Al-Ahmadi K, See L, Heppenstall A, et al. Calibration of a fuzzy cellular automata model of urban dynamics in Saudi Arabia. Ecological Complexity. 2009.6:80–101.

[12] Zhang J, Zhou YK, Li RQ,et al. Accuracy assessments and uncertainty analysis of spatially explicit modeling for land use/cover change and urbanization: A case in Beijing metropolitan area. Science China Earth Sciences. 2010.53(2):173-180。

发明内容

针对现有技术存在的不足,本发明提供了一种基于云模型元胞自动机的城市扩展预测方法,该方法能更真实、准确地模拟城市扩展。

为了解决以上技术问题,本发明采用如下的技术方案:

一种基于云模型元胞自动机的城市扩展预测方法,包括以下步骤:

1)栅格化城市土地利用图,每一栅格对应一个元胞;

2)对各元胞转换概率变量及非城市元胞距各空间变量的空间距离采用概念集“近、中、远”表达,并采用云模型定量表达概念集,实现定性到定量的表达,所述的空间变量为城市中心、镇中心、工业中心和主干道;所述的元胞转换概率变量值在0~1之间,并等分为三个区间,分别对应“高、中、低”概念;

3)构造规则发生器,本步骤进一步包括以下子步骤:

3-1构建云推理引擎,即定性规则库,采用“If A,then B”形式表达,其中,A为各非城市元胞距各空间变量的空间距离所属的概念集;B为在空间变量影响下的元胞转换概率概念集;

3-2根据所构建的定性规则库,采用多维云模型构造对应于If部分的前件规则发生器,并采用二维云模型构造对应于then部分的后件规则发生器;

4)根据判断各非城市元胞所属的概念,通过规则发生器得到空间变量影响下的元胞转换概率,本步骤进一步包括以下子步骤:

4-1将各非城市元胞距各空间变量的空间距离值输入至前件规则发生器,采用极大判定法判断各空间距离值属于概念集{远、中、近}中的某个概念,并激活前件规则发生器中对应的定性规则库,通过前件规则发生器推理引擎,输出中间隶属度值;

4-2将上述中间隶属度值作为后件规则发生器的输入值,激活对应的后件规则发生器,通过后件规则发生器推理引擎,得到各非城市元胞对于某一概念的最终隶属度值,即,非城市元胞城市化的转换概率;

5)考虑邻域、宏观因素,制定非城市元胞的行为转变规则,并根据非城市元胞的行为转变规则计算每个元胞的最终元胞转换概率,若元胞转换概率大于预设阈值,则该非城市元胞所代表的土地城市化;否则,该非城市元胞状态不变;所述预设阈值根据元胞转换概率值分布范围和分布特征确定。

上述步骤2)进一步包括以下子步骤:

2-1获取各非城市元胞距各空间变量的空间距离值,根据各非城市元胞距各空间变量的空间距离构建空间距离矩阵                                                ,其中,取1、2、3、4,分别为各非城市元胞距城市中心、镇中心、工业中心、主干道的空间距离矩阵; 

2-2将各非城市元胞距各空间变量的空间距离值离散化、概念化;

2-3从空间距离矩阵中各随机抽取20%的元胞;

2-4 对步骤2-3所抽取的各元胞距各空间变量的空间距离值进行K均值聚类分析,将空间距离矩阵的聚类中心分为三类,将空间距离矩阵对应的论域记为,根据聚类中心将划分为三个离散区间:,其中,

2-5将空间距离表达为三个概念粒度,即远、中、近;将由空间变量决定的元胞转换概率记为,其对应的论域记为,将抽象为概念集{高、中、低};

2-6利用云模型表示空间距离矩阵的概念集合,其中,分别表示基本概念“远”、“中”、“近”;期望值,分别对应论域的中心值,这里,期望值;熵,根据“3δ”原则,为论域中心值,为论域最小值;超熵决定着云模型的离散程度,控制着模型输入与输出的不确定性关系;

2-7 根据步骤1-5得到云模型的数学特征值、熵的数学公式:

上述步骤2-1中利用ArcGIS工具获取各元胞距各空间变量的空间距离值。

上述步骤2-3中采用随机分层法从空间距离矩阵中各随机抽取20%的元胞。

上述步骤3-2中采用多维云模型构建前件规则发生器进一步包括以下步骤:

3-2-1假设变量构成的多维云模型的数学特征为分别为变量的期望值、熵、超熵,其中,变量指各非城市元胞距各空间向量的空间距离值;

3-2-2以为期望值,为方差,生成维正态随机数

3-2-3以为期望值,以的绝对值为方差,生成具有稳定性的随机数

3-2-4计算输出值,产生一个云滴

3-2-5重复步骤2-1-2 ~2-1-4,至产生N个云滴,N可取任意值,作为优选,N取2000。 

上述步骤5)中制定元胞的行为转换规则进一步包括以下子步骤:

5-1采用四方网格对元胞空间进行划分,元胞空间中元胞邻域采用半径R=4的圆型邻居;

5-2对在空间变量驱动下的元胞转换概率修正为,其中,为半径R=4的圆型邻域对中心元胞的局部转换概率,表示中心元胞受宏观约束的转换概率,用0、1表示,0表示不能发展或发展阻力较大,1表示具有较大的发展潜力或无发展阻力,若有一项取值为0,则

上述可表示为:,其中,为中心元胞的邻域内第个邻居元胞的状态;为判断是否为城市元胞,若为真,则返值1,若为假,则返值0;为邻域内元胞总数。

 

与现有技术相比,本发明具有以下优点和有益效果:

本发明方法构建了不确定性与模拟结果精度之间的定量模型,并将上述模型应用于城市扩展预测中,有助于更真实、准确地模拟城市扩展。

附图说明

    图1为实施例中空间距离的云模型;

图2为多维云模型;

图3为实施例中采用规则发生器计算隶属度的过程示意图;

图4为本实施例中值与对应预测图所计算的Kappa系数的关系图;

图5为不同超熵He值下的结果模拟图对比;

图6为超熵He值与新增元胞数量之间的关系;

图7为采用本发明方法预测的2005年和2007年江夏区城市扩张图及实际图的对比,(a)为2005年江夏区城市扩展实际图,(b)为采用本发明方法预测的2005年江夏区城市扩展图;(c)为2007年江夏区城市扩展实际图,(d)为采用本发明方法预测的2007年江夏区城市扩展图;

图8为2010年江夏区城市扩展预测图;

图9为2020年江夏区城市扩展预测图;

图10为2002-2007年三种不确定CA模型的模拟结果对比,(a)为采用本发明方法得到的2005年江夏区城市扩展预测图;(b)为采用模糊CA法得到的2005年江夏区城市扩展预测图;(c)为采用Monte Carlo法得到的2005年江夏区城市扩展预测图;(d)为采用本发明方法得到的2007年江夏区城市扩展预测图;(e)为采用模糊CA法得到的2007年江夏区城市扩展预测图;(f)为采用Monte Carlo法得到的2007年江夏区城市扩展预测图。

具体实施方式

下面将结合附图和具体实施例对本发明技术方案的优势做进一步说明。

本发明所提出的一种基于云模型元胞自动机的城市扩展预测方法,包括以下步骤:

一、栅格化城市土地利用图,每一个栅格对应一个元胞。

二、对各元胞转换概率变量(该变量值在0~1之间,该变量一开始并不需要求取,知道其值的变化区间[0, 1],等分为三个区间即可,分别对应高中低即可)及非城市元胞距各空间变量的空间距离采用概念集“近、中、远”表达,并采用云模型定量表达概念集,实现定性到定量的表达树,所述的空间变量为城市中心、镇中心、工业中心和主干道;所述的元胞转换概率变量值在0~1之间,并等分为三个区间,分别对应“高、中、低”概念;

本步骤的一种具体实施方式如下:

2-1利用ArcGIS或Matlab工具获取各非城市元胞距各空间变量的空间距离值,根据各非城市元胞距各空间变量的空间距离构建空间距离矩阵,其中,取1、2、3、4,分别为各非城市元胞距城市中心、镇中心、工业中心、主干道的空间距离矩阵; 

2-2将各非城市元胞距各空间变量的空间距离值离散化、概念化;

2-3采用随机分层法从空间距离矩阵中各随机抽取20%的元胞;

2-4对步骤2-3所抽取的各元胞距各空间变量的空间距离值进行K均值聚类分析,将空间距离矩阵的聚类中心分为三类,将空间距离矩阵对应的论域记为,根据聚类中心将划分为三个离散区间:,其中,

2-5根据人类的认知习惯,可将距离表达为三个概念粒度,即远、中、近;同理,将由空间变量决定的元胞转换概率记为,其对应的论域记为,将抽象为概念集{高、中、低};

2-6利用云模型表示空间距离矩阵的概念集合,其中,分别表示基本概念“远”、“中”、“近”;期望值,分别对应论域的中心值,这里,期望值;熵,根据“3δ”原则,为论域中心值,为论域最小值;超熵决定着云模型的离散程度,控制着模型输入与输出的不确定性关系;

2-7 根据步骤1-5得到云模型的数学特征值、熵的数学公式:

2-6 为了处理定性概念中存在的不确定性,包括模糊性和随机性,这里采用云模型来表示各个定向概念,并用云模型构造一种具有不确定性的泛概念树。定性概念的集合可定义为由论域上的基本概念组成,即利用云模型表示空间距离矩阵的概念集合,其中,分别表示基本概念“远”、“中”、“近”;期望值,分别对应论域的中心值,这里,期望值;熵,根据“3δ”原则,为论域中心值,为论域最小值;超熵决定着云模型的离散程度,控制着模型输入与输出的不确定性关系;

2-7 根据步骤1-5得到云模型的数学特征的数学公式:

按照上述云模型的数学特征的数学公式可生成空间距离的云模型,如图1所示,从图中可知,用云模型表达的各个概念之间的区分是软性的,允许一定的交叠。这与模糊理论有相似之处,但是由云模型建立的论域输入与输出之间为一对多的映射关系,从而有力结合了模糊性和随机性,完整表达了城市扩展中驱动因素与元胞转化概率之间的必然性和不确定性关系。

三、构造规则发生器,本步骤进一步包括以下子步骤:

3)构造规则发生器,本步骤进一步包括以下子步骤:

3-1构建云推理引擎,即定性规则库,采用“If A,then B”形式表达,其中,A为各非城市元胞距各空间变量的空间距离所属的概念集;B为在空间变量影响下的元胞转换概率概念集;由于在城市扩展研究问题中,有四个空间变量,即城市中心、镇中心、工业中心和主干道,四个空间变量对应四个概念集,因此采用多规则知识库,具体为If A1 and A2 … and A4, then B,其中,表示第个空间变量的定性状态,即{远,中,近}中的一个,;B表示元胞转换概率的定性状态,即{高、中、低}中的一个。

3-2根据所构建的定性规则库,构建对应于If部分的前件规则发生器和对应于then部分的后件规则发生器。对于If部分,由于有4个空间变量,所以选择多维云模型构建前件规则发生器,可以避免各空间变量赋权的复杂问题;对于then部分,由于只有元胞转换概率一个变量,所以选择二维云模型构建后件规则发生器。图2所示为以非城市元胞与城市中心距离和与工业中心距离为例,生成的多维云模型。

本步骤中采用多维云模型构建前件规则器的一种具体实施方式如下:

3-2-1假设变量构成的多维云模型的数学特征为分别为变量的期望值、熵、超熵;

3-2-2以为期望值,为方差,生成维正态随机数

3-2-3以为期望值,以的绝对值为方差,生成具有稳定性的随机数

3-2-4计算输出值,产生一个云滴

3-2-5重复步骤2-1-2 ~2-1-4,至产生N个云滴。 

四、根据判断各非城市元胞所属的概念,通过规则发生器得到空间变量影响下的元胞转换概率,本步骤进一步包括以下子步骤:

4-1将各非城市元胞距各空间变量的空间距离值输入至前件规则发生器,采用极大判定法判断各空间距离值属于概念集{远、中、近}中的某个概念,并激活前件规则发生器中对应的定性规则库,通过前件规则发生器推理引擎,输出与输入对应的确定度,即,中间隶属度值;

4-2将步骤4-1得到的确定度作为后件规则发生器的输入值,激活对应的后件规则发生器,通过后件规则发生器推理引擎,得到随机值为在一定区间范围内波动的随机值,即,各元胞对于某一概念的最终隶属度值,即为非城市元胞城市化的转换概率。

图3所示为采用规则发生器计算隶属度的过程示意图。

五、考虑邻域、宏观因素,制定非城市元胞的行为转变规则,并根据非城市元胞的行为转变规则计算每个元胞的最终元胞转换概率,若元胞转换概率大于上预设阈值,则该非城市元胞所代表的土地城市化;否则,该非城市元胞状态不变;所述预设阈值根据元胞转换概率值分布范围和分布特征确定:

本步骤的一种具体实施方式如下:

5-1为了便于计算机表达,采用四方网格对元胞空间进行划分;为了克服传统元胞邻域的局限性,元胞空间中元胞邻域采用半径R=4的圆型邻居;

5-2除了考虑空间变量的驱动因素外,转换规则的构建也需考虑邻域、宏观因素的综合影响。将在空间变量驱动下的元胞转换概率修正为,其中,为半径R=4的圆型邻域对中心元胞的局部转换概率,表示中心元胞受宏观约束的转换概率,用0、1表示,0表示不能发展或发展阻力较大,1表示具有较大的发展潜力或无发展阻力。若规划用途为基本农田、自然保护区等,则取0,否则取1;若有山体、水体等自然障碍存在,则取0,否则取1。若有一项取值为0,则=0。

上述可表示为:,其中,为中心元胞的邻域内第个邻居元胞的状态;为判断是否为城市元胞,若为真,则返值1,若为假,则返值0;为邻域内元胞总数。

 

下面将以武汉市江夏区2002年的土地利用现状图为基础,采用本发明方法来预测江夏区在2002-2005年、2005-2007年两个时间段的城市扩展情况。

一、采用本发明方法预测江夏区在2002-2005年、2005-2007年两个时间段的城市扩展情况:

在实施本方法前,需进行云模型参数的设置,参数的设置直接影响到预测结果的好坏。前件规则发生器为五维云模型,江夏区在2002-2005年、2005-2007年两个时间段的前件规则发生器和后件规则发生器的数学特征如表1所示。超熵用于控制不确定性的传播,为降低复杂度,本实施例均采用统一的值,其取值区间为[0,1],在区间[0,1]内每隔0.05对进行取值,共有21次取值,即,分别预测2005年江夏区城市扩展情景,并采用Kappa系数对预测图和实际图进行精度评价。Kappa系数是点对点评价中常用的方法,其计算公式如下:,其中,为观测精度估计;为期望精度估计。图4所示为本实施例中值与对应预测图所计算的Kappa系数的关系。

表1 江夏区云模型参数设置表

从图4中易看出,随着值的增大,Kappa系数的大体变化趋势为先上升后下降。总体而言,以为临界值,当,Kappa系数值均低于Kappa系数的平均值,而当时,Kappa均大于均值,且在时Kappa系数达到极大点。这说明了模拟图的精度受值的影响,由于值的大小决定着模型的不确定性程度,即模糊性和随机性的程度,所以这说明不确定性是影响模型精度的主要因素,不确定性程度的提高有助于提高模拟精度,因此在模型中考虑不确定性因素的影响是十分必要的。

从不同值对应的模拟结果来看,不确定性主要存在于城市发展边缘,图5为本实施例不同值下的结果模拟图对比。随着值的上升,不确定性程度也逐渐增大,城市边缘元胞发展的模糊性和随机性逐渐增加,表现为城市边缘元胞越难发展为城市元胞。这与黎夏的研究结论,即CA模拟结果的不确定性主要体现在城市的边缘,是相符的。随着值的增加,即不确定性程度增加,模拟图的新增城市元胞数量大体呈下降的趋势,如图6所示。

时,Kappa系数最高,说明该值对应下的预测精度最佳。采用此值,分别模拟2005年和2007年江夏区的城市发展状况,如图7所示,Kappa系数分别为94.16%和85.44%,预测精度高,说明采用本发明方法可以很好的预测城市扩展。

 

以2005年城市发展现状图为基点,模型运行N个周期后与2007年实际城市发展图的匹配精度最高,则推算出模型每运行一个周期代表的实际年数为N/2。采用Cloud-CA进行预测时,以2007年城市发展现状图为基点,以模拟2007年Cloud-CA的参数进行设置,假定模型的运行周期代表的实际年数为N/2。模型运行若干周期后,得到2010年和2020年的城市发展预测图,分别如图8和图9所示。

二、本发明方法和蒙特卡洛方法(Monte Carlo)用于预测江夏区城市扩展结果的对比

    从现有研究看,CA模型中不确定性的表达始终不够完整,不能有机地整体考虑模糊性和随机性问题,全面表达不确定性的两个方面,即复杂事物及模型均具有模糊性和随机性。Cloud-CA有力地集成了概念的模糊性和随机性,完整表达了城市扩展中不确定性问题;模糊CA主要反映了城市扩展中不确定性的模糊性;采用Monte Carlo的CA主要反映了城市扩展中不确定性的随机性。

从图10所示的两组模拟结果对比可知,在同等参数设置的条件下,Cloud-CA的不确定性程度高于模糊元胞自动机,且主要表现在城市扩展的边缘区域。原因主要是Cloud-CA集成了模糊元胞自动机表达的模糊性和Monte Carlo方法的随机性,模糊性和随机性的结合增加了模型的不确定性表达,故不确定性程度增加,新增城市元胞数量比模糊CA有所减少。而不确定性主要集中在城市边缘这一现象是与云模型的原理密不可分的。将云模型引入CA中,元胞的不确定性程度与距离变量(如与城市中心、镇中心、工业中心、主干道等的距离)存在负相关关系,即距离变量越小,说明越靠近城市中心、镇中心、工业中心、主干道等城市用地,该元胞的地理区位优势越明显,多维云模型会产生高的隶属度值。根据云模型的云滴分布状况,当输入的隶属度越高时,输出值(即元胞转换概率)的云滴随机性越低,不确定性程度降低,反之亦然。因此,云模型能有效地表达城市扩展的不确定性分布特征。

而对比Cloud-CA与Monte Carlo方法,易发现Monte Carlo方法的不确定性程度高于Cloud-CA,不确定性的分布不仅在城市边缘也分布在城市内部,造成城市被“蜂窝洞”填满。虽然Monte Carlo方法被广泛运用,用以反映随机因素对城市模拟的影响。但是从效果上看,Monte Carlo方法属于无约束条件的干扰,城市内部各处受干扰的概率一致,造成城市内部“蜂窝状”,不符合城市土地集约利用需要和城市紧凑式发展的现状。Cloud-CA将干扰概率与空间距离约束条件结合起来,区别了城市内部和城市边缘的不确定性分布规律,完整表达了不确定性的内涵。采用Kappa系数来衡量模拟结果的精度,三种CA模型的模拟精度如表2所示,对比可知Cloud-CA的模拟精度均高于模糊CA和Monte Carlo方法。

表2 三种CA模型2002-2007年模拟结果精度对比

Kappa系数云模型元胞自动机模糊元胞自动机蒙特卡洛方法2005年94.16%93.65%82.16%2007年85.44%84.12%80.66%

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号