度量对比敏感度阈值的自适应优化方法研究

王丽萍，陆佳峰，邱飞岳，陈 宏

1(浙江工业大学 计算机科学与技术学院，杭州 310023)

2(浙江工业大学 教育科学与技术学院，杭州 310023)

E-mail：690589058@qq.com

1 引 言

对比敏感度函数(contrast sensitivity function)描述对比度阈值(contrast threshold)的倒数随空间频率的变化，是视觉科学中衡量人眼对物体分辨边界清晰程度的工具.对比敏感度作为衡量人类视觉的重要指标之一，不仅可以用于预测视觉神经病理学，还可用于其他视觉感知等研究[1，2].

对比敏感度阈值的检测早在1760年Bouguer就使用2根不同距离的蜡烛通过照亮屏幕和不透明的杆模拟得到[3].随后由于印刷的发展出现了条栅图法和字母表法，条栅图法引入调制传递函数，依据灰度不同调整曲线的变化，制成宽窄和明暗不同的条栅图，然后通过调整条栅的空间频率和对比度来测量对比敏感度阈值.但是存在被试接受度低、猜测率较高、耗时过长和固定的空间频率和对比度等问题.字母表法也是引入调制传递函数，不同的是它将条栅改为大小和明暗不同的字母或数字.相比条栅图法，字母表法更易被识别，被试接受度高，但是仍然存在耗时过长和固定空间频率和对比度的缺陷.

随着计算机编程和硬件的不断发展，可以通过计算机产生不同空间频率和对比度的刺激.同时由于心理测量函数(psychometric function)描述了被试者随着刺激强度变化而发生知觉变化的反应，成为测量对比敏感度阈值的重要工具.随着近几十年的研究，如何更好的拟合心理测量函数一直是研究的热点，并且有源源不断的研究者投身其中[4-7].研究者们研究心理测量函数时更多的是预先假定一个数学模型，然后估计这个模型的参数，以便更高效地测量对比敏感度阈值.

早期的恒定刺激法[8]，实验者选定刺激范围，并在此范围内随机的选择刺激，在一次实验中，刺激可能会被重复选择，实验完后根据一般线性法进行拟合.一般线性法是对实验数据进行反向累计高斯分布变化，然后引入线性回归进行分析.但是一般线性法难以对心理测量函数的另外两个参数(猜测率和失误率)进行估计[9]，因而限制了它的发展.随后便出现了最大似然法[10]和贝叶斯法[11].贝叶斯法的出现解决了一般线性法的上述问题，成为一种参数估计的常用办法.贝叶斯法主要是把先验的约束和实验中获得的信息相联系，按一定准则整合了观测到的数据和对模型已有的认知，得出后验信息，然后根据得到的后验信息估计未知的参数.

由于拟合方法的不断发展和计算机的普及，自适应测量方法的提出改善了采样低下的问题，加快阈值的估计，因此被广泛应用[12-15]，如变阶梯法[16]和贝叶斯自适应法.

变阶梯法主要使刺激可以分布在阈值附近，其原理是设计一个上升和下降的固定值，当被试者回答正确时，下降的固定值会减去1，然后再次将重复的刺激呈现给被试者，直到下降的固定值为0时，刺激的强度将会减小.同样，若被试者回答错误，上升的固定值会减去1，直到上升的固定值为0时，刺激的强度将会增加.但是由于变阶梯法能选择的上升和下降的固定值有限，例如2进1退、3进1退或者4进1退，导致能选择的阈值收敛的点有限.

基于贝叶斯法的一系列自适应方法相继被提出，如Watson 和 Pelli提出了QUEST法[17]，得到了在心理测量函数阈值处的有效采样，提高对比敏感度阈值测量的速度.随后，King-Smith等人提出了ZEST法[18]，计算了后验的平均值来估计阈值.其后，Kontsevich 和 Tyler提出了psi法[19]，它不仅估计了心理测量函数的阈值，还对斜率也进行了采样估计.Prins则在psi法的基础上，将心理测量函数的另外两个参数(猜测率和失误率)加以优化，提高阈值估计的精度，提出了psi-marginal法[20].

