張海勇，陳立軍，周　朋

(1．海軍大連艦艇學(xué)院通信系，遼寧大連116018; 2．中國人民解放軍91404部隊，河北秦皇島066000)

基于改進BP神經(jīng)網(wǎng)絡(luò)的GPS信號反演蒸發(fā)波導(dǎo)研究

張海勇1，陳立軍1，周朋2

(1．海軍大連艦艇學(xué)院通信系，遼寧大連116018; 2．中國人民解放軍91404部隊，河北秦皇島066000)

摘要:海上蒸發(fā)波導(dǎo)的存在，對艦載通信、雷達等電子裝備產(chǎn)生嚴重影響，蒸發(fā)波導(dǎo)的準確預(yù)測是保證艦載電子裝備對波導(dǎo)有效利用的前提。研究基于神經(jīng)網(wǎng)絡(luò)利用GPS信號接收功率反演蒸發(fā)波導(dǎo)的算法流程，提出采用貝葉斯正則化訓(xùn)練方法改進的BP神經(jīng)網(wǎng)絡(luò)的方法原理和計算步驟，實現(xiàn)了對蒸發(fā)波導(dǎo)的仿真反演。仿真結(jié)果表明，貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)能夠?qū)崟r高效反演蒸發(fā)波導(dǎo)高度和大氣修正折射指數(shù)剖面，且在存在干擾的條件下反演結(jié)果更為真實可信。

關(guān)鍵詞:GPS信號;貝葉斯正則化; BP神經(jīng)網(wǎng)絡(luò);蒸發(fā)波導(dǎo)反演

Inversion for evaporation duct using GPS signals based on improved BP neural network

ZHANG Hai-yong1，CHEN Li-jun1，ZHOU Peng2
(1．Department of Communication，Dalian Navy Academy，Dalian 116018．China; 2．No．91404 Unit of PLA，Qinhuangdao 066000，China)

Abstract:The existence of evaporation duct has severely influenced the electronic equipments，such as shipboard communications and radar，and so on．The accurate predication of evaporation duct is the presupposition of efficient using of duct by shipboard electronic equipments．This paper studied the inversion for evaporation duct using the received power of GPS signals which based on neural network，raised the principle and figure of improved BP neural network which is modified by Bayesian regularization algorithm．Stimulated the inversion for evaporation duct．The stimulation result shows: improved BP neural network can inverse the height of evaporation duct and the atmospheric modified reflection index on time and sufficiently，the inversion result under interference is more reliable．

Key words:GPS signals; bayesian regularization algorithm; BP neural network; inversion for evaporation duct

0　引言

目前，國內(nèi)外對海上蒸發(fā)波導(dǎo)預(yù)測的研究主要集中在2個方面:一種是傳統(tǒng)的以莫寧—奧布霍夫(Monin－Obukhov)相似理論為基礎(chǔ)，基于氣象水文參數(shù)預(yù)測蒸發(fā)波導(dǎo)的方法［1－6］;另一種是利用雷達或衛(wèi)星數(shù)據(jù)反演蒸發(fā)波導(dǎo)的方法［7－9］。與傳統(tǒng)的基于氣象水文參數(shù)的蒸發(fā)波導(dǎo)預(yù)測方法相比，蒸發(fā)波導(dǎo)反演方法具有輸入?yún)⒘繂我弧⒁撰@取、反演速度快、時效性高等優(yōu)點。與基于雷達?；夭ǚ囱菡舭l(fā)波導(dǎo)模型相比，基于GPS信號的反演方法不需要主動發(fā)射探測信號，只需要一個被動的海面信號接收機即可，設(shè)備簡單、易行，并具有非常好的隱蔽性，適合軍事通信的需求。人工神經(jīng)網(wǎng)絡(luò)是模擬大腦神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和功能建立的一種信息處理系統(tǒng)［10］，在反演的優(yōu)化算法中，神經(jīng)網(wǎng)絡(luò)方法具有優(yōu)良的非線性映射功能，在解決復(fù)雜的非線性問題上具有一定的優(yōu)勢，因此，本文采用貝葉斯正則化訓(xùn)

練方法的BP神經(jīng)網(wǎng)絡(luò)利用GPS信號接收功率反演蒸發(fā)波導(dǎo)。

