首页> 中国专利> 基于种子集的半监督权重核模糊聚类的图像分割方法

基于种子集的半监督权重核模糊聚类的图像分割方法

摘要

本发明公开了一种基于种子集的半监督权重核模糊聚类的图像分割方法,其实现步骤为:(1)选择图像;(2)提取待分割图像纹理特征;(3)产生聚类对象数据矩阵;(4)初始化聚类中心、隶属度和核参数;(5)计算点密度函数值;(6)更新聚类中心、隶属度和核参数;(7)计算目标函数值;(8)判断是否满足终止条件;如果是,执行(9),否则,执行(6);(9)产生分割图像。本发明提取图像每个像素的纹理特征,用基于种子集的半监督权重核模糊聚类的图像分割方法对该纹理特征进行划分,提高了图像分割的稳定性,获得更加准确的分割结果。

著录项

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2016-07-06

    授权

    授权

  • 2014-01-15

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

    实质审查的生效

  • 2013-12-18

    公开

    公开

说明书

技术领域

本发明属于图像处理技术领域,更进一步涉及图像分割技术领域的一种基于 种子集的半监督权重核模糊聚类的图像分割方法。本发明可用于对纹理图像、自 然图像和SAR图像进行分割,以达到目标识别的目的。

背景技术

近年来,图像分割一直是图像处理领域的一个热门研究方向。现有的图像分 割方法根据先验知识的不同主要分为有监督分割方法,半监督分割方法和无监督 分割方法,半监督图像分割方法是近年来提出的一种新方法,该方法只需要少量 的先验知识,因而在适应度上优于有监督分割方法,在分类精度上优于无监督分 割方法。

半监督聚类主要包括基于约束对的方法和基于种子集的方法。从分割结果的 角度看,图像分割的过程就是给每个像素赋予一个类标号,该类标号反映像素在 分割结果中所属的类别,因此,只要找到这些特征的标号,就能实现对像素的分 类,从而得到图像分割的结果。而传统的图像分割技术对纹理图像中的噪声敏感, 易造成过分割现象。纹理作为图像的一个重要特征,在计算机视觉和图像处理中 有重要应用,比如早期的癌细胞识别和遥感图像中军事和民用目标的识别。

华中科技大学在其申请的专利“一种遗传模糊聚类的图像分割方法”(专利 申请号200910273517.7,公开号CN102622761A)中公开了一种利用遗传模糊聚 类分配像素标号的图像分割方法。该方法在聚类过程中,加入聚类中心间距惩罚 措施,虽然能有效分割噪声干扰严重且待分割目标较小的图像,获得正确的聚类 中心,但该方法存在的不足是,分割结果依赖于像素的空间分布,如果该空间分 布的边界是线性不可分的,以及类分布为非高斯分布或类分布为非椭圆分布的样 本不能更好的聚类,导致该分割方法鲁棒性不强,降低了多次分割运行的平均准 确度。

Huaxiang Zhang和Jing Lu在文章“Semi-supervised fuzzy clustering:A  kernel-based approach”(Knowledge-Based Systems,2009,pp.477-481)中提出 半监督核模糊聚类方法,该方法使用种子集来初始化聚类中心,并使用硬化分来 初始化标记样本的隶属度矩阵,同时随机初始化为标记样本的隶属度矩阵,然后 使用折中的方法来计算目标函数值,该方法能避免无监督聚类方法随机初始化聚 类中心陷入局部局部最优的缺点,但是该方法存在的不足是,它对数据集有等划 分趋势,所以对于团状、每类样本数相差较大的数据集都不是最佳聚类方法。

发明内容

本发明的目的在于克服上述已有技术的不足,提出一种基于种子集的半监督 权重核模糊聚类的图像分割方法。本发明提取图像每个像素的小波纹理特征,用 基于核模糊聚类方法和点密度权重的思想对该小波纹理特征向量矩阵进行聚类, 进而对像素进行类划分,达到图像分割的目的。

实现本发明目的的思路是:首先,从所选待分割图像中提取小波特征信息进 行前期处理以产生聚类对象数据矩阵;然后,在聚类过程中结合核模糊聚类方法 和点密度权重的思想寻找最佳的目标函数值;最后,对分割后的每一个类标号, 从灰度值范围[0,255]中任意选择一个整数作为该类标号对应聚类对象数据的灰 度值,实现对图像的分割。

