吳 夢劉宏偉 王 旭

(西安電子科技大學(xué)雷達(dá)信號處理國家重點(diǎn)實(shí)驗(yàn)室 西安 710071)

一種循環(huán)迭代的MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法

吳 夢*劉宏偉 王 旭

(西安電子科技大學(xué)雷達(dá)信號處理國家重點(diǎn)實(shí)驗(yàn)室 西安 710071)

現(xiàn)有多輸入多輸出(MIMO)雷達(dá)發(fā)射方向圖設(shè)計(jì)是通過優(yōu)化信號相關(guān)矩陣來逼近期望方向圖，可采用凸優(yōu)化方法求解，但其計(jì)算量較大，不利于工程實(shí)現(xiàn)。針對上述問題，該文提出一種循環(huán)迭代的MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法。該方法以加權(quán)最小二乘為準(zhǔn)則，通過引入輔助變量，將對信號相關(guān)矩陣的優(yōu)化問題轉(zhuǎn)化為關(guān)于其Hermite平方根的二次優(yōu)化問題，再以循環(huán)迭代的方式進(jìn)行求解。對于均勻線陣，當(dāng)采用均勻加權(quán)且離散化方位角在歸一化空間頻率域均勻采樣時(shí)，可采用快速傅里葉變換(FFT)的方式進(jìn)行求解，進(jìn)一步提高計(jì)算效率。仿真結(jié)果表明，該方法所得發(fā)射方向圖可以很好地逼近期望方向圖，且具有較高的實(shí)時(shí)性。

多輸入多輸出雷達(dá)；發(fā)射方向圖設(shè)計(jì)；循環(huán)迭代；快速傅里葉變換

1 引言

多輸入多輸出(Multiple Input Multiple Output, MIMO)雷達(dá)作為一種新體制雷達(dá)，現(xiàn)已成為雷達(dá)技術(shù)領(lǐng)域的研究熱點(diǎn)[116]-。根據(jù)天線位置的空間分布，

MIMO雷達(dá)可分為分布式MIMO雷達(dá)[1]和集中式

MIMO雷達(dá)[2]。對于分布式MIMO雷達(dá)，其天線間距較大，各個(gè)天線對目標(biāo)的觀測視角不同，且目標(biāo)回波之間相互獨(dú)立，可以有效克服目標(biāo)的閃爍效應(yīng)，提高雷達(dá)的檢測性能。對于集中式MIMO雷達(dá)，其布陣方式與相控陣?yán)走_(dá)相同，但各個(gè)發(fā)射天線發(fā)射不同的信號，能夠獲得更高的角度分辨率、更好的參數(shù)辨別能力和雜波抑制能力[7]。集中式MIMO雷達(dá)的發(fā)射信號可以是正交的[79]-或是部分相關(guān)的[1012]-，此時(shí)發(fā)射信號相關(guān)矩陣的秩大于等于1，使得集中式MIMO雷達(dá)具有更多的自由度，因此可以通過設(shè)計(jì)發(fā)射信號相關(guān)矩陣來靈活地設(shè)計(jì)滿足實(shí)際需要的發(fā)射方向圖[1216]-。本文主要對集中式MIMO雷達(dá)(以下簡稱為MIMO雷達(dá))的發(fā)射方向圖設(shè)計(jì)進(jìn)行研究。