但是，psi-marginal法在测量对比敏感度的阈值时，由于需要测量5个甚至更多的空间频率，导致完成一次实验需要耗费较多的时间.故本文提出了一种基于扰动算子的psi-marginal自适应优化刺激强度算法(psi-marginal/d).该算法可以根据被试者的响应，在预期下一次呈现刺激强度的基础上，设计扰动算子，使呈现的刺激强度更快地朝阈值方向逼近，提高收敛速度，加快对比敏感度阈值的测量.

2 相关工作

2.1 测量工具

拟合心理测量函数是数据建模中通用方法的一种变体.建模数据一般分为3个步骤，首先是选择一个模型，并且调整参数以最小化误差度量或者损失函数，其次是得出参数的误差估计值，最后评估模型和数据之间的拟合度.

典型的心理测量函数被定义为一个S型曲线，它的横坐标由刺激强度给出，纵坐标表示被试者的回答正确的可能性，如图1所示.常用的一般形式如公式(1)所示：

Ψ(x；α，β，γ，δ)=γ+(1-γ-δ)F(x；α，β)

(1)

图1 S型心理测量函数Fig.1 S-type psychometric function

从图1可以发现，在较低刺激强度下，被试者回答正确的可能性较低，是因为被试者在刺激强度较低时的回答具有随机性；随着刺激强度的不断增加，在阈值附近处于快速上升阶段；当刺激强度增加到一定值时，又趋于稳定值.由此可见，心理测量函数可以被4个参数所定义：快速上升阶段的刺激强度(阈值α)、阈值α处曲线的坡度(斜率β)、在较低刺激强度下的渐近线(猜测率γ)和在较高刺激强度下的渐近线(失误率δ)[9，21].常用的S型模型有逻辑斯特模型(Logistic function)、耿贝尔模型(Gumbel function)和基于韦伯分布的韦伯模型(Weibull function)等.

2.2 贝叶斯方法

贝叶斯方法源于Thomas Bayes生前为解决一个“逆概”问题而撰写的一篇文章.贝叶斯法是一种常用的参数估计方法[22-24].假设先验概率为p(B)，其包含了对参数B的先验认知.然后通过观测数据A，将其整合进来，先验概率p(B)将被更新，得到后验概率p(B|A)，从而完成贝叶斯拟合.其贝叶斯原理如公式(2)所示：

p(B|A)=p(A|B)p(B)p(A)

(2)

其中数据A的概率p(A)由参数B的条件概率p(A|B)和参数B的概率p(B)乘积之和得出，如公式(3)所示：

p(A)=∑mj=1p(A|Bj)p(Bj)

(3)

其中m为参数B的最大范围，即∑mj=1p(Bj)的值为1.

2.3 相关算法

up/down算法需要设定初始值、单次上升或下降的间隔值和两个固定值.这两个固定值分别为上升值和下降值，被试回答正确，下降值减1；回答错误，上升值减1，当下降值为0时，呈现的刺激强度下降一个设定的下降间隔值，同时下降值恢复原值；当上升值为0时，呈现的刺激强度上升一个设定的上升间隔值，同时上升值恢复原值，重复上述过程，直到达到实验指定的最大代数.

QUEST算法设定了心理测量函数的模型和固定的先验概率，通过计算最小化二项式标准差和心理测量函数斜率之比的平方作为下一次呈现的刺激强度，在得到被试的反馈后通过贝叶斯法则估计得到后验概率，通过最大化后验概率来估计心理测量函数的参数，重复上述过程直到达到最大代数.

Psi-marginal算法是对psi算法的改进，基本的psi算法主要包括计算被试对下一次刺激的反馈概率，然后通过贝叶斯规则进行后验估计，得到使期望熵达到最小值时的刺激强度，将其作为下一次呈现的刺激，随后在得到被试对呈现刺激的反应后计算获得参数估计，重复上述过程直到实验代数达到指定的最大代数.