为实现上述目的,本发明具体实现步骤包括如下:

(1)选择图像:

1a)从纹理图像库中下载多幅纹理图像,从所下载多幅纹理图像中任选一 幅图像作为待分割图像;

1b)从纹理图像库中下载与待分割图像对应的参考图像。

(2)提取待分割图像纹理特征:

2a)在待分割图像中,以待提取特征的像素点为中心,选取一个大小为16 ×16的窗口,将该窗口作为子图像块;

2b)利用小波分解公式,得到表示子图像块的纹理特征的10维小波特征向 量矩阵;

2c)将10维小波特征向量矩阵合并,得到待分割图像的10维小波特征向 量矩阵。

(3)产生聚类对象数据矩阵:

采用线性公式,将待分割图像的10维小波特征向量矩阵映射到闭区间[-1,1] 内,得到聚类对象数据矩阵。

(4)获得各变量的初始值:

4a)从聚类对象数据矩阵中,随机选择10%的聚类对象数据作为标记聚类对 象数据,从标记聚类对象数据中,找出类标号相同的聚类对象数据,分别求每个 类标号相同的聚类对象数据的平均值,将所求的平均值作为初始聚类中心值;

4b)采用隶属度方法,初始化隶属度矩阵;

4c)按照下式,获得点密度函数值:

Wi=Σj=1,jiN1||xi-xj||/Σj=1NΣj=1,jiN1||xi-xj||

其中,Wi表示聚类对象数据矩阵中第i个聚类对象数据的点密度函数值,∑ 表示求和操作,||·||表示求欧氏距离操作,xi表示聚类对象数据矩阵中第i个聚类 对象数据,xj表示聚类对象数据矩阵中第j个聚类对象数据,i=1,...,N, j=1,...,N,N表示聚类对象数据矩阵中聚类对象数据的个数;

4d)按照下式,获得初始核参数:

σ0=(1N-1Σi=1N(ai-ai/N)2)1/2

其中,σ0表示控制函数径向作用范围宽度的初始核参数,N表示聚类对象 数据矩阵中聚类对象数据的个数,∑表示求和操作,ai表示聚类对象数据xi到 聚类对象数据矩阵中所有聚类对象数据均值的欧氏距离,xi表示聚类对象数据矩阵中第i个聚类对象数据,||·||表示求欧氏距离操作;

4e)按照下式,计算核函数矩阵:

k(xi,xj)=exp(-||xi-xj||2/2σ02)

其中,k(xi,xj)表示核函数矩阵,xi表示聚类数据对象矩阵中的第i个聚类 对象数据,xj表示聚类数据对象矩阵中第j个聚类对象数据,exp表示取指数操 作,||·||表示求欧式距离操作,σ0表示控制函数径向作用范围宽度的初始核参数。

(5)划分聚类对象数据矩阵:

5a)利用核函数矩阵更新隶属度矩阵,从更新后的隶属度矩阵中找出每个聚 类对象数椐对应的最大隶属度值,并给每个聚类对象数据标上最大隶属度所在的 类标号;

5b)采用聚类中心优化公式,更新聚类中心值;

5c)采用迭代更新公式,更新核参数;

5d)按照下式,获得目标函数值:

Js=wlΣm=1lΣk=1KWlJl+wuΣt=1uΣk=1KWuJu

其中,Js表示第s次的目标函数值,s表示迭代次数,s=1,..,100;wl表示 聚类对象数据矩阵中标记聚类对象数据的折中系数,wl=u/(u+l);u表示聚类 对象数据矩阵中未标记聚类对象数据的个数,l表示聚类对象数据矩阵中标记聚 类对象数据的个数;∑表示求和操作;m表示聚类对象数据矩阵中第m个未标 记聚类对象数据,m=1,...,l;k表示聚类对象数据矩阵中的第k类,k=1,…,K,K 表示聚类对象数据的类别数;Wl表示聚类对象数据矩阵中l个标记聚类对象数据 点密度函数值;Jl表示聚类对象数据矩阵中标记聚类对象数据的目标函数值;wu表示聚类对象数据矩阵中未标记聚类对象数据的折中系数,wu=l/(u+l);t表 示聚类对象数据矩阵中第t个未标记聚类对象数据,t=1,...,u;Wu表示聚类对象 数据矩阵中u个未标记聚类对象数据的点密度函数值;Ju表示聚类对象数据矩阵 中未标记聚类对象数据的目标函数值;

