陳純毅，熊琴琴，于海洋，姚海峰，宋佳雪，婁巖

（1.長春理工大學(xué)計算機(jī)科學(xué)技術(shù)學(xué)院，吉林 長春 130022；2.長春理工大學(xué)空地激光通信技術(shù)國防重點學(xué)科實驗室，吉林 長春 130022）

1 引言

光波的軌道角動量（OAM,orbital angular momentum）橫向空間模為經(jīng)典和量子光信息科學(xué)中的信息編碼提供了一種額外的自由度，近年來受到國內(nèi)外研究者的廣泛關(guān)注[1-3]。在經(jīng)典光通信中，人們使用OAM 來實現(xiàn)空間模分復(fù)用，以便提高光通信系統(tǒng)容量[1-2]。在光量子密鑰分發(fā)領(lǐng)域，OAM被用于實現(xiàn)高維量子態(tài)，以便把超過1 bit 的信息編碼在一個光子中，從而增加每個光子的密鑰生成率[3]。從是否需要光纜的角度可將光通信分成光纖通信和自由空間光通信2 種形式。在這2 種形式中，前者穩(wěn)定性好，但需要預(yù)先在通信雙方之間鋪設(shè)光纜，實施成本高；后者不需要鋪設(shè)光纜，靈活機(jī)動，實施成本低，但容易受到外界環(huán)境的干擾[4]。綜合比較來說，對于難以鋪設(shè)光纜或者對靈活性要求很高的應(yīng)用場景，自由空間光通信是一個較好的方案。

大氣湍流是自由空間光通信中最重要的外界環(huán)境干擾因素之一[4-5]。為了理解自由空間光通信中大氣湍流的影響，國內(nèi)外學(xué)者對受大氣湍流擾動的光信號的統(tǒng)計特性做了大量研究。目前，學(xué)者已經(jīng)提出了多種描述受大氣湍流擾動的光信號起伏概率分布模型，例如對數(shù)正態(tài)分布、K分布和Gamma-Gamma 分布等[6-7]。然而，這些傳統(tǒng)分布模型在建立時沒有考慮OAM 橫向空間模接收的特殊性。大氣湍流對OAM 光通信的影響表現(xiàn)為信號光功率在不同OAM 橫向空間模之間的隨機(jī)再分布。Chen 等[8]發(fā)現(xiàn)，上述光信號起伏概率分布模型不太適合用來描述因大氣湍流擾動而導(dǎo)致的OAM接收光信號起伏，因此提出了雙Johnson SB分布模型。模擬實驗表明該分布模型能夠很好地描述受大氣湍流擾動的 OAM 接收光信號起伏統(tǒng)計特性。

在用蒙特卡洛方法來分析 OAM 光通信系統(tǒng)性能時，需要模擬生成受大氣湍流擾動的OAM 接收光信號的時域隨機(jī)樣本。受大氣湍流擾動的光信號起伏具有時間相關(guān)性，屬于時域自相關(guān)信號。模擬生成受大氣湍流擾動的OAM 接收光信號的關(guān)鍵是，要使生成的隨機(jī)光信號時域樣本在統(tǒng)計上服從雙Johnson SB分布，同時具有指定的時間自相關(guān)。本文將滿足這2 個條件的隨機(jī)光信號稱為雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號。文獻(xiàn)[9-10]基于隨機(jī)微分方程（SDE,stochastic differential equation）分別生成了具有指數(shù)自相關(guān)函數(shù)的Gamma-Gamma 分布和Johnson SB分布隨機(jī)數(shù)序列。但目前還沒有生成雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號時域樣本的報道。文獻(xiàn)[11]利用自回歸（AR,autoregression）模型生成了具有指定時間自相關(guān)的復(fù)高斯隨機(jī)樣本。文獻(xiàn)[12]利用逆變換法將自相關(guān)高斯隨機(jī)樣本變換為服從特定統(tǒng)計分布的隨機(jī)樣本。然而，采用逆變換法將自相關(guān)高斯樣本轉(zhuǎn)換為雙 Johnson SB分布樣本的過程較復(fù)雜。文獻(xiàn)[13]提出的秩匹配技術(shù)通過排序操作將自相關(guān)瑞利序列的時間自相關(guān)傳遞給服從Nakagami 分布的白噪聲隨機(jī)數(shù)序列，避免了逆變換法中對目標(biāo)分布的求逆運算，對生成服從雙Johnson SB分布的時間自相關(guān)隨機(jī)光信號的時域樣本具有參考價值。

