1. 程式人生 > >DVL聲學多普勒速度儀技術剖析(理論模擬篇)

DVL聲學多普勒速度儀技術剖析(理論模擬篇)

     聲學多普勒速度儀(DVL)是一種測量相對於水底速度的聲納裝置。早年DVL一開始被設計主要是用於船舶的定位導航中,但是近年來隨著水下機器人領域,尤其是民用領域的興起(據有關機構測算這一市場規模已達100億美元),DVL的作用顯得愈發重要,因為水下機器人想要在水下實現自主導航有三樣東西必不可少,一個是姿態測量單元,一個是速度測量單元,一個是位置測量單元。姿態測量單元通常是通過MEMS慣導元件或者光纖慣導元件等測量得到的,而MEMS通常價格非常便宜,非常容易獲得,因此它可以被大量廣泛的用在水下機器人中;位置測量單元可以由GPS、USBL等裝置構成,獲取的難度相對容易,成本也不算太高;最後一個速度測量單元目前最有效的手段就是通過DVL來獲取,但其成本高昂,售價一般在十幾萬以上甚至更高,因此在目前民用的水下機器人上很少有它的身影。目前水下機器人市場還沒有完全爆發,主要原因就是相關的感測器不夠便宜,許多關鍵部件的成本沒有降下來,國內外在這一領域還沒有出現獨角獸。筆者正是意識到了這一點,所以才下定決心來研究DVL,希望有志同道合的朋友能和筆者一起將DVL的成本降下來,讓水下機器人市場徹底爆發的這一天來的更早一些。

1.多普勒測速原理