1　GPS衛(wèi)星系統(tǒng)簡介及反演步驟

全球定位系統(tǒng)主要包括GPS衛(wèi)星、地面監(jiān)控系統(tǒng)以及GPS信號接收機三大部分。該系統(tǒng)設(shè)計的最初目的是用于目標(biāo)的導(dǎo)航與定位，隨著研究的深入，在環(huán)境與目標(biāo)的遙感方面得到成功應(yīng)用，在對流層蒸發(fā)波導(dǎo)反演中使用的GPS海面散射信號來自于低仰角的GPS衛(wèi)星［11］。利用GPS信號反演蒸發(fā)波導(dǎo)是對波導(dǎo)高度的反演，因為蒸發(fā)波導(dǎo)大氣修正折射指數(shù)模型是單參數(shù)模型，只要反演出波導(dǎo)高度即可以得到大氣修正折射指數(shù)剖面，反演步驟如下:

1)蒸發(fā)波導(dǎo)參數(shù)建模。即構(gòu)建合適的蒸發(fā)波導(dǎo)剖面模型，確定折射指數(shù)參數(shù)的維數(shù)。在構(gòu)建模型時，認為折射指數(shù)剖面隨距離的變化緩慢，可假定在水平距離上不變。本文采用P－J模式給出的單參數(shù)蒸發(fā)波導(dǎo)模型，該模型描述蒸發(fā)波導(dǎo)剖面準確，應(yīng)用廣泛，其表達式如下:式中: z0為海面粗糙度高度，通常取z0=1.5×10－4m; M(z0)為z0處的大氣修正折射指數(shù); d為蒸發(fā)波導(dǎo)高度。

2)正向計算GPS信號接收功率。構(gòu)建合適的GPS信號電波傳播模式，根據(jù)步驟1中確定的波導(dǎo)傳播模型，利用拋物線方程方法計算接收機接收到的GPS功率，在傳播距離(r1，r2，…，rn)處進行離散，得到離散功率Pr。

利用拋物線方程可以有效計算電磁波在大氣環(huán)境中的傳播特征，求解基于GPS信號的拋物線方程過程中需如下初始場及邊界條件:

①初始場計算公式:

式中: E0為入射波的幅值; k為波數(shù); z為垂直高度。

②上邊界條件:處理上邊界條件時，在每一步計算得到的場分布乘以如下窗函數(shù)，保證上邊界區(qū)域內(nèi)場分布緩慢。

③下邊界條件:求解拋物線方程時選取下邊界(z = 0 )是平滑的理想導(dǎo)體平面，即下邊界條件如下:

3)訓(xùn)練神經(jīng)網(wǎng)絡(luò)。將一部分離散功率Pr作為輸入數(shù)據(jù)，所對應(yīng)的蒸發(fā)波導(dǎo)高度作為目標(biāo)數(shù)據(jù)對神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練，得到訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)。

4)反演波導(dǎo)高度。將另一部分離散功率Pr作為輸入數(shù)據(jù)代入已訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)中，反演輸出蒸發(fā)波導(dǎo)高度。

5)誤差檢驗。對反演高度與實際高度進行誤差分析，以檢驗反演結(jié)果與實際的符合程度，如有必要，修改模型重新反演。

2　改進的BP神經(jīng)網(wǎng)絡(luò)

用反向傳播算法進行訓(xùn)練的多層前行網(wǎng)絡(luò)稱為反向傳播神經(jīng)網(wǎng)絡(luò)，簡稱BP網(wǎng)絡(luò)，BP網(wǎng)絡(luò)是一種具有3層或3層以上神經(jīng)元的前向、有導(dǎo)師學(xué)習(xí)方式訓(xùn)練的神經(jīng)網(wǎng)絡(luò)。網(wǎng)絡(luò)結(jié)構(gòu)如圖1所示，網(wǎng)絡(luò)的第1層為輸入層，最后1層為輸出層，中間的若干層稱為隱含層，相鄰層之間通過權(quán)值連接。

圖1　BP網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.1　BP network structure

