位置: IT常识 - 正文

DEMATEL-ISM模型的Python实现——方法介绍以及代码复现(dematel模型有什么缺点)

编辑:rootadmin
DEMATEL-ISM模型的Python实现——方法介绍以及代码复现 DEMATEL-ISM模型的Python实现——方法介绍以及代码复现前言DEMATEL-ISM分析方法方法简介步骤明确系统要素确定直接影响矩阵规范影响矩阵计算综合影响矩阵计算各个要素的影响度、被影响度、中心度和原因度绘制因果图确定整体影响矩阵确定可达矩阵划分层级绘制因素之间的递阶层次结构实例与代码DEMATEL-ISM计算原因-结果图节点度图ISM图后记前言

推荐整理分享DEMATEL-ISM模型的Python实现——方法介绍以及代码复现(dematel模型有什么缺点),希望有所帮助,仅作参考,欢迎阅读内容。

文章相关热门搜索词:ism模型步骤,dematel模型有什么缺点,dea-sbm模型,ism模型步骤,dea-sbm模型,m/m/s模型,dematel-ism模型的缺点,m-shel模型,内容如对您有帮助,希望把文章链接给更多的朋友!

本文源于笔者的《系统工程》课程的小组作业,笔者尝试运用DEMATEL-ISM方法来进行分析,建模求解,但在网络上并没有找到相应的,特别是集合DEMATEL-ISM方法的代码。因此自己码了DEMATEL-ISM模型的Python代码,并作为第一个博客发布~

参考文献中,笔者主要参考了李广利等1的研究,本文也将依此论文进行方法解读和代码复现。

网上查找资料的过程中,笔者发现了一个MCDA方法的python代码库2,里面有很多多准则决策分析模型的相关代码,其中就有DEMATEL的代码,笔者做了一定的参考。

DEMATEL-ISM分析方法方法简介

DEMATEL(Decision Making Trial and Evaluation Laboratory),全称为“决策试验和评价实验法”,是一种运用图论与矩阵工具进行系统要素分析的方法,通过分析系统中各要素之间的逻辑关系与直接影响关系,可以判断要素之间关系的有无及其强弱评价。ISM(Interpretative Structural Modelling)法全称为“解释结构模型”,其特点是把复杂的系统分解为若干子系统(要素),通过代数运算将系统构造成一个多级递阶的结构模型。

DEMATEL 模型可利用矩阵运算求出因素间的因果关系和影响强度,通过可视化因素间的因果关 系,得以揭示复杂问题中的关键影响因素及影响程度;但该方法无法有效识别系统中因素的层级结 构。ISM 法则通过分析构成系统的各子系统( 因素或要素) 之间的直接二元相关关系,基于布尔代数运算等,构造多级递阶有向拓扑图,但无法确定要素对系统的影响程度。

将两种方法结合,可以识别系统中关键要素及其影响程度,并构建要素的层级结构。DEMATEL-ISM方法的过程如下:

步骤明确系统要素

明确分析系统所构成的要素,将构成系统的要素标记为x1x_1x1​, x2x_2x2​, x3x_3x3​, …\ldots…, xnx_nxn​。

确定直接影响矩阵

采用专家打分法,比较xix_ixi​对xjx_jxj​的影响,由于因素与自身比较为没有影响,直接影响矩阵的对角线值为0。通过比较得到直接影响矩阵AAA。 A=[x12⋯x1nx21⋯x2n⋮⋮⋱⋮xm1xm2…]\begin{align} A= \begin{bmatrix} 0&x_{12}&\cdots&x_{1n}\\ x_{21}&0&\cdots&x_{2n}\\ \vdots &\vdots&\ddots &\vdots \\ x_{m1}&x_{m2}&\dots &0 \end{bmatrix} \end{align}A=​0x21​⋮xm1​​x12​0⋮xm2​​⋯⋯⋱…​x1n​x2n​⋮0​​​​ 式中因素xij(i=1,2,…,m;j=1,2,…,n;i≠j)x_{ij}(i=1,2, \ldots, m;j=1,2, \ldots,n;i\neq j)xij​(i=1,2,…,m;j=1,2,…,n;i=j)表示因素xix_ixi​对xjx_jxj​的直接影响。i=ji=ji=j时,xij=x_{ij}=0xij​=0。