1.1多普勒效應

   多普勒效應是由於聲源與觀測者之間存在著相對運動,使觀測者聽到的聲音訊率不同於聲源所發出聲音訊率的現象。這是由奧地利物理學家Christian Doppler 於十九世紀在聲學領域首先發現的物理學原理。

   典型的多普勒效應發生於聲源和觀測者之間。在聲源靜止條件下,當觀測者向著聲源運動時,他收到的聲波頻率高於他相對聲源靜止時收到的聲波頻率;當觀測者遠離聲源而去時,他收到的聲波頻率低於他相對聲源靜止時收到的聲波頻率。聲源於觀測者之間的相對運動速度越大,則所產生的聲波頻率變化也就越大。有關多普勒效應的研究不僅在聲學領域可以看到,利用多普勒效應進行的研究還可擴充套件到電磁學,光學等其它領域。

     下面推導任意方向上的多普勒效應公式。如圖 1 所示,S 為聲源,D為觀測者,\overrightarrow{S{S}'}為聲源的運動向量,\overrightarrow{D{D}'}為觀測者的運動向量,這時從聲源 S 到觀測者 D 傳播方向的向量 SD 隨時在變。設聲源在時刻t=t0 和t0+dt 的位置分別為S 和 S′ ,相位分別為\varphi{\varphi }' ,其中相位的增量為d\varphi ={\varphi }'-\varphi =2\pi f_{t}dt。在t0時刻,由聲源 S 發出的相位為 \varphi的聲波傳播到 D,觀測者的位置在 D ;t0+dt 時刻,聲源到達 S′ ,由於觀測者也在運動,所以,聲源 S′ 發出的相位為\varphi的聲波只能到達觀測者的新位置D′,觀測者從D走到D′所用的時間dt′和他感受到的頻率f_{r},與聲源的dtf_{t}是不一樣的。在D接收的頻率f_{t},在D′接收的頻率f_{r},所接收聲波的相位變化了d\varphi =2\pi f_{r}dt^{'}

,由於相位增量是一樣的,即d\varphi =2\pi f_{t}dt=2\pi f_{r}dt^{'}

                                                                                          f_{r}/f_{t} = dt/dt^{'}                                                                                                   (1)

     相位\varphi從聲源 S 傳播到觀測者的位置 D 的時刻為 t=t_{0}+\left \| \vec{SD} \right \|/c;而相位\varphi +d\varphi由聲源 S′ 傳播到觀測者的位置 D′ 的時刻為t=t_{0}+dt+\left \| \vec{{S}'{D}'} \right \|/c。二者之差即為 

                                                                              d{t}'=dt+\left \| \vec{{S}'{D}'} \right \|/c-\left \| \vec{SD} \right \|/c                                                                               (2)

如圖 1 所示,從 {S}'{D}'作SD的垂線,令相應的垂足分別為S_{0}D_{0}。由於S_{0}D_{0}{S}'{D}'的長度相差高階無窮小量,可認為二者大小近似相等,於是 

                                                          \left \| \vec{SD} \right \|-\left \| \vec{{S}'{D}'} \right \|\approx \left \| \vec{SD} \right \|-\left \| \vec{S_{0}D_{0}} \right \|=\left \| \vec{SS_{0}} \right \|+\left \| \vec{DD_{0}} \right \|

                                                                                                                          =\left \| \vec{S{S}'} \right \|\cos \theta +\left \| \vec{D{D}'} \right \|\cos \psi                                        (3)

式中,θ 是 \vec{S{S}'}\vec{SD}之間的夾角,ψ 則是\vec{S{D}'}\vec{DS}之間的夾角。

   由於聲波在水體介質中的傳播速度遠大於聲源和觀測者的運動速度,可認為聲源和觀測者近似作勻速直線運動,因此\left \| \vec{S{S}'} \right \|=v_{t}dt\left \| \vec{D{D}'} \right \|=v_{r}d{t}'。由式(2)和式(2-3)解得 

                                                                       d{t}'=dt-v_{t}dt\cos \theta /c-v_{r}d{t}'\cos \psi /c                                                                              (4)

由此得到

                                                                                 \frac{f_{r}}{f_{t}}=\frac{dt}{d{t}'}=\frac{c+v_{r}\cos \psi }{c-v_{t}\cos \theta }                                                                                            (5)

這便是多普勒效應的普遍公式。

1.2多普勒測速公式

不難看出,當θ 、ψ 都等於0 時,式(5)可以過渡到共線情形

                                                                                           f_{r}=\frac{c+v_{r}}{c-v_{t}}f_{t}                                                                                                   (6)

由式(6)可知,當聲源與觀測者相互靠近時,f_{r}> f_{t};當聲源與觀測者相互遠離時,f_{r}< f_{t}

       對於DVL而言,收發合制換能器。換能器首先作為聲源,水底作為觀測者,於是根據是(6)水底接收的聲波頻率為:

                                                                                            {f}'=\frac{c+v_{r}}{c-v_{t}}f_{t}                                                                                                  (7)

之後水底作為聲源,而換能器則作為觀測者,於是換能器接收的聲波頻率為:                                                                                              f_{r}=\frac{c+v_{t}}{c-v_{r}}{f}'                                                                                                  (8)

整合式(7)與式(8),得到:

                                                                                       f_{r}=\frac{(c+v_{r})(c+v_{t})}{(c-v_{t})(c-v_{r})}f_{t}                                                                                       (9)

換能器相對靜止、即當v_{t}=0時,v_{r}就是水底相對換能器的徑向流速。於是,式(9)經過整理得到:

                                                                            f_{r}=\frac{c+v_{r}}{c-v_{r}}f_{t}=\frac{1+2v_{r}/c+(v_{r}/c)^{2}}{1-(v_{r}/c)^{2}}f_{t}                                                                     (10)

由此經過整理可以得到通用的徑向測流基本公式

                                                                                     v_{r}=\frac{f_{r}-f_{t}}{f_{r}+f_{t}}c=\frac{f_{d}}{2f_{t}+f_{d}}c                                                                                 (11)

其中f_{d}=f_{r}-f_{t}為多普勒頻移。由於水中聲速 c 約為1500m/s,遠大於水體流速v_{r},於是一般可以忽略式(10)中v_{r}/c的二次項,再經過整理就可以得到近似的多普勒徑向測流公式:

                                                                                            v_{r}=\frac{c}{2f_{t}}(f_{r}-f_{t})=\frac{\lambda }{2}f_{d}                                                                             (12)

式(2-12)中,\lambda =c/f_{t}為換能器發射的訊號波長,f_{d}為收發訊號之間的多普勒頻移。由此通過測量f_{d}就可以得到換能器相對徑向水底速度v_{r},並可以進一步得到水平與垂直方向的速度。

2.DVL換能器陣型與座標變換

2.1Janus陣型結構

本節介紹多普勒測流常用的四波束Janus陣型結構。具有該換能器陣型結構的DVL可以同時發射和同時接收四個窄波束聲波,再經處理就可以得到相對於聲學換能器的徑向流速資訊。在此定義DVL的右手座標系D 。D 的原點位於 Janus 陣型結構的等效中心O,它的三個軸x_{d}y_{d}z_{d}分別指向測流系統所規定的前、右和下的方向。DVL的四波束Janus陣型結構是由四個換能器構成,各換能器收發聲波的軸線與測流系統座標系z軸之間的夾角為α,如圖 2 所示。

                                                          圖2.四波束Janus陣型結構

      另外,各換能器收發聲波的軸線在基陣座標系 xoy 平面的投影則形成兩條相互垂直的直線,且波束 1 的投影與x 軸之間的夾角為 β,如圖3 所示。通常測流系統成“X”字形式安裝時的 β 等於 45°,成“十”字形式安裝時的 β 等於 0°。配有四波束 Janus 陣型結構的DVL在工作時,換能器同時發射的窄波束聲波,沿著各自換能器的軸線方向傳播。在傳播的過程中,每路聲波被該路水體散射,會有散射回波被各換能器接收。DVL換能器和水體之間的徑向運動使得收發聲波頻率之間存在多普勒頻移。通過測量多普勒頻移就可利用式(12)解算出徑向流速v_{1}v_{2}v_{3}v_{4}

                                                      圖3.四波束Janus陣型結構

2.2基陣座標系下的速度

      設水底在D下的速度向量為\mathbf{V}_{d}=[v_{dx},v_{dy},v_{dz}]^{T}。同時假設DVL所測水底是均勻水平的。配有四波束 Janus 陣型結構的DVL利用任意三個聲波所測流速就可以轉換為D下的流速向量\mathbf{V}_{d}。即使因為某一波束受到干擾不能正常測量時系統也可以工作,因此提高了測量的可靠性。以第 1、2 和 3 波束測得的流速列向量 \mathbf{V}_{b}=[v_{1},v_{2},v_{3}]^{T}為例,得到\mathbf{V}_{b}\mathbf{V}_{d}的關係為:

                                                                                              \mathbf{V}_{b}=\mathbf{A}\mathbf{V}_{d}                                                                                                   (13)

                                                                                      \mathbf{V}_{d}=(\mathbf{A}^{T}\mathbf{A})^{-1}\mathbf{A}^{T}\mathbf{V}_{b}                                                                                        (14)

其中

                                                                     \mathbf{A}=\begin{bmatrix} \sin \alpha \cos \beta &\sin \alpha \sin \beta &\cos \alpha \\ -\sin \alpha \cos \beta&\sin \alpha \cos \beta &\cos \alpha \\ -\sin \alpha \cos \beta &-\sin \alpha \sin \beta & \cos \alpha \end{bmatrix}                                                                       (15)

                                                \left ( \mathbf{A}^{T}\mathbf{A} \right )^{-1}\mathbf{A}^{T}=\begin{bmatrix} 1/\sin \alpha \cos \beta &-1/\sin \alpha \sin \beta &-1/\sin \alpha \cos \beta \\ 1/\sin \alpha \sin \beta&1/\sin \alpha \cos \beta &-1/\sin \alpha\sin \beta \\ 1/\cos \alpha &1/\cos \alpha & 1/\cos \alpha \end{bmatrix}                                                   (16)

      實際上,所測水底很少是完全均勻的,即DVL在同一深度上所測相對水底的速度大小與方向往往不同。水底的不均勻性會引入不同程度的速度固有誤差。四波束 Janus 陣型結構較三波束 Convex 陣型結構可以多獲得第 4 波束的冗餘資訊,因此可以驗證水底的均勻性情況。為了量化水底不均勻性造成的影響,可以引入流速固有誤差公式,這是評估資料質量是否有效的重要因素。具體為:

                                                                                         v_{error}=v_{1}-v_{2}+v_{3}-v_{4}                                                                              (17)

通過分析可知,只要所測水底均勻,則無論系統怎樣搖擺v_{error}都基本趨近於零,且四路回波認為是有效資料;如果所測流場不均勻,也要看v_{error}來確定不均勻的程度,用於判斷四路回波資料的有效性。

      四波束 Janus 陣型結構較三波束 Convex 陣型結構還有一個優勢,就是可以更有效地減小由DVL基陣搖擺引入的流速測量誤差,這是通過四個波束所測流速轉化為D 下的流速向量\mathbf{V}_{d}來實現的。設水底的四個徑向速度構成的列向量為\mathbf{V_{b}}=[v_{1},v_{2},v_{3},v_{4}]^{T},由於水底在D下的流速列向量為\mathbf{V_{d}}=[v_{dx},v_{dy},v_{dz}]^{T},於是得到\mathbf{V_{b}}\mathbf{V_{d}}之間的關係為:

                                                                                                     \mathbf{V_{b}}=\mathbf{BV_{d}}                                                                                         (18)

                                                                                              \mathbf{V_{d}}=(\mathbf{B^{T}B})^{-1}\mathbf{B^{T}V_{b}}                                                                              (19)

式(18)和式(19)中的轉換矩陣分別為:

                                                                           \mathbf{B}=\left ( \begin{matrix} \sin \alpha \cos \beta &\sin \alpha \sin \beta & \cos \alpha \\ -\sin \alpha \sin \beta &\sin \alpha \cos\beta & \cos \alpha\\ -\sin \alpha \cos \beta &-\sin \alpha \sin \beta & \cos \alpha\\ \sin \alpha \sin \beta &-\sin \alpha \cos\beta & \cos \alpha \end{matrix} \right )                                                              (20)

                                          \mathbf{(B^{T}B)^{-1}B^{T}}=\left ( \begin{matrix} 1/\sin\alpha\cos\beta &-1/\sin\alpha\sin\beta &-1/\sin\alpha\cos\beta &1/\sin\alpha\sin\beta \\ 1/\sin\alpha\sin\beta &1/\sin\alpha\cos\beta &-1/\sin\alpha\sin\beta &-1/\sin\alpha\cos\beta \\ 1/\cos \alpha & 1/\cos \alpha & 1/\cos \alpha & 1/\cos \alpha \end{matrix} \right )                            (21)

由式(14)和式(19)的關係,最後就可以將測得的徑向速度\mathbf{V_{b}}轉換為測速系統座標系D下的速度\mathbf{V_{d}}

2.3大地座標系下的速度

       相對於聲學換能器的徑向流速和基陣座標系下的流速往往並不是使用者想要得到的最終流速資訊。在實際應用過程中,大地座標系下的流速資訊往往具有更好的實用價值,因此需要將所測換能器徑向流速轉換成基陣座標系下的流速後,再進一步轉換為大地座標系下的流速,而這就涉及到了兩個座標系下的流速座標轉換問題。

                                                              圖4. 座標轉換

    在此定義大地的右手座標系G。大地座標系G的原點位於測速系統的等效中心O,它的三個軸\mathbf{x}_{g}\mathbf{y}_{g}\mathbf{z}_{g}分別指向北、東和垂直向下。接下來,要將基陣座標系D下的速度轉換為大地座標系G下的速度。如圖4 所示設某一右手直角座標系xyz就是基陣座標系D,首先繞oz軸順時針旋轉一個角度ψ,則相應的方向餘弦矩陣\mathbf{\Psi}為:

                                                                              \mathbf{\Psi }=\left ( \begin{matrix} \cos \psi & \sin \psi &0 \\ -\sin \psi & \cos \psi &0 \\ 0 & 0 & 1 \end{matrix} \right )                                                                                    (22)

接下來繞oy軸順時針旋轉一個角度θ,則相應的方向餘弦矩陣\mathbf{\Theta }為:

                                                                               \mathbf{\Theta} =\left ( \begin{matrix} \cos \theta & 0 &-\sin \theta \\ 0 & 1 & 0\\ \sin \theta & 0 & \cos \theta \end{matrix} \right )                                                                                    (23)

最後繞ox 軸順時針旋轉一個角度\varphi,則相應的方向餘弦矩陣\mathbf{\Phi }為: 

                                                                               \mathbf{\Phi }=\left ( \begin{matrix} 1 & 0 &0 \\ 0 & \cos \varphi & \sin \varphi \\ 0 & -\sin \varphi & \cos \varphi \end{matrix} \right )                                                                                        (24)

此時ψ、θ與\varphi為三個尤拉角,他們分別稱為航向角、縱搖角和橫搖角,這是通過DVL內嵌的姿態感測器提供的。

       設經過以上三次旋轉後的新右手直角座標系\mathbf{{x}'{y}'{z}'}為大地座標系G,則\mathbf{xyz}\mathbf{{x}'{y}'{z}'}的座標轉換矩陣\mathbf{R}為:

                                                                                              \mathbf{R}=\mathbf{\Psi \Theta \Phi }                                                                                                  (25)

設速度向量在D內座標為\mathbf{V_{d}},在G內座標為\mathbf{V_{g}},則有如下關係成立: 

                                                                                              \mathbf{V_{g}}=\mathbf{RV_{d}}                                                                                                  (26)

                                                                                     \mathbf{V_{d}}=\mathbf{R^{-1}V_{g}}=\mathbf{R^{T}V_{g}}                                                                                     (27)

由式(26)的關係,就可以將基陣座標系D下的流速轉換為大地座標系G下的流速。最後由式(18)、(19)與式(26)、(27),得到\mathbf{V_{b}}\mathbf{V_{g}}的通用關係: 

                                                                                    \mathbf{V_{g}}=\mathbf{R(B^{T}B)^{-1}B^{T}V_{b}}                                                                                     (28)

                                                                                  \mathbf{V_{b}}=\mathbf{BR^{-1}V_{g}}=\mathbf{BR^{T}V_{g}}                                                                                  (29)

由式(28)與式(29),就可以將波束所測徑向速度直接轉換為大地座標系下的速度。

3.DVL的訊號處理

3.1頻率估計演算法簡述

     如何從多普勒測速回波中提取流速資訊是一個重要的問題。該問題涉及到一定空間範圍內平均多普勒頻移的測量。目前,常用的解決方法是將該問題簡化為單頻復正弦訊號加上高斯白噪聲的頻率估計問題。 

     基於此觀點,有很多可以利用的演算法來估計頻率。其中大部分是基於最大似然(ML)的方法。實際的 ML 估計器需要定位週期圖頻譜譜峰的位置,其中,Palmer 利用了 FFT 來估計譜峰位置,但效能並不理想。Rife 和 Boorstyn給出了最大似然估計器的公式,此時的該估計是無偏的且在信噪比(SNR)達到門限值以上時該估計的均方誤差可以達到克拉美-羅下限(CLRB)。需要注意,以上基於頻域的估計演算法計算時較為複雜,目前並不適合DVL快速靈活的測量要求。

     減小頻域估計器計算量常用的方法是基於時域相位的估計思想。通過將加性高斯噪聲近似為高斯相位噪聲,Tretter給出了一種時域估計器。該近似在高 SNR 條件下是有效的,且作為一種最小二乘估計器其效能等價於 ML的效能。由於需要解決相位模糊問題,該估計器在低 SNR 時效能很難滿足。Kay給出了一個延遲取樣間隔的相位差分估計器,該估計器在高 SNR 條件下可以達到 CRLB。由於改變了提取相位與求和運算的順序,該估計器可以被認為是由原來的一種相位平均估計器(PA)轉變為一種線性預測估計器(LP)。以 Kay 的工作為基礎,不同延遲和不同權重的時域頻率估計器相繼出現。這些估計器在效能上比 Kay 的估計器有所提高,但這是以限制頻率估計範圍或增加計算量為前提的。Brown 和 Wang,肖揚燦分別給出了各自的迴圈頻率估計器。但這些迴圈形式的頻率估計器效能的提高也是以進一步增加計算複雜性為代價的。

      還有其它頻率估計的典型演算法,包括過零檢測法,自適應法等。過零檢測法是在一個過零點處開始以非常高的時鐘脈衝計數,來確定 N 個訊號週期所需的時間,從而進一步估計頻率。該演算法的運算量小且實現邏輯簡單,但當信噪比較低時精度不能令人滿意。自適應法是用最小均方演算法(LMS)自適應調整基於自迴歸 AR 模型的 LP 的係數,然後通過對頻率軸的掃描根據 LP譜峰值來確定訊號的頻率。自適應法可以實現高精度和連續調節的窄帶頻率估計,但該演算法需要一定的自適應時間。由過零檢測法和自適應法的以上特點可以看出,它們無法應用到寬帶回波的多普勒頻移估計中。

   由波束散射模型可以看出,任意時刻的回波都是全部發射訊號在對應某空間範圍內的響應。這一情況使得回波與發射訊號之間差異很大,特別是在發射訊號形式較為複雜的時候。由於訊號形式的不同情況,選擇流速估計器的主要原則是在保證一定測速效能的基礎上,考慮快速靈活的測量要求。可以看到,基於時域相位思想的復自相關演算法因其快速靈活可控的特點而成為了流速估計的一個合適選擇。

3.2復自相關演算法

    復自相關演算法的主要思想是確定兩段回波訊號之間的幅值和相位關係,從而確定兩段回波訊號之間的頻率。復自相關演算法適合於不同訊號下的流速的測量,因而是應用較為廣泛的方法。從功率譜的觀點出發,多普勒頻移測量的問題就轉化為確定觀測訊號功率譜密度的一階矩。在此令回波的複數觀測訊號為:

                                                                                 x(t)=s(t)+n(t), 0<t<T                                                                                (30)

     這裡s(t)是含有多普勒資訊的複數觀測訊號,n(t)為零均值獨立平穩復高斯噪聲。x(t) 對 應 的 協 方 差 函 數 與 功 率 譜 密 度 可 以 分 別 表 示 為 :R(\tau )=R_{s}(\tau )+R_{n}(\tau ),S(f)=S_{s}(f)+S_{n}(f)。於是主要問題就是估計多普勒頻率均值 \mu_{s} (f),它可以用s(t)功率譜的一階矩表示為

                                                                                      \mu_{s} (f)=\frac{\int_{-\infty }^{\infty }fS_{s}(f)df}{\int_{-\infty }^{\infty }S_{s}(f)df}                                                                                      (31)

    首先,將 x(t),s(t)與n(t)的自相關函式表示成極座標形式

                                                                                     R(\tau )=\left | R(\tau ) \right |\exp(j\phi (\tau ))                                                                                  (32)

                                                            R_{s}(\tau )=\left | R_{s}(\tau ) \right |\exp(j\phi_{s} (\tau )),R_{n}(\tau )=\left | R_{n}(\tau ) \right |\exp(j\phi_{n} (\tau ))                                                   (33)

由於假定噪聲是非零獨立平穩高斯的,所以有

                                                                                    R_{n}(\tau )=\left\{\begin{matrix} R_{n}(0), & \tau =0\\ 0, & \tau \neq 0 \end{matrix}\right.                                                                                    (34)因此就有R(\tau )=R_{s}(\tau ),\tau \neq 0。這說明\tau \neq 0時對R(\tau )的計算可得到R_{s}(\tau ),即可以用R(\tau )來代替R_{s}(\tau )

其次,對式(33)的R_{s}(\tau )求導,有

                                                              {R_{s}}'(\tau )=\frac{dR_{s}(\tau )}{d\tau }=[{\left |R_{s}(\tau ) \right |}'+j\left | R_{s}(\tau ) \right |{\phi _{s}}'(\tau )]\exp (j\phi _{s}(\tau ))                                                (35)

由於\left | R_{s}(\tau ) \right |為偶函式,{\left | R_{s}(\tau ) \right |}'為奇函式,所以{\left | R_{s}(0 ) \right |}'=0。因此

                                                                                       {R_{s}}'(0)=jR_{s}(0){\phi _{s}}'(0)                                                                                       (36)

另外根據維納-辛欽定理

                                                                             R_{s}(\tau )=\int_{-\infty }^{\infty }S_{s}(f)\exp (j2\pi \tau f)df                                                                            (37)

                                                                             S_{s}(f)=\int_{-\infty }^{\infty }S_{s}(\tau)\exp (-j2\pi f\tau)d\tau                                                                         (38)

對式(37)的R_{s}(\tau )求導,有

                                                                          {R_{s}}'(\tau )=j2\pi \int_{-\infty }^{\infty }fS_{s}(f)\exp (j2\pi \tau f)df                                                                     (39)

\tau =0時,式(37)與式(39)式變為

                                                                                 {R_{s}}'(0)=j2\pi \int_{-\infty }^{\infty }fS_{s}(f)df                                                                                   (40)

                                                                                       R_{s}(0)=\int_{-\infty }^{\infty }S_{s}(f)df                                                                                       (41)

最後,將式(36)、式(40)與式(41)帶入到式(31),得到

                                                                    \mu_{s} (f)=\frac{\int_{-\infty }^{\infty }fS_{s}(f)df}{\int_{-\infty }^{\infty }S_{s}(f)df}=\frac{-j{R_{s}}'(0)}{2\pi R_{s}(0)}=\frac{1}{2\pi }{\phi _{s}}'(0)                                                             (42)

 由式(42)可以看出,多普勒頻率均值\mu_{s} (f)可以由回波複數觀測訊號自相關函式相位的導數來表示。於是進一步可以看出,對於線性相位情況,就有如下關係成立

                                                                                         \mu_{s} (f)=\frac{\phi (\tau )}{2\pi \tau }                                                                                                   (43)

式(43)就是得到的多普勒頻率均值\mu_{s} (f)的復自相關演算法的估計公式,其中\phi (\tau )是自相關函式的相位。由復自相關函式及其相位的公式

                                                        R(\tau )=\frac{1}{T-\left | \tau \right |}\int_{0}^{T-\left | \tau \right |}x(t+\tau )x^{*}(t)dt,\phi (\tau )=\arctan \frac{Im{R(\tau )}}{Re{R(\tau )}}                                               (44)

對於複數觀測訊號x(t)=x_{r}(t)+jx_{i}(t),\phi (\tau )可以表示為

                                                             \phi(\tau )=\arctan \frac{\int_{0}^{T-\left | \tau \right |}x_{r}(t)x_{i}(t+\tau )-x_{i}(t)x_{r}(t+\tau )dt}{\int_{0}^{T-\left | \tau \right |}x_{r}(t)x_{r}(t+\tau )+x_{i}(t)x_{i}(t+\tau )dt}                                                          (45)

以上的式(45)等價於 Kay