倪智宇，劉金國，吳志剛

(1. 中國科學(xué)院沈陽自動化研究所機器人學(xué)國家重點實驗室，沈陽 110016； 2. 大連理工大學(xué)航空航天學(xué)院，大連 116024)

0 引 言

航天器系統(tǒng)模態(tài)參數(shù)的辨識實驗?zāi)壳耙呀?jīng)在一些在軌空間結(jié)構(gòu)上開展并獲得成功，例如伽利略探測器、哈勃空間望遠鏡以及工程實驗衛(wèi)星(Engineering Test Satellite，ETS)系列等[1-3]，但這些在軌辨識實驗均是基于定常系統(tǒng)模型進行的，比如ETS-VI衛(wèi)星就在實驗中將太陽能帆板固定在某一個角度，以防止它的旋轉(zhuǎn)對整星模態(tài)參數(shù)辨識造成影響。從實際在軌運行的角度考慮，系統(tǒng)模態(tài)參數(shù)可能由于航天器的對接與分離、太陽能帆板的轉(zhuǎn)動、大型天線的展開或者對目標(biāo)物體的捕獲等情況發(fā)生改變[4]。例如美國的土壤濕度主/被動觀測衛(wèi)星，攜帶的大型可展開式天線在工作時會始終以4 s/r的速度進行旋轉(zhuǎn)[5]，而日本ETS-VIII衛(wèi)星帆板的轉(zhuǎn)動則可以導(dǎo)致結(jié)構(gòu)參數(shù)的變化幅度高達25%[6]?？梢娺@種由于結(jié)構(gòu)構(gòu)型改變而對系統(tǒng)模態(tài)參數(shù)造成的時變影響是不能夠輕易忽略的，因此有必要考慮航天器運行過程中的時變模態(tài)參數(shù)辨識問題。

針對航天器的模態(tài)參數(shù)辨識問題，Juang等人[3]采用特征系統(tǒng)實現(xiàn)算法(Eigensystem Realization Algorithm，ERA)，使用控制力矩輸入信號和陀螺儀量測信號的數(shù)據(jù)，對哈勃空間望遠鏡進行了結(jié)構(gòu)頻率和阻尼等模態(tài)參數(shù)的在軌辨識。黎康等人[7]僅利用被噪聲污染的輸入信號，基于子空間算法對撓性航天器的模態(tài)參數(shù)進行了辮識。Kasai等人[8]利用航天器姿態(tài)調(diào)整過程中的星體姿態(tài)信號和撓性附件振動響應(yīng)信號，在ETS-VIII衛(wèi)星上成功進行了系統(tǒng)模態(tài)參數(shù)的在軌辨識實驗。文獻[9]提出了一種基于視覺測量的太陽翼模態(tài)參數(shù)辨識方法，通過相機測量得到測點處的振動位移響應(yīng)，然后采用ERA算法辨識太陽翼的模態(tài)參數(shù)。在現(xiàn)有的航天器模態(tài)參數(shù)辨識方法中，ERA和子空間方法是兩種最為常用的時域辨識方法[10]。這兩種方法的核心思想是通過構(gòu)建相應(yīng)的Hankel矩陣，利用奇異值分解(Singular Value Decomposition，SVD)計算得到系統(tǒng)的狀態(tài)空間模型和對應(yīng)的模態(tài)參數(shù)。這兩種方法已經(jīng)被工程人員證明是成熟可靠的辨識手段，但該類方法主要針對的是定常系統(tǒng)，對于時變系統(tǒng)的辨識能力極為有限。此外，對于航天器系統(tǒng)的控制問題而言，一些控制方法需要實時地獲得系統(tǒng)的時變模態(tài)參數(shù)，以便于及時實現(xiàn)控制器參數(shù)的調(diào)整。因此有關(guān)航天器時變模態(tài)參數(shù)的辨識問題，尋找一種計算效率較快、更加適合于在線辨識的方法是有必要的[11]。