5e)按照下式,获得目标函数值的差:

J=|Js-Js-1|

其中,J表示目标函数值的差,|·|表示取绝对值操作,Js表示第s次的目标 函数值,s表示迭代次数,s=1,..,100,Js-1表示第s-1次的目标函数值;

5f)判断目标函数值的差是否小于阈值10-5,如果是,则执行步骤5g),否 则,执行步骤5a);

5g)输出聚类对象数据的聚类的类标号。

(6)产生分割图像:

对聚类的每一个类标号,从灰度值范围[0,255]中任意选择一个整数作为该类 标号对应聚类对象数据的灰度值,产生分割图像。

本发明与现有技术相比具有以下优点:

第一,本发明在初始化步骤中,利用样本间的方差来初始化核参数,避免了 现有技术手动调整的不合适核参数对方法聚类性能的影响,从而使得本发明方法 的准确度大大提高,因此可以获得更加准确的图像分割结果。

第二,本发明在图像分割的聚类过程中,采用了点密度权重,使得那些团状、 每类样本数相差较大的数据集能得到正确的划分,从而使得本发明可以获得更加 准确的图像分割结果。

第三,本发明在图像分割的前期处理中,利用小波分解提取每个样本数据的 10维特征,为后面的聚类过程提供了更多的图像细节信息,解决了现有技术在 分割复杂图像的情况下由于过平滑而丢失过多局部信息的缺点,从而使得本发明 提高了识别图像中不显著目标的能力。

附图说明

图1是本发明的流程图;

图2是本发明与现有技术在两类纹理图像上的分割结果对比图;

图3是本发明与现有技术在五类纹理图像上的分割结果对比图;

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1,本发明的具体步骤如下:

步骤1,选择图像。

从纹理图像库中下载多幅纹理图像,从所下载多幅纹理图像中任选一幅图像 作为待分割图像。

从纹理图像库中下载与待分割图像对应的参考图像。

步骤2,提取待分割图像纹理特征。

在待分割图像中,以待提取特征的像素点为中心,选取一个大小为16×16 的窗口,将该窗口作为子图像块。

子图像块的纹理特征的10维小波特征向量矩阵的获得公式如下:

e=116×16Σi=116Σj=116|coef(a,b)|

其中,e表示10维小波特征向量矩阵,∑表示求和操作,|·|表示取绝对值 操作,表示子图像块第a行第b列的系数值,a=1,...,16,b=1,...,16。

将10维小波特征向量矩阵合并,得到待分割图像的10维小波特征向量矩阵。

步骤3,产生聚类对象数据矩阵。

将待分割图像的10维小波特征向量矩阵映射到闭区间[-1,1]内,得到聚类对 象数据矩阵,映射公式如下:

b0=(aij-minj)/(maxj-minj)

其中,b0表示聚类对象数据矩阵,aij表示10维小波特征向量矩阵中第i行 第j列的元素值,minj、maxj分别表示10维小波特征向量矩阵中第j列的最小 值和最大值,i=1,…,N,j=1,…,10,N表示聚类对象数据矩阵中聚类对象数据的 个数。

步骤4,获得各变量的初始值。

从聚类对象数据矩阵中,随机选择10%的聚类对象数据作为标记聚类对象数 据,从标记聚类对象数据中,找出类标号相同的聚类对象数据,分别求每个类标 号相同的聚类对象数据的平均值,将所求的平均值作为初始聚类中心值。

隶属度矩阵的初始化步骤如下:

第一步,按照硬化分的方法,将标记聚类对象数据的隶属度矩阵初始化。

第二步,采用下式,对未标记聚类对象数据的隶属度进行初始化:

Σk=1Kμkt=1

其中,∑表示求和操作,μkt表示第t个未标记聚类对象数据属于第k类的 隶属度值,μkt∈[0,1],k=1,…,K,K表示聚类对象数据的类别数,t=1,…,u,u 表示聚类对象数据矩阵中未标记聚类对象数据的个数。

