基于GVF力的拓扑保持分割模型及其对偶算法

沈梦洁 潘振宽 宋金涛 魏伟波

摘要：针对自排斥Snake模型对于狭窄图像区域作用力不足，传统加性算子分裂方法计算复杂，内存用量也会随着图像大小的增加而迅速增长等问题，提出了在原模型的基础上增加梯度矢量流有向力场，以加快轮廓线在图像狭窄区域的演化速度，并为改进的模型设计快速对偶算法以简化算法设计，提高求解效率。数值实验表明，改进模型及算法在计算效率方面较经典模型及算法有较大提高。

关键词：自排斥Snake模型;拓扑保持分割;对偶算法;梯度矢量流;变分法

中图分类号：TP391 文献标志码：A

图像分割是把图像分成各具特性的区域并提取感兴趣目标的技术和过程，在医学影像[1]、遥感影像[2]、深度学习[3]等多个领域得到了广泛应用。在某些特定情况下，需要在图像分割的基础上保持分割目标的拓扑结构来实现更好的分割效果，比如在医学图像中，大脑皮层表面的分割结果必须与真实世界大脑皮层结构一致[4]，再比如具有噪声或遮挡的对象的分割，可以使用拓扑约束来防止过度分割或欠分割等。学者在这一领域提出了许多模型及方法，如基于拓扑先验的最小切割/最大流算法[5];具有可变视场拓扑约束的自动磁共振脊髓分割方法[6];能在多标签图像分割中保持拓扑的方法[7];约束轮廓演化的拓扑保持方法[8];非局部的拓扑保持分割模型和基于非局部形状描述符的分割模型[9-10];使用形状的Beltrami表示进行拓扑保持图像分割的方法[11];基于超弹性正则化的拓扑保持三维图像分割模型[12]等。为使测地线活动轮廓模型[13]具有拓扑结构的适应性，模型求解时往往采用变分水平集方法[14]。变分水平集方法由于具有自适应复杂拓扑结构变化、二维和三维图像分割模型表达一致、数值计算方法稳定以及具有集成多模型信息能力等特点，在图像分割领域被广泛应用[15]。变分水平集方法可以方便地追踪物体的拓扑结构改变，但对于有粘连的物体分割效果不理想。为了在轮廓演化过程中保持拓扑结构，通常在变分公式中添加一个约束项，以防止轮廓自相交，即合并或分裂，如隐式水平集框架中的简单点检测方案，该方案通过引入算法，在每次迭代时都监视水平集函数的符号变化并防止水平集函数改变符号不简单的网格点[16];或将拓扑保持问题重铸为水平集的形状优化问题，使用了符号距离函数，避免了演化轮廓的窄带重叠[17];利用基于结能量最小化的方法来实现同样的效果[18];或在引入非局部正则化项时使用了类似的思想，并应用于遥感图像中细长物体的跟踪[19];而后出现了自排斥Snake模型（Self-repelling Snake，SRS）[20]，该模型使用隐式水平集方法表示曲线，并在经典测地活动轮廓线模型中添加了非局部排斥项[13]。SRS模型对于拓扑保持分割效果较好，但也存在一些问题，如对于狭窄图像区域作用力不足，梯度下降方程的推导复杂，需要使用迎风差分离散格式，内存用量也会随着图像大小的增加而迅速增长。文献[21]使用了Split Bregman算法求解原SRS模型，相对于SRS模型，减少了内存用量，缩短了迭代时间。本文针对狭窄图像区域作用力不足，设计一种加入梯度矢量流有向力场[22]的模型，通过引入对偶变量，利用凸优化模型的快速对偶算法[23]，使用投影方法进行约束，并对其迭代求解，与原模型及方法在多幅图像上开展对比实验，通过对比分割效果、迭代次数以及运行时间，验证了所提模型及方法的有效性和高效性。

1 相关研究工作基础

1.1 SRS模型

SRS模型是在测地活动轮廓线模型[10]基础上加入气球力项和排斥力项，使用隐式水平集方法来表示曲面，其中分割轮廓隐式表示为符号距离函数[24]的零水平线，模型的定义为：令fx：Ω→R为标量值图像，x∈Ω，Ω为图像域，Ω为其边界，标准边缘检测函数gx