神經(jīng)網(wǎng)絡(luò)的神經(jīng)元是一個多輸入、單輸出且具有非線性特性的單元，ωji為連接權(quán)值，θj為神經(jīng)元閾值，用bj表示，φ為神經(jīng)元激勵函數(shù)，BP網(wǎng)絡(luò)的神經(jīng)元激活函數(shù)包括雙曲正切S型函數(shù)、對數(shù)S型函數(shù)、線性函數(shù)、階梯函數(shù)和斜坡函數(shù)等，連接權(quán)值和閾值的初值一般在(－1，1)內(nèi)隨機取值。因此神經(jīng)元輸入、輸出關(guān)系數(shù)學(xué)表達式如下:

BP網(wǎng)絡(luò)的學(xué)習(xí)過程主要分為2個階段:第1階段是輸入已知的學(xué)習(xí)樣本，通過設(shè)置的網(wǎng)絡(luò)結(jié)構(gòu)和前一次迭代的權(quán)值和閾值，計算出各個神經(jīng)元的輸出; 第2階段是對權(quán)值和閾值進行修改，從最后一層向前

計算各權(quán)值和閾值對總誤差的影響，據(jù)此對各權(quán)值和閾值進行修改。2個階段反復(fù)進行，直到收斂為止。

雖然BP網(wǎng)絡(luò)理論依據(jù)充分，過程嚴謹并得到廣泛的應(yīng)用，但自身也存在一定的局限性，例如網(wǎng)絡(luò)收斂速度慢，易于陷入局部極小點等問題，考慮到本文對蒸發(fā)波導(dǎo)的預(yù)測具有實時性或近實時性，且預(yù)測高度需滿足一定的精確度，因此本文利用在BP網(wǎng)絡(luò)基礎(chǔ)上采用貝葉斯正則化訓(xùn)練方法的改進的BP網(wǎng)絡(luò)。

貝葉斯正則化神經(jīng)網(wǎng)絡(luò)的一個最主要特點是增強網(wǎng)絡(luò)的泛化能力。正則化是指為提高神經(jīng)網(wǎng)絡(luò)的泛化能力而限制權(quán)、閾值規(guī)模，通過修正神經(jīng)網(wǎng)絡(luò)的訓(xùn)練性能函數(shù)來提高其泛化能力［12］。

在一般情況下，神經(jīng)網(wǎng)絡(luò)的誤差為:

式中: a為樣本對個數(shù); ei為目標(biāo)輸出與實際輸出的差值，正則化理論在誤差基礎(chǔ)上加入一個約束項，代表平滑性約束，這時目標(biāo)函數(shù)為:

式中: Eω為網(wǎng)絡(luò)權(quán)重的平方和;α和β為參數(shù)，用于控制權(quán)值及閾值的分布形式。在引入新的性能指標(biāo)函數(shù)后，正則化方法在確保網(wǎng)絡(luò)誤差盡可能小的情況下盡量減少網(wǎng)絡(luò)權(quán)值，有助于提高神經(jīng)網(wǎng)絡(luò)的泛化能力。

貝葉斯理論用來確定正則化參數(shù)α和β為大小，根據(jù)貝葉斯理論α和β的后驗概率分布為:

式中: M為隱層神經(jīng)元的個數(shù); D為樣本數(shù)據(jù)集。

可以求得:

式中:γ為有效的網(wǎng)絡(luò)參數(shù)數(shù)目，γ=N－2αtr(H)－1; H為目標(biāo)函數(shù)的Hessian矩陣; N為網(wǎng)絡(luò)的所有參數(shù)的數(shù)目。γ的取值范圍為0～N，如果γ非常接近N，則說明網(wǎng)絡(luò)的隱含層神經(jīng)元個數(shù)己經(jīng)足夠大而可以正確表達網(wǎng)絡(luò)輸入?yún)?shù)和網(wǎng)絡(luò)目標(biāo)值之間的關(guān)系，當(dāng)繼續(xù)增大隱含層神經(jīng)元個數(shù)時，且γ保持不變，則表明已獲得了最優(yōu)的網(wǎng)絡(luò)結(jié)構(gòu)及網(wǎng)絡(luò)有效參數(shù)γ［13］。采用貝葉斯正則化訓(xùn)練方法的BP神經(jīng)網(wǎng)絡(luò)主要步驟如下:

1)設(shè)置網(wǎng)絡(luò)結(jié)構(gòu)，初始設(shè)定神經(jīng)網(wǎng)絡(luò)的隱含層數(shù)和連接權(quán)值，初始化網(wǎng)絡(luò)參數(shù)α和β，為α和β賦予先驗初值。

2)利用標(biāo)準的BP算法訓(xùn)練神經(jīng)網(wǎng)絡(luò)，使總誤差最小，并利用Hessian矩陣計算有效網(wǎng)絡(luò)參數(shù)的數(shù)目γ。

3)利用式(9)計算優(yōu)化正則化參數(shù)α和β。

4)對不同的參數(shù)初值重復(fù)以上3步，直到網(wǎng)絡(luò)收斂為止。

3　蒸發(fā)波導(dǎo)的反演

根據(jù)蒸發(fā)波導(dǎo)的反演步驟貝葉斯正則化BP神經(jīng)網(wǎng)絡(luò)方法原理，研究利用GPS信號接收功率對蒸發(fā)波導(dǎo)進行反演。通過拋物線方程的正向出傳播模型可計算得到不同距離處GPS信號接收功率，GPS衛(wèi)星系統(tǒng)參數(shù)如表1所示。

蒸發(fā)波導(dǎo)大氣修正折射指數(shù)模型是單參數(shù)模型，反演出波導(dǎo)高度即可以得到大氣修正折射指數(shù)剖面，因此，本文將不同蒸發(fā)波導(dǎo)高度的GPS信號接收功率與其對應(yīng)的波導(dǎo)高度作為神經(jīng)網(wǎng)絡(luò)訓(xùn)練的輸入數(shù)據(jù)和目標(biāo)數(shù)據(jù)。蒸發(fā)波導(dǎo)高度一般不超過40 m，目標(biāo)數(shù)據(jù)選取高度從2～40 m，每隔2 m取一個波導(dǎo)高度，共20個采樣高度;利用拋物線方程計算接收天線高度15 m處的GPS信號接收功率，輸入數(shù)據(jù)選取距離從31～60 km以0.5 km為間隔的各點的接收功率，共59個采樣頻率。因此輸入為59×20的不同距離和波導(dǎo)高度處接收功率矩陣，輸出為1× 20的不同波導(dǎo)高度矩陣，以此輸入輸出來訓(xùn)練神經(jīng)網(wǎng)絡(luò)，神經(jīng)網(wǎng)絡(luò)參數(shù)如表2所示。

表2　神經(jīng)網(wǎng)絡(luò)參數(shù)Tab.2 Neural network parameters

訓(xùn)練神經(jīng)網(wǎng)絡(luò)共用時1'36″，網(wǎng)絡(luò)實際訓(xùn)練次數(shù)為994次，即在訓(xùn)練994次之后網(wǎng)絡(luò)已滿足精度要求，網(wǎng)絡(luò)訓(xùn)練過程均方差下降曲線如圖2所示。

圖2　Tainbr函數(shù)BP網(wǎng)絡(luò)訓(xùn)練過程曲線Fig.2　The training process curve by trainbr function of BP network

以波導(dǎo)高度為15 m和25 m的GPS信號接收功率作為檢驗數(shù)據(jù)來驗證反演的準確性，圖3給出了波導(dǎo)高度為15 m和25 m時的GPS信號接收功率隨距離的變化情況。

圖3　GPS信號接收功率隨傳播距離的變化曲線Fig.3　The curve of GPS signal reception power with the propagation distance

圖3給出了在沒有干擾和含有干擾情況下不同波導(dǎo)高度條件下GPS信號接收功率隨傳播距離的變化曲線，雖然含有干擾的情況下GPS信號接收功率會發(fā)生振蕩，對反演精度帶來一定的影響，但其與實際情況更為相符。分別對上述情況的接收功率代入訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)進行反演，反演結(jié)果如表3所示。

實際高度/m無干擾時反演高度/m無干擾時相對誤差/%有干擾時反演高度/m有干擾時相對誤差/% 15　15.0268　0.18　14.5230　3.18 25　24.9911　0.04　25.9920　3.97