点密度函数值的获得公式如下:

Wi=Σj=1,jiN1||xi-xj||/Σj=1NΣj=1,jiN1||xi-xj||

其中,Wi表示聚类对象数据矩阵中第i个聚类对象数据的点密度函数值,∑ 表示求和操作,||·||表示求欧氏距离操作,xi表示聚类对象数据矩阵中第i个聚类 对象数据,xj表示聚类对象数据矩阵中第j个聚类对象数据,i=1,…,N, j=1,…,N,N表示聚类对象数据矩阵中聚类对象数据的个数。

初始核参数的获得公式如下:

σ0=(1N-1Σi=1N(ai-ai/N)2)1/2

其中,σ0表示控制函数径向作用范围宽度的初始核参数,N表示聚类对象 数据矩阵中聚类对象数据的个数,∑表示求和操作,ai表示聚类对象数据xi到 聚类对象数据矩阵中所有聚类对象数据均值的欧氏距离,xi表示聚类对象数据矩阵中第i个聚类对象数据,||·||表示求欧氏距离操作。

核函数矩阵的计算公式如下:

k(xi,xj)=exp(-||xi-xj||2/2σ02)

其中,k(xi,xj)表示核函数矩阵,xi表示聚类数据对象矩阵中的第i个聚类 对象数据,xj表示聚类数据对象矩阵中第j个聚类对象数据,exp表示取指数操 作,||·||表示求欧式距离操作,σ0表示控制函数径向作用范围宽度的初始核参数。

步骤5,划分聚类对象数据矩阵。

5a)利用核函数矩阵更新隶属度矩阵,从更新后的隶属度矩阵中找出每个聚 类对象数据对应的最大隶属度值,并给每个聚类对象数据标上最大隶属度所在的 类标号。

隶属度矩阵的更新步骤如下:

第一步,按照下式,获得更新后的标记聚类对象数据的隶属度:

μmn(l)=(Σk=1K(k(xm,xm)-2k(xm,vn)+k(vn,vn)k(xm,xm)-2k(xm,vk)+k(vk,vk))1/m-1)-1

其中,μmn(l)表示更新后的标记聚类对象数据的隶属度,∑表示求和操作, k(xm,xm)表示xm和自身的核函数的值,xm表示聚类对象数据矩阵中第m个标记 聚类对象数据,m=1,…,l,l表示聚类对象数据矩阵中标记聚类对象数据的个数, k(xm,vn)表示xm和vn的核函数的值,vn表示聚类对象数据矩阵的第n个聚类中 心,n=1,…,K,K表示聚类对象数据矩阵中聚类对象数据的类别数,k(vn,vn)表 示vn和自身的核函数的值,k(xm,vk)表示xm和vk的核函数的值,vk表示聚类对 象数据矩阵的第k个聚类中心,k=1,…,K,k(vk,vk)表示vk和自身的核函数的 值。

第二步,按照下式,获得更新后的未标记聚类对象数据的隶属度:

μtn(u)=μtnk(xt,vn)Σk=1Kμtnk(xt,vk)

其中,μtn(u)表示更新后的第t个未标记聚类对象数据到第n个聚类中心的 隶属度,μtn表示更新前的第t个未聚类对象数据到第n个聚类中心的距离, k(xt,vn)表示xt和vn的核函数的值,xt表示第t个未标记聚类对象数据,t=1,…,u, u表示聚类对象数据矩阵中未标记聚类对象数据的个数,vn表示聚类对象数据矩 阵的第n个聚类中心,n=1,…,K,K表示聚类对象数据的类别数,∑表示求和 操作,k(xt,vk)表示xt和vk的核函数的值,vk表示聚类对象数据矩阵的第k个聚 类中心,k=1,…,K。

5b)聚类中心的更新公式如下:

vk=Σi=1NWiμkim0φ(xi)Σi=1NWiμkim0

