H264编码- 码率控制 RQ 模型参数推导过程以及JM代码分析

摘要：本文主要介绍H264码率控制过程中，RQ模型参数更新推导过程，并结合JM19.0代码分析其功能实现。

H264码率控制中比较重要的一个模型是RQ模型，不管是图像级码率控制还是基本单元码率控制都会用到。

RQ模型主要作用是根据目标比特数以及预测MAD，预测得到当前编码单元（帧/宏块）的量化步长QStep，得到量化步长后，才可以进行后续的率失真优化RDO获取最佳编码方式。

RQ模型参数初始值是固定值，肯定不能满足所有的编码场景，因此，在编码过程中需要根据已编码的数据不断更新参数，从而使 RQ 模型更加精准。

1. RQ模型

RQ模型公式如下：
$$Target=c1\ast \frac{MAD}{QStep}+c2\ast \frac{MAD}{QSte{p}^2}$$
其中，$Target$ 表示预测得到的目标比特数，
$MAD$表示编码后数据的平均误差，由于这个时候还没有开始编码，所以这个值也是通过已编码帧的MAD预测当前帧得到的。

$c1、c2$就是我们要不断更新的参数。

2. 参数更新过程

使用线性回归算法来更新RQ模型的参数。更新过程由下面三步组成（这部分代码实现参考JM的updateRCModel函数）：

1. 选择数据集

模型的准确度依赖于做线性回归的数据集大小及其质量。一般来说，选择的数据量越多，越能反映出平均情况，但难以反映出实时情况。例如在场景切换时，需要根据最近的情况及时更新模型。

在滑动窗口内的已编码视频帧的信息$（QSte{p}_i、MA{D}_i、bit{s}_i）$构成了一组数据集。每编码一个帧或基本编码单元，会把该次编码生成的数据样本放入滑动窗，作为数据集的一部分。

2. 计算模型参数

根据滑动窗内数据集（QStep、MAD、编码bit数），使用线性回归法更新RQ模型参数。计算过程在JM的RCModelEstimator函数
使用线性回归方法，统计数据集中所有样本，使得真实编码比特数和RQ模型预测目标比特数最接近。

3. 移除滑动窗内误差太大的数据样本，并重新计算RQ模型参数。

使用上一步骤计算得到的RQ模型参数计算每个样本的均方差，均方差超过阈值K的样本则认为该数据是错误数据，应该从数据集剔除。然后再根据步骤2重新计算一次RQ模型参数。

3. 线性回归推导

参数更新通过线性回归方法实现.
可得到滑动窗内所有样本平方差代价函数如下：
$$E=\sum _{i=1}^n(c_1\ast MA{D}_i/QSte{p}_i+c_2\ast MA{D}_i/QSte{p}_i^2-Targe{t}_i{)}^2$$

令$x=QStep$，表示实际编码使用的量化步长；
$y=bits$，表示最终编码后的真实比特数。

代价函数可表示为
$$E(c_1,c_2)=\sum _{i=1}^n(c_1\ast MA{D}_i/x_i+c_2\ast MA{D}_i/x_i^2-y_i{)}^2$$
我们要做的实际上是求得$c_1,c_2$，使$E$最小。

$分别对c_1和c_2求偏导，可得E的最小值$
$$\frac{\partial E}{\partial c_1}=\frac{\partial {\sum }_{i=1}^n(c_1\ast MA{D}_i/x_i+c_2\ast MA{D}_i/x_i^2-y_i{)}^2}{\partial c_1}=0$$

$$\sum _{i=1}^n\frac{\partial (c_1\ast MA{D}_i/x_i+c_2\ast MA{D}_i/x_i^2-y_i{)}^2}{\partial c_1}=0$$

$$2\ast \sum _{i=1}^n(c_1\ast MA{D}_i/x_i+c_2\ast x_i/x_i^2-y_i)\ast \frac{MA{D}_i}{x_i}=0$$

$$\sum _{i=1}^n(c_1+\frac{c_2}{x_i}-y_i\ast x_i/MA{D}_i)=0$$

同理对$c_2$的偏导可得方程：
$$\sum _{i=1}^n(c_1/x_i+c_2/x_i^2-y_i/MA{D}_i)=0$$