目前，較為經(jīng)典的發(fā)射方向圖設(shè)計(jì)方法是由文獻(xiàn)[12]針對MIMO雷達(dá)發(fā)射方向圖匹配問題提出的半正定二次規(guī)劃方法(Semidefinite Quadratic Programming, SQP)，該算法根據(jù)給定的發(fā)射方向圖在約束各個(gè)發(fā)射天線發(fā)射功率相同的條件下，以最小二乘為準(zhǔn)則得到全局最優(yōu)的發(fā)射信號相關(guān)矩陣。文獻(xiàn)[15]通過修正已有發(fā)射信號相關(guān)矩陣的非對角線元素來獲得低旁瓣方向圖。文獻(xiàn)[16]則提出了兩種基于凸優(yōu)化的MIMO雷達(dá)發(fā)射信號相關(guān)矩陣設(shè)計(jì)算法，可以有效逼近期望發(fā)射方向圖，其中一種方法使得綜合出的方向圖在全角域范圍內(nèi)與期望方向圖匹配，另一種方法則在保證主瓣不失真的情況下獲得低于閾值的旁瓣電平。雖然這些算法的優(yōu)化問題都是凸優(yōu)化問題，可以直接使用凸優(yōu)化工具包(如CVX[17])進(jìn)行求解，但其計(jì)算復(fù)雜度高，不利于實(shí)時(shí)處理。而且在實(shí)際應(yīng)用中，動(dòng)目標(biāo)的方位、距離和徑向速度時(shí)刻在發(fā)生變化，為提高雷達(dá)資源的利用率，人們希望發(fā)射的電磁能量分配也要自適應(yīng)于它們的變化。因此，發(fā)射方向圖設(shè)計(jì)方法不僅要考慮其逼近期望方向圖的性能，還要具有較高的計(jì)算效率。

針對發(fā)射方向圖設(shè)計(jì)中計(jì)算量較大的問題，本文提出了一種循環(huán)迭代的MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法。該方法基于加權(quán)最小二乘準(zhǔn)則建立方向圖匹配模型，首先將對發(fā)射信號相關(guān)矩陣的求解轉(zhuǎn)化為對其Hermite平方根的求解，再通過引入輔助變量，將關(guān)于Hermite平方根的四次優(yōu)化問題近似地轉(zhuǎn)化為二次優(yōu)化問題，最后以循環(huán)迭代的方式進(jìn)行求解。特別地，對于均勻線陣，當(dāng)采用均勻加權(quán)且離散化方位角在歸一化空間頻率域均勻采樣時(shí)，該方法在迭代過程中還可以利用快速傅里葉變換(FFT)進(jìn)一步提高運(yùn)算效率。由于該算法的主要步驟由FFT完成，運(yùn)算量較小，工程應(yīng)用性較強(qiáng)。實(shí)驗(yàn)結(jié)果表明該算法不但能夠有效逼近期望發(fā)射方向圖，而且具有較高的實(shí)時(shí)性。

2 問題描述

假設(shè)MIMO雷達(dá)系統(tǒng)的發(fā)射陣列是由M個(gè)天線構(gòu)成的均勻線陣，陣元間距為d，如圖1所示。各天線發(fā)射窄帶的相位編碼信號，則第l個(gè)碼元時(shí)刻陣列發(fā)射的信號矢量為

其中，xm(l)為第m個(gè)天線在第l個(gè)碼元時(shí)刻的發(fā)射信號，L為發(fā)射信號的碼元長度或子脈沖數(shù)。不失一般性，以陣列左邊第1個(gè)天線為相位參考點(diǎn)，忽略傳播衰減，則子脈沖x(l)在遠(yuǎn)場θ方位處的基帶信號為

其中，λ為發(fā)射信號波長，a(θ)=[1,ej2πdsinθ/λ,…, ej2π(M-1)dsinθ/λ]T為發(fā)射導(dǎo)向矢量?？沼蛐盘杫(θ,l)在一個(gè)脈沖內(nèi)的平均功率為

圖1 MIMO雷達(dá)陣列示意圖

P(θ)反映了雷達(dá)發(fā)射能量在空間上的分布情況，稱作發(fā)射方向圖。從式(3)可以看出，對于確定的發(fā)射陣列結(jié)構(gòu)，發(fā)射方向圖的形狀主要取決于發(fā)射信號的相關(guān)矩陣R。通過對R的優(yōu)化，可以靈活地對發(fā)射方向圖進(jìn)行設(shè)計(jì)，以滿足不同場景下的需求，提高雷達(dá)系統(tǒng)資源的利用率。目前MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)多采用方向圖匹配模型[12]，即對相關(guān)矩陣R進(jìn)行優(yōu)化使其方向圖P(θ)在一定準(zhǔn)則下逼近期望的發(fā)射方向圖Pd(θ)，通常以加權(quán)最小二乘為準(zhǔn)則建立優(yōu)化模型，可以表示為

3 快速發(fā)射方向圖設(shè)計(jì)方法

3.1 循環(huán)迭代法

