[發(fā)明專(zhuān)利]一種全球電離層數(shù)據(jù)同化和預(yù)報(bào)方法

有效

• 申請(qǐng)?zhí)枺?/td>: 202011304313.8
• 申請(qǐng)日: 2020-11-19
• 公開(kāi)（公告）號(hào): CN112649899B
• 公開(kāi)（公告）日: 2023-01-24
• 發(fā)明（設(shè)計(jì)）人: 歐明;陳亮;於曉;許娜;陳龍江;吳家燕;熊雯;馮健;甄衛(wèi)民
• 申請(qǐng)（專(zhuān)利權(quán)）人: 中國(guó)電波傳播研究所（中國(guó)電子科技集團(tuán)公司第二十二研究所）
• 主分類(lèi)號(hào): G01W1/10
• 分類(lèi)號(hào): G01W1/10;G01S19/14;G01S19/37
• 代理公司: 青島博雅知識(shí)產(chǎn)權(quán)代理事務(wù)所(普通合伙) 37317
• 代理人: 封代臣

1.一種全球電離層數(shù)據(jù)同化和預(yù)報(bào)方法，其特征在于，包括如下步驟：

步驟1，獲取全球分布的地基GNSS接收機(jī)、LEO衛(wèi)星掩星、垂測(cè)儀、衛(wèi)星信標(biāo)接收機(jī)和小型光度計(jì)探測(cè)數(shù)據(jù)，分別提取電離層總電子含量TEC和電子密度的觀(guān)測(cè)數(shù)據(jù)，按照設(shè)定的時(shí)間窗口將數(shù)據(jù)存儲(chǔ)在一個(gè)指定的文件內(nèi)，該文件包括：數(shù)據(jù)類(lèi)型編號(hào)、觀(guān)測(cè)時(shí)刻、測(cè)站經(jīng)緯高坐標(biāo)、衛(wèi)星測(cè)站經(jīng)緯高坐標(biāo)、TEC或電子密度測(cè)量值；

步驟2，獲取電離層數(shù)據(jù)同化和預(yù)報(bào)需要的日地空間環(huán)境參量數(shù)據(jù)，包括太陽(yáng)輻射通量和地磁指數(shù)，將數(shù)據(jù)同化所對(duì)應(yīng)的太陽(yáng)輻射通量、地磁指數(shù)及模型運(yùn)行所需要的配置參量，包括運(yùn)行模式設(shè)定、輸入的觀(guān)測(cè)數(shù)據(jù)存儲(chǔ)路徑、輸出數(shù)據(jù)的存儲(chǔ)和網(wǎng)格分辨率設(shè)定，一并存儲(chǔ)到指定的配置文件內(nèi)；

步驟3，構(gòu)建電離層數(shù)據(jù)同化觀(guān)測(cè)算子，將電離層TEC和電子密度數(shù)據(jù)按照“線(xiàn)”型和“點(diǎn)”型觀(guān)測(cè)數(shù)據(jù)進(jìn)行分類(lèi)；

對(duì)于“點(diǎn)”型電子密度觀(guān)測(cè)數(shù)據(jù)，插值方法如下：

其中：Ne表示電子密度，K表示級(jí)數(shù)展開(kāi)的級(jí)數(shù)，r表示地理位置，h_k(r)為基函數(shù)，a_k為插值系數(shù)，基函數(shù)計(jì)算方法如下：

對(duì)于“線(xiàn)”型TEC觀(guān)測(cè)數(shù)據(jù)，根據(jù)接收機(jī)與衛(wèi)星間的幾何位置關(guān)系，TEC和同化估計(jì)的狀態(tài)參量之間的關(guān)系為：

其中：STEC表示電離層總電子含量，N_e表示電離層電子密度值，積分上下限分別為衛(wèi)星和接收機(jī)位置，因此幾何轉(zhuǎn)換矩陣計(jì)算方法如下：

其中ΔR_in表示第i個(gè)TEC觀(guān)測(cè)值在第n個(gè)網(wǎng)格內(nèi)的射線(xiàn)截距，i為T(mén)EC觀(guān)測(cè)值的序號(hào)，N為電子密度網(wǎng)格總數(shù)，K為展開(kāi)的總級(jí)數(shù)；

遍歷同化窗口內(nèi)所有電離層電子密度和TEC數(shù)據(jù)，根據(jù)式(1)和式(3)同化觀(guān)測(cè)方程可構(gòu)建為：

y＝Hx (4)

其中：y為電離層電子密度和TEC數(shù)據(jù)組成的矢量，矩陣H為觀(guān)測(cè)數(shù)據(jù)與同化狀態(tài)參量“電子密度”對(duì)應(yīng)的幾何轉(zhuǎn)換矩陣，x為待估計(jì)的各網(wǎng)格點(diǎn)的電子密度；

將矩陣H按照行優(yōu)先稀疏矩陣CSR格式存儲(chǔ)，即只存儲(chǔ)H矩陣中非零元素所在的起始行的序號(hào)、非零元素所在的列編號(hào)及非零元素本身；

步驟4，構(gòu)建觀(guān)測(cè)誤差協(xié)方差矩陣，根據(jù)不同觀(guān)測(cè)數(shù)據(jù)間觀(guān)測(cè)誤差的不同，計(jì)算數(shù)據(jù)同化的觀(guān)測(cè)誤差協(xié)方差矩陣R，將觀(guān)測(cè)誤差協(xié)方差矩陣R用對(duì)角矩陣表示：

