首发于SIGAI
理解变分自动编码器

理解变分自动编码器

其它机器学习、深度学习算法的全面系统讲解可以阅读《机器学习-原理、算法与应用》,清华大学出版社,雷明著,由SIGAI公众号作者倾力打造。

本文PDF全文链接:理解变分自动编码器

今天的文章用深入浅出的语言和形式为大家介绍变分自动编码器(VAE)的基本原理,以帮助初学者入门,真正理解这一较为晦涩的模型。还是那种熟悉的风格和味道!读懂本文需要读者理解KL散度包括正态分布之间的KL散度计算公式、KL散度的非负性(涉及到变分法的基本概念),蒙特卡洛算法等基本知识,自动编码的知识。

本文是《机器学习-原理,算法与应用》(2019年8月底出版,《机器学习与应用》的第二版)以及《机器学习与应用》雷明著(2019年1月出版),清华大学出版社,这两本书的补充,限于篇幅,没有在书中详细介绍VAE。如果想系统的学好机器学习与深度学习,它们是不错的入门教材,对于想在这一领域提高与进阶的读者,也是不错的选择。

数据生成问题

数据生成可要解决与模式识别相反的问题,其目标是根据语义信息生成样本数据。语义信息可以是类别标签,也可以是其他抽象信息如笔画、风格等。数据生成模型以这些语义信息作为输入,输出是符合概率条件并具有随机性的样本数据。符合概率分布条件可以形象化的理解为生成的样本数据要“像”真实的样本数据。


数据生成模型以生成图像、声音、文字等数据为目标,生成的数据服从某种未知的概率分布。以图像生成为例,假设要生成狗,汉堡,风景等图像。算法输出向量X,该向量由图像的所有像素拼接而成。各类样本X都服从各自的概率分布,定义了对所有像素值的约束。例如对于狗来说,如果图像看上去要像真实的狗,则每个位置处的像素值必须满足某些约束条件,且各个位置处的像素值之间存在相关性。算法生成的样本要有较高的概率值,即像真的样本。图1为典型的生成模型-生成对抗网络所生成的逼真图像。


图1生成对抗网络生成的逼真图像


用概率分布变换生成数据


为了保证生成的样本具有随机性,生成算法通常都要借助于随机数。分布变换是生成随机数的一种常用手段,它将服从简单分布的随机数映射成服从复杂分布的随机数。计算机里可以直接生成的是均匀分布的随机数,借助它即对它进行变换,可以生成各种类型分布的随机数。


简单的分布变换可以人工设计。以生成正态分布的随机数为例,广为使用的Box-Muller算法将均匀分布的随机数映射成正态分布的随机数。假设随机变量 u_{1}u_{2} 服从 [0,1] 内的均匀分布,则随机数 z_{1}z_{2}


\begin{aligned} {z_{1}=\sqrt{-2 \ln x_{1}} \cos 2 \pi x_{2}} \\ {z_{2}=\sqrt{-2 \ln x_{1}} \sin 2 \pi x_{2}} \end{aligned}



相互独立并且服从正态分布 N(0,1) 。借助于均匀分布随机数,通过上面的变换就可以得到正态分布的随机数。


一个重要结论是:以服从正态分布的随机数作为输入,对它们进行映射,可以生成任意分布的随机数。下面用一个简单的例子进行说明。假设要构造位于圆环上的2D随机向量,如果\mathbf{Z} 是2D正态分布的随机数,则


\begin{aligned} g(\mathbf{z})=\frac{\mathbf{z}}{20}+\frac{\mathbf{z}}{\|\mathbf{z}\|} \end{aligned}


即为目标概率分布的随机数,如图2所示。


图2将正态分布随机数映射为圆环上的分布


复杂数据的生成同样可通过分布变换实现。假设输入随机向量 \mathbf{z} 服从概率分布 p(\mathbf{z}) ,此分布的类型一般已知,称为隐变量,典型的是正态分布和均匀分布。通过函数 \mathbf{g} 对此随机向量进行映射,得到输出向量


\begin{aligned} \mathbf{x}=g(\mathbf{z}) \end{aligned}