可得矩阵运算方程
$$\left[\begin{array}{cc}{\sum }_{i=1}^{n}1& {\sum }_{i=1}^{n}1/x_i\\ {\sum }_{i=1}^{n}1/x_i& {\sum }_{i=1}^{n}2/x_i^2\end{array}\right]\ast \left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]=\left[\begin{array}{c}{\sum }_{i=1}^n{y}_i\ast x_i/MA{D}_i\\ {\sum }_{i=1}^n{y}_i/MA{D}_i\end{array}\right]$$
求解以上方程可得$c_1,c_2$

JM代码中求解过程在函数RCModelEstimator中，其中用到的几个变量如下：
$a_{00}={\sum }_{i=1}^{n}1,a_{01}={\sum }_{i=1}^{n}1/x_i$
$a_{10}={\sum }_{i=1}^{n}1/x_i,a_{11}={\sum }_{i=1}^{n}2/x_i^2$
$b_0={\sum }_{i=1}^n{y}_i\ast x_i/MA{D}_i,b_1={\sum }_{i=1}^n{y}_i/MA{D}_i$

void RCModelEstimator (VideoParameters *p_Vid, InputParameters *p_Inp, RCQuadratic *p_quad, int n_windowSize, Boolean *m_rgRejected)
{
  int n_realSize = n_windowSize;
  int i;
  double oneSampleQ = 0;
  double a00 = 0.0, a01 = 0.0, a10 = 0.0, a11 = 0.0, b0 = 0.0, b1 = 0.0;
  double MatrixValue;
  Boolean estimateX2 = FALSE;

// 获取滑动窗内有效样本个数
  for (i = 0; i < n_windowSize; i++)
  {// find the number of samples which are not rejected
    if (m_rgRejected[i])
      n_realSize--;
  }

  // default RD model estimation results，m_X1、m_X2对应公式里的c1、c2
  p_quad->m_X1 = p_quad->m_X2 = 0.0;

  for (i = 0; i < n_windowSize; i++)
  {
    if (!m_rgRejected[i])
      oneSampleQ = p_quad->m_rgQp[i]; // m_rgQp存放的是QStep值
  }
  for (i = 0; i < n_windowSize; i++)
  {// if all non-rejected Q are the same, take 1st order model
    if ((p_quad->m_rgQp[i] != oneSampleQ) && !m_rgRejected[i])
      estimateX2 = TRUE;
    if (!m_rgRejected[i])
      p_quad->m_X1 += (p_quad->m_rgQp[i] * p_quad->m_rgRp[i]) / n_realSize;
  }

  // take 2nd order model to estimate X1 and X2
  /*
  求解方程,也就是上面推导的矩阵公式
  | a00  a01 |  |  m_x1 | = | b0 |
  | a10  a11 |  |  m_x2 |    | b1 |
  */
  if ((n_realSize >= 1) && estimateX2)
  {
    for (i = 0; i < n_windowSize; i++)
    {
      if (!m_rgRejected[i])
      {
        a00  = a00 + 1.0;
        a01 += 1.0 / p_quad->m_rgQp[i];
        a10  = a01;
        a11 += 1.0 / (p_quad->m_rgQp[i] * p_quad->m_rgQp[i]);
        b0  += p_quad->m_rgQp[i] * p_quad->m_rgRp[i];
        b1  += p_quad->m_rgRp[i];
      }
    }
    // solve the equation of AX = B
    MatrixValue=a00*a11-a01*a10;
    if(fabs(MatrixValue) > 0.000001)
    {
      p_quad->m_X1 = (b0 * a11 - b1 * a01) / MatrixValue;
      p_quad->m_X2 = (b1 * a00 - b0 * a10) / MatrixValue;
    }
    else
    {
      p_quad->m_X1 = b0 / a00;
      p_quad->m_X2 = 0.0;
    }
  }
  if( p_Vid->type == P_SLICE || (p_Inp->RCUpdateMode == RC_MODE_1 && (p_Vid->number != 0)) )
  {
    p_quad->Pm_X1 = p_quad->m_X1;
    p_quad->Pm_X2 = p_quad->m_X2;
  }
}