為克服基于SVD的子空間辨識方法計算效率較低的問題，Houtzager等人[12]提出了一種名叫基于遞推預(yù)測器的子空間辨識(Recursive Predictor-Based Subspace Identification，RPBSID)方法。這種遞推方法通過自適應(yīng)濾波，利用外生向量自回歸(Vector Autoregressive with Exogenous，VARX)預(yù)測器來提供漸進一致估計，從而通過求解最小二乘線性問題計算狀態(tài)空間模型，進而辨識得到系統(tǒng)的時變模態(tài)參數(shù)。這類遞推方法已經(jīng)在一些機械和航空結(jié)構(gòu)的辨識實驗中得到驗證，但是還尚未見到其在航天器參數(shù)辨識方面的應(yīng)用研究。

本論文針對RPBSID方法在估計系統(tǒng)狀態(tài)量時需要構(gòu)建廣義Hankel矩陣而導(dǎo)致數(shù)據(jù)量較大的問題，提出了一種改進后的RPBSID方法，并將其用于航天器時變模態(tài)參數(shù)的辨識中。這種改進方法通過引入仿射投影方法(Affine Projection Algorithm，APA)對狀態(tài)量進行遞推求解，不需要逐個時刻構(gòu)建高維數(shù)的廣義Hankel矩陣和進行矩陣求逆計算，從而提高了方法的計算效率。在數(shù)值算例中，建立帶有大型撓性附件的衛(wèi)星動力學(xué)模型，分別對線性變化、突然改變和周期變化三種情況下的系統(tǒng)模態(tài)頻率進行辨識，辨識結(jié)果證明了這種遞推方法能夠有效地辨識得到航天器的時變模態(tài)頻率值。通過與傳統(tǒng)RPBSID方法計算結(jié)果的比較，證明本文提出的這種改進方法不僅能夠保證計算的精度，而且具有更高的計算效率。

1 剛-撓耦合航天器的建模描述

對于常見的攜帶有太陽能帆板或者可展開式天線等附件的航天器結(jié)構(gòu)，通常將其簡化為中心剛體帶有N個撓性附件的剛-撓耦合模型，因此這類航天器的動力學(xué)方程可以描述為[13]：

(1)

(2)

航天器在軌運行中，如果航天器的總質(zhì)量或者結(jié)構(gòu)構(gòu)型發(fā)生變化，則會導(dǎo)致式(1)和(2)中的轉(zhuǎn)動慣量矩陣J和耦合影響矩陣F是時變的。在這種情況下，我們首先定義一個狀態(tài)矢量δ為：

(3)

那么，如果矩陣J和F是時變的，則式(1)和(2)可以表述為：

(4)

(5)

式中:I為單位陣。選擇航天器的姿態(tài)信號ψ和附件的振動位移信號s作為m×1維的輸出y(t)：

(6)

那么量測方程可以寫為：

(7)

式中：

而C為m×n維的輸出矩陣，Φ為相應(yīng)的模態(tài)矩陣。

2 改進后的RPBSID方法辨識時變模態(tài)參數(shù)

2.1 Markov參數(shù)的遞推估計

將式(5)和(7)進行離散化并考慮噪聲對系統(tǒng)的影響，可以進一步將其寫為如下的新息形式：

xk+1=Akxk+Bkuk+Kkek

(8)

yk=Ckxk+ek

(9)

式中：下標(biāo)k為離散后的時間點，ek表示新息白噪聲序列，而矩陣Kk表示卡爾曼增益矩陣。那么對式(8)和(9)做進一步的變換，可以寫為：

(10)

yk=Ckxk+ek

(11)

定義一個如下形式的VARX預(yù)測器，則k時刻的輸出yk可以表述為：

(12)

yk=Ξkφk

(13)

式中：

那么利用自適應(yīng)濾波技術(shù)，可以得到各個采樣時刻的矩陣Ξk的最小二乘遞推形式為[12]：

(14)

(15)

式中：β1為遺忘因子，需要保證0<β1≤1；Zk為定義的臨時變量矩陣，通常初值給為Z0=I。通過式(14)的遞推計算，可以得到各個時刻的Markov參數(shù)矩陣Ξk的值，接下來將利用該參數(shù)矩陣Ξk遞推估計系統(tǒng)的狀態(tài)量xk。

2.2 基于仿射投影算法的系統(tǒng)狀態(tài)量估計

在傳統(tǒng)子空間方法中，狀態(tài)量xk的計算通?；赟VD或QR分解；而原RPBSID方法在計算狀態(tài)量時則需要分別構(gòu)建各時刻的廣義Hankel矩陣并對其中的參數(shù)矩陣進行求逆處理。上述方法往往均需要較大的計算量，為了減少辨識過程中的計算復(fù)雜度，本文這里擬采用遞推形式來估計各個采樣時刻的狀態(tài)量xk。