gx=11+ρGσfs（1）

其中，s=1或2，ρ是缩放参数，Gσ表示图像的高斯卷积，标准差为σ，对象边界C用水平集函数的零水平集表示

C=x∈Ωx=0（2）

水平集函数被定义符号距离函数

x=-dx，C  x inside C0      x∈Cd（x，C）   x outside C （3）

其中，dx，C是点x与轮廓C之间的欧氏距离。作为一个符号距离函数，满足以下约束条件

=1（4）

为了表示图像区域和轮廓，引入了Heaviside函数H和Dirac函数δ，由于原始 Heaviside 函数是不连续的，因此不可微，所以采用以下正则化方案

Hε=121+ε+1πsinπε ≤ε1            >ε0            <-ε （5）

δε=12ε1+cosπε    ≤ε0            >ε （6）

SRS模型的能量泛函為

E（）=γ∫ΩgxHεxdx+α∫Ωgx1-Hεxdx-

β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy（7）

其中，等式右边第一项为等高线的测地线长度[13]，改写为γ∫Ωgxxδεxdx;第二项为气球力项，能够推动分割轮廓越过弱边缘;第三项为排斥力项，在分割轮廓线较近时能够限制拓扑结构的改变。γ，α，β为平衡三个条件的惩罚参数，e-x-y2d2是用来测量两个点x和y的接近度，距离越远的点排斥力越小，hεx和hεy表示x点和y点周围的窄带（l为实数）

hεx=Hεx+l1-Hεx-l（8）

hεy=Hεy+l1-Hεy-l（9）

SRS模型的能量极值问题

minE=γ∫Ωgxxδεxdx+α∫Ωgx1-Hεxdx-      β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy s.t.=1 （10）

应用变分方法解以上三个能量项，演化方程为

t=δεxγ·gx+αgx+4βd2hεx∫Ωe-x-y2d2x-y·yhεydy x∈Ωx，0=0x     t=0=0         x∈Ωs.t.=1 （11）

关于=1的约束，采用以下动态初始化方案

t+signx-1=0x，0=x （12）

式（12）是一个典型的Hamilton-Jacobi方程，可以通过迎风差分格式[24]进行离散求解。

1.2 SRS模型的AOS（Additive Operator Splitting）方法

文献[20]中应用的是传统AOS方法。用半点差分格式和谐波平均近似对式（11）中的第一个式子等式右边的第一项进行离散，后两项采用迎风方案。通过将图像的行和列分别连接起来，构造了两个半隐式格式

1-2τAx1kvk+1=k+τT2k+T3k（13）

1-2τAx2kwk+1=k+τT2k+T3k（14）

其中，Ax1和Ax2为两个串联矩阵，v和w为中间变量，T2和T3是演化方程中的第一个式子等式右边的第二项和第三项的迎风离散化，对于每个All∈x1，x2

Alijk=2γοkiοkigi+οkjgj     j∈Nli-Σm∈Nl（i）2γοkiοkigi+οkmgm j=i0              else （15）

其中，i，j为图像中的两点，Nli是矩阵Al中i的邻域

οkj=i+1，j-i-1，j22+i，j+1-i，j-122（16）

最后，k+1计算为

k+1=12vk+1+wk+1（17）

总体上说，AOS方法解决了SRS模型的计算问题，但是引入了新的参数，且算法计算复杂，内存用量也会随着图像大小的增加而迅速增长。

2 ISRS（Improved Self-repelling Snake）模型及对偶算法

2.1 对偶算法

通过使用对偶算法[23]，可以简化算法设计，提高求解效率。对偶算法的主要解决思路是引入对偶变量

γ∫ΩgxHεxdx=max：≤γg（x）∫ΩHεx·dx（18）

对原变量和对偶变量进行快速迭代计算，得到SRS模型的交替优化形式为（k为迭代次数）

k+1，k+1=argmin：=1max：≤γgx∫ΩHεx·dx+α∫Ωgx1-Hεxdx-

β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy  （19）