本文基于AR 模型和秩匹配法，提出一種模擬生成雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號的方法。所提方法分為2 個步驟：首先用AR 模型生成具有指定時間自相關(guān)的高斯隨機(jī)樣本；然后生成一個服從雙Johnson SB分布的白噪聲隨機(jī)數(shù)序列，并對該白噪聲隨機(jī)數(shù)序列進(jìn)行適當(dāng)?shù)闹匦屡帕?，以匹配指定的時間自相關(guān)。其中，重排的順序由自相關(guān)高斯過程的統(tǒng)計秩決定。根據(jù)指定的指數(shù)形式的時間自相關(guān)函數(shù)，本文利用所提方法生成了雙Johnson SB分布湍流擾動OAM 光通信接收信號，并將生成信號的自相關(guān)系數(shù)與理論自相關(guān)系數(shù)的均方誤差作為生成信號的自相關(guān)性匹配精度，由此分析了不同的雙Johnson SB分布模型參數(shù)對OAM接收光信號時域樣本的自相關(guān)性匹配精度的影響。

2 OAM 光通信湍流擾動接收信號模擬模型

生成自相關(guān)高斯隨機(jī)數(shù)序列的常用方法為采用線性濾波，高斯白噪聲隨機(jī)數(shù)序列經(jīng)過線性濾波后可獲得指定的自相關(guān)，并且服從其原有的高斯分布[11]。當(dāng)目標(biāo)分布為非高斯分布時，僅進(jìn)行線性濾波不再可行，例如，對指數(shù)序列進(jìn)行線性濾波時，輸出的自相關(guān)隨機(jī)數(shù)序列將不再服從指數(shù)分布。對于雙Johnson SB分布來說，這個結(jié)論同樣適用。為了對雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號進(jìn)行蒙特卡洛模擬，本文提出如圖1 所示的原理框架。首先使用AR 模型生成具有指定自相關(guān)的高斯隨機(jī)信號，將其作為輔助集；然后利用蒙特卡洛接受?拒絕采樣法[15]生成服從雙Johnson SB分布的白噪聲隨機(jī)信號，即樣本集；最后通過秩匹配將樣本集按照輔助集的統(tǒng)計秩進(jìn)行重新排列，在不改變分布特性的情況下使樣本集具有指定的自相關(guān)。

2.1 基于AR 模型的自相關(guān)高斯過程

AR 模型本質(zhì)上是一種線性濾波器，能生成具有指定自相關(guān)的高斯隨機(jī)過程。由文獻(xiàn)[9]可知，指數(shù)形式的歸一化時間自相關(guān)函數(shù)如式(1)所示。

其中，τ表示時延，τc表示時間自相關(guān)長度?；贏R模型生成的自相關(guān)高斯過程稱為AR高斯過程，一個p階AR高斯過程通過時域遞歸[11]表示為

其中，ak(k=1,2,…,p)表示AR 模型系數(shù)，p表示AR 模型的階數(shù)；w(n)表示高斯白噪聲信號，其均值為0、方差為。AR 模型的參數(shù)包括系數(shù)ak和高斯白噪聲的方差，它們和指定的自相關(guān)序列R(n)間的非線性關(guān)系用Yule-Walker方程[14]可表示為

其中，R表示由R(n)構(gòu)成的自相關(guān)矩陣，R(n)表示式(1)對應(yīng)的自相關(guān)值序列，具體為

其中，Δt表示采樣時間間隔；對于本文的p階AR高斯過程，n滿足?p≤n≤p的條件。

本文利用Levinson-Durbin 遞推算法來快速地求解式(3)，得到AR 模型的參數(shù)。該算法的初始值[14]為

對于k=2,3,…,p，遞推寫為[14]

其中，

其中，ρk為k階預(yù)測誤差方差。遞推完成得到參數(shù)集{a1(1),ρ1},{a2(1),a2(2),ρ2},…,{ap(1),ap(2),…,ap(p),ρp}。最后的p階參數(shù)集就是Yule-Walker方程的解，即ap(i)=ai,i=1,2,…,p，。對式(4)的自相關(guān)序列進(jìn)行一次Levinson-Durbin 遞推后，式(9)中的預(yù)測誤差方差便不再隨p的增加而變化，由此確定AR 模型的階數(shù)p為1。

由式(2)可以看出，若要得到x(n)，除了AR 模型的參數(shù)外，還需要計算前p個平穩(wěn)輸入序列元素的初始值[11]。對于本文的一階AR 模型，其初始條件x(0)是均值為0、方差為1 的高斯變量。最后將AR 模型的參數(shù)和初始條件代入式(2)，就得到自相關(guān)高斯過程。

圖1 本文的原理框架

2.2 自相關(guān)雙Johnson SB 分布接收光信號