规范影响矩阵

归一化原始关系矩阵得到规范影响矩阵。归一化方法有很多种,这里可以采用行最大值法进行归一化,即将矩阵AAA的每一行求和,在这些值中取最大值,将矩阵AAA中元素除以最大值,得到规范直接影响矩阵BBB。 B=xijmax(∑j=1nxij)\begin{align} B=\frac{x_{ij}}{max(\sum\limits_{j=1}^n x_{ij})} \end{align}B=max(j=1∑n​xij​)xij​​​​

计算综合影响矩阵

综合系统矩阵体现系统中各个元素间的影响的综合效应,规范影响矩阵不断自乘后矩阵的所有值会趋近于0,即lim⁡k→∞Bk=\lim\limits_{k \to \infty} B^k=0k→∞lim​Bk=0。计算综合影响矩阵时,得到 T=(B+B2+⋯+Bk)=∑k=1∞Bk=B(I−B)−1\begin{align} T=(B+B^2+\dots +B^k)=\sum\limits_{k=1}^\infty B^k=B(I-B)^{-1} \end{align}T=(B+B2+⋯+Bk)=k=1∑∞​Bk=B(I−B)−1​​ 式中为III单位矩阵。

计算各个要素的影响度、被影响度、中心度和原因度

影响度指矩阵TTT中各行值和,表示各行要素对其他所有要素的综合影响值,记作DiD_iDi​,有: Di=∑j=1nxij,(i=1,2,…,n)\begin{align} D_i=\sum\limits_{j=1}^n x_{ij},(i=1,2,\dots,n) \end{align}Di​=j=1∑n​xij​,(i=1,2,…,n)​​ 被影响度指矩阵TTT中各列值和,表示各列要素对其他所有要素的综合影响值,记作CiC_iCi​,有: Ci=∑j=1nxji,(i=1,2,…,n)\begin{align} C_i=\sum\limits_{j=1}^n x_{ji},(i=1,2,\dots,n) \end{align}Ci​=j=1∑n​xji​,(i=1,2,…,n)​​ 中心度表示因素在评价体系中的位置及其所起作用的大小,某要素的中心度为其影响度于被影响度之和,记作MiM_iMi​,有: Mi=Di+Ci\begin{align} M_i=D_i+C_i \end{align}Mi​=Di​+Ci​​​ 原因度由某要素的影响度和被影响度相减得到,记作,有: Ri=Di−Ci\begin{align} R_i=D_i-C_i \end{align}Ri​=Di​−Ci​​​ 当原因度大于0时,表示该要素对其他要素的影响程度大,称其为原因要素,当原因度小于0时,该要素为结果要素。

绘制因果图DEMATEL-ISM模型的Python实现——方法介绍以及代码复现(dematel模型有什么缺点)

将中心度MiM_iMi​为横坐标,将原因度RiR_iRi​作为纵坐标,绘制因果关系图。该图可以直观体现因果关系。

确定整体影响矩阵

整体影响矩阵由综合影响矩阵和单位矩阵相加得到。

确定可达矩阵