交替优化方法通过固定其他变量，从而求解一个变量的极值问题，式（19）可转化为以下子优化问题

k+1=argmax E1=Ek，s.t.≤γgx （20）

k+1=argmin E2=E，k+1（21）

首先，固定变量k，由式（19）得到关于k+1的能量泛函极值问题

k+1=max≤γgxEk，=max≤γgx∫ΩHεkx·xdx （22）

求解式（22），得到的梯度下降方程

t=-Hεkxs.t.≤γgx （23）

得

x=kx+τkpHεkx（24）

求对偶变量的解，得

k+1x=kx+τkpHεkxs.t. k+1≤γgx （25）

由式（25）得的投影公式

k+1xγgxk+1xmaxk+1（x），γgx（26）

固定k+1，求解k+1，由式（19）得到关于的能量泛函极值问题

k+1=min：=1E，k+1=min：=1∫ΩHεx·k+1dx+α∫Ωgx1-Hεxdx-

β∫Ω∫Ωe-x-y2d2x·yhεxhεydxdy （27）

对式（27）求偏导，得

E+εηε=∫Ω·k+1δεxηxdx-α∫Ωgxδεxηxdx-

4βd2∫Ωhεx∫Ωe-x-y2d2x-y·yhεydyηxdx（28）

其中，ε為一个极小增量，η为一个与性质相同的任意函数。

通过梯度降方法来求解，得

t=αgx-·k+1δεkx+   4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy=1 （29）

整理，得

t=αgx-·k+1δεkx+   4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydyt+signx-1=0 （30）

2.2 GVF模型

本文提出的方法添加了GVF（Gradient Vector Flow）模型[22]，该模型可以有效地克服模型不能收敛到轮廓凹陷的问题，且该模型是基于外力矢量的能量函数构造的，利用该能量函数可以产生约束轮廓向边界演化的力。其中GVF力的力场为：x，y=ux，y，vx，y，能量泛函定义为

E=α∫ 2dx+∫f 2-f2dx（31）

其中，f为与图像数据有关的外部能量项，f为外部作用力。α为惩罚参数用来平衡能量项中的第一项与第二项，α越大，该力场就越平滑。在图像边缘附近f2较大，外部能量项（第二项）起主要作用，此时趋近于f。而在灰度比较均匀的区域，式中第二项为零，只有第一项起作用，此时具有平滑的作用。

令Qk=αgx-·k+1，用代替δεx，并将GVF力·加入到式（30）当中，得

t=Qk+·+4βd2hεkx∫Ωe-x-y2d2x-y·kyhεkydy（32）

2.3 离散化及迭代方案

对初始化

k + 1，0i，j = ki，j（33）

令

Qki，j=αgi，j--·k+1i，j（34）

使用迎风差分格式离散求解，得

k+1，li，j+=max-x1k+1，li，j，02+min+x1k+1，li，j，02+max-x2k+1，li，j，02+min+x2k+1，li，j，02（35）

k+1，li，j-=min-x1k+1，li，j，02+max+x1k+1，li，j，02+min-x2k+1，li，j，02+max+x2k+1，li，j，02（36）

其中

+i，j=i+1，j-i，j（37）

-i，j=i，j-i-1，j（38）

式（32）的求解结果为

k+1，l+1i，j-k+1，li，jτ=maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+

minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+

4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy（39）

整理，得

k+1，l+1i，j=k+1，li，j+τ＼[maxQki，j，0k+1，li，j-+minQki，j，0k+1，li，j++maxF1，0+x1k+1，li，j+

minF1，0-x1k+1，li，j+maxF2，0+x2k+1，li，j+minF2，0-x2k+1，li，j+

4βd2hεk+1，li，jx∫Ωe-x-y2d2x-y·k+1，li，jyhεk+1，li，jydy＼]（40）

ISRS模型的对偶算法为：

1）初始化。根据式（1）计算gx，初始化0x为符号距离函数并设置0=0;

2）设置惩罚参数的值;

3）设置时间步长和迭代次数;

4）迭代过程：

for i =0，1，2，…，step

compute k+1 from （26）;

for k=0，1，2，…，K

for s=0，1，2，…，S