向量 \mathbf{x} 服从真实样本的概率分布p _{r}(\mathbf{x}) 。如果已知要生成的概率分布 p_{r}(\mathbf{x}) ,借助于概率论中的逆变换算法,可以人工显式地构造出分布变换来生成服从概率分布的随机数。但一般的数据生成问题用这种方法存在以下问题:


1.对于图像、声音生成等问题,样本所服从的概率分布 p_{r}(\mathbf{x}) 是未知的。我们只有一批服从这种未知概率分布的样本,需要从这组样本进行学习,生成与它们相似但又不完全相同的样本。此时无法直接得到此概率分布的表达式,因此无法手工设计出针对它的分布变换。


2.即使样本的概率分布 p_{r}(\mathbf{x}) 已知,当 \mathbf{x} 维数很高且所服从的概率分布非常复杂时人工设计概率分布变换的做法将面临困难。


学习出一个概率分布变换函数


可以用机器学习的方法解决数据生成问题,分布变换函数 \mathbf{g} 通过训练确定。这体现了机器学习中的一贯套路:一切皆可学习,从分类函数,回归函数,到距离度量学习,强化学习的策略函数,只要给我样本,我就可以解决问题!

给定一组样本 \mathbf{x}_{i}, i=1, \dots, l ,它们采样自某种概率分布 p_{r}(\mathbf{x}) 。训练目标是模型生成的样本所服从的分布与真实的样本分布一致,即看上去像又不完全相同。直观来看,生成模型生成的样本与真实样本“相似”。这一过程如图3所示。


图3通过学习得到概率变换函数


这里的映射函数是通过学习得到而非人工设计,在深度生成模型中用神经网络表示。主流的深度生成模型,如变分自动编码器、生成对抗网络均采用了这种思路。问题的关键是:


1.如何判断模型所生成的样本与真实的样本分布一致。

2.如何在训练过程中迫使映射函数生成的样本逐步趋向于真实的样本分布。

本文将要讲述的变分自动编码器使用变分推断和神经网络作为工具来解决此问题。


变分自动编码器


变分自动编码器(Variational Auto-Encoder,简称VAE)由Kingma等人提出[1],是对复杂的概率分布进行无监督学习的典型方法。VAE是变分推断与神经网络相结合的产物。整个系统遵循自动编码器的结构,由编码器和解码器构成。在训练时,编码器将训练样本映射成隐变量所服从的概率分布的参数,然后从此概率分布进行采样得到隐变量,解码器则将隐变量映射回样本变量,即进行重构。这种方法在标准自动编码器的基础上加入了随机性,从而保证可以输出带有随机性的数据。由于用神经网络进行实现,因此可以用随机梯度下降法训练。



使用隐变量



生成模型对概率分布p(x)建模,随机向量x是高维空间的样本,所有样本的集合称为样本空间。对于图像,每个样本点是上千维(由图像所有像素拼接而成)的向量,生成模型的任务是刻画像素之间的依赖关系。如位置相近的像素有相似的颜色,它们被组织成各种物体。如果确定了概率分布函数,给定样本点x,可计算出p(x)的值即对样本进行概率评估。看上去像真实图像的x有更大的概率值,看上去像随机噪声的x的概率值则很小。


但现实应用中更关心的是如何生成与训练样本集中的样本相似但又不完全相同的新样本,而非对已知样本进行概率评估。算法用训练集中的样本数据进行学习,然后生成新的、没有见过的样本,即获得服从某种潜在且未知的概率分布 p_{gt}(x) 的样本x。算法需要学习出一个概率模型p(x),它尽可能与$ p_{gt}(x) 相似,从这个概率分布可以采用出样本。


直接生成样本x或得到函数p(x)是困难的,通常要借助隐变量。假设要生成手写数字的图像,如果算法在给每个像素赋值之前先确定要写0-9之中的哪个数字,则问题变得更为容易。在这里数字的类别即为隐变量。模型在画图之前,先从0-9之中随机地选择一个数字值,然后确保所有的笔画都符合这个数字的要求。这种做法如图4所示。


图 4 通过隐变量生成数字图像