相较psi算法，psi-marginal算法在后验分布中添加了失误率，防止当固定失误率和实际实验中出现的失误率不一致时出现较大阈值和斜率估计的偏差[25].同时通过减少边缘阈值和斜率(或者仅边缘阈值)后验分布的期望熵来选择刺激.前者优化阈值和斜率的估计，并且允许将失误率作为目标，后者仅优化阈值的估计，但是如果它有助于阈值参数的估计，则允许将斜率和失误率作为目标.

3 psi-marginal/d算法

3.1 设计扰动算子

在一些以阈值为主要关注点的心理测量实验(如对比敏感度)中，由于实验要在多个空间频率[26]上进行，导致完成一次测量需要耗费较多的时间.

另外如果在将心理测量函数准确的定位在强度轴上后，斜率就成为主要关注的对象.而在这一阶段，熵的分布将呈现出具有两个局部最小值的形状，并且随着试验的进行，呈现的刺激将在两个最小值之间交替.随着试验进一步迭代，最终这两个最小值会位于两个较相近的概率水平相对应的强度处.

为了能够使刺激强度更早地出现在阈值附近，以提高收敛速度，本文提出了一种基于扰动算子的自适应调整刺激强度的方法，该方法会在下一次呈现地刺激强度上添加扰动，扰动算子如公式(4)所示.这种方法可以根据当前被试者的反馈，在使期望熵达到最小值时的最小刺激强度上动态调整刺激强度，使其更快地向阈值方向靠近.

dt=μ·range(x)t·12πe-(xt+1-x-)22

(4)

公式(4)中的μ为扰动系数，范围为(0，1)，其值的大小将直接影响偏移量.range(x)为刺激强度的范围，其值设定是为了试验前期可以充分带动偏移量，防止dt在试验前期就指向零.分母的t为当前实验代数，随着试验的进行，扰动算子将逐渐趋于零，刺激强度最终在阈值附近波动.因为试验测量的是人的心理测量函数，符合正态分布，所以扰动算子基于正态分布构建，其最高点的横坐标为刺激强度的均值.

3.2 算法框架

本文提出的psi-marginal/d算法步骤如算法1所示.

算法1.基于扰动算子的psi-marginal自适应优化刺激强度算法

输入：测量函数Ψλ(x)，最大迭代数N；

输出：测量函数的4个参数λ

1.初始化先验概率密度p0(λ)和条件概率的查找表p(r|λ，x)

2.计算下一次得到反馈后的概率pt(r|x)

3.通过贝叶斯法则估计每个测量函数的概率pt(λ|x，r)

4.估计测量函数在整个参数空间内概率密度函数的熵Ht(x，r)

5.估计每个刺激的期望熵E[Ht(x)]

6.找出使期望熵达到最小时的刺激强度xt+1

7.添加扰动算子后的刺激强度x′t+1

8.将x′t+1作为刺激后获得的反馈计算得到后验分布pt+1(λ)

9.根据步骤8的后验分布得出新的参数估计λt+1

10.如果当前迭代数等于最大迭代数N，停止并输出参数λ；否则跳转到步骤2

算法首先定义测量函数Ψλ(x)，λ是S型测量函数的4个参数(α，β，γ，δ)的向量，x是刺激强度变量.本文的测量函数用到耿贝尔模型，如公式(5)所示.

Ψ(x；α，β，γ，δ)=γ+(1-γ-δ)(1-e-eβ(x-α))

(5)

在每次实验开始前，需要初始化两个函数，一个是心理测量函数先验概率密度p0(λ)，由λ的4个参数的标准正态分布相乘得到，另一个是条件概率的查找表p(r|λ，x)，r是被试者的反应.查找表由公式(6)和公式(7)得出.

p(success|λ，x)=Ψλ(x)

(6)

p(failure|λ，x)=1-Ψλ(x)

(7)

步骤2中计算下一次试验呈现刺激强度x后得到的反应r的概率pt(r|x)可以由公式(8)得出.

pt(r|x)=∑λp(r|λ，x)pt(λ)

(8)