calculate k+1，s+1 from （40）

end for s

end for k

end for i

3 数值实验与分析

本部分针对本文提出的ISRS模型的对偶算法进行分割实验，并与经典SRS模型的AOS方法进行比较，对比两者的拓扑保持的分割效果及计算效率。为了客观比较实验结果，对每个实验中的两种方法使用相同的初始化位置。实验中用黑色线条表示分割区域边界。实验环境：CPU为Intel（R）Core（TM）i5-8500T，内存8G，操作系统为Windows10，软件平台为Matlab2018b。

图1为SRS模型的AOS方法分割合成手部图像的效果图，图2为ISRS模型的对偶算法分割合成手部图像的效果图。

图1 SRS模型的AOS方法分割合成手部图像

（a）初始轮廓线;（b）迭代200次效果;（c）迭代500次效果;（d）迭代900次效果

图2 ISRS模型的对偶算法分割合成手部图像

（a）初始轮廓线;（b）迭代20次效果;（c）迭代40次效果;（d）迭代60次效果;（e）迭代100次效果;（f）迭代200次效果

图1参数分别为1、0.5、1、2×2、-0.2、-0.2、1，其中h、τ、l均根据文献[20]中参数设置;window表示搜索框的大小，α，β，ε均根据文献[20]中参数做了适量调整，图像大小127×150像素。圖2参数γ、α、β、ε、τ、τ2、τ3、l、d、window分别为1、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，其中，γ为Snake模型的参数，用来保持正常的分割作用;α为气球力项参数，能使模型在平滑区域加速分割;β为排斥力项参数，用于保持0水平集的拓扑结构;ε为域值，让水平线函数保持0-1的形态，从而使函数变为分段常值函数，便于计算和保持良好的形态;τ为以及k迭代过程中的时间步长;τ2为s迭代过程中的时间步长;τ3为迭代过程中的时间步长;以上三项均根据时间步长基本原则进行取值，为长期实验所取得的经验值;l为拓扑保持距离，l越大，拓扑保持距离越远;d为高斯卷积核的标准差， window为搜索框大小;图像大小300×354像素。

对比图1图2两组图像，可以看出两组运行结果基本一致，但本文改进的模型及方法相比原来的模型及AOS方法迭代次数有了大幅度地减少。本文改进的模型及方法在迭代到第100步时便达到收敛条件，完成分割，并且拓扑结构得到了保持，而传统的模型及方法900步才完成分割。

图3为SRS模型的AOS方法分割人为噪声图像的效果图，图4为ISRS模型的对偶算法分割人为噪声图像的效果图。图3参数分别为1、0.5、1、2×2、-0.2、-0.2、1，图像大小200×138像素。图4参数γ、α、β、ε、τ、τ2、τ3、l、d、window分别为8、10、0.01、0.1、0.5、0.2、0.125、1、4、3×3，图像大小400×275像素。

图3 SRS模型的AOS方法分割人为噪声图像

（a）初始轮廓线;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果

图4 ISRS模型的对偶算法分割人为噪声图像

（a）初始轮廓线;（b）迭代10次效果;（c）迭代20次效果;（d）迭代40次效果;（e）迭代80次效果;（f）迭代200次效果

对比图3图4，两种方法的运行结果都较好地保持了图像本身的拓扑结构，但本文改进的模型及方法在迭代到第80步时便达到收敛條件，完成分割，并且拓扑结构得到了保持，而传统模型及AOS方法在第800步才完成分割。

图5为SRS模型的AOS方法分割脑部切片1的效果图，图6为ISRS模型的对偶算法分割脑部切片1的效果图。图5参数h、τ、l、window、α、β、ε分别为1、0.5、1、2×2、-0.2、-0.2、1，图像大小125×125像素。图6参数γ、α、β、ε、τ、τ2、τ3、l、d、window分别为1、50、0.01、1、0.3、0.12、0.125、1、4、5×5，图像大小300×300像素。

图5 SRS模型的AOS方法分割脑部切片1

（a）初始轮廓线;（b）迭代200次效果;（c）迭代400次效果;（d）迭代700次效果

图6 ISRS模型的对偶算法分割脑部切片1