z之所以被称为隐变量,是因为如果只给定一张由模型生成的数字图像,我们并不知道这变量z是取哪个值的时候生成了此图像,也就是说该变量的值无法直接观察到。如果要知道这个变量的值,必须进行推断,可以采用计算机视觉技术,这就是图像分类问题。


对于样本集中的每个样本点x,至少存在一个或多个隐变量取值可以导致模型生成类似于它的样本。对于高维空间中的隐变量向量z,它采样自一个概率分布p(z),有某一确定性的函数族g (z;θ) ,其参数为θ,将隐变量映射到样本空间


z\rightarrow x


即将服从简单分布的隐变量映射为服从复杂分布的样本变量。在这里θ和映射函数是确定的,而z是随机变量,因此g(z;θ)是样本空间中的一个随机变量。训练时的目标是优化参数θ从而确定映射函数。目标是根据概率分布p(z)采样出一个隐变量值z,使得g(z;θ)以很高的概率像数据集中的样本x。这里的映射函数是随机变量的变换函数而非概率密度函数。


由于要使得生产的样本概率最大化,因此训练时采用最大似然估计,最大化样本集中的样本x在生成过程中的概率值


p(\mathbf{x})=\int_{\mathbf{z}} p(\mathbf{x}, \mathbf{z} ; \boldsymbol{\theta}) d \mathbf{z}=\int_{\mathbf{z}} p(\mathbf{x} | \mathbf{z} ; \boldsymbol{\theta}) p(\mathbf{z}) d \mathbf{z}(1)


在这里g(z;θ)被条件概率p(x z;θ)所取代,此概率分布显式的表达了x和z之间的概率依赖关系。式1基于全概率公式,对隐变量求积分。其含义是在所有各种可能的隐变量取值情况下,样本x出现的概率(严格来说是概率密度值)有多大。该概率值也可以看做是条件概率的数学期望



\mathrm{E}_{\mathbf{z} \sim p(\mathbf{z})}[p(\mathbf{x} | \mathbf{z} ; \mathbf{\theta})]



如果生成模型可能产生训练样本集中的样本,则它可能会产生与其类似的样本。在VAE中,上面的条件概率一般用正态分布进行建模


p(\mathbf{x} | \mathbf{z} ; \boldsymbol{\theta})=N(\mathbf{x} | g(\mathbf{z} ; \boldsymbol{\theta}), \mathbf{\Sigma})(2)


即给定隐变量值的条件下,样本值服从正态分布,其均值为映射函数的输出值g(z;θ),协方差矩阵通常为单位矩阵I与方差的乘积 σ^2 ,其中σ是人工指定的超参数。p(z)是已知的先验分布,通常设定为正态分布。


编码器-解码器结构



根据式1,VAE训练时的目标是近似地优化概率值p(x)。计算此概率值需要解决以下几个问题:怎样定义隐变量z,如何计算对z的积分。


第一个问题是如何选择隐变量z以捕获数据中的隐含信息。以生成数字图像为例,模型在绘制数字图像之前要做的隐决策非常复杂。不仅要选择绘制哪个数字,还要决定数字的倾角、笔画宽度、风格特征等。这些特征之间可能还存在相关性,手工设计特征显然不现实。在VAE中并不需要人工设计z的每一维,只是假设z服从某一概率分布,如N(0,1)。根据之前的结论,这种做法是可行的。


只要在训练时学习到一个函数,将各个分量独立、正态分布的随机数z映射成模型想要的任意的隐变量,问题即可解决。根据式2,如果g(z;θ)是一个多层神经网络,则该神经网络前面的层将正态分布的随机数映射为隐变量,后面几层将隐变量映射为样本向量。


接下来需要解决的问题是对于训练样本集如何最大化p(x)。假设隐变量服从标准正态分布


p(\mathbf{z})=N(\mathbf{z} | \mathbf{0}, \mathbf{I})


使用蒙特卡洛算法计算式1所定义的数学期望。对于单个训练样本x,先采样出服从概率分布p(z)的隐变量值{z1,...zn},然后计算数学期望值


p(\mathbf{x}) \approx \frac{1}{n} \sum_{i} p\left(\mathbf{x} | \mathbf{z}_{i}\right)(3)