為了便于計(jì)算并保證發(fā)射信號相關(guān)矩陣R的半正定性，采用的分解為

在式(6)的約束條件下對矩陣U進(jìn)行優(yōu)化，再由式(5)得到矩陣R，可以保證獲得的R滿足式(4)中的兩個(gè)約束，即半正定約束和等功率約束。因此，將求解R的問題式(4)轉(zhuǎn)化為求解U的優(yōu)化問題：

由于式(7)中目標(biāo)函數(shù)是關(guān)于變量umn的四次函數(shù)，對U進(jìn)行優(yōu)化比較困難，因此參考文獻(xiàn)[9]中的處理方式，引入輔助變量∈?M×1(k=1,2,…,K)，將式(7)近似地轉(zhuǎn)化為關(guān)于umn的二次優(yōu)化模型：

該問題的解為

由式(10)得到K個(gè)sk并組合成矩陣形式，可以表示為S=[s1,s2,…,sK]T。

在矩陣ΛAT的左乘作用下，U的行向量之間相互耦合，故無法像式(9)一樣進(jìn)行求解。因此為獲得，繼續(xù)采用循環(huán)交替的方式對其進(jìn)行求解。令B=ΛAT有

其解為

按照式(15)和式(17)依次計(jì)算iZ和iu，并循環(huán)重復(fù)直至式(13)中的目標(biāo)函數(shù)滿足終止條件。

根據(jù)以上思路，對式(8)的求解可轉(zhuǎn)化為對式(9)、式(11)和式(13)的循環(huán)迭代求解，具體算法流程如表1所示。

表1 發(fā)射方向圖設(shè)計(jì)算法流程

將關(guān)于發(fā)射信號相關(guān)矩陣R的優(yōu)化問題轉(zhuǎn)化為對其Hermite平方根矩陣U的二次優(yōu)化問題，并采用循環(huán)迭代的方式進(jìn)行求解，避免了凸優(yōu)化工具包的使用，以及由于矩陣半正定約束引入的復(fù)雜數(shù)值計(jì)算，因此有利于該發(fā)射方向圖設(shè)計(jì)方法在實(shí)際工程中的應(yīng)用。

3.2 迭代FFT方法

對于陣元間距為半波間距(d=λ/2)的均勻線陣，空域合成的基帶信號式(2)可以表示為

其中，f=dsinθ/λ=sinθ/2為歸一化空間頻率。如果對f在區(qū)間[-0.5,0.5]進(jìn)行Ns點(diǎn)均勻采樣，即

式(18)相當(dāng)于對發(fā)射信號()lx進(jìn)行快速傅里葉逆變換(IFFT),sN就是IFFT點(diǎn)數(shù)。這意味著空域合成信號和發(fā)射信號之間存在離散傅里葉逆變換關(guān)系[18,19]。

當(dāng)方位角按照式(19)進(jìn)行離散化時(shí)，方位角kθ處的導(dǎo)向矢量可以表示為

此時(shí)，導(dǎo)向矢量矩陣A=[a(θ-Ns/2),a(θ-Ns/2+1),…, a(θNs/2-1)]。定義IFFT矩陣：

上一小節(jié)考慮了加權(quán)最小二乘準(zhǔn)則下的發(fā)射方向圖設(shè)計(jì)，由于矩陣ΛAT的左乘作用，矩陣U的行向量umT之間存在耦合，無法有效求解。對于均勻加權(quán)的情況，即Λ=IK，其中IK為K維單位陣，當(dāng)離散方位角按照式(19)進(jìn)行采樣時(shí)，式(13)中的目標(biāo)函數(shù)可以表示為

因此對于均勻加權(quán)且離散方位角為式(19)的情況，發(fā)射方向圖設(shè)計(jì)的算法流程表1中，步驟4可以按照式(23)進(jìn)行更新。在發(fā)射方向圖設(shè)計(jì)算法中，主要計(jì)算代價(jià)為式(10)和式(23)的計(jì)算，由于矩陣F為IFFT矩陣，導(dǎo)向矢量矩陣A為矩陣F的前M行，因此式(10)中的UTa(θk)及式(23)中V的計(jì)算可以采用IFFT和FFT操作，有效提高計(jì)算效率。

4 仿真實(shí)驗(yàn)