（a）初始轮廓线;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果

对比图5图6，本文改进的模型及方法同原模型的AOS方法相比能更好地进入到凹槽内，并且在迭代次数上也有了大幅度地减少。本文改进的模型及方法在迭代到第50步时便达到收敛条件，完成分割，并且拓扑结构得到了保持，而传统模型及AOS方法在第700步才完成分割。

图7为SRS模型的AOS方法分割脑部切片2的效果图，图8为ISRS模型的对偶算法分割脑部切片2的效果图。图7参数h、τ、l、window、α、β、ε分别为1、0.5、1、2×2、-0.2、-0.2、1，图像大小125×125像素。图8参数γ、α、β、ε、τ、τ2、τ3、l、d、window分别为1、50、0.01、10、0.3、0.12、0.125、1、4、5×5，图像大小300×300像素。

（a）初始轮廓线;（b）迭代200次效果;（c）迭代400次效果;（d）迭代800次效果

（a）初始轮廓线;（b）迭代10次效果;（c）迭代20次效果;（d）迭代30次效果;（e）迭代50次效果;（f）迭代100次效果

对比图7图8，两种方法都保持了对象本身的拓扑结构，本文改进的模型及方法在迭代到第50步时便达到收敛条件，完成分割，并且拓扑结构得到了保持，而传统模型及AOS方法在第800步才完成分割。

表1给出了以上四个对比实验的迭代次数、运算时间以及运算时间加速比的具体数值。可知，在有参数影响的情况下，本文也能较好地保持分割对象的拓扑结构，并在保持对象拓扑结构的基础上大幅度减少了算法的迭代次数，模型的收敛速度也得到了大幅度提高。

4 结论

为了进一步提高拓扑保持分割模型的计算效率，本文提出了ISRS模型及其对偶算法，简化了计算过程，提高了运算速度;GVF有向力场明显加快了轮廓线在图像狭窄区域的演化速度。对比实验表明，本文模型及其对偶算法的解类似于经典SRS模型及其方法的解，可以有效防止分割轮廓的合并或分裂，保持原来的拓扑结构，并且在计算效率方面较经典模型与算法有较大提高。下一步工作可以尝试用高阶模型来忽略分割过程中的微小细节，保持边缘的纹理特征。

参考文献

[1]国凯，王国栋，潘振宽，等.CT影像中的主动脉血管钙化检测[J].青岛大学学报（自然科学版），2013，26（1）：50-54.

[2]田佑仕，黄宝香，姜涛，等.融合图像恢复的遥感数据变分分割方法研究[J].青岛大学学报（自然科学版），2020，33（4）：1-8.

[3]吴则举，王嘉琦，焦翠娟，等.基于深度学习网络U-Net的轮胎带束层分割算法研究[J].青岛大学学报（自然科学版），2019，32（4）：22-29+35.

[4]MACDONALD D， KABANI N， AVIS D， et al. Automated 3-D extraction of inner and outer surfaces of cerebral cortex from MRI[J]. NeuroImage， 2000， 12 （3）：340-356.

[5]ZENG Y， SAMARAS D， CHEN W， et al. Topology cuts： A novel min-cut/max-flow algorithm for topology preserving segmentation in N-D images[J]. Computer Vision and Image Understanding， 2008， 112 （1）：81-90.

[6]CHEN M， CARASS A， OH J， et al. Automatic magnetic resonance spinal cord segmentation with topology constraints for variable fields of view[J]. NeuroImage， 2013，83： 1051-1062.

[7]WAGGONER J， ZHOU Y J， SIMMONS J， et al. Topology-preserving multi-label image segmentation[C]// IEEE Winter Conference on Applications of Computer Vision （WACV 2015）， Waikoloa， 2015： 1084-1091.

[8]LIU G Q， LI H F， YANG L. A topology preserving method of evolving contours based on sparsity constraint for object segmentation[J]. IEEE Access， 2017， 5：19971-19982.

[9]DEBROUX N， OZERE， S， LE GUYADER C. A non-local topology-preserving segmentation-guided registration model[J]. Journal of Mathematical Imaging and Vision， 2017， 59 （3）：432-455.