若z的维数很高,则需要大量的样本才能计算出准确的概率值。事实上,对于大多数的z取值,其p(x 丨z)接近于0,即这些隐变量取值不太可能会生成该样本数据。由于这些z值对p(x)的贡献接近于0,在计算p(x)时可偏向于采样那些产生过x的z值。由此引入一个概率分布q(z丨x) ,其输入为x,输出为那些很可能产生过x的隐变量z值的概率分布,即近似后验概率p(z丨x)。q(z丨x)所定义的隐变量空间要比隐变量的原始空间p(z)小。因此可以用Ez~q p(x丨z)近似p(x),而计算这个数学期望则相对简单。不直接使用p(z丨x)是因为该分布的形式未知且通常很复杂,不便于计算。我们可以将q(z丨x)设定为简单的概率分布,从而简化问题的处理难度。


VAE在训练时需要找到一个概率分布q(z丨x)来近似p(z丨x),同时需要最大化概率值p(x)。衡量两个概率分布之间的差异可使用KL散度,如果对这一概念不清楚,请查阅相关材料。后验概率及其近似分布的KL散度为


D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})]=\mathrm{E}_{\mathbf{z} \sim q}[\ln q(\mathbf{z} | \mathbf{x})-\ln p(\mathbf{z} | \mathbf{x})](4)


根据贝叶斯公式


p(\mathbf{z} | \mathbf{x})=\frac{p(\mathbf{x} | \mathbf{z}) p(\mathbf{z})}{p(\mathbf{x})}(5)


将式5代入式4,由于p(x)与z无关,因此


\mathrm{E}_{\mathbf{z} \sim q}[\ln p(\mathbf{x})]=\ln p(\mathbf{x})


从而可以得到


\begin{aligned} D[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})] &=\mathrm{E}_{z-q}\left[\ln q(\mathbf{z} | \mathbf{x})-\ln \frac{p(\mathbf{x} | \mathbf{z}) p(\mathbf{z})}{p(\mathbf{x})}\right] \\ &=\mathrm{E}_{z-q}[\ln q(\mathbf{z} | \mathbf{x})-\ln p(\mathbf{x} | \mathbf{z})-\ln p(\mathbf{z})]+\ln p(\mathbf{x}) \end{aligned}


对上式进行变形,可以得到


\ln p(\mathbf{x})-D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})]=\mathrm{E}_{z-q}[\ln p(\mathbf{x}, \mathbf{z})]-D_{\mathrm{KL}}[q(\mathbf{z} |) \| p(\mathbf{z})](6)


在这里令变分下界函数为


L=\mathrm{E}_{z-q}[\ln p(\mathbf{x} | \mathbf{z})]-D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z})]


由于KL散度非负,因此是对数似然函数的变分下界。式6是VAE的核心,左侧需要最大化lnp(x)同时最小化


D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})]


此KL散度可以看做是误差项。q(z丨x)产生z,根据这个可以复现出。只要函数的容量足够,此误差项的值非常小。


式6左侧的项不易优化,而最大化左侧等价于最大化右侧。右侧的项即变分下界函数容易优化,因为q(z),q(z丨x)以及p(x丨z)均被限定为类型已知的概率分布,通常为正态分布。因此优化下降函数的问题为优化这些概率分布的参数问题。


变分下界函数定义了一个编码器-解码器结构。q(z丨x)充当编码器的角色,将x编码为z。更准确地说,给定一个x,输出其对应的隐变量的概率分布。需要注意的是编码器的输出为隐变量的均值和方差,而非隐变量本身。p(x丨z)充当解码器,从z重构出x。


VAE的系统结构图下图5所示。编码器由神经网络实现,其输入为样本向量,输出为该样本向量的隐变量所服从的正态分布的均值与方差。即编码器的输出为专属于输入样本的隐变量的概率分布。接下来根据该概率分布进行采样,得到隐变量值,然后由解码器进行映射,重构出样本向量。解码器同样由神经网络实现。


图-5 VEA 的系统结构


与标准的自动编码器相比,这里最大的不同是编码器输出的是概率分布的参数,而不是直接的隐编码向量,且在重构时进行了随机采样,从而注入了随机性。而标准的自动编码器只能原样重构出输入样本数据。