然后根据因刺激强度x而产生的反应r，通过贝叶斯法则估计每个测量函数的后验概率pt(λ|x，r)，其可以用于估计测量函数在整个参数空间内的概率密度函数的熵Ht(x，r)，由公式(9)和公式(10)得到.

pt(λ|x，r)=pt(λ)p(r|λ，x)∑λp(r|λ，x)pt(λ)

(9)

Ht(x，r)=-∑λpt(λ|x，r)log(pt(λ|x，r))

(10)

之后根据步骤2的概率pt(r|x)和步骤4的熵Ht(x，r)估计每个试验刺激的期望熵E[Ht(x)]，如公式(11)所示.随后找出最小化期望熵的刺激强度，添加如公式(4)的扰动算子作为刺激强度得到被试的反应后，得出后验分布pt+1(λ)，如公式(12)-公式(14)所示.

E[Ht(x)]=Ht(x，success)pt(success|x)+
Ht(x，failure)pt(failure|x)

(11)

xt+1=argminxE[Ht(x)]

(12)

x′t+1=xt+1+dt，r=failure
xt+1-dt，r=success

(13)

pt+1(λ)=pt(λ|x′t+1，rt+1)

(14)

其中公式(14)表示为从公式(9)中得到满足刺激强度和被试反馈条件的分布，并且将当前得到的后验分布作为下一次迭代的先验分布.

最后基于后验概率分布得出一个新的参数估计λt+1，如公式(15)所示.

λt+1=∑λλpt+1(λ)

(15)

跳转到步骤2继续实验，直到实验代数达到最大的迭代数，最终得到心理测量函数的4个参数.

4 实验结果和分析

4.1 实验参数设置

为了得到真实准确的测试数据，同时防止当同一个人在进行多次实验时出现不确定性因素，导致对比实验的结论不具有说服性，本次实验的数据由某三甲医院眼科中心在OPTEC6500视功能检测仪上测量得出.OPTEC6500视功能检测仪是基于ANSI公布的CSF测试的统一规范标准制造的[27]，测量的数据具有高准确性.测得的数据为10个人在5个空间频率上的共50个对比度阈值，空间频率的单位为cycle/deg，具体的阈值如表1所示.

表1 5个空间频率的对比度阈值Table 1 Contrast thresholds for five spatial frequencies

实验通过测量上述对比敏感度的5个空间频率上的对比度阈值来比较up/down自适应算法[28]、QUEST自适应算法[17]、psi-marginal自适应算法[20]和psi-marginal/d自适应算法这4种算法.up/down自适应算法的初始值设为刺激范围的均值，步长设为刺激的最小取值步长.psi-marginal/d自适应算法的扰动系数设为0.3.4种对比算法在10名被试者的5个空间频率上分别运行100代.

本次实验中的3种自适应算法已经被整合到Palamedes Toolbox工具箱[29]内.实验平台为Intel(R)Core(TM)i5-8250U CPU @ 1.60GHz 1.80GHz的64位Windows10系统，本文所有对比实验均在MATLAB R2014b版本下运行.

4.2 实验结果及分析

在运行完所有的自适应算法后，将被试1在同一空间频率下4种对比算法呈现的刺激强度进行整合，如图2所示.图2(a)-图2(e)这5个子图分别展示了5种空间频率下刺激强度的变化.每幅图的横坐标代表实验运行代数，规定每次实验的运行代数为100代，纵坐标表示刺激的强度，图中的4种不同的折线描述了4种自适应算法在不同空间频率下的刺激强度.

图2 5个空间频率刺激强度图Fig.2 Five spatial frequency stimulus intensity map

从图2中可以发现，up/down算法呈现的刺激强度波动呈直线下降趋势，直到在阈值附近后才平稳收敛，但是收敛的代数也有较大差异，如图2(e)中，就和前面4幅图的收敛代数差距较大.而QUEST算法前期浮动较大，在阈值上下来回波动，较长代数之后才逐渐收敛到阈值.psi-marginal算法和psi-marginal/d算法都是在前期快速下降后开始有小幅度波动，然后较早开始收敛于阈值.