蒸發(fā)波導(dǎo)高度反演結(jié)果表3所示，可以看出，在未加干擾的情況下，波導(dǎo)高度在15 m和25 m時，反演結(jié)果與實際結(jié)果相差不大，相對誤差小于1%，反演精度較高;而在存在干擾的情況下，蒸發(fā)波導(dǎo)高度反演結(jié)果與實際值差距較大一些，相對誤差在3%左右，但在可接受范圍內(nèi)，考慮到實際海洋氣象環(huán)境下噪聲等干擾是不可避免的，因此，在有干擾情況下的反演結(jié)果更可信。利用神經(jīng)網(wǎng)絡(luò)反演得到的大氣修正折射指數(shù)廓線如圖4所示，可以看出，在未加干擾情況下大氣修正折射指數(shù)廓線與實際廓線幾乎重合，而存在干擾情況下的反演結(jié)果更為真實。

圖4　蒸發(fā)波導(dǎo)大氣修正折射指數(shù)廓線反演結(jié)果Fig.4　The inversion results of evaporation duct atmospheric refraction corrcction index profile

為進一步驗證利用貝葉斯正則化BP網(wǎng)絡(luò)反演結(jié)果的可信性，在其他條件不變的情況下，改變GPS衛(wèi)星系統(tǒng)頻率參數(shù)，以9 GHz頻率的GPS信號接收功率來訓(xùn)練和反演蒸發(fā)波導(dǎo)，以波導(dǎo)高度23 m 和33 m時的接收功率來驗證反演結(jié)果的準確性，反演結(jié)果如表4所示。

表4　反演高度計算結(jié)果Tab.4　The results of the inversion height

圖5　蒸發(fā)波導(dǎo)大氣修正折射指數(shù)廓線反演結(jié)果Fig.5　The inversion results of evaporation duct atmospheric refraction corrcction index profile

從表4和圖5可看出，在GPS衛(wèi)星系統(tǒng)頻率不同的情況下，反演結(jié)果趨勢相同，即在存在干擾條件下反演結(jié)果與實際情況更相符，可信度更高;同時，貝葉斯正則化BP網(wǎng)絡(luò)對不同頻率的GPS信號功率的反演具有較強的適應(yīng)性。

綜上所述，雖然在存在干擾的情況下利用貝葉斯正則化BP網(wǎng)絡(luò)反演蒸發(fā)波導(dǎo)的反演結(jié)果的準確性較差于不存在干擾情況下的反演結(jié)果，但存在干擾情況下的反演結(jié)果與實際環(huán)境下的預(yù)測結(jié)果更為接近，并且反演結(jié)果誤差均在可接受的范圍內(nèi)，說明利用貝葉斯正則化BP網(wǎng)絡(luò)反演蒸發(fā)波導(dǎo)具備可應(yīng)用性;從反演時間來看，每次神經(jīng)網(wǎng)絡(luò)的訓(xùn)練時間均在一分鐘左右，與傳統(tǒng)的基于水文氣象環(huán)境預(yù)測方法的大量數(shù)據(jù)測量和計算相比，時間縮短程度較大，反演速度具有近實時性特點。因此貝葉斯正則化BP網(wǎng)絡(luò)是一種適合艦載裝備實時高效反演蒸發(fā)波導(dǎo)的方法。

4　結(jié)語

隨著國內(nèi)外對蒸發(fā)波導(dǎo)研究的深入，蒸發(fā)波導(dǎo)的反演問題逐步成為對蒸發(fā)波導(dǎo)預(yù)測研究的主要方向，本文基于海洋大氣環(huán)境中蒸發(fā)波導(dǎo)存在特征，提出基于貝葉斯正則化訓(xùn)練方法的BP神經(jīng)網(wǎng)絡(luò)利用GPS信號接收功率反演蒸發(fā)波導(dǎo)的方法，仿真分析了蒸發(fā)波導(dǎo)條件下GPS信號接收功率隨傳播距離變化特征，分別在加入干擾和未加干擾的情況下利用檢驗數(shù)據(jù)對該算法進行了驗證和蒸發(fā)波導(dǎo)的反演。反演結(jié)果表明，貝葉斯正則化BP網(wǎng)絡(luò)能夠反演出蒸發(fā)波導(dǎo)高度和大氣修正折射指數(shù)廓線，在存在干擾情況下反演結(jié)果誤差較大，但仍滿足精度要求，且與實際情況更為相符，結(jié)果可信度更高。