训练算法


式6左侧在最大化lnp(x)的同时最小化


D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})]


p(z丨x)无法得到解析解,通常是复杂的概率分布。方程左侧的第二项致力于q(z丨x)使得靠近p(z丨x)。由于使用了高容量的函数,因此可以确信q(z丨x)能够逼近p(z丨x)。在理想情况下,这会导致方程左侧的KL散度值为0,变成直接优化lnp(x)。这样避免了直接处理难以计算的p(z丨x),而是用一个易于计算的q(z丨x)近似代替它。


现在的问题是如何计算式6的右端。首先需要确定q(z丨x)的形式,通常的选择是使用多维正态分布,即


q(\mathbf{z} | \mathbf{x})=N(\mathbf{z} | \boldsymbol{\mu}(\mathbf{x} ; \boldsymbol{\theta}), \mathbf{\Sigma}(\mathbf{x} ; \boldsymbol{\theta}))


在这里μ 和Σ 是正态分布的均值向量和协方差矩阵,它们都是确定性的函数,输入值为x ,参数为θ 。这两个函数通过训练得到,由神经网络实现,是编码器的输出值。一般将协方差矩阵限定为对角矩阵。


方程6右侧的第二项


D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z} | \mathbf{x})]


是两个正态分布之间的KL散度,称为KL损失,需要最小化该值。根据两个正态分布的KL散度的计算公式,其值为


D_{\mathrm{KL}}\left[N\left(\boldsymbol{\mu}_{0}, \boldsymbol{\Sigma}_{0}\right) \| N\left(\boldsymbol{\mu}_{1}, \boldsymbol{\Sigma}_{1}\right)\right]=\frac{1}{2}\left(\operatorname{tr}\left(\boldsymbol{\Sigma}_{1}^{-1} \boldsymbol{\Sigma}_{0}\right)+\left(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}\right)^{\mathrm{T}} \boldsymbol{\Sigma}_{1}^{-1}\left(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{0}\right)-k+\ln \frac{\left|\boldsymbol{\Sigma}_{1}\right|}{\left|\boldsymbol{\Sigma}_{0}\right|}\right)


如果第二个正态分布为标准正态分布,即迫使编码器输出的概率分布靠近标准正态分布,则上式可以简化为


D_{\mathrm{KL}}[N(\boldsymbol{\mu}, \mathbf{\Sigma}) \| N(\mathbf{0}, \mathbf{I})]=\frac{1}{2}\left(\operatorname{tr}(\mathbf{\Sigma}(\mathbf{x}))+(\boldsymbol{\mu}(\mathbf{x}))^{\mathrm{T}} \boldsymbol{\mu}(\mathbf{x})-k+\ln \boldsymbol{\Sigma}(\mathbf{x})\right)


此时可以根据编码器的输出值直接计算得到。式21-5右侧的第一项称为重构损失,需要使用采样来估计数学期望


\mathrm{E}_{\mathbf{z} \sim q}[\ln p(\mathbf{x} | \mathbf{z})]


的值,要得到精确的估计值需要将大量的z送入g进行计算,成本太高。可以使用随机梯度下降法,每次只用一个z,将p(x丨z)作为此数学期望的近似。用训练样本集中的所有样本进行梯度下降,即可优化式6定义的目标函数。对于训练集中的所有样本,训练时的优化目标为最大化下面的值


\begin{array}{c}{\mathrm{E}_{\mathrm{x}-\mathrm{p}}\left[\ln p(\mathrm{x})-D_{\mathrm{KL}}\left[q(\mathrm{z} | \mathrm{x}) \|_{p}(\mathrm{z} | \mathrm{x})\right]\right]=} \\ {\quad \mathrm{E}_{\mathrm{x}-\mathrm{p}}\left[\mathrm{E}_{z-q}[\ln p(\mathrm{x} | \mathrm{z})]-D_{\mathrm{kL}}[q(\mathrm{z} | \mathrm{x}) \| p(\mathrm{z})]\right]}\end{array}