[10] DEBROUX N， LE GUYADER C. A joint segmentation/registration model based on a nonlocal characterization of weighted total variation and nonlocal shape descriptors[J]. SIAM Journal on Imaging Sciences， 2018， 11 （2）：957-990.

[11] CHAN H L， YAN S， LUI L M， et al. Topology-preserving image segmentation by beltrami representation of shapes[J]. Journal of Mathematical Imaging and Vision， 2018， 60 （3）：401-421.

[12] ZHANG D P， LUI L M. Topology-Preserving 3D image segmentation based on hyperelastic regularization[J]. Journal of Scientific Computing， 2021， 87（3）： 74.

[13] CASELLES V， KIMMEL R， SAPIRO G. Geodesic active contour. IJCV[J]. International Journal of Computer Vision， 1997， 22（1）：61-79.

[14] ZHAO H， CHAN T， MERRIMAN B， et al. A variational level set approach to multiphase motion[J]. Journal of Computational Physics， 1996， 127（1）： 179-195.

[15] 郭振波.基于變分水平集方法的多相图像分割研究[D]. 青岛： 中国海洋大学， 2008.

[16] HAN X， XU C Y， PRINCE J L. A topology preserving level set method for geometric deformable models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2003， 25（6）：755-768.

[17] ALEXANDROV O， SANTOSA F. A topology-preserving level set method for shape optimization[J]. Journal of Computational Physics， 2005， 204（1）：121-130.

[18] SUNDARAMOORTHI G， YEZZI A. Global regularizing flows with topology preservation for active contours and polygons[J]. IEEE Transactions on Image Processing， 2007， 16（3）：803-812.

[19] ROCHERY M， JERMYN I H， ZERUBIA J. Higher order active contours[J]. International Journal of Computer Vision， 2006， 69（1）：27-42.

[20] GUYADER C L， VESE L A. Self-repelling snakes for topology-preserving segmentation models[J]. IEEE Transactions on Image Processing， 2008， 17（5）：767-779.

[21] PAN H， SONG J， LIU W， et al. Using the split bregman algorithm to solve the self-repelling snake model＼[J＼]. Computer Vision and Pattern Recognition， 2021＼[2021-08-15＼]https：//doi.org/10.1007/s10851-021-01065-9.

[22] XU C， PRINCE J L. Snakes， shapes， and gradient vector flow[J]. IEEE Transactions on Image Processing， 1998：7（3）：359-369.

[23] CHAMBOLLE A. An algorithm for total variation minimization and applications[J]. Journal of Mathematical Imaging and Vision， 2004， 20（1-2）： 89-97.

[24] OSHER S， SETHIAN J A. Fronts propagating with curvature-dependent speed： Algorithms based on Hamilton-Jacobi formulations＼[J＼]. Journal of Computational Physics， 1988， 79 （1）： 12-49.

Topology Preserving Segmentation Model Based on

GVF Force and Its Dual Algorithm

SHEN Meng-jie， PAN Zhen-kuan， SONG Jin-tao， WEI Wei-bo

（College of Computer Science & Technology， Qingdao University， Qingdao 266071， China）

Abstract：

In view of the insufficient force of the self-repelling Snake model on narrow image regions， the calculation of traditional additive operator splitting method is complicated， and the memory usage will also increase rapidly with the increase of the image size. It is proposed that the Gradient Vector Flow directed force field based on the original model was added to accelerate the evolution speed of the contour line in the narrow area of the image， and a fast Dual algorithm is designed for the improved model to simplify the algorithm design and improve the efficiency of the solution. Numerous numerical experiments show that the proposed improved model and algorithm have a greater improvement in computational efficiency than the classic model and algorithm.

Keywords：

self-repelling snake model; topology-preserving segmentation; dual algorithm; gradient vector flow; variational method

收稿日期：2021-09-10

基金項目：

国家自然科学基金（批准号：61772294，11472144） 资助;山东省联合基金（批准号：ZR2019LZH002） 资助。

通信作者：

潘振宽，男，博士，教授，主要研究方向为计算机视觉、图像处理，E-mail： zkpan@126.com