對于式(10)-(11)，如果暫不考慮量測噪聲ek，那么從k-p到k時刻的輸出信號{yk-p,yk-p+1,…,yk}可以分別寫為：

yk-p=Ck-pxk-p

(16)

(17)

(18)

則式(16)-(18)可以用矩陣形式表示為：

(19)

式中：

τk-p=Γk-pxk-p

(20)

那么對于k-1時刻，式(20)可以表示為如下的最小二乘形式：

τk-1=Γk-1xk-1

(21)

其中

(22)

對于式(21)，這里引入APA算法[15]，從而得到各個時刻狀態(tài)量xk的遞推形式為：

(23)

Φk-1=τk-1-Γk-1xk-1

(24)

式中：μ為遞推過程中的仿射投影因子，一般應(yīng)滿足0<μ?1。在式(23)的遞推迭代中，矩陣Γk-1的初值Γ0應(yīng)保證為滿秩矩陣，如果沒有其它先驗知識，則Γ0一般可以選擇為：

(25)

在式(23)中，在遞推得到xk后，矩陣Γk的更新同樣可以通過如下的遞推最小二乘算法得到：

(26)

(27)

Γk=Γk-1+(τk-Γk-1xk)Wk

(28)

式中：β2為遺忘因子，而矩陣Lk的初值選取方式和Γ0相似。根據(jù)式(23)-(28)，可以遞推計算得到各個時刻的狀態(tài)量xk。

2.3 時變狀態(tài)空間模型參數(shù)的遞推估計

在得到系統(tǒng)的狀態(tài)量xk后，則可以遞推計算各個時刻的狀態(tài)空間模型參數(shù)[Ak,Bk,Ck]。首先將式(8)和(9)重寫為如下形式：

(29)

(30)

(31)

(32)

(33)

(34)

式中：矩陣Δk的初值Δ0的選取需要滿足Δ0=(1/α1)I,α1>0。而式(30)中的新息噪音序列ek可以通過下式得到：

(35)

(36)

(37)

辨識得到系統(tǒng)矩陣Ak后，可以利用偽模態(tài)分析理論計算得到航天器的時變模態(tài)參數(shù)。對系統(tǒng)矩陣進行特征根分解計算：

(38)

Φk=CkΛk

(39)

2.4 改進的RPBSID方法與原算法計算復(fù)雜度比較

改進后的RPBSID方法的計算步驟可以簡要地總結(jié)如下：

步驟1.通過式(14)和(15)，求解系統(tǒng)的Markov參數(shù)矩陣Ξk的值；

步驟2.根據(jù)式(22)-(24)以及(26)-(28)，計算得到狀態(tài)量xk；

步驟4.利用式(38)和(39)計算得到航天器的時變模態(tài)參數(shù)。

本文對于原方法的改進主要集中在對狀態(tài)量xk的遞推估計，即步驟2。在遞推過程中，由于參數(shù)p、r和m通常是選定好的，因此總的計算復(fù)雜度主要由系統(tǒng)的未知階次數(shù)目n來決定，通過算法復(fù)雜度的估計可知，系統(tǒng)的最高計算復(fù)雜度為O(n2)。而原方法由于需要計算求逆估計狀態(tài)量xk，因此其計算復(fù)雜度可以達到O(n3)。原方法計算狀態(tài)量的內(nèi)容請參見文獻[12]，這里限于篇幅不再詳述。此外，由于改進的遞推方法不需要構(gòu)造高維數(shù)的Hankel矩陣，從而減少了所需的數(shù)據(jù)量，這種計算優(yōu)勢在系統(tǒng)階次較高時將更加明顯，接下來將通過數(shù)值計算來證明。

3 仿真校驗

在數(shù)值仿真中，選擇帶有大型撓性附件的ETS-VIII衛(wèi)星作為建模對象，考慮當(dāng)系統(tǒng)的模態(tài)參數(shù)值分別發(fā)生線性改變、突變和周期變化的情況，利用改進的遞推方法辨識系統(tǒng)的時變模態(tài)頻率值。

3.1 ETS-VIII衛(wèi)星模型描述