其中i，j表示觀(guān)測(cè)點(diǎn)i和觀(guān)測(cè)點(diǎn)j，d為觀(guān)測(cè)的TEC或電子密度值，R_ij為觀(guān)測(cè)點(diǎn)i與觀(guān)測(cè)點(diǎn)j間的誤差協(xié)方差值，α是比例系數(shù)；

將矩陣R按照CSR格式存儲(chǔ)，即只存儲(chǔ)R矩陣中非零元素所在的起始行的序號(hào)、非零元素所在的列編號(hào)及非零元素本身；

步驟5，構(gòu)建背景場(chǎng)誤差協(xié)方差矩陣，確定電離層相關(guān)距離；

電離層垂直方向相關(guān)距離的計(jì)算方法如下：

分別表示i和j點(diǎn)在垂直方向的電離層相關(guān)長(zhǎng)度，表示i和j點(diǎn)之間的垂直方向電離層相關(guān)距離；

電離層水平相關(guān)距離計(jì)算方法如下：

其中，和L_λ分別表示緯向和經(jīng)向的相關(guān)長(zhǎng)度，θ表示兩個(gè)網(wǎng)格點(diǎn)之間的方位角，表示i和j點(diǎn)之間的水平方向電離層相關(guān)距離；表示i點(diǎn)緯向的相關(guān)長(zhǎng)度；表示j點(diǎn)緯向的相關(guān)長(zhǎng)度；表示i點(diǎn)經(jīng)向的相關(guān)長(zhǎng)度；表示j點(diǎn)經(jīng)向的相關(guān)長(zhǎng)度；

其中，L_z表示垂直方向上的相關(guān)長(zhǎng)度，表示磁緯度，z表示高度，γ是地方時(shí)LT的函數(shù)，表示為：

在存儲(chǔ)數(shù)據(jù)同化背景場(chǎng)誤差協(xié)方差矩陣P時(shí)，設(shè)定的限帶方法是任意兩點(diǎn)間的距離超過(guò)1200km，其背景協(xié)方差為零值，根據(jù)(6)和(7)式，數(shù)據(jù)同化背景場(chǎng)誤差協(xié)方差矩陣P具體表示為：

其中：分別代表背景模型在i點(diǎn)和j點(diǎn)的背景值，e為自然指數(shù)，z_ij代表第i點(diǎn)和第j點(diǎn)在高度上的距離，r_ij代表第i點(diǎn)和第j點(diǎn)的直線(xiàn)距離，表示電離層在高度方向的相關(guān)距離，表示電離層在水平方向的相關(guān)距離，β表示模式誤差與模式值之間的比例系數(shù)，g_ij代表第i點(diǎn)和第j點(diǎn)在水平方向上的大圓距離，計(jì)算方法如下：

式中，是第i點(diǎn)的緯度和經(jīng)度，是第j點(diǎn)的緯度和經(jīng)度；

將矩陣P按照CSR格式存儲(chǔ)，即只存儲(chǔ)P矩陣中非零元素所在的起始行的序號(hào)、非零元素所在的列編號(hào)及非零元素本身；

步驟6，讀取步驟2設(shè)定的配置文件參量，讀取太陽(yáng)輻射通量和地磁指數(shù)，根據(jù)設(shè)定的數(shù)據(jù)同化網(wǎng)格分辨率，運(yùn)行電離層背景模型，計(jì)算背景場(chǎng)電離層電子密度分布x_b；

步驟7，采用限帶卡爾曼濾波算法進(jìn)行地基和天基TEC或電子密度觀(guān)測(cè)數(shù)據(jù)的同化，模型使用的卡爾曼濾波方程如下：

x_a＝x_b+PH^T(R+HPH^T)^-1(y-Hx_b) (14)

其中x_a表示同化后電子密度分析場(chǎng)；

對(duì)矩陣作以下變換：

(R+HPH^T)^-1(y-Hx_b)＝T (15)

(R+HPH^T)T＝(y-Hx_b) (16)

采用廣義最小殘差法求解方程，使用不完全下三角矩陣—上三角矩陣LU分解算法進(jìn)行預(yù)處理，具體方法如下：

C^-1[(R+HPH^T)T-(y-Hx_b)]＝0 (17)

其中：C是預(yù)條件矩陣，同化給出電離層TEC和電子密度現(xiàn)報(bào)結(jié)果后，將相關(guān)現(xiàn)報(bào)結(jié)果存儲(chǔ)到指定的輸出路徑下；

步驟8，利用高斯—馬爾科夫預(yù)報(bào)方法對(duì)同化的電子密度進(jìn)行預(yù)報(bào)，具體公式為：

其中，t表示預(yù)報(bào)時(shí)間，表示t+1時(shí)刻的電離層電子密度背景場(chǎng)，表示t時(shí)刻的電離層電子密度背景場(chǎng)，表示同化后t+1時(shí)刻的電子密度分析場(chǎng)，K為增益矩陣，L代表轉(zhuǎn)換矩陣，其表現(xiàn)形式為對(duì)角矩陣，具體計(jì)算方法為：

其中：i，j表示網(wǎng)格點(diǎn)i和網(wǎng)格點(diǎn)j，ΔT為提前預(yù)報(bào)的時(shí)間間隔，τ表示電離層的時(shí)間相關(guān)尺度，取5h，磁暴期間取值縮短為2h，同化給出電離層TEC和電子密度預(yù)報(bào)結(jié)果后，將相關(guān)現(xiàn)報(bào)結(jié)果存儲(chǔ)到指定的輸出路徑下。

2.根據(jù)權(quán)利要求1所述的全球電離層數(shù)據(jù)同化和預(yù)報(bào)方法，其特征在于：步驟(2)中的運(yùn)行模式為區(qū)域或全球。