下面通過仿真實(shí)驗(yàn)驗(yàn)證本文方法的有效性和可行性。假設(shè)MIMO雷達(dá)系統(tǒng)的發(fā)射陣列為半波間距的均勻線陣，各陣元發(fā)射功率相同且均為1。表1中本文算法的終止閾值分別為ε1=10-2, ε2=10-2。將本文方法與文獻(xiàn)[12]中的半正定二次規(guī)劃方法(SQP)進(jìn)行對比。同時(shí)，采用均方誤差(Mean-Squared Error, MSE)來衡量本文方法對發(fā)射方向圖的優(yōu)化質(zhì)量，其定義為

4.1 非均勻加權(quán)情況

設(shè)方位角采樣間隔為1?，陣元數(shù)16M=，分別給出在采用非均勻加權(quán)時(shí)本文所提循環(huán)迭代方法單次實(shí)驗(yàn)的方向圖綜合結(jié)果，及多次實(shí)驗(yàn)的均方誤差隨迭代次數(shù)的變化曲線。獨(dú)立實(shí)驗(yàn)次數(shù)為100，每次實(shí)驗(yàn)的初始值隨機(jī)給定。

首先考慮期望發(fā)射方向圖為指向0?，波束寬度為60?的寬波束，即

圖2中虛線為期望發(fā)射方向圖，實(shí)線為本文方法所得方向圖，其中權(quán)重wk在區(qū)間[-33?,-29?]∪[29?,33?]為0，在其他方位為1；點(diǎn)線為SQP方法采用非均勻加權(quán)時(shí)進(jìn)行發(fā)射方向圖設(shè)計(jì)的結(jié)果，權(quán)重與本文方法對應(yīng)的權(quán)重相同；點(diǎn)劃線為SQP方法采用均勻加權(quán)所得方向圖，權(quán)重為1。從圖2可以看出，采用非均勻加權(quán)方式可以獲得比均勻加權(quán)更低的旁瓣，且本文方法的性能與SQP方法接近。從圖3可以看出多次實(shí)驗(yàn)下所提循環(huán)迭代方法的均方誤差曲線下降趨勢相同，都隨著迭代次數(shù)的增加逐漸減小并趨于平穩(wěn)。

其次考慮期望發(fā)射方向圖為分別指向-40?, 0?和40?，波束寬度為20?的多波束，即

圖4中實(shí)線為本文方法所得方向圖，權(quán)重在區(qū)間[-54?,-49?]∪[-31?,-27?]∪[-12?,-9?]和其對稱區(qū)間的值為0，在其他方位為1。SQP方法采用非均勻加權(quán)時(shí)的權(quán)重選擇和本文方法相同。由圖4可以看出本文方法在采用非均勻加權(quán)時(shí)能夠獲得比SQP方法更低的旁瓣。圖5為循環(huán)迭代方法的均方誤差隨迭代次數(shù)的變化曲線，說明隨著迭代次數(shù)的增加，所設(shè)計(jì)方向圖逐漸逼近期望方向圖。

圖2 寬波束發(fā)射方向圖

圖3 均方誤差隨迭代次數(shù)的變化曲線

圖4 多波束發(fā)射方向圖

4.2 均勻加權(quán)情況

當(dāng)采用均勻加權(quán)，按照式(19)對方位角進(jìn)行采樣時(shí)，本文方法在循環(huán)迭代過程中可以進(jìn)行FFT/IFFT操作，以提高計(jì)算效率。設(shè)采樣點(diǎn)數(shù)Ns=256。為測試本文方法與SQP方法設(shè)計(jì)發(fā)射方向圖的計(jì)算時(shí)間，本文以2.8 GHz主頻和32 GB內(nèi)存的PC機(jī)作為計(jì)算平臺，采用Matlab語言編寫程序，用凸優(yōu)化工具包CVX[17]求解SQP方法。計(jì)算時(shí)間是由100次Monte Carlo實(shí)驗(yàn)平均所得。