其中D为训练样本集,需要对训练集中的所有样本计算期望值。实现时在每次迭代时从训练样本集中选取一个样本x,经过编码器计算得到q(z丨x),然后从这个概率分布采样出一个z,计算目标函数的梯度值然后用梯度下降法计算下式的极大值


\ln p(\mathbf{x} | \mathbf{z})-D_{\mathrm{KL}}[q(\mathbf{z} | \mathbf{x}) \| p(\mathbf{z})]


上式做了简化,用lnp(x丨z)替代了


\mathrm{E}_{\mathbf{z} \sim q}[\ln p(\mathbf{x} | \mathbf{z})]


,去掉了对q(z丨x)的依赖。因为假设p(x丨z)为正态分布,因此极大化等价于极小化欧氏距离,下面给出证明。


由于p(x丨z)各个分量的方差相等,因此-lnp(x丨z)与x和g (z;θ)之间的欧氏距离的平方成正比,解码器可以用欧氏距离损失,即计算重构值g (z;θ)与x之间的误差。由于


\begin{aligned} p(\mathbf{x} | \mathbf{z}) &=\frac{1}{(2 \pi)^{\frac{d}{2}}|\mathbf{\Sigma}|^{\frac{1}{2}}} \exp \left(-\frac{1}{2}(\mathbf{x}-g(\mathbf{z} ; \mathbf{\theta}))^{\mathrm{T}} \mathbf{\Sigma}^{-1}(\mathbf{x}-g(\mathbf{z} ; \boldsymbol{\theta}))\right) \\ &=\frac{1}{(2 \pi)^{\frac{d}{2}} \sigma^{d}} \exp \left(-\frac{1}{2 \sigma^{2}}(\mathbf{x}-g(\mathbf{z}, \boldsymbol{\theta}))^{\mathrm{T}}(\mathbf{x}-g(\mathbf{z} \boldsymbol{\theta}))\right) \end{aligned}


因此有


\begin{aligned}-\ln p(\mathbf{x} | \mathbf{z}) &=-\ln \frac{1}{(2 \pi) \sigma^{d}}-\ln \left(\exp \left(-\frac{1}{2 \sigma^{2}}(\mathbf{x}-g(\mathbf{z} ; \boldsymbol{\theta}))^{\mathrm{T}}(\mathbf{x}-g(\mathbf{z}, \boldsymbol{\theta}))\right)\right) \\ &=-\frac{d}{2} \ln (2 \pi)-d \ln \sigma+\frac{1}{2 \sigma^{2}}(\mathbf{x}-g(\mathbf{z} \boldsymbol{\theta}))^{\mathrm{T}}(\mathbf{x}-g(\mathbf{z} \boldsymbol{\theta})) \end{aligned}


上式前两项为常数,第三项即为欧氏距离损失函数。因此极小化-lnp(x丨z)等价于极小化此欧氏距离。这里需要设定超参数σ的值,合理的设置该值存在困难。




还有一个问题未解决,正向传播时需要对q(z丨x)进行采样,根据此概率密度函数采样出一个z作为解码器的输入,采样操作不可导,因此无法进行反向传播。解决这个问题的方法是重参数化技术,将采样操作作为神经网络的输入,而不是网络的某一层。给定μ(x)和Σ(x),即q(z丨x)的均值和协方差,为了从N (μ(x), Σ(x))采样出样本,可以先从N(0,I)采样出一个ε,然后计算z每个分量的值


\mu_{i}+\sigma_{i} \varepsilon_{i}


即可得到想要的随机样本。采用这种方式之后损失函数对均值和方差均变为可导,可用标准的反向传播算法计算编码器和解码器各层参数的梯度值。


预测时用训练的模型生成样本,此时只需要解码器而无需再使用编码器。具体做法是从N(0,I)采样出一个随机ε,然后送入解码器进行计算,即可得到样本。需要注意的是,这时的的随机噪声不再用均值和方差进行变换,因为编码器输出的概率分布已经接近标准正态分布,这是训练时的目标之一。


参考文献

[1] Kingma D P, Welling M. Auto-Encoding Variational Bayes[J]. international conference on learning representations, 2014.

发布于 2019-08-30

文章被以下专栏收录