WebRTC視訊接收緩衝區基於KalmanFilter的延遲模型
在WebRTC的視訊處理流水線中，接收端緩衝區JitterBuffer是關鍵的組成部分：它負責RTP資料包亂序重排和組幀，RTP丟包重傳，請求重傳關鍵幀，估算緩衝區延遲等功能。其中緩衝區延遲JitterDelay對視訊流的單向延遲有重要影響，很大程度上決定著應用的實時性。本文不打算全面分析接收端緩衝區的實現細節，只針對緩衝區延遲JitterDelay的計算這一議題進行深入分析。</br>
1 接收端延遲的組成
</br> WebRTC視訊接收端延遲包括三部分：緩衝區延遲JitterDelay，解碼延遲DecodeDelay和渲染延遲RenderDelay。其中DecodeDelay和RenderDelay相對比較穩定，而JitterDelay受傳送端位元速率和網路狀況影響較大。JitterDelay也是造成接收端延遲的最大因素。</br>
緩衝區延遲由兩部分延遲構成：傳輸大尺寸視訊幀造成位元速率burst引起的延遲和網路噪聲引起的延遲。WebRTC採用卡爾曼濾波Kalman Filter估算網路傳輸速率和網路排隊延遲，進而確定緩衝區延遲。本文在理論學習卡爾曼濾波的基礎上，結合WebRTC原始碼，詳細深入分析緩衝區延遲的計算和更新過程。</br>
2 卡爾曼濾波理論學習
</br> 本節主要參考引用文獻[1][2][3]，以此為入門資料學習卡爾曼濾波的基本原理。簡單來說，卡爾曼濾波器是一個最優化自迴歸資料處理演算法，文獻[2]概括卡爾曼濾波的基本思想：</br>
“本質上來講，濾波就是一個訊號處理與變換(去除或減弱不想要的成分，增強所需成分)的過程。卡爾曼濾波屬於一種軟體濾波方法，其基本思想是：以最小均方誤差為最佳估計準則，採用訊號與噪聲的狀態空間模型，利用前一時刻的估計值和當前時刻的觀測值來更新對狀態變數的估計，求出當前時刻的估計值。演算法根據建立的系統方程和觀測方程對需要處理的訊號做出滿足最小均方誤差的估計。”</br>
下面以文獻[1]為基礎，簡要分析卡爾曼濾波的基本過程。</br>
2.1 建立系統數學模型
</br> 首先，我們先要引入一個離散控制過程的系統。該系統可用一個線性隨機微分方程(Linear Stochastic Difference equation)來描述：</br>
X(K) = AX(K-1) + BU(k) + W(k) (2.1.1)
再加上系統的測量值：</br>
Z(k) = HX(k) + V(k) (2.1.2)
上兩式子中，X(k)是k時刻的系統狀態，U(k)是k時刻對系統的控制量。A和B是系統引數，對於多模型系統，他們為矩陣。Z(k)是k時刻的測量值，H是測量系統的引數，對於多測量系統，H為矩陣。W(k)和V(k)分別表示過程噪聲和測量噪聲。他們被假設成高斯白噪聲，他們的協方差分別是Q，R。</br>
2.2 卡爾曼濾波過程
</br> 首先我們要利用系統的過程模型，來預測下一狀態的系統狀態。假設現在的系統狀態是k，根據系統的過程模型，可以基於系統的上一狀態而預測出現在狀態：</br>
X(k|k-1)=AX(k-1|k-1) + BU(k) (2.2.1)
式(2.2.1)中，X(k|k-1)是利用上一狀態預測的結果，X(k-1|k-1)是上一狀態最優的結果，U(k)為現在狀態的控制量。到現在為止，我們的系統結果已經更新。可是，對應於X(k|k-1)的誤差協方差還沒更新。我們用P表示誤差協方差：</br>
P(k|k-1)=AP(k-1|k-1)A’+ Q (2.2.2)
式(2.2.2)中，P(k|k-1)是X(k|k-1)對應的誤差協方差，P(k-1|k-1)是X(k-1|k-1)對應的誤差協方差，A’表示A的轉置矩陣，Q是系統過程的過程噪聲協方差。式子1，2就是卡爾曼濾波器5個公式當中的前兩個，也就是對系統的預測。</br>
現在我們有現在狀態的預測值，然後我們再收集現在狀態的測量值。結合預測值和測量值，我們可以得到現在狀態(k)的最優化估算值X(k|k)：</br>
X(k|k) = X(k|k-1) + Kg(k)[Z(k) - HX(k|k-1)] (2.2.3)
其中Kg為卡爾曼增益(Kalman Gain)：</br>
Kg(k)= P(k|k-1) H’/ (H P(k|k-1) H’+ R) (2.2.4)
到現在為止，我們已經得到了k狀態下最優的估算值X(k|k)。但是為使卡爾曼濾波器不斷的執行下去直到系統過程結束，我們還要更新k狀態下X(k|k)的誤差協方差：</br>
P(k|k)=(I-Kg(k) H) P(k|k-1) (2.2.5)
其中I 為單位矩陣。當系統進入k+1狀態時，P(k|k)就是式子(2)的P(k-1|k-1)。這樣，演算法就可以自迴歸的運算下去。</br>
上述式子1～5是卡爾曼濾波的核心演算法所在，包括預測值計算，誤差協方差計算，最優值估算，卡爾曼增益計算，誤差協方差更新等五個重要步驟。</br>
3 WebRTC中JitterDelay的計算過程
</br> 本節結合WebRTC原始碼，分析視訊接收端緩衝區延遲JitterDelay的計算過程。</br>
3.1 JitterDelay的計算公式
</br> 由第一節分析可知，JitterDelay由兩部分延遲造成：傳輸大幀引起的延遲和網路噪聲引起的延遲。其計算公式如下：</br>
JitterDelay = theta[0] * (MaxFS – AvgFS) + [noiseStdDevs * sqrt(varNoise) – noiseStdDevOffset] (3.1.1)
其中theta[0]是通道傳輸速率的倒數，MaxFS是自會話開始以來所收到的最大幀大小，AvgFS表示平均幀大小。noiseStdDevs表示噪聲係數2.33，varNoise表示噪聲方差，noiseStdDevOffset是噪聲扣除常數30。</br>
解碼執行緒從緩衝區獲取一幀視訊資料進行解碼之前，會根據公式3.1.1計算當前幀的緩衝區延遲，然後再加上解碼延遲和渲染延遲，得到當前幀的預期渲染結束時間。然後根據當前時刻，確定當前幀在解碼之前需要等待的時間，以保證視訊渲染的平滑性。 </br>
3.2 JitterDelay的更新過程
</br> 解碼執行緒從緩衝區獲取一幀視訊資料進行解碼之前，會根據當前幀的大小、時間戳和當前本地時刻，更新緩衝區本地狀態，包括最大幀大小、平均幀大小、噪聲平均值、通道傳輸速率、網路排隊延時等引數。 更新過程如圖1所示：</br>
圖1 視訊接收端緩衝區狀態更新過程.png
首先，確定幀間相對延遲frameDelay和幀間大小差值deltaFS：</br>
frameDelay = t(i) – t(i-1) – (T(i) – T(i-1)) (3.2.1)
deltaFS = frameSize – prevFrameSize (3.2.2)
其中t(i)表示當前幀本地接收時刻，t(i-1)表示上一幀本地接收時刻；T(i)表示當前幀的時間戳，T(i-1)表示上一幀的時間戳。frameDelay表示相鄰兩幀的相對延遲，frameDelay > 0表示幀i相對於幀i-1在路上花費的時間更長。frameSize和prevFrameSize分別表示當前幀和上一幀的大小。</br>
然後計算幀大小的平均值和方差：</br>
avgFrameSize = phi * avgFrameSize + (1-phi) * frameSize (3.2.3)
varFrameSize = phi * varFrameSize + (1-phi) * (frameSize – avgFramesize)^2 (3.2.4)
接下來計算延遲殘差(反映網路噪聲的大小)，並據此計算噪聲均值和方差：</br>
residual = frameDelay – (theta[0] * deltaFSBytes + theta[1]) (3.2.5)
avgNoise = alpha*avgNoise + (1-alpha)*residual (3.2.6)
varNoise = alpha*varNoise + (1 – alpha)*(residual – avgNoise)^2 (3.2.7)
alpha = pow(399/400, 30/fps) (3.2.8)
其中alpha表示概率係數，受幀率fps影響：當fps變低時，alpha會變小，表示當前噪聲變大，噪聲方差受當前噪聲影響更大一些。實時幀率越接近30 fps越好。avgNoise表示自開始以來的平均噪聲。theta[0]表示通道傳輸速率，theta[1]表示網路排隊延遲。</br>
最後一步，卡爾曼濾波器根據當前幀間延遲和幀間大小差值，動態調整通道傳輸速率theta[0]和網路排隊延遲theta[1]，具體過程在下一節進行詳細討論。</br>
至此，JitterDelay的一次更新過程結束。當下一幀資料到來時，使用本次更新結果計算JitterDelay，並再次執行更新過程。</br>
4 卡爾曼濾波更新緩衝區狀態
</br> 本節結合WebRTC原始碼，深入分析卡爾曼濾波更新更新通道傳輸速率theta[0]和網路排隊時延theta[1]的過程，這也是緩衝區延遲估計的核心演算法所在。</br>
4.1 數學符號定義
</br>
X_bar： 向量X
X_hat： 向量X的估計值
X(i)： 向量X的第i個分量
[x y z]：包含元素xyz的行向量
X_bar^T：向量X_bar的轉置向量
E{X}： 隨機變數X的期望
4.2 理論推導過程
</br>
d(i) = t(i) – t(i-1) – (T(i) – T(i-1)) (4.2.1)
t(i)表示當前幀本地接收時刻，t(i-1)表示上一幀本地接收時刻；T(i)表示當前幀的傳送端取樣時刻即時間戳，T(i-1)表示上一幀的取樣時刻。d(i) > 0表示幀i相對於幀i-1在路上花費的時間更長。式子4.2.1是系統的測量值。</br>
d(i) = dL(i) / C(i) + m(i) + v(i) (4.2.2)
dL(i)表示幀i和幀i-1的長度之差，C(i)表示通道傳輸速率，m(i)表示幀i的網路排隊時延，v(i)表示測量噪聲，其協方差為R。其中[1/C(i) m(i)]是我們要求的目標值，即通道傳輸速率和網路排隊時延。式子4.2.2是系統的預測值。</br>
theta_bar(i) = [1/C(i) m(i)]^T 為幀i狀態列向量；
h_bar(i) = [dL(i) 1]^T 為幀間長度差列向量；
theta_bar(i+1) = theta_bar(i) + u_bar(i)； u_bar(i)表示過程噪聲；
Q(i) = E{u_bar(i) * u_bar(i)^T} 表示過程噪聲協方差矩陣；
diag(Q(i)) = [10^-13 10^-3]^T Q(i)是對角陣；
估計過程：</br>
theta_hat(i) = [1/C_hat(i) m_hat(i)]^T； // 目標估計值
z(i) = d(i) – h_bar(i)^T * theta_hat(i-1)
= d(i) – (dL(i)/C_hat(i-1) + m_hat(i-1)) (4.2.3)
上述式子的意義是，用上一時刻估計值估算本時刻的時間消耗：</br>
dL(i)/C_hat(i-1) + m_hat(i-1)
然後用當前觀測值d(i)和估算值求出殘差z(i)。</br>
theta_hat(i) = theta_hat(i-1) + z(i) * k_bar(i) (4.2.4)
上述式子表示i時刻和i-1時刻的狀態迭代關係：i-1時刻的狀態加上i時刻的殘差z(i)和i時刻卡爾曼濾波的乘積，其中i時刻的kalman gain計算公式如下：</br>
k_bar(i)=[(E(i-1)+Q(i))*h_bar(i)]/[var_v_hat(i)+h_bar(i)^T*(E(i-1)+Q(i))*h_bar(i)] (4.2.5)
E(i) = (I–k_bar(i)*h_bar(i)^T)*[E(i-1)+Q(i)] (4.2.6)
var_noise(i) = max(alpha*var_noise(i-1)+(1–alpha)*z(i)^2, 1) (4.2.7)
alpha = pow(399 / 400, 30 / fps) (4.2.8)
var_v_hat(i)=(300*exp[-fabs(deltaFS)/maxFrameSize)+1]*sqrt(varNoise) (4.2.9)
其中I是2*2單位陣，E(i)是誤差協方差(它是卡爾曼濾波迭代的關鍵引數)。var_noise是噪聲方差，var_v_hat是噪聲標準差指數過濾平均後取值，其實就是幀i的測量噪聲協方差R，它對最終計算出的卡爾曼增益有較大影響：var_v_hat(i)較大時，卡爾曼增益較小，表示本次測量中噪聲較大，最終估計值更靠近上次估計值而較少受本次殘差的影響。</br>
4.3 WebRTC程式碼實現
</br> 下面虛擬碼是對WebRTC中VCMJitterEstimator::KalmanEstimateChannel函式的精簡概括，描述了卡爾曼濾波的程式碼實現。</br>
void VCMJitterEstimator::KalmanEstimateChannel(frameDelayMS, deltaFSBytes) {
// 計算誤差協方差和過程噪聲協方差的和：E = E + Q；
_thetaCov[0][0] += _Qcov[0][0];
_thetaCov[0][1] += _Qcov[0][1];
_thetaCov[1][0] += _Qcov[1][0];
_thetaCov[1][1] += _Qcov[1][1];
// 計算卡爾曼增益:
// K = E*h'/(sigma + h*E*h')
// h = [deltaFS 1], Eh = E*h'
// hEh_sigma = h*E*h' + sigma
Eh[0] = _thetaCov[0][0] * deltaFSBytes + _thetaCov[0][1];
Eh[1] = _thetaCov[1][0] * deltaFSBytes + _thetaCov[1][1];
// sigma為測量噪聲標準差的指數平均濾波結果，即測量噪聲協方差R。
double sigma = (300.0 * exp(-fabs(static_cast<double>(deltaFSBytes)) /
(1e0 * _maxFrameSize)) +1) * sqrt(varNoise);
hEh_sigma = deltaFSBytes * Eh[0] + Eh[1] + sigma;
kalmanGain[0] = Eh[0] / hEh_sigma;
kalmanGain[1] = Eh[1] / hEh_sigma;
// 計算殘差，獲得最優估計值
measureRes = frameDelayMS - (deltaFSBytes * _theta[0] + _theta[1]);
_theta[0] += kalmanGain[0] * measureRes;
_theta[1] += kalmanGain[1] * measureRes;
// 更新誤差協方差：E = (I - K*h)*E；
t00 = _thetaCov[0][0]; t01 = _thetaCov[0][1];
_thetaCov[0][0] = (1 - kalmanGain[0] * deltaFSBytes) * t00 -
kalmanGain[0] * _thetaCov[1][0];
_thetaCov[0][1] = (1 - kalmanGain[0] * deltaFSBytes) * t01 -
kalmanGain[0] * _thetaCov[1][1];
_thetaCov[1][0] = _thetaCov[1][0] * (1 - kalmanGain[1]) -
kalmanGain[1] * deltaFSBytes * t00;
_thetaCov[1][1] = _thetaCov[1][1] * (1 - kalmanGain[1]) -
kalmanGain[1] * deltaFSBytes * t01;
}
至此，卡爾曼濾波估計通道傳送速率和網路排隊時延的一次迭代完成。</br>
5 總結
</br> 本文在理論學習卡爾曼濾波基本原理的基礎上，結合WebRTC原始碼，深入分析了WebRTC視訊接收端緩衝區延遲的計算方法和更新過程。通過本文，更進一步加深對WebRTC的學習和了解。</br>
參考文獻
</br> [1] 卡爾曼濾波的原理說明 http://blog.sciencenet.cn/blog-1009877-784428.html [2] 卡爾曼濾波的基本原理及應用[J]. 軟體導刊, Vol.8 No.11, Nov. 2009. [3] An Introduction to the Kalman Filter [J]. University of North Carolina at Chapel Hill Department of Computer Science Chapel Hill, NC 27599-3175 [4] A Google Congestion Control Algorithm for Real-Time Communicationhttps://tools.ietf.org/html/draft-alvestrand-rmcat-congestion-03 [5] Analysis and Design of the Google Congestion Control for Web Real-time Communication[J]. MMSys ’16 Article No.13
作者：weizhenwei 連結：https://www.jianshu.com/p/bb34995c549a 來源：簡書 簡書著作權歸作者所有，任何形式的轉載都請聯絡作者獲得授權並註明出處。