當(dāng)陣元數(shù)較少時(shí)，即16M=時(shí)對比兩種方法所得發(fā)射方向圖，結(jié)果示于圖6。圖6(a)中期望發(fā)射方向圖為指向0?波束寬度為40?的寬波束，與SQP方法所得方向圖相比，采用本文所提迭代FFT方法獲得的方向圖具有更低的旁瓣。圖6(b)中期望發(fā)射方向圖為分別指向40-?, 0?和40?的多波束，其各波束主瓣由傳統(tǒng)相控陣?yán)走_(dá)的主波束(加30 dB切比雪夫權(quán))構(gòu)成，可以看出兩種方法所設(shè)計(jì)方向圖在主瓣區(qū)域與期望方向圖很接近。從表2中給出的計(jì)算時(shí)間對比結(jié)果可以看出，由于采用FFT/IFFT操作，迭代FFT方法具有更高的計(jì)算效率。

為了進(jìn)一步說明迭代FFT方法在陣元數(shù)較大時(shí)的有效性，圖7給出了陣元數(shù)40M=時(shí)的發(fā)射方向圖設(shè)計(jì)結(jié)果，期望發(fā)射方向圖與圖6中的相同。從圖7的對比結(jié)果可以看出，采用迭代FFT方法設(shè)計(jì)的方向圖依然能夠獲得更低的旁瓣，且由表3可以看出，在陣元數(shù)較大情況下迭代FFT方法的求解速度比SQP方法快很多。

表2 陣元數(shù)M=16時(shí)計(jì)算時(shí)間(s)對比

表3 陣元數(shù)M=40時(shí)計(jì)算時(shí)間(s)對比

5 結(jié)束語

圖5 均方誤差隨迭代次數(shù)的變化曲線

圖6 陣元數(shù)M=16時(shí)發(fā)射方向圖對比

圖7 陣元數(shù)M=40時(shí)發(fā)射方向圖對比

本文提出了一種循環(huán)迭代的MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法，該方法以加權(quán)最小二乘為準(zhǔn)則建立方向圖匹配模型，通過引入輔助變量，將對信號相關(guān)矩陣的優(yōu)化問題轉(zhuǎn)化為關(guān)于其Hermite平方根的二次優(yōu)化問題，所得相關(guān)矩陣不但具有半正定性且滿足各陣元發(fā)射功率相同。尤其是對于均勻線陣，當(dāng)采用均勻加權(quán)且離散化方位角在歸一化空間頻率域均勻采樣時(shí)，迭代過程可以采用FFT/IFFT操作，進(jìn)一步提高計(jì)算效率。通過仿真實(shí)驗(yàn)可以看出，本文方法綜合出的發(fā)射方向圖可以有效逼近期望方向圖，與已有方法相比，不但能夠獲得更低的旁瓣，而且計(jì)算時(shí)間大大減少。由于所提方法求解過程無需利用凸優(yōu)化工具包，實(shí)時(shí)性較高，利于工程實(shí)現(xiàn)。

[1] Haimovich A M, Blum R S, and Cimini L J. MIMO radar with widely separated antennas[J]. IEEE Signal Processing Magazine, 2008, 25(1): 116-129.

[2] Li J and Stoica P. MIMO radar with colocated antennas[J]. IEEE Signal Processing Magazine, 2007, 24(5): 106-114.

[3] 羅濤, 關(guān)永峰, 劉宏偉, 等. 平面陣MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法[J]. 電子與信息學(xué)報(bào), 2013, 35(11): 2707-2713.

[4] Aittomaki T and Koivunen V. Low-complexity method for transmit beamforming in MIMO radars[C]. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Honolulu, Hi, USA, 2007, 2: 305-308.

[5] Ahmed S and Alouini M S. MIMO radar transmit beampattern design without synthesising the covariance matrix [J]. IEEE Transactions on Signal Processing, 2014, 62(9): 5326-5337.

[6] 王旭, 糾博, 劉宏偉, 等. 一種基于先驗(yàn)信息的MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)方法[J]. 電子與信息學(xué)報(bào), 2013, 35(12): 2802-2808.

[7] Stoica P, Li J, and Zhu X. Waveform synthesis for diversity-based transmit beampattern design[J]. IEEE Transactions on Signal Processing, 2008, 56(6): 2593-2598.

[8] Stoica P, He H, and Li J. New algorithms for designing unimodular sequences with good correlation properties[J]. IEEE Transactions on Signal Processing, 2009, 57(4): 1415-1425.