由自相關(guān)高斯過程生成自相關(guān)非高斯過程的標(biāo)準(zhǔn)方法是非線性逆變換法[12]。在這種方法中，目標(biāo)分布的累積分布函數(shù)（CDF,cumulative distribution function）非常重要，因為在構(gòu)建合適的非線性變換和高斯過程的自相關(guān)函數(shù)時都需要目標(biāo)CDF的反函數(shù)。對雙Johnson SB分布來說，其概率密度函數(shù)比較復(fù)雜，不便于獲得其CDF 的反函數(shù)，因此本文不考慮使用該方法而是利用統(tǒng)計秩來匹配指定自相關(guān)[13]。

統(tǒng)計秩是按指定（遞減或遞增）順序排列的樣本列表中一個值的序數(shù)[13]，本文選擇升序排序。假設(shè)樣本長度為5，對應(yīng)的統(tǒng)計秩向量如圖2 所示，其中，第i個元素表示排列中第i個原始樣本的統(tǒng)計秩。對于平穩(wěn)隨機(jī)過程，樣本的自相關(guān)函數(shù)近似于其統(tǒng)計秩的自相關(guān)函數(shù)[13]。因此，AR 高斯過程和雙Johnson SB分布隨機(jī)過程具有式(10)所示關(guān)系。

圖2 統(tǒng)計秩向量原理

其中，RGau和分別表示AR 高斯過程及其統(tǒng)計秩的自相關(guān)函數(shù)；RJoh和分別表示雙Johnson SB分布隨機(jī)過程及其統(tǒng)計秩的自相關(guān)函數(shù)。若雙Johnson SB分布隨機(jī)過程的統(tǒng)計秩與AR 高斯過程的統(tǒng)計秩相同，則二者統(tǒng)計秩的自相關(guān)函數(shù)相同，為

由式(10)和式(11)可得

由此可知，將服從雙Johnson SB分布的樣本按照AR 高斯過程的統(tǒng)計秩進(jìn)行重排，可以得到與AR高斯過程秩匹配的自相關(guān)雙Johnson SB分布隨機(jī)過程，實現(xiàn)過程如圖3 所示，具體步驟如下。

步驟1生成服從雙Johnson SB分布的白噪聲隨機(jī)信號。雙Johnson SB分布的概率密度函數(shù)（PDF,probability density function）為[8]

其中，4 個模型參數(shù)滿足條件δ1＞0，δ2＞0，γ1≤γ2；pz(·)表示Johnson SB分布的PDF。為了避免數(shù)值擬合雙Johnson SB分布CDF 的反函數(shù)，本文采用蒙特卡洛接受?拒絕采樣法[15]生成服從雙Johnson SB分布的白噪聲隨機(jī)信號，計算方法 如下。

1)分別生成在[c,d]和[0,1]范圍內(nèi)均勻分布的隨機(jī)變量u和v。

圖3 雙Johnson SB 分布湍流擾動OAM 光通信自相關(guān)接收信號的秩匹配流程

3)重復(fù)1)和2)，直到生成N個被接受樣本，其中，N表示樣本長度。

步驟2找到自相關(guān)高斯隨機(jī)信號x(n)的統(tǒng)計秩向量。由圖2 可以看出，原始樣本的統(tǒng)計秩向量與其進(jìn)行遞增排序后的索引向量相同。因此，對x(n)升序排序，得到x(Jn),n=0,…,N?1，其中索引Jn使x(J0)≤x(J1)≤…≤x(JN?1)，Jn即為x(n)的統(tǒng)計秩向量。

步驟3根據(jù)Jn對雙Johnson SB分布白噪聲隨機(jī)信號z(n)進(jìn)行重新排列，使z(n)的最小值放置在與x(n) 最小值相同的位置，z(n)的第二小值位于x(n) 第二小值的位置，以此類推。首先將z(n)升序排序，使z(I0)≤z(I1)≤…≤z(IN?1)，有序值記為z(In),n=0,…,N?1；然后重排z(n)，即

最后得到的重排信號z(n) 就具有指定的時間自相關(guān)并且服從雙Johnson SB分布。

3 仿真實驗與討論

通過模擬仿真實驗來驗證圖1 所示的原理框架，所有仿真的樣本長度N=105、采樣時間間隔Δt=0.1 ms。為了更準(zhǔn)確地描述雙Johnson SB分布OAM 光通信接收信號的自相關(guān)性，所有仿真實驗重復(fù)100 次，并根據(jù)這100 次仿真得到的數(shù)據(jù)來計算自相關(guān)系數(shù)。