ETS-VIII衛(wèi)星攜帶一對大型可展開式天線反射器和一對太陽能帆板，是一個典型的剛-撓耦合航天器結(jié)構(gòu)。對該衛(wèi)星模型做如下的簡化：將所攜帶的大型天線反射器和衛(wèi)星主體部分分別簡化為平面桁架和長方體結(jié)構(gòu)，太陽能帆板和天線反射器通過連桿鉸接于衛(wèi)星中心體上，并且認(rèn)為是由均質(zhì)材料所組成。由于衛(wèi)星在地球同步軌道上運行，因此這里忽略重力梯度的影響。全局坐標(biāo)系的原點建立在整星的質(zhì)心上，而附件坐標(biāo)系的原點則建立在中心體與附件的鉸接點上。簡化后的衛(wèi)星模型如圖1所示，其中符號{s1,s2,a1,a2}分別表示該衛(wèi)星的兩副帆板s1/s2和兩副天線反射器a1/a2。

圖1 簡化后的ETS-VIII模型Fig.1 Simplified ETS-VIII model

在仿真中，太陽能帆板s1/s2和天線a1/a2的頻率可以通過有限元建模得到，而附件的阻尼比ζs1=ζs2=ζa1=ζa2=0.01。在建立動力學(xué)方程的過程中，選取每個附件的前3階頻率，因此該衛(wèi)星模型的系統(tǒng)階次為n=(3+4×3)×2=30。

3.2 時變模態(tài)頻率的辨識

仿真條件中給定噪聲的信噪比(Signal Noise Ratio，SNR)為40 dB，而辨識方法中的仿真參數(shù)如下：分塊矩陣的參數(shù)p=15，遞推系數(shù)α1,2=0.9、μ=0.05，遺忘因子β1,2,3,4=0.98，系統(tǒng)的采樣時間Δt=0.1 s。分別考慮以下三種工況的時變頻率辨識問題：

工況1：模態(tài)參數(shù)線性變化的情況

在該工況中，假設(shè)衛(wèi)星的推進劑消耗導(dǎo)致結(jié)構(gòu)整體質(zhì)量發(fā)生線性改變的情況。當(dāng)然，實際工作中，不同于運載火箭在發(fā)射階段燃料迅速減少的情況，在衛(wèi)星在軌過程中，推進劑的消耗是比較小而且非連續(xù)的。但在本算例中，為直觀地表現(xiàn)出這種變化，假設(shè)在運行的過程中，衛(wèi)星所攜帶的推進劑是較為快速且連續(xù)性的隨著時間而不斷消耗，基于這種假設(shè)，衛(wèi)星中心剛體的質(zhì)量(定義為mr)有著如下的線性變化關(guān)系：

mr2=mr1-ρtp

(40)

式中：mr1和mr2分別為衛(wèi)星中心剛體在推進劑消耗前和消耗后的質(zhì)量，而ρ為推進劑的消耗率，tp則是運行的時間。在仿真中給定mr1=2400 kg。給定推進劑的消耗率ρ=0.2 kg/s。由于系統(tǒng)質(zhì)量的不斷減少，從而導(dǎo)致某些階次的系統(tǒng)頻率發(fā)生時變變化。利用本文提出的改進RPBSID方法對這種工況下的結(jié)構(gòu)模態(tài)參數(shù)進行辨識。從動力學(xué)方程(4)中可以看出，由于該方程包含有衛(wèi)星的姿態(tài)運動方程(1)和各個附件的振動方程(2)，而式(1)中三個方向的姿態(tài)運動對應(yīng)的模態(tài)參數(shù)值為零，所以在動力學(xué)方程(1)中，衛(wèi)星的前三階頻率表征的是系統(tǒng)的剛體模態(tài)從而對應(yīng)的值為零。此外，在動力學(xué)方程(1)-(4)中，由于轉(zhuǎn)動慣量矩陣J(t)中的質(zhì)量參數(shù)mr發(fā)生了變化(具體的轉(zhuǎn)動慣量計算表達式請參見文獻[6])，那么從質(zhì)量和矢徑矩陣的乘積中可以看出，轉(zhuǎn)動慣量矩陣中只有部分元素發(fā)生了變化，從而導(dǎo)致只有部分階次的模態(tài)參數(shù)隨著質(zhì)量變化而發(fā)生改變。對于本算例建立的衛(wèi)星動力學(xué)模型，在前15階系統(tǒng)頻率中，只有第5、7、9和11階的模態(tài)參數(shù)隨著質(zhì)量變化而改變，而其它階次的頻率、阻尼和振型參數(shù)仍為常值。

