ORB-SLAM内的卡方检验
- 1. 概念
- 2. 卡方检验的基本思想
- 3. 卡方检测示例
- 4. ORB-SLAM2中卡方检测剔除外点的策略
- 4.1 示例，卡方检验计算置信度得分: CheckFundamental()、CheckHomography()
Reference:
- 卡方检验(Chi-square test/Chi-Square Goodness-of-Fit Test)
- 卡方检验详解分析与实例
1. 概念
-
卡方值: χ 2 \chi^2 χ2 值表示观察值与理论值之间的偏离程度。计算这种偏离程度的基本思路如下：
- 设
O
O
O 代表某个类别的
观察频数， E E E 代表基于某个假设 H 0 H_0 H0 计算出的期望频数， O O O 与 E E E 之差称为残差； - 残差可以表示某一个类别观察值和理论值的偏离程度，但如果将残差简单相加以表示各类别观察频数与期望频数的差别，则有一定的不足之处。因为残差有正有负，相加后会彼此抵消，综合仍然为 0 0 0，为此可以将残差平方后求和；
- 另一方面，残差大小是一个相对概念，相对于期望频数为 10 10 10 时，期望频数为 20 20 20 的残差非常大，但相对于期望频数为 1000 1000 1000 时 20 20 20 的残差就很小了。考虑到这一点，人们又将残差平方除以期望频数再求和，以估计观察频数与期望频数的差别。
进行上述操作之后，就得到了常用的 χ 2 \chi^2 χ2 统计量，其计算公式为：
χ 2 = ∑ ( O − E ) 2 E = ∑ i = 1 k ( O i − E i ) 2 E i = ∑ i = 1 k ( O i − n p i ) 2 n p i ( i = 1 , 2 , 3 , . . . , k ) \chi^{2}=\sum{\frac{(O-E)^{2}}{E}}=\sum_{i=1}^{k}{\frac{(O_{i}-E_{i})^{2}}{E_{i}}}=\sum_{i=1}^{k}{\frac{(O_{i}-n p_{i})^{2}}{n p_{i}}}\quad{\mathrm{(i=1,2,3,...,k)}} χ2=∑E(O−E)2=i=1∑kEi(Oi−Ei)2=i=1∑knpi(Oi−npi)2(i=1,2,3,...,k)其中， O i O_i Oi 为 i i i 水平的观察频数， E i E_i Ei 为 i i i 水平的期望频数， n n n 为总频数， p i p_i pi 为 i i i 水平的期望概率， i i i 水平的期望频数 E i E_i Ei 等于总频数 n × i n\times i n×i 水平的期望概率 p i p_i pi， k k k 为单元格数。当 n n n 比加大时， χ 2 \chi^2 χ2 统计量近似服从 k − 1 k-1 k−1(计算 E i E_i Ei 时用到的参数个数)个自由度的卡方分布。由卡方的计算公式可知，当观察频数与期望频数完全一致时， χ 2 \chi^2 χ2 值为 0 0 0；观察频数与期望频数越接近，两者之间的差异越小， χ 2 \chi^2 χ2 值越小；反之，观察频数与期望频数差别越大，两者之间的差异越大， χ 2 \chi^2 χ2 值越大。换言之，大的 χ 2 \chi^2 χ2 值表明观察频数远离期望频数，即表明远离假设。小的 χ 2 \chi^2 χ2 值表明观察频数接近期望频数，接近假设。因此， χ 2 \chi^2 χ2 是观察频数与期望频数之间距离的一种度量指标，也是假设成立与否的度量指标。如果 χ 2 \chi^2 χ2 值小，研究者就倾向于不拒绝 H 0 H_0 H0；如果 χ 2 \chi^2 χ2 值大，就倾向于拒绝 H 0 H_0 H0。至于 χ 2 \chi^2 χ2 在每个具体研究中究竟要达到什么程度才能拒绝 H 0 H_0 H0，则要借助于卡方分布求出所对应的 P P P 来确定，这就引出了显著性水平的概念。
- 设
O
O
O 代表某个类别的
-
显著性水平: 显著性水平是估计总体参数落在某一区间内可能犯错的概率，用 α \alpha α 表示，即当原假设为正确时，却把它滤掉了的概率或风险，通常取 α = 0.05 \alpha=0.05 α=0.05 或 α = 0.01 \alpha=0.01 α=0.01。这表明当做出接受原假设的决定时，其正确的可能性(概率)为 95 % 95\% 95% 或 99 % 99\% 99%。
2. 卡方检验的基本思想
卡方检验是以 χ 2 \chi^2 χ2 分布为基础的一种常用假设检验方法，它的无效建设 H 0 H_0 H0 是：观察频数与期望频数没有差别。
该检验的基本思想是：首先假设 H 0 H_0 H0 成立，基于此前提计算出 χ 2 \chi^2 χ2 的值，它表示观察值与理论值之间的偏离程度。根据 χ 2 \chi^2 χ2 分布及自由度可以确定在 H 0 H_0 H0 假设成立的情况下获得当前统计量及更极端情况的概率 P P P。如果 P P P 值很小，说明观察值与理论值偏离程度太大，应当拒绝无效假设，表示比较资料之间有显著差异；否则就不能拒绝无效假设，尚不能认为样本所代表的实际情况和理论假设有差别。
根据自由度和显著性水平查询检验统计量临界值，其中纵轴为自由度，横轴为显著性水平阈值。比如检查单应阵的函数，点到点的重投影误差自由度为
2
2
2，在显著性水平为
0.05
0.05
0.05 时通过查表得阈值为
5.99
5.99
5.99：
3. 卡方检测示例
为了验证肺癌与吸烟的关系，我们得到如下数据：
| 是否肺病患者 | 吸烟 | 不吸烟 | 合计 | 吸烟比例 |
|---|---|---|---|---|
| 是 | 158 | 169 | 327 | 48% |
| 否 | 82 | 311 | 393 | 20% |
| 合计 | 240 | 480 | 720 | 33% |
首先假设吸烟与肺癌两者之间没有关系，我们计算期望值(因为上面吸烟的比例为33%，因此在吸烟与肺癌没有关系的时候，肺癌患者吸烟与不吸烟的比例应该是33%，没有得肺癌的吸烟与不吸烟的比例也应该是33%)：
| 是否肺病患者 | 吸烟 | 不吸烟 | 合计 | 吸烟比例 |
|---|---|---|---|---|
| 是 | 109 | 218 | 327 | 33% |
| 否 | 131 | 262 | 393 | 33% |
| 合计 | 240 | 480 | 720 | 33% |
带入卡方计算公式：
χ
2
=
(
158
−
109
)
2
109
+
(
169
−
218
)
2
218
+
(
82
−
131
)
2
131
+
(
311
−
262
)
2
262
=
60.53
\chi^2=\frac{(158-109)^2}{109}+\frac{(169-218)^2}{218}+\frac{(82-131)^2}{131}+\frac{(311-262)^2}{262} =60.53
χ2=109(158−109)2+218(169−218)2+131(82−131)2+262(311−262)2=60.53自由度的计算方法，可以简单抽象成：(行数-1)*(列数-1)，所以这里的自由度为
1
1
1。
通过查表可得自由度为 1 1 1 时，显著性水平为 0.05 0.05 0.05，当卡方值小于 3.84 3.84 3.84 时，可以接受原假设，即变量之间没有相关性。卡方值越小，不相关的概率越大。现在卡方值远大于 3.84 3.84 3.84，说明两者不相关的概率很小，即抽烟与肺癌有关。
4. ORB-SLAM2中卡方检测剔除外点的策略
就特征点法的视觉SLAM而言，一般会计算重投影误差。具体而言，记
u
\mathbf{u}
u 为特征点的
2
D
2D
2D 位置，
u
ˉ
\bar{\mathbf{u}}
uˉ 为由地图点投影到图像上的
2
D
2D
2D 位置。重投影误差为：
e
=
u
−
u
ˉ
\mathbf{e}=\mathbf{u}-\mathbf{\bar{u}}
e=u−uˉ重投影误差服从高斯分布：
e
∼
N
(
0
,
Σ
)
\mathbf{e}\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma})
e∼N(0,Σ)其中，协方差
Σ
\mathbf{\Sigma}
Σ 一般根据特征点提取的金字塔层级决定。具体的，记提取ORB特征时，图像金字塔的每层缩小尺度为
s
s
s(ORB-SLAM中为1.2)。在ORB-SLAM中假设第
0
0
0 层的标准差为
p
p
p 个pixel(ORB-SLAM中设为了1个pixel)；那么，一个在金字塔第
n
n
n 层提取的特征的重投影误差的协方差为：
Σ
=
(
s
n
×
p
)
2
[
1
0
0
1
]
\boldsymbol{\Sigma}=(s^n\times p)^2\begin{bmatrix}1&0\\ 0&1\end{bmatrix}
Σ=(sn×p)2[1001]式
e
=
u
−
u
ˉ
\mathbf{e}=\mathbf{u}-\mathbf{\bar{u}}
e=u−uˉ 中的误差是一个
2
2
2 维向量，这里阈值不好设置，就把它变成一个标量，计算向量的内积
r
r
r(向量元素的平方和)。但是，这又有一个问题，不同金字塔层的特征点都用同一个阈值，是不是不合理呢？于是，在计算内积的时候，利用协方差进行加权(协方差表达了不确定度)。金字塔层级越高特征点提取精度越低，协方差
Σ
\mathbf{\Sigma}
Σ 越大，那么就有了：
r
=
e
T
Σ
−
1
e
r=\mathbf{e}^{T}\mathbf{\Sigma}^{-1}\mathbf{e}
r=eTΣ−1e利用协方差加权，起到了归一化的作用，具体如上式，可以变为：
r
=
(
Σ
−
1
2
e
)
T
(
Σ
−
1
2
e
)
r=(\boldsymbol{\Sigma}^{-\frac{1}{2}}\mathbf{e})^T(\boldsymbol{\Sigma}^{-\frac{1}{2}}\mathbf{e})
r=(Σ−21e)T(Σ−21e)而：
(
Σ
−
1
2
e
)
∼
N
(
0
,
I
)
(\mathbf\Sigma^{-\frac{1}{2}}\mathbf e)\sim\mathcal N(\mathbf0,\mathbf I)
(Σ−21e)∼N(0,I)为多维标准正态分布，也就是说不同金字塔层提取的特征，计算的重投影误差都被归一化了，或者说去量纲化了，那么，我们只用一个阈值就可以了。
4.1 示例，卡方检验计算置信度得分: CheckFundamental()、CheckHomography()
根据重投影误差构造统计量 χ 2 \chi^2 χ2，其值越大，观察结果和期望结果之间的差别越显著，某次计算越可能用到了外点：
float Initializer::CheckHomography(const cv::Mat &H21, const cv::Mat &H12, vector<bool> &vbMatchesInliers, float sigma) {
const int N = mvMatches12.size();
//取出单应矩阵H各位上的值
const float h11 = H21.at<float>(0, 0);
const float h12 = H21.at<float>(0, 1);
const float h13 = H21.at<float>(0, 2);
const float h21 = H21.at<float>(1, 0);
const float h22 = H21.at<float>(1, 1);
const float h23 = H21.at<float>(1, 2);
const float h31 = H21.at<float>(2, 0);
const float h32 = H21.at<float>(2, 1);
const float h33 = H21.at<float>(2, 2);
const float h11inv = H12.at<float>(0, 0);
const float h12inv = H12.at<float>(0, 1);
const float h13inv = H12.at<float>(0, 2);
const float h21inv = H12.at<float>(1, 0);
const float h22inv = H12.at<float>(1, 1);
const float h23inv = H12.at<float>(1, 2);
const float h31inv = H12.at<float>(2, 0);
const float h32inv = H12.at<float>(2, 1);
const float h33inv = H12.at<float>(2, 2);
vbMatchesInliers.resize(N);//标记是否是内点
float score = 0;//置信度得分
const float th = 5.991;//<---自由度为2(u,v),显著性水平为0.05的卡方分布对应的临界阈值
const float invSigmaSquare = 1.0 / (sigma * sigma);//信息矩阵,方差平方的倒数
for (int i = 0; i < N; i++) {//双向投影,计算加权投影误差
bool bIn = true;
//step1. 提取特征点对
const cv::KeyPoint &kp1 = mvKeys1[mvMatches12[i].first];
const cv::KeyPoint &kp2 = mvKeys2[mvMatches12[i].second];
const float u1 = kp1.pt.x;
const float v1 = kp1.pt.y;
const float u2 = kp2.pt.x;
const float v2 = kp2.pt.y;
// Reprojection error in first image
// x2in1 = H12*x2
//step2. 计算img2到img1的重投影误差
const float w2in1inv = 1.0 / (h31inv * u2 + h32inv * v2 + h33inv);
const float u2in1 = (h11inv * u2 + h12inv * v2 + h13inv) * w2in1inv;
const float v2in1 = (h21inv * u2 + h22inv * v2 + h23inv) * w2in1inv;
const float squareDist1 = (u1 - u2in1) * (u1 - u2in1) + (v1 - v2in1) * (v1 - v2in1);
const float chiSquare1 = squareDist1 * invSigmaSquare;
//step3. 离群点标记上,非离群点累加计算得分
if (chiSquare1 > th)
bIn = false;
else
score += th - chiSquare1;
// Reprojection error in second image
// x1in2 = H21*x1
//step4. 计算img1到img2的重投影误差
const float w1in2inv = 1.0 / (h31 * u1 + h32 * v1 + h33);
const float u1in2 = (h11 * u1 + h12 * v1 + h13) * w1in2inv;
const float v1in2 = (h21 * u1 + h22 * v1 + h23) * w1in2inv;
const float squareDist2 = (u2 - u1in2) * (u2 - u1in2) + (v2 - v1in2) * (v2 - v1in2);
const float chiSquare2 = squareDist2 * invSigmaSquare;
//step5. 离群点标记上,非离群点累加计算得分
if (chiSquare2 > th)
bIn = false;
else
score += th - chiSquare2;
if (bIn)
vbMatchesInliers[i] = true;
else
vbMatchesInliers[i] = false;
}
return score;
}