[9] He H, Stoica P, and Li J. Designing unimodular sequence sets with good correlations—including an application to MIMO radar[J]. IEEE Transactions on Signal Processing, 2009, 57(11): 4391-4405.

[10] Ahmed S, Thompson J S, Petillot Y R, et al.. Finite alphabet constant-envelop waveform design for MIMO radar[J]. IEEE Transactions on Signal Processing, 2011, 59(11): 2278-2289.

[11] Fuhrmann D R and Antonio G S. Transmit beamforming for MIMO radar systems using signal cross-correlation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(1): 171-186.

[12] Stoica P, Li J and Xie Y. On probing signal design for MIMO radar[J]. IEEE Transactions on Signal Processing, 2007, 55(8): 4151-4161.

[13] Hua G and Abeysekera S. MIMO radar transmit beampattern design with ripple and transition band control[J]. IEEE Transactions on Signal Processing, 2013, 61(11): 2963-2974.

[14] Ahmed S, Thompson J S, Petillot Y R, et al.. Unconstrained synthesis of covariance matrix for MIMO radar transmit beampattern[J]. IEEE Transactions on Signal Processing, 2011, 59(8): 3837-3849.

[15] 羅濤, 關(guān)永峰, 劉宏偉, 等. 低旁瓣MIMO雷達(dá)發(fā)射方向圖設(shè)計(jì)[J]. 電子與信息學(xué)報(bào), 2013, 35(12): 2815-2822.

[16] Gong P, Shao Z, Tu G, et al.. Transmit beampattern design based on convex optimization for MIMO radar systems[J]. Signal Processing, 2014, 94: 195-201.

[17] Grant M and Boyd S. CVX: Matlab software for disciplined convex programming, version 2.0 beta. http://cvxr.com/cvx, 2013, 7.

[18] Keizer W P M N. Fast low-sidelobe synthesis fore large planar array antennas utilizing successive fast Fourier transforms of the array factor[J]. IEEE Transactions on Antennas and Propagation, 2007, 55(3): 715-722.

[19] Keizer W P M N. Low-sidelobe pattern synthesis using iterative Fourier techniques coded in MATLAB[EM Programmer’s Notebook][J]. IEEE Antennas and Propagation Magazine, 2009, 51(2): 137-150.

吳 夢： 女，1985年生，博士生，研究方向?yàn)镸IMO雷達(dá)、自適應(yīng)信號處理.

劉宏偉： 男，1971年生，博士，教授，博士生導(dǎo)師，研究方向?yàn)槔走_(dá)系統(tǒng)、雷達(dá)信號處理、雷達(dá)自動(dòng)目標(biāo)識別等.

王 旭： 男，1987年生，博士生，研究方向?yàn)镸IMO雷達(dá)、認(rèn)知雷達(dá)發(fā)射波形設(shè)計(jì).

A Cyclic Iterative Method for MIMO Radar Transmit Beampattern Design

Wu Meng Liu Hong-wei Wang Xu
(National Laboratory of Radar Signal Processing, Xidian University, Xi’an 710071, China)

Transmit beampattern design for Multiple Input Multiple Output (MIMO) radar is obtained by synthesizing the signal covariance matrix, which can be achieved by convex optimization approaches. Due to the high computational complexity, these approaches are not easy for practical implementation. In this paper, a cyclic iterative method for MIMO radar transmit beampattern design is proposed. Based on the weighted least square criterion, signal covariance matrix can be obtained by optimizing the quadratic cost function with respect to its Hermite square root. For the uniform linear array, especially when the sampling grid in normalized spatial frequency is uniform and the weights for grid points are the same, FFT can be used to further reduce the computational complexity of the proposed algorithm. Simulation results show that the resulting beampattern matches the desired pattern closely and the proposed algorithm is efficient for the real time application.

MIMO radar; Transmit beampattern design; Cyclic iterative; FFT

TN958

A

1009-5896(2015)02-0322-06

10.11999/JEIT141043

2014-08-04收到，2014-12-01改回

國家自然科學(xué)基金(61401329, 61201285)資助課題

*通信作者：吳夢 mengwu121@163.com