利用式(38)-(39)，計算得到系統(tǒng)的時變頻率、阻尼比和模態(tài)振型。原RPBSID方法和改進后的遞推方法對于這四階時變頻率的辨識結(jié)果如圖2所示。值得注意的是，由于這類遞推方法需要一個迭代逼近的過程，所以初始時刻的辨識值是有較大誤差的，而圖2只給出了系統(tǒng)在基本達到穩(wěn)定跟蹤后的辨識結(jié)果。表1則給出了這兩種方法對于時變頻率值和阻尼比的平均相對誤差結(jié)果，其中RPBSID和改進后的遞推方法均需要一個迭代過程來逼近模型的原始值，因此計算的是50 s-120 s區(qū)間內(nèi)已達到穩(wěn)定跟蹤時的頻率和阻尼比的平均相對誤差值。計算結(jié)果證明這兩種方法均能夠有效地辨識線性時變系統(tǒng)的頻率值，對于該算例，模態(tài)頻率的最大誤差不超過3%，但阻尼比的辨識結(jié)果則存在著較大的誤差。這是由于阻尼分析的機理比較復(fù)雜，而且阻尼比還容易受到頻率或者剛度的影響，特別是當(dāng)含有噪聲干擾時，噪聲對于阻尼比辨識結(jié)果的影響要遠大于對頻率的影響，因此阻尼比往往難于準(zhǔn)確辨識得到。

圖2 兩種方法辨識得到的系統(tǒng)時變頻率值Fig.2 Identification results of system time-varying frequencies with two methods

表1 兩種方法所得到的頻率和阻尼比的平均相對誤差Table 1 Average relative errors of frequencies and damping ratios for two methods

對于模態(tài)振型的辨識結(jié)果，這里通過引入模態(tài)置信準(zhǔn)則(Modal Assurance Criteria，MAC)，來判斷和比較原振型向量與辨識得到的振型之間的相似度。模態(tài)振型向量的MAC計算公式如下：

(41)

表2 兩種方法所得到的時變模態(tài)振型的MAC值(t=80 s)Table 2 MAC values of mode shape for two methods (t=80 s)

工況2：模態(tài)參數(shù)突然變化的情況

在工況2中，模擬空間對接的情況，考慮衛(wèi)星的中心剛體在Z軸(偏航軸)的兩端分別與兩個長方體結(jié)構(gòu)對接的情況，如圖3所示。在本仿真中，我們假設(shè)與其分別對接的兩個長方體結(jié)構(gòu)的尺寸均為2.35 m×2.45 m×3 m，質(zhì)量均為600 kg，而衛(wèi)星中心體的質(zhì)量隨時間做如下的變化：

(42)

由于這種對接會導(dǎo)致衛(wèi)星中心剛體的尺寸和質(zhì)量發(fā)生改變，從而影響到系統(tǒng)的轉(zhuǎn)動慣量值，導(dǎo)致衛(wèi)星結(jié)構(gòu)的某些階次的模態(tài)參數(shù)值發(fā)生突變。假設(shè)該對接過程并不會對衛(wèi)星的整體質(zhì)心位置造成影響，當(dāng)然這只是一種比較理想的情況，在實際的空間對接情況中，由于對接的結(jié)構(gòu)分布的不同，航天器的質(zhì)心將發(fā)生改變，從而導(dǎo)致航天器動力學(xué)方程中的參數(shù)將更加復(fù)雜，因此這里質(zhì)心不變的情況只是為了簡化建模步驟和后續(xù)討論而作的假設(shè)。那么原RPBSID和改進后的遞推方法對于發(fā)生突變的第5、7、9和11階頻率值的辨識結(jié)果如圖4所示。

圖4 兩種方法辨識得到的系統(tǒng)時變頻率值Fig.4 Identification results of system time-varying frequencies with two methods