圖4(a)為τc=10 ms、δ1=3、γ1=1.1、δ2=5、γ2=2.4的雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號，其中，雙Johnson SB分布的4 個模型參數(shù)的取值源自文獻(xiàn)[8]。圖4(b)為對應(yīng)仿真信號的自相關(guān)系數(shù)隨時延變化的曲線，與式(1)的理論歸一化自相關(guān)函數(shù)基本一致。為了進(jìn)一步分析仿真信號的自相關(guān)性是否與理論值匹配，本文以仿真信號的自相關(guān)系數(shù)與理論自相關(guān)系數(shù)的均方誤差（MSE,mean square error）作為自相關(guān)性匹配精度，即

圖4 雙Johnson SB 分布OAM 光通信接收信號及其自相關(guān)函數(shù)示例

為了驗證所提方法的有效性，根據(jù)文獻(xiàn)[8]提供的查找表，另外設(shè)置了6 組雙Johnson SB分布的模型參數(shù)進(jìn)行仿真實驗，如表1 所示。當(dāng)τc=10ms 時，根據(jù)表1 的參數(shù)仿真生成信號的概率密度如圖5所示。OAM 功率傳輸系數(shù)為接收的特定OAM 模中的功率與特定OAM 模的發(fā)射功率之比[8]。本文通過生成的雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)信號對應(yīng)的功率傳輸系數(shù)擾動來表征OAM光信號功率接收時域擾動。

表1 6 組雙Johnson SB 分布模型參數(shù)值

圖5 根據(jù)表1 的參數(shù)仿真生成信號的樣本概率密度與理論概率密度的對比(τ c=10ms)

圖6 表示當(dāng)τc為10 ms 和20 ms 時，根據(jù)表1的參數(shù)仿真生成信號的自相關(guān)系數(shù)與理論自相關(guān)系數(shù)。從圖6 可以看出，仿真信號的自相關(guān)性基本匹配理論自相關(guān)性。當(dāng)雙Johnson SB分布的形狀與高斯分布比較相似時，誤差會減小。通過式(16)的自相關(guān)性匹配精度對此進(jìn)行詳細(xì)分析，M的取值如圖4 所示。

表1 中模型參數(shù)對應(yīng)的匹配精度如表2 所示。從表2 可以看出，τc=10 ms 時的匹配精度與τc=20 ms 時的精度基本上接近。其中第4 組和第6 組參數(shù)的匹配精度最好，結(jié)合圖5 的概率密度曲線可以看出，這是由于第4 組和第6 組模型參數(shù)對應(yīng)的概率密度曲線形狀類似于高斯分布。值得注意的是，AR 模型的輸出是均值為0、方差為1 的自相關(guān)高斯過程。進(jìn)一步分析可知，當(dāng)進(jìn)行秩匹配的2 種分布相似時，本文所提方法的仿真精度更高。

為分析雙Johnson SB分布的模型參數(shù)值對仿真信號的自相關(guān)性的影響，對表2 進(jìn)一步研究：將第1 組參數(shù)與第6 組參數(shù)的自相關(guān)性匹配精度進(jìn)行比較可以看出，隨著δ1和δ2的增加，接近或者大于2 時，自相關(guān)性匹配更好；根據(jù)第6 組參數(shù)的匹配精度優(yōu)于第4 組參數(shù)的精度和第2 組參數(shù)的匹配精度優(yōu)于第5 組參數(shù)的精度可以發(fā)現(xiàn)，當(dāng)越小時，自相關(guān)性匹配越好；通過對比第3 組與第5 組的匹配精度發(fā)現(xiàn)，當(dāng)?shù)谋嚷式咏? 時，匹配精度更高。對于表2 中的6 組雙JohnsonSB分布的模型參數(shù)，對應(yīng)的自相關(guān)性匹配精度基本保持在10?3以下，再次說明了所提的雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號蒙特卡洛模擬方法的有效性。該方法生成的信號可用于研究大氣湍流對OAM 光通信的影響。

表2 表1 的6 組模型參數(shù)的自相關(guān)性匹配精度

4 結(jié)束語

本文提出一種基于AR 模型和秩匹配的雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號蒙特卡洛模擬方法。由于AR 模型的參數(shù)計算簡單，秩匹配法僅涉及排序，因此本文所提方法具有復(fù)雜度低、便于工程應(yīng)用的優(yōu)點。通過模擬仿真實驗發(fā)現(xiàn)，模擬生成的信號在統(tǒng)計上滿足雙Johnson SB分布且具有指定的時間自相關(guān)，進(jìn)而證明了本文所提方法的可靠性。本文所提方法生成的雙Johnson SB分布湍流擾動OAM 光通信自相關(guān)接收信號，可以為用蒙特卡洛方法分析大氣OAM 光通信系統(tǒng)性能提供支持。

圖6 不同條件下的仿真生成信號自相關(guān)函數(shù)與理論自相關(guān)函數(shù)的對比