为了方便比较，设定连续3次呈现刺激强度和阈值差距为5%以内的试验代数为收敛代数，将同一空间频率下的同种算法的10名被试去除收敛代数的最大值和最小值，然后求均值同时四舍五入得出10名被试在不同空间频率下的平均收敛代数，如表2所示.从表2可以发现，up/down算法大约需要50代才会收敛，但是18空间频率下，收敛代数只要24.QUEST算法则相对平稳，大约34代便可以收敛到阈值.psi-marginal算法普遍在25代开始收敛，但是在3空间频率下，收敛代数却为38.而psi-marginal/d算法则最为稳定，约20代便收敛到阈值，且相比于psi-marginal算法，在5个空间频率下提高的平均收敛速度约为27.75%，其中在3空间频率下提高了约47.73%的收敛速度.

表2 4种对比算法平均收敛代数Table 2 Average convergence iteration times of the four comparison algorithms

为了更直观地比较，将表2的5种空间频率作为横坐标，收敛代数作为纵坐标，绘制了4种对比算法在5种空间频率下的平均收敛代数折线图，如图3所示，图中的4种不同的折线表示4种自适应算法.

图3 4种对比算法在5种空间频率下的平均收敛代数Fig.3 Average convergence iteration times of four comparison algorithms at five spatial frequencies

由图3可知，QUEST算法、psi-marginal算法和psi-marginal/d算法总体上比up/down算法能够更早地收敛到阈值，并且在40代试验前基本收敛.但是当空间频率为18时，up/down算法的平均收敛代数只有24左右，这是由于up/down算法的收敛完全依赖初始值的设定如图4所示，不同初始值收敛至同一阈值，所需的代数完全不同.由此可知，因为实验up/down算法的初始值设为刺激范围的均值，而18空间频率的阈值都很靠近此均值，故平均收敛代数较低.

图4 不同初始值下up/down算法刺激强度呈现图Fig.4 Up/down algorithm stimulus intensity presentation diagram under different initial values

综上所述，up/down算法较为不稳定，其收敛代数极度依赖初始值的设定，一旦设定的初始值离阈值较远，收敛速度就会下降.

QUEST算法则相对来说较为稳定，其5种空间频率的平均收敛代数都在34左右，虽然在3空间频率上的平均收敛代数比psi-marginal算法略少，但是其整体的平均收敛代数较高.

由图3可以发现，本文提出的psi-marginal/d算法比psi-marginal算法除3空间频率外普遍提高了约20%的收敛速度，psi-marginal算法普遍要25代试验才开始收敛，而psi-marginal/d算法平均只要20代就开始收敛，特别是在3空间频率下可以发现psi-marginal算法平均收敛代数约为38，而psi-marginal/d算法的平均收敛代数仍然为20左右，提高了接近47%的收敛速度.并且从表1的3空间频率可以发现，3空间频率的阈值普遍较小，较为靠近刺激范围的边缘.

实验表明，up/down算法不稳定，收敛速度依赖初始值的设定，QUEST算法较为稳定，但收敛速度相对来说较慢，psi-marginal算法整体收敛速度尚可，但是当阈值在刺激范围边缘时，收敛速度不理想，而本文所提出的psi-marginal/d算法可以普遍加快收敛，提高算法的收敛速度，特别是当阈值靠近刺激范围的边缘时，收敛速度提升较大.

5 结束语

本文提出了一种基于扰动算子的psi-marginal自适应优化刺激强度算法，该算法通过设计扰动算子，加快对阈值的求解.不管阈值是在刺激的边缘还是中间，psi-marginal/d算法都可以加快求解阈值的收敛速度.特别是当阈值在刺激范围的边缘时，psi-marginal/d算法收敛速度的提升较大.QUEST算法和up/down算法虽然最终都可以收敛到阈值，但是试验耗时过长，而本文提出的psi-marginal/d算法在准确求解阈值时，有较好的收敛速度，同时若阈值在刺激范围的边缘时，依旧保持较好的收敛速度.