從圖4的結(jié)果中可以看出，與文獻[16]中對于參數(shù)突變系統(tǒng)的辨識結(jié)論相一致，基于預(yù)測器的這兩類辨識方法在頻率值發(fā)生突變的時間點處需要一個較長的迭代過程來使得遞推結(jié)果趨于穩(wěn)定，但是對于突變后的穩(wěn)定跟蹤階段還是能夠有效地跟蹤的。如果不考慮兩個突變階段的較大誤差，兩種方法所得到的各階時變頻率和阻尼比的平均相對誤差值如表3所示，結(jié)果證明這種改進的遞推方法可以較為準(zhǔn)確的辨識突變系統(tǒng)的頻率值，但辨識的阻尼比結(jié)果與原始值相比仍然存在著較大的誤差。

表3 兩種方法所得到的頻率和阻尼比的平均相對誤差Table 3 Average relative errors of frequencies and damping ratios for two methods

工況3：模態(tài)參數(shù)周期變化的情況

圖5 兩種方法對于系統(tǒng)時變頻率的辨識結(jié)果Fig.5 Identification results of system time-varying frequencies for two methods

表4 兩種方法對于第7和第9階時變頻率值的平均相對誤差Table 4 Average relative errors of the 7th and 9th time-varying frequencies for two methods

在算例的最后，我們通過MATLAB軟件的仿真結(jié)果，對比這兩種方法的計算時間，以驗證2.4節(jié)計算復(fù)雜度的理論結(jié)果。為了更方便的比對計算結(jié)果，我們分別統(tǒng)計這兩種遞推方法在系統(tǒng)模型階次n分別為24、30、36和48時的計算時間(每組進行30次蒙特卡羅實驗)，則這兩種方法的仿真計算效率結(jié)果如圖7所示。從圖7中可以看出，改進的遞推方法與原方法相比，在應(yīng)用于系統(tǒng)階次較高的模型辨識時其計算優(yōu)勢更加明顯，這是和上文關(guān)于該算法計算復(fù)雜度的理論分析結(jié)果相符合的。由于實際航天器對象的系統(tǒng)模型階次往往選取得比較高，并考慮到在線辨識與控制參數(shù)實時修正的需求，這種計算效率的改進將具有重要的意義。

圖6 兩種方法對于系統(tǒng)時變頻率的辨識結(jié)果Fig.6 Identification results of system time-varying frequencies for two methods

表5 兩種方法對于第7和第9階時變頻率值的平均相對誤差Table 5 Average relative errors of the 7th and 9th time-varying frequencies for two methods

圖7 兩種方法所用的平均計算時間Fig.7 Average computation time of two methods

4 結(jié) 論

本文提出了一種改進的RPBSID方法并應(yīng)用于航天器模態(tài)參數(shù)的辨識。不同于傳統(tǒng)的子空間或RPBSID等參數(shù)辨識方法，這種改進方法利用仿射投影和最小二乘技術(shù)來實現(xiàn)狀態(tài)量的遞推辨識，從而避免了傳統(tǒng)方法在計算系統(tǒng)狀態(tài)量的過程中，由于奇異值分解或者構(gòu)建廣義Hankel矩陣所導(dǎo)致的數(shù)據(jù)量較大的問題，提高了方法的計算效率。仿真結(jié)果證明了該方法能夠有效得到系統(tǒng)的時變模態(tài)參數(shù)，與傳統(tǒng)的子空間方法相比具有更高的計算效率。

由于這種方法是基于RPBSID算法而提出，因此還存在著一些未來需要解決的問題：(1)該方法需要一定的先驗知識，例如參數(shù)p的選取往往需要提前知道系統(tǒng)的階次n。而且這類方法需要一定的迭代次數(shù)來逐步逼近于實際值，例如圖2中初始時間所示的情況，從而達到參數(shù)值的穩(wěn)定跟蹤；(2)這種方法主要適用于慢時變系統(tǒng)，例如工況1的情況；而對于變化較快的系統(tǒng)，這種遞推方法的跟蹤效果將會顯著變差。三種工況的驗證結(jié)果也能反映出這種方法需要一定的反應(yīng)時間來達到參數(shù)跟蹤的目的，因此對于慢變系統(tǒng)的跟蹤效果較好，通過RPBSID方法的原始作者所提供的MATLAB源程序算例和對應(yīng)的文獻[12,16]所給出的結(jié)論也證明了這一點。因此，如何提高這種方法對于快變系統(tǒng)的辨識精度，是未來需要繼續(xù)探索和深入研究的問題。