确定可达矩阵时需要引入一个阈值λ\lambdaλ以剔除因素之间影响程度较小的关系,从而明确层级结构的划分。引入阈值λ\lambdaλ对整体影响矩阵EEE进行处理,可以得到可达矩阵FFF。则有: fij={1eij≥λ(i,j=1,2,…,n)eij<λ(i,j=1,2,…,n)\begin{equation} f_{ij}= \begin{cases} 1 & & {e_{ij}\geq\lambda(i,j=1,2,\dots,n)}\\ 0 & & {e_{ij}<\lambda(i,j=1,2,\dots,n)} \end{cases} \end{equation}fij​={10​​eij​≥λ(i,j=1,2,…,n)eij​<λ(i,j=1,2,…,n)​​​

划分层级

由可达矩阵FFF的第iii行上值为1的列对应的因素求得可达集R(x1)={xi∣Fij=1}R(x_1)=\{x_i|F_{ij} =1\}R(x1​)={xi​∣Fij​=1},表示从因素xix_ixi​出发可以达到的全部因素的合集,由可达矩阵FFF的第iii列上值为1的行对应的因素求得先行集S(x1)={xi∣Fji=1}S(x_1)=\{x_i|F_{ji} =1\}S(x1​)={xi​∣Fji​=1},表示可以达到因素的全部因素的合集。如果R(x1)R(x_1)R(x1​)和S(x1)S(x_1)S(x1​)满足R(x1)∩S(x1)=R(x1)R(x_1)\cap S(x_1)=R(x_1)R(x1​)∩S(x1​)=R(x1​),则表示中对应的元素均能在中找到前因,将该元素称为高层级的元素。然后从可达矩阵中去除对应的行和列,再从矩阵中抽取最高级的因素,不断重复该过程,直到所有的行和列均被去除。

绘制因素之间的递阶层次结构

根据去除因素的顺序,绘制系统要素间多级递阶有向拓扑图。

实例与代码DEMATEL-ISM计算

这里笔者直接给出一组数据,10项因素,直接影响矩阵AAA如下。

下面是DEMETAL-ISM计算的Python代码:

#导入所需库import numpy as npimport pandas as pd#直接影响矩阵AA = np.array([[0, 2, 0, 2, 3, 2, 3, 2, 3, 2], [0, 0, 0, 1, 3, 2, 0, 1, 1, 2], [0, 1, 0, 0, 2, 1, 0, 0, 0, 1], [0, 2, 0, 0, 3, 2, 0, 0, 0, 1], [1, 3, 0, 2, 0, 1, 1, 0, 1, 1], [0, 2, 0, 2, 1, 0, 1, 1, 0, 2], [2, 2, 0, 0, 2, 1, 0, 1, 2, 1], [0, 2, 0, 0, 1, 1, 1, 0, 1, 2], [3, 2, 0, 1, 2, 2, 2, 2, 0, 1], [0, 2, 0, 0, 1, 2, 1, 1, 2, 0]])#行和最大值归一化,得到规范影响矩阵Brow_sum = np.sum(A, axis = 1)max_sum = np.max(row_sum)B = A/max_sum#综合影响矩阵TT = np.matmul (B, np.linalg.inv(np.identity(A.shape[0]) - B))#计算影响度D,被影响度C,中心度M,原因度RD = np.sum(T, axis = 1)C = np.sum(T, axis = 0)M = D + CR = D - C#判断输出结果要素和原因要素causal_factors = ''result_factors = ''for i in range(len(R)): if R[i] > 0: causal_factors = causal_factors + 'x' + str(i + 1) + ' ' elif R[i] < 0: result_factors = result_factors + 'x' + str(i + 1) + ' 'print("原因要素:" + causal_factors)print("结果要素:" + result_factors)#使用相应的lambda值计算可达矩阵F,并输出节点度排序E = T + np.identity(A.shape[0])lbd = 0.20 # 相应的lambda值F = Efor i in range(A.shape[0]): for j in range(A.shape[1]): if E[i, j] > lbd: F[i,j] = 1 elif E[i, j] < lbd: F[i, j] = 0node_degree = np.sum(F, axis = 0) + np.sum(F, axis = 1)print(sorted(node_degree, reverse=True))

上面代码计算出的一些矩阵和值比较重要,如综合影响矩阵TTT,被影响度CCC、影响度DDD、原因度MMM、中心度RRR等,数据在下面的画图也有用的,笔者在此没有输出,在做研究的过程中可以用pandas库将相应的值输出到excel表中。 λ\lambdaλ的取值也很关键,具体取何值下文有说明,这里取0.20,得到最后的可达矩阵如下图所示。

原因-结果图

接下来运用上文计算的结果,绘制原因-结果图,这张图可以看到各个因子之间互相影响的大小,以及对研究内容影响的大小,如下图所示。 第1类为强原因因子集(第Ⅰ区),这类因子对研究对象的形成具有非常显著的影响,且对其他结果型因子有较大的影响;第2类为弱原因因子集(第Ⅱ区),这类因子对研究对象的形成也具有重要影响,对其他结果型因子也有一定的影响;第3类为弱结果因子集(第Ⅲ区),这类因子是其他原因型因子综合作用的结果,对研究对象的形成具有一定的影响;第4类为强结果因子集(第Ⅳ区),这类因子也是其他原因型因子综合作用的结果,但是对研究对象的形成具有非常重要的影响。DEMATEL计算所得因子中心度越大,其影响程度也就越大,因此,需要重点关注第Ⅰ区的强原因因子集及第Ⅳ区的强结果因子集。

下面是原因-结果图的代码,为了图形美观,笔者代码水平也很有限,所以做了大量的调整,大家可以看着自己来。

import matplotlib as mplimport matplotlib.pyplot as plt#设置图形格式config = { "font.family": 'serif', "font.size": 14, "mathtext.fontset": 'stix', "font.serif": ['SimSun'], 'axes.unicode_minus': False }mpl.rcParams.update(config)#x轴为上文计算出的中心度x = [2.9456079466124003, 3.243636237536544, 0.5565493379272154, 1.9429001271455832, 3.205862930468748, 2.627384244237508, 2.4431692291267364, 1.891879956077375, 3.017054245988732, 2.609503379763571]y = [1.5128414543332274, -0.9939990519336557, 0.5565493379272154, -0.1661886145842798, -0.8173828786366384, -0.6568575842419181, 0.3753242725311403, -0.015046149922503416, 0.6625117023194962, -0.4577524877920831]#因子名factors_name = [r'$x_1$', r'$x_2$', r'$x_3$', r'$x_4$', r'$x_5$', r'$x_6$', r'$x_7$', r'$x_8$', r'$x_9$', r'$x_{10}$']#画散点图,并增加相应名称、线段和调整大小和位置plt.scatter(x, y, s=3, c='k')plt.xlabel('中心度')plt.ylabel('原因度')for i in range(len(x)): if i == 0: plt.text(x[i]+0.025, y[i]-0.085, factors_name[i], fontsize=17) elif i == 3: plt.text(x[i]+0.025, y[i]-0.1, factors_name[i], fontsize=17) elif i == 6: plt.text(x[i]-0.18, y[i]-0.1, factors_name[i], fontsize=17) elif i == 7: plt.text(x[i]-0.15, y[i]-0.12, factors_name[i], fontsize=17) else: plt.text(x[i]+0.025, y[i]+0.025, factors_name[i], fontsize=17)plt.vlines(sum(x)/len(x), -1.55, 1.55, colors='k', linestyles='dashed')plt.hlines(0, sum(x)/len(x)-2, sum(x)/len(x)+2, colors='k', linestyles='dashed')plt.xlim(sum(x)/len(x)-2, sum(x)/len(x)+2)plt.ylim(-1.55, 1.55)plt.text(sum(x)/len(x)+2-0.2, 1.4-0.1, 'Ⅰ')plt.text(sum(x)/len(x)-2+0.1, 1.4-0.1, 'Ⅱ')plt.text(sum(x)/len(x)-2+0.1, -1.55+0.1, 'Ⅲ')plt.text(sum(x)/len(x)+2-0.2, -1.55+0.1, 'Ⅳ')plt.show()节点度图

可达矩阵中每个因素所在行与所在列的和称为该因素的节点度,将各阈值对应的节点度由大到小排列,可得到不同阈值下疫情下大学生焦虑情绪形成因子的节点度衰减散点图。 为了得到可达矩阵,需要引入阈值剔除影响程度较小的关系,λ\lambdaλ通常通过经验取值比较。这里分别取0.18、0.20、0.24和0.27 相应代码如下:

import matplotlib as mplimport matplotlib.pyplot as plt#设置图片格式config = { "font.family": 'serif', "font.size": 14, "mathtext.fontset": 'stix', "font.serif": ['SimSun'], 'axes.unicode_minus': False }mpl.rcParams.update(config)#x轴为递增的整数,y轴分别为不同lambda值时候的节点度,节点度在上文代码中可以算出x = list(range(1, 11))y_1 = [13.0, 12.0, 12.0, 11.0, 10.0, 10.0, 9.0, 7.0, 6.0, 2.0] # lbd = 0.15y_2 = [12.0, 11.0, 11.0, 8.0, 7.0, 7.0, 6.0, 5.0, 5.0, 2.0] # lbd = 0.18y_3 = [11.0, 8.0, 8.0, 7.0, 5.0, 4.0, 4.0, 4.0, 3.0, 2.0] # lbd = 0.20y_4 = [8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 3.0, 2.0, 2.0, 2.0] # lbd = 0.24#画图plt.plot(x, y_1, marker='o', linestyle='-', label='λ=0.15', linewidth=1, color='black')plt.plot(x, y_2, marker='x', linestyle='-', label='λ=0.18', linewidth=1, color='black')plt.plot(x, y_3, marker='*', linestyle='-', label='λ=0.20', linewidth=1, color='black')plt.plot(x, y_4, marker='D', linestyle='-', label='λ=0.24', linewidth=1, color='black')plt.legend()plt.xlabel("序号")plt.ylabel("节点度")plt.show()ISM图

通过最后的步骤可以做出系统要素间多级递阶有向拓扑图,这里就不展示了,根据步骤即刻得到之间的关系并作图。

后记

DEMATEL-ISM模型,更多只是矩阵间的相互计算,想要得出结果还是比较容易的。 它可以用来分析因子间的相互关系,做出阶梯式的关系模型图,但是也有主观性强等的缺点。 因为笔者水平有限,所以分析和代码可能不尽人意,多多包含~


李广利,严一知,刘文琦,陈耀光,吴泽玉.基于DEMATEL-ISM的矿工不安全情绪形成因子研究[J]…中国安全科学学报,2021,31(07):30-37.DOI:10.16265/j.cnki.issn1003-3033.2021.07.005. ↩︎

A python library for MCDA methods. ↩︎

本文链接地址:https://www.jiuchutong.com/zhishi/299645.html 转载请保留说明!

上一篇:HTML+CSS简单漫画网页设计成品 蜡笔小新3页 大学生个人HTML网页制作作品(html动画教程)

下一篇:【中秋快乐】如何用three.js实现我的太空遐想3D网页(中秋快乐的祝福语怎么说)

  • 浅谈影响网站权重因素及怎样增加网站权重(浅谈影响网站权重的因素)

    浅谈影响网站权重因素及怎样增加网站权重(浅谈影响网站权重的因素)

  • 红米10xpro搭载的处理器是什么(红米10xpro缺点)

    红米10xpro搭载的处理器是什么(红米10xpro缺点)

  • 为什么群聊突然消失了(为什么群聊突然发不了信息)

    为什么群聊突然消失了(为什么群聊突然发不了信息)

  • 算法一般都可以用哪几种控制结构组合而成(算法一般都可以由顺序循环嵌套的控制结构组合而成)

    算法一般都可以用哪几种控制结构组合而成(算法一般都可以由顺序循环嵌套的控制结构组合而成)

  • qq里的视频保存在手机是存在哪了(qq里的视频保存不到手机)

    qq里的视频保存在手机是存在哪了(qq里的视频保存不到手机)

  • 红米k30充电器多少瓦(红米K30充电器多少w原装)

    红米k30充电器多少瓦(红米K30充电器多少w原装)

  • 小红书直播能看回放吗(小红书直播能看见在线吗)

    小红书直播能看回放吗(小红书直播能看见在线吗)

  • 华为手机自动拨号怎么解决(华为手机自动拨打紧急电话怎么办)

    华为手机自动拨号怎么解决(华为手机自动拨打紧急电话怎么办)

  • 华为售后换机条件(华为售后换机保修怎么算)

    华为售后换机条件(华为售后换机保修怎么算)

  • 腾讯视频能不能看春晚直播(腾讯视频能不能投屏电视)

    腾讯视频能不能看春晚直播(腾讯视频能不能投屏电视)

  • vivo手机挂电话挂不了(vivo手机挂电话怎么有声音)

    vivo手机挂电话挂不了(vivo手机挂电话怎么有声音)

  • word文档生成目录步骤(word文档生成目录时章节内容串位)

    word文档生成目录步骤(word文档生成目录时章节内容串位)

  • realme x2怎么新建联系人(realmex添加桌面小部件)

    realme x2怎么新建联系人(realmex添加桌面小部件)

  • 华为路况目的地怎么更改(华为手机的路况在哪里设置)

    华为路况目的地怎么更改(华为手机的路况在哪里设置)

  • 淘宝直播号怎么换头像(淘宝直播号怎么样可以养起来)

    淘宝直播号怎么换头像(淘宝直播号怎么样可以养起来)

  • 苹果有单独显示器吗(苹果单独显示器在哪买)

    苹果有单独显示器吗(苹果单独显示器在哪买)

  • 消防控制室的布置要求(消防控制室的布置有哪些要求?)

    消防控制室的布置要求(消防控制室的布置有哪些要求?)

  • 苹果x录音怎么没声音(苹果x录音怎么剪辑)

    苹果x录音怎么没声音(苹果x录音怎么剪辑)

  • 华为mate9怎么设置快捷键(华为Mate9怎么设置至第二空间)

    华为mate9怎么设置快捷键(华为Mate9怎么设置至第二空间)

  • 摄像头离线怎么恢复(摄像头离线怎么重新连接wifi)

    摄像头离线怎么恢复(摄像头离线怎么重新连接wifi)

  • jovi是什么功能(vivo手机的jovi是什么意思)

    jovi是什么功能(vivo手机的jovi是什么意思)

  • 电脑屏幕突然黑屏但主机仍在运作的原因?(电脑屏幕突然黑屏怎么回事)

    电脑屏幕突然黑屏但主机仍在运作的原因?(电脑屏幕突然黑屏怎么回事)

  • js中var如何声明作用域(js中var的用法)

    js中var如何声明作用域(js中var的用法)

  • 外籍人士如何在中国办理延长签证
  • 原材料的采购成本包括运费吗
  • 银行电子承兑到期多久时间之内可以兑现
  • 用友财务软件怎么设置二级科目
  • 工厂宿舍的水电安装规范图
  • 可以先注销银行账户吗
  • 专票当成普票入账怎么调整
  • 售后回租产生的费用
  • 机关和事业单位哪个好
  • 广告位租赁交印花税吗
  • 小规模纳税人是简易计税吗
  • 增值税进项和销项税的计算方法
  • 开增值税专票交多少税费
  • 专票开错了可以作废吗
  • 上传失败显示网络不佳
  • 旅游业务开什么发票
  • 代建项目增值税税率
  • 汇兑损失的原因
  • 企业的应收帐款增长率超过销售收入增长率是正常现象
  • 暂估应付款借方
  • 水利基金免税额度是多少
  • 购买的样品怎么做账
  • 小规模纳税人每季度超过30万交税
  • 员工离职再入职要重新签订合同吗
  • 新注册的公司用不用开公户
  • 电子承兑到期提示付款选择线上还是线下
  • 普通增值税发票左上角的二维码可以扫出来吗?
  • 外贸出口换汇成本怎么计算
  • 收到棚户区改造资金
  • 无偿受赠房产出售税费
  • 给予减免
  • php web3
  • php rediscluster
  • 作废的发票对方钱能打出来吗
  • 长期股权投资投资收益计算
  • rgbt目标跟踪
  • 客户多付款入哪个科目
  • 姆科马齐国家公园
  • php 收邮件
  • 低调低调
  • 魔改apk
  • sass和less的区别使用
  • php脚本加密
  • python子类init
  • phpcms怎么样
  • 企业向个人借款合同模板
  • 预收账款可以用吗
  • 可供出售金融资产名词解释
  • 法人和税务可以是一个人吗
  • 企业所得税特殊事项
  • 有限合伙企业收入怎么做账
  • 溢价发行债券取得的收益应如何处理
  • 存货报废的会计分录怎么做
  • 为什么要把工资说的好高呢
  • 企业研制新产品
  • 厂家给的促销费可以退吗
  • 销售费用包括哪些内容?其明细科目有哪些?
  • 生产费用总分类核算的程序包括
  • sql多表关联修改操作
  • linux系统关机重启命令
  • Fedora25踏上Wayland的路上讲解
  • winXP运行内存查看
  • windows8.1备份
  • linux安装软件可以用哪些方式
  • centos获取权限
  • 如何给win7系统升级
  • 近期win7蓝屏
  • Tree、Unformat、Vsafe命令的区别与使用说明
  • unity3D游戏开发
  • [置顶]马粥街残酷史
  • 批处理应用实例
  • 批处理计数
  • 检查shell脚本有无语法错误
  • node.js ts
  • jq获取form表单数据
  • js动态改变网页标题
  • javascript学习指南
  • 国税电子版
  • 个体工商户怎样交税
  • 江西各市财政收入和支出
  • 免责声明:网站部分图片文字素材来源于网络,如有侵权,请及时告知,我们会第一时间删除,谢谢! 邮箱:opceo@qq.com

    鄂ICP备2023003026号

    网站地图: 企业信息 工商信息 财税知识 网络常识 编程技术

    友情链接: 武汉网站建设