其中,vk表示优化后的第k个聚类中心,k=1,…,K,K表示聚类对象数据 的类别数,Σ表示求和操作,Wi表示聚类对象数据矩阵中第i个聚类对象数据 的点密度函数值,i=1,...,N,N表示聚类对象数据矩阵中聚类对象数据的个数, 表示聚类对象数据矩阵中第i个聚类对象数据隶属于第k类的隶属度值,m0表示权重指数,m0=2,φ(xi)表示聚类对象数据矩阵中第i个聚类对象数据对 应的核空间中的聚类对象数据,xi表示聚类对象数据矩阵中第i个聚类对象数据。

5c)核参数的更新公式如下:

σ=σ0-ησJs-1σ0

其中,σ表示更新后控制函数径向作用范围宽度的核参数,σ0表示更新前 控制函数径向作用范围宽度的核参数,ησ表示控制函数径向作用范围宽度参数 的学习率,Js-1表示第s-1次的目标函数值,s表示迭代次数,s=1,..,100。

5d)目标函数值的获得步骤如下:

第一步,按照下式,获得标记聚类对象数据的目标函数值:

Jl=(μkmm0-μkm,om0)(||φ(xm)-φ(vk)||2-||φ(xm)-φ(vk,o)||2)

其中,Jl表示聚类对象数据矩阵中标记聚类对象数据的目标函数值,表 示更新后的标记对象数据m属于第k类的隶属度,k表示聚类对象数据的第k 类,k=1,…,K,K表示聚类对象数据的类别数,m表示第m个标记聚类对象数 据,m=1,…,l,l表示聚类对象数据矩阵中标记聚类对象数据的个数,m0表示权 重指数,等于常数2,表示初始化的标记对象数据m属于第k类的隶属度, xm表示第m个聚类对象数据,vk表示更新后的第k个聚类中心,vk,o表示初始 化的第k个聚类中心。

第二步,按照下式,获得未标记聚类对象数据的目标函数值:

Ju=μktm0(||φ(xt)-φ(vk)||)2

其中,Ju表示聚类对象数据矩阵中未标记聚类对象数据的目标函数值,表示聚类对象数据矩阵中第t个未标记数据属于第k类的隶属度,k表示聚类对 象数据的第k类,k=1,...,K,K表示聚类对象数据的类别数,t表示第t个未标 记聚类对象数据,t=1,...,u,u表示聚类对象数据矩阵中未标记聚类对象数据的 个数,xt表示聚类对象数据矩阵中第t个未标记聚类对象数据,vk表示更新后的 第k个聚类中心。

第三步,按照下式,获得目标函数值:

Js=wlΣm=1lΣk=1KWlJl+wuΣt=1uΣk=1KWuJu

其中,Js表示第s次的目标函数值,s表示迭代次数,s=1,..,100;wl表示 聚类对象数据矩阵中标记聚类对象数据的折中系数,wl=u/(u+l);u表示聚类 对象数据矩阵中未标记聚类对象数据的个数,l表示聚类对象数据矩阵中标记聚 类对象数据的个数;Σ表示求和操作;m表示聚类对象数据矩阵中第m个未标 记聚类对象数据,m=1,…,l;k表示聚类对象数据矩阵中的第k类,k=1,…,K,K 表示聚类对象数据的类别数;Wt表示聚类对象数据矩阵中l个标记聚类对象数据 点密度函数值;Jl表示聚类对象数据矩阵中标记聚类对象数据的目标函数值;wu表示聚类对象数据矩阵中未标记聚类对象数据的折中系数,wu=l/(u+l);t表 示聚类对象数据矩阵中第t个未标记聚类对象数据,t=1,…,u;Wu表示聚类对象 数据矩阵中u个未标记聚类对象数据的点密度函数值;Ju表示聚类对象数据矩阵 中未标记聚类对象数据的目标函数值。

5e)目标函数值的差的获得公式如下:

J=|Js-Js-1|

其中,J表示目标函数值的差,|·|表示取绝对值操作,Js表示第s次的目标 函数值,s表示迭代次数,s=1,..,100,Js-1表示第s-1次的目标函数值。

5f)判断目标函数值的差是否小于阈值10-5,如果是,则执行步骤5g),否 则,执行步骤5a);

5g)输出聚类对象数据的聚类的类标号。

步骤6,产生分割图像。

对聚类的每一个类标号,从灰度值范围[0,255]中任意选择一个整数作为该类 标号对应聚类对象数据的灰度值,产生分割图像。

本发明的效果可通过以下仿真进一步说明:

1.仿真实验环境与参数设置:

仿真实验环境为:MATLAB7.8.0(R2009a),Hewlett-Packard2.80GHz,32.0GB 内存,Windows XP Professional平台。

仿真实验参数设置为:实验中随机选择10%的标记聚类对象数据,聚类精确 率是10次仿真实验结果的平均值。

2.仿真实验内容:

图2为本发明仿真实验中使用的两类纹理图像和分割结果图,此纹理图像 是从纹理图像库下载的,有两种类标,图像大小为128×128像素。其中,图2 (a)为待分割纹理图像,图2(b)为待分割图像对应的参考图,图2(c)为本 发明的分割结果图,图2(d)为现有技术中的核模糊聚类方法的分割结果图, 图2(e)为现有技术中固定核参数的半监督核模糊聚类方法的分割结果图,图2 (f)为现有技术中调整核参数的半监督核模糊聚类方法的分割结果图。

图3为本发明仿真实验中使用的三类纹理图像和分割结果图,此纹理图像 是从纹理图像库下载的,有五种类标,图像大小为128×128像素。其中,图3 (a)为待分割纹理图像,图3(b)为待分割图像对应的参考图,图3(c)为本 发明的分割结果图,图3(d)为现有技术中的核模糊聚类方法的分割结果图, 图3(e)为现有技术中固定核参数的半监督核模糊聚类方法的分割结果图,图3 (f)为现有技术中调整核参数的半监督核模糊聚类方法的分割结果图。

3.仿真实验结果分析:

通过上述两幅纹理图像的仿真实验以及利用本发明方法和现有技术的分割 结果对比图,可以看出本发明方法在不同纹理图像中都可以获得更加准确的分割 结果。

两类纹理图像的仿真结果如图2所示,本发明仿真采用的纹理图像有两种不 同的纹理特征区域。由图2看出,虽然图2(c)中K均值方法和图2(d)中核K均 值方法的分割结果能把平滑区域分割出来,但边界区域的分割效果并不理想,而 且平滑区域也产生了一部分错分割的点,因此会损失一部分的边缘和细节特征, 而本发明除了对平滑区域取得了理想的分割效果外,同样较好的分割了纹理图像 中的边界区域,而且也使边界区域的分割更平滑获得了更准确的分割结果。

五类纹理图像的仿真结果如图3所示,本发明仿真采用的纹理图像有五种不 同的纹理特征区域。由图3看出,虽然图3(c)中K均值方法和图3(d)中核K均 值方法的分割结果能把平滑区域分割出来,但边界区域的分割效果并不理想,因 此会损失一部分的边缘和细节特征,而本发明除了对平滑区域取得了理想的分割 效果外,同样较好的分割了纹理图像中的边界区域,获得了更准确的分割结果。

用不同算法在相同图像上的聚类精确率作为分割结果的定量评价指标,如果 算法的聚类精确率越高,那么表示算法的分割能力越强。上述三种现有技术和本 发明方法在不同纹理图像上的聚类精确率被列在表1中,表中同时给出了不同算 法在不同图像上的运行时间对比。

表1现有与本发明方法在不同纹理图像上的聚类结果对比

从表1中可以看出,对于不同形状、不同类别的纹理图像,现有的模糊核聚 类方法、固定核参数的半监督模糊核聚类方法和调整核参数的半监督模糊核聚类 方法都有一定的分割效果,它们的聚类精确率都可以达到90%以上,但是相比本 发明方法的聚类精确率,它们的分割效果还是相对差一些,这充分说明了本发明 方法在对图像细节和边缘信息的分割上可以取得较好的分割结果,也正体现了点 密度权重和选择合适初始核参数的优越性。

综合分析实验,我们可以得出结论:本发明最明显的优点在于使用点密度权 重对聚类中心位置进行调整,使其更接近实际的聚类中心,达到正确分类的目的, 这就使得那些每类样本数相差较大的数据集得到正确的划分;同时使用了方差来 计算初始核参数,使得免于手动调整的初始核参数对最终的聚类性能不会产生很 大的影响,因此可以提高算法的聚类性能。

去获取专利,查看全文>

相似文献

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

客服邮箱:kefu@zhangqiaokeyan.com

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

  • 服务号