參考文獻:

［1］JESKE H．The state of radar-range prediction over sea．in Tropospheric Radio Wave Propagation．PartⅡ，1971: AGARD－CP－70．

［2］PAULUS R A．Specification for evaporation duct height calculations［D］．Naval Ocean Systems Center，1989．

［3］LUE Musson-Genon，GAUTHIER S，BRUTH E．A simple method to determine evaporation duet height in the sea surface boundary layer［J］．Radio Science，1992，27(5) :635－644．

［4］BABIN S M，YOUNG G S，CARTON J A．A new model of the oceanic evaporation duct［J］．Journal of Applied Meteorology，1997，36:193－204．

［5］劉成國，黃際英，江長蔭，等．偽折射率和相似理論計算海上蒸發(fā)波導(dǎo)剖面［J］．電子學(xué)報，2001，29(7) :970－972．LIU Cheng-guo，HUANG Ji-ying，JIANG Chang-yin，et al．Modeling evaporation duct over sea with pseudo-refractivity and similarity Theory［J］．Acta Electronica Sinica，2001，29(7):970－972．

［6］李云波，張永剛，唐海川，等．基于海氣通量算法的海上蒸發(fā)波導(dǎo)診斷模型［J］．應(yīng)用氣象學(xué)報，2009，20(5) : 628 －633．

［7］Yardim Caglar，Peter Gerstoft，WilliamS．Hodgkiss．Estimation of Radio from Radar Clutter Using Bayesian Monte Carlo Analysis．IEEE Trans，Antennas Propagate，2006，54(4) :1318 －1327．

［8］成印河，何宜軍，趙振維，等．利用AMSR－E衛(wèi)星數(shù)據(jù)反演蒸發(fā)波導(dǎo)高度的神經(jīng)網(wǎng)絡(luò)方法［J］．海洋技術(shù)，2008，27(4) :63－66．CHENG Yin-he，HE Yi-jun，ZHAO Zhen-wei，et al．BP Neural networks method for retrieval of ocean evaporation duct height from AMSR-E［J］．Ocean Technology，2008，27 (4) :63－66．

［9］王波，吳振森，趙振維，等．基于蟻群算法的需達海雜波反演蒸發(fā)波導(dǎo)研究［J］．電波科學(xué)學(xué)報，2009，24(4) :598－604．WANG Bo，WU Zhen-sen，ZHAO Zhen-wei，et al．Retrieving evaporation duct from radar sea clutter based on ant colony algorithm［J］．Chinese Journal of Radio Sience， 2009，24(4) :598－604．

［10］張治國．人工神經(jīng)網(wǎng)絡(luò)及其在地學(xué)中的應(yīng)用研究［D］．吉林:吉林大學(xué)，2006．ZHANG Zhi-guo．Study on artificial neural networks and their applications in geoscience［D］．Jilin: Jilin Univeersity，2006．

［11］張金鵬．海上對流層波導(dǎo)的雷達海雜波/GPS信號反演方法研究［D］．西安:西安電子科技大學(xué)，2012．ZHANG Jin-peng．Methods of retrieving tropospheric ducts above ocean surface using radar seaclutter and GPS signals ［D］．Xi'an: Xi'an Electronic and Enginecring University，2012．

［12］ZHANG Xin-tao．WANG Yu，GAO Yong-sheng．Estimation of non-statistical uncertainty using fuzzy-set theory［M］．Meas．Se i．Technol，2000:430－435．

［13］李鵬．Bayesian算法在某運動軌跡測量儀中的應(yīng)用研究［D］．哈爾濱:哈爾濱工程大學(xué)，2008．LI Peng．The application and research of Bayesian algorithm in a trajectory measuring meter［D］．Harbin:Harbin Engineering University，2008．

作者簡介:張海勇(1966－)，男，博士后，教授，研究方向為通信與信息系統(tǒng)。

基金項目:國家自然科學(xué)基金資助項目(11374001)

收稿日期:2015－01－04;修回日期: 2015－01－26

文章編號:1672－7649(2015) 07－0073－05doi:10.3404/j.issn.1672－7649.2015.07.017

中圖分類號:TP183

文獻標(biāo)識碼:A