国产日韩欧美一区二区三区三州_亚洲少妇熟女av_久久久久亚洲av国产精品_波多野结衣网站一区二区_亚洲欧美色片在线91_国产亚洲精品精品国产优播av_日本一区二区三区波多野结衣 _久久国产av不卡

?

旋翼振動監(jiān)測信號的演化分析與損傷跟蹤方法*

2015-06-10 00:34:46朱旭程
振動、測試與診斷 2015年1期
關鍵詞:相空間槳葉維數

朱旭程

(海軍航空工程學院飛行器工程系 煙臺,264001)

?

旋翼振動監(jiān)測信號的演化分析與損傷跟蹤方法*

朱旭程

(海軍航空工程學院飛行器工程系 煙臺,264001)

為獲得一種對損傷敏感而對正常狀態(tài)擾動不敏感的旋翼損傷觀測信號、解決直升機旋翼在飛行狀態(tài)下的損傷監(jiān)測問題,研究了旋翼在重構相空間中的損傷演化特性并提出一種新的損傷跟蹤方法。該方法首先應用嵌入技術將旋翼氣彈有限元模型產生的振動監(jiān)測信號重構到維數更高的相空間中,采用非線性Volterra級數建立旋翼基準狀態(tài)預測模型,以基準狀態(tài)預測結果與實測相軌跡之件的差異作為狀態(tài)預測殘差,在多個鄰域內對其統(tǒng)計平均形成一個損傷觀測特征向量;然后,應用奇異譜分解方法從損傷觀測特征向量的時間序列中提取出槳葉損傷演化的維度和趨勢信息,采用雙指數平滑方法建立損傷演化趨勢預測模型并估計出槳葉損傷故障剩余壽命。采用槳葉損傷模型和旋翼氣彈模型仿真數據驗證了所提方法的可行性和有效性。結果表明:該方法能充分利用監(jiān)測信號的非線性特性在高維相空間中重建系統(tǒng)動力學本質,以不同的時間尺度來觀測系統(tǒng)的演化特性;具有損傷模式自動識別能力,可用于難以事先確定系統(tǒng)損傷演化模型或維數信息的場合。

直升機;旋翼;健康監(jiān)測;故障診斷;相空間

引言

旋翼是直升機飛行的關鍵動力部件,在槳葉周期變距和氣動交變載荷作用下槳葉、槳轂、減擺器和變距拉桿等部件會出現磨損和疲勞失效故障。能夠準確監(jiān)測旋翼的損傷狀況,對于提高直升機飛行安全性極具現實意義。文獻[1]對旋翼損傷監(jiān)測問題開展了研究,提出了通過槳葉相對位移和機體振動識別旋翼故障的方案。Ganguli等[2-3]分析了典型損傷對旋翼氣彈響應信號的影響,運用神經網絡、模糊推理等人工智能技術建立了一些旋翼故障診斷算法。目前,旋翼在飛行條件下的損傷監(jiān)測技術還存在瓶頸[4]:在高速旋轉的旋翼上加裝各種數據采集系統(tǒng)尚存在一些技術上的限制;此外,由于槳距時變、氣彈耦合和非定常尾流等因素影響,傳統(tǒng)損傷跟蹤方法難以獲得一種對損傷敏感而對正常狀態(tài)擾動不敏感的損傷觀測信號,有效的旋翼損傷跟蹤方法依然缺乏。

筆者應用混沌信號嵌入分析理論的相空間重構技術,將旋翼狀態(tài)監(jiān)測數據變換到更高維數的狀態(tài)空間內、以不同的時間尺度來研究旋翼系統(tǒng)的損傷演化特性。采用非線性Volterra級數建立旋翼基準狀態(tài)預測模型,用基準狀態(tài)預測結果與實測相軌跡之間的差異作為狀態(tài)預測殘差,發(fā)展一種新的旋翼損傷跟蹤方法。

1 旋翼在重構相空間的損傷演化特性

1.1 旋翼氣彈有限元模型

為研究損傷對槳葉位移、槳轂力等旋翼狀態(tài)監(jiān)測信號的影響,需要建立旋翼氣彈有限元模型。由Hamilton能量變分原理

(1)

其中:Π為系統(tǒng)總能;T為動能;U為應變能;W為外力虛功。

氣動載荷采用葉素理論建立,旋翼誘導速度采用White-Blake入流模型[5]計算。

采用有限元方法離散上述能量積分,進一步采用分部積分法建立旋翼氣彈有限元方程

(2)

其中:M,C,K,F分別表示慣性矩陣、阻尼矩陣、剛度矩陣和載荷向量,在定常平飛狀態(tài)下它們?yōu)榉轿唤堑闹芷诤瘮怠?/p>

旋翼氣彈有限元方程可采用時域有限元方法化為代數方程求解,或者將其轉換為非線性狀態(tài)方程采用數值積分方法求解。

旋翼的故障模式種類多樣,不同類型的故障對系統(tǒng)參數會有不同影響。采用復合材料槳葉剛度疲勞損傷演化模型說明本研究方法的有效性。

(3)

槳葉損傷變量隨使用壽命的擴展如圖1所示。

圖1 槳葉剛度損傷變量擴展曲線Fig.1 Damage accumulation curve of blade stiffness

目前,直升機健康和使用監(jiān)測系統(tǒng)(health and usage monitor system,簡稱HUMS)能夠提供的旋翼狀態(tài)監(jiān)測信號主要有二類:a. 槳葉的揮舞、擺動和扭轉位移,可通過在機身上安裝激光多普勒測振儀實時監(jiān)測[6];b. 槳轂力和力矩,可通過機身上的應變計或振動傳感器結合飛行參數得到。下面以槳尖位移、槳轂力和力矩作為旋翼典型狀態(tài)監(jiān)測信號。

將旋翼氣彈有限元模型和槳葉損傷演化模型結合,分析槳葉損傷對槳尖位移和槳轂載荷的影響。筆者以BO-105直升機旋翼為參考設置實例旋翼的基本參數:槳葉數為4,旋翼半徑R為4.94 m,旋翼轉速Ω為40 rad/s,槳葉密度m為6.46 kg/m,槳葉弦長為0.055R,弦向剛度為0.026 8mΩ2R4,揮向剛度為0.010 8mΩ2R4,扭轉剛度為0.006 15mΩ2R4,洛克數為5.2,槳葉升力線斜率為6.0,型阻系數為0.002,誘導阻力系數為0.2,俯仰力矩系數為0.00。飛機在前進比0.3,旋翼拉力系數0.004 9的條件下定常平飛。旋翼氣彈模型計算程序采用Matlab軟件建立(如圖2),槳葉在健康(t/tF=0)和損傷(t/tF=1/3)狀態(tài)下仿真產生的槳尖位移、槳根力和力矩如圖3所示。圖中:w為槳尖揮舞位移;v為槳尖擺動位移;Fz為槳轂力;Mx為槳轂力矩;n為樣本長度。

圖2 旋翼氣彈模型計算流程Fig.2 Flowchart of rotor simulation procedure

圖3 健康和損傷旋翼槳尖位移和槳轂載荷Fig.3 Health and damaged blade tip deflections and Hub loads

1.2 旋翼狀態(tài)監(jiān)測信號的相空間重構

旋翼損傷監(jiān)測的一個基本問題是找到一種對損傷敏感而對系統(tǒng)狀態(tài)變化、干擾噪聲不敏感的觀測信號。然而槳葉損傷與氣彈耦合現象會引起氣動載荷沿槳葉結構的分布發(fā)生變化[7],隨著槳葉損傷程度的擴展,在揮舞、擺動、扭轉方向的一些固有頻率互相穿越與交叉,使得槳尖位移和槳轂載荷的峰值等特征信號在單調性關系上出現突變,因而難以由這些特征信號辨識出損傷的變化趨勢。源自混沌理論的相空間重構技術可以在不同維數的空間內、以不同的時間尺度來研究演化系統(tǒng)的動力學特性,為非線性振動系統(tǒng)的損傷跟蹤提供了一條新途徑。只要時間尺度和空間維數選擇合適,可通過相軌跡的移動變化定量跟蹤系統(tǒng)損傷的演化過程[8]。將某監(jiān)測信號的時序數據{y(i):i=1,2,…,M} 進行時間坐標延遲,構造一狀態(tài)點列 {x(n):n=1,2,…,N}。其中:x(n) = [y1(n),y2(n), …,ym(n)]T;yk(n) =y(n+(k-1)τ);m為嵌入維數;τ為延遲時間。該點列描述了系統(tǒng)在m維狀態(tài)空間中的相軌跡,由它構成的空間稱為相空間,這個過程稱為相空間重構。根據嵌入理論[9],如果重構系統(tǒng)與原動力學系統(tǒng)的動態(tài)響應可限制在相同維數的吸引子流形上,重構系統(tǒng)與原系統(tǒng)在拓撲關系或微分意義上屬于同一等價(微分同胚),具有相同的動力學性質和幾何性質,系統(tǒng)的許多特性如分維、K熵、李指數等均相同。

圖4 平均互信息法和假近鄰法求重構參數流程Fig.4 Flowchart for time delay and embedding dimension based on average mutual information and FNN percentage

圖5 相軌跡在三維空間中的投影Fig.5 3D views of reconstructed phase space trajectories

為直觀表明旋翼在相空間中的損傷演化過程,將帶槳葉損傷BO-105旋翼的槳尖位移和槳轂載荷數據進行相空間重構。在相空間重構時,延遲時間和嵌入維數是兩個主要參數,采用平均互信息法和假近鄰法來確定[10],計算流程如圖4所示。求得延遲時間為30,嵌入維數為4,信號相軌跡在三維空間中的投影如圖5所示??梢姡硐到y(tǒng)的穩(wěn)態(tài)相軌跡(吸引子)在重構空間中完全展開,通過相軌跡隨槳葉損傷擴展而出現的漂移和變形等幾何變化可直觀地了解和跟蹤旋翼損傷演化過程。

2 在重構相空間跟蹤旋翼損傷的方法

筆者所提的旋翼損傷跟蹤方法首先在多個鄰域內對狀態(tài)預測殘差統(tǒng)計平均形成一個損傷觀測特征向量,然后應用奇異譜分解方法從損傷觀測特征向量的時間序列中提取出槳葉損傷演化的維度和趨勢信息,采用雙指數平滑(double exponential smoothing,簡稱DES)方法[11]建立損傷演化趨勢預測模型并估計出槳葉損傷故障剩余壽命。程序主要包括損傷觀測信號生成模塊和損傷狀態(tài)優(yōu)化估計模塊(如圖6)。損傷觀測信號生成模塊事先根據健康旋翼在相空間內的運動軌跡建立基準狀態(tài)預測模型,運行時利用現場采集數據和預測模型產生短時狀態(tài)預測殘差并統(tǒng)計生產損傷觀測特征向量。損傷分析和剩余壽命預測模塊采用奇異譜分解方法獲得旋翼槳葉損傷維數和演化趨勢信息,利用損傷演化趨勢預測模型估計槳葉剩余壽命。

圖6 旋翼損傷監(jiān)測程序原理圖Fig.6 Rotor health monitoring system framework

2.1 采用Volterra預測模型生成損傷觀測信號

在旋翼系統(tǒng)狀態(tài)空間中,用損傷旋翼和健康旋翼軌跡之間的短時偏差作為系統(tǒng)的損傷觀測信號。設旋翼在某任意時刻的狀態(tài)為xn,損傷為φ,經過短時間τ后旋翼狀態(tài)變?yōu)閤(xn,τ,φ)。若按健康系統(tǒng)(φ=φR)估計其狀態(tài),則應為x(xn,τ,φR),兩狀態(tài)之間的殘差為

S(xn,τ)(φ-φR)+O(ε)

(4)

短時狀態(tài)殘差利用系統(tǒng)的動力學模型去除信號中因運行時間和系統(tǒng)狀態(tài)的正常轉移而造成的變化,保留了損傷造成的變化,是一種更好的損傷觀測信號。若能忽略高階小量的影響,式(4)可近似為

(5)

(6)

其中:k5(x0;φ)=(k1+2k2+2k3+k4)/6;k1=F(x0;φ);k2=F(x0+0.5τk1;φ);k3=F(x0+0.5τk2;φ);k4=F(x0+τk3;φ)。

應用式(6)計算得到短時狀態(tài)殘差后,就可研究靈敏度S(xn,τ)隨系統(tǒng)初始狀態(tài)和損傷大小的變化情況,具體分析高階小量的影響。下面分析實例旋翼槳尖揮舞位移的短時殘差靈敏度隨初始狀態(tài)和槳葉剛度損傷的變化情況。首先,在健康旋翼的狀態(tài)空間中隨機選擇16個點作為計算初始狀態(tài),在槳葉壽命周期(0

本文提出了(N,m)維修模型,在充分考慮不同部件的維修優(yōu)先級不同的條件下,對高維修優(yōu)先級部件,采用修復性維修以及預防性維修;對低維修優(yōu)先級部件,采用成組維修。不同于傳統(tǒng)的依賴于單一團隊的集中式維修模式,本文提出了異構多部件系統(tǒng)的競爭式維修模式。以全周期的最低維修費率為目標,通過計算不同部件的失效次序統(tǒng)計量,對參數進行優(yōu)化,進而得到最優(yōu)維修策略。通過數值算例驗證了維修策略的有效性,并對參數敏感度進行了分析。

圖7 短時狀態(tài)殘差與壽命和損傷的關系Fig.7 Curves of short time state error vs. life time and damage variable

(7)

其中:y為傳感器輸出序列;x為相點;xR為預測器輸出的參考點;ex為狀態(tài)預測殘差;ek為損傷觀測信號。Pha為相空間重構算子;m,τ分別為嵌入維數和時延參數;Bi為相空間鄰域(i=1,2,…,M);Ni為鄰域內相點數目;Vol為Volterra預測器;p為預測器階數,Volterra預測器的濾波系可根據健康旋翼的相軌跡采用自適應算法確定[12]。

該方法的計算流程如圖8所示。根據健康旋翼在相空間內的運動軌跡,采用非線性Volterra級數建立基準狀態(tài)短時預測模型。監(jiān)測時將旋翼的狀態(tài)采集數據運用延遲嵌入方法重構出實測相軌跡,采用Volterra預測模型計算出參考軌跡,將參考軌跡和實測相軌跡之差統(tǒng)計平均生成損傷觀測信號。由于Volterra預測模型在不同區(qū)域內的精度變化,單個短時狀態(tài)殘差靈敏度具有較大的不確定性。為提高損傷檢測的魯棒性并增加損傷模式的跟蹤數目,將參考相空間劃分成靈敏度大體相等且互不重疊的幾個區(qū)域Bi(i=1,2,…,M),狀態(tài)預測殘差對多種損傷類型的非線性特性可通過其在不同相空間區(qū)域內的靈敏度差異來體現,用各個區(qū)域內所有殘差的平均值來消除靈敏度矩陣的隨機誤差,將一條記錄中各區(qū)域內殘差均值組合形成一個維數足夠高的損傷觀測特征向量ek。

圖8 損傷觀測信號的計算流程Fig.8 Flow chart of damage measurement prediction

應用該方法對實例旋翼在不同損傷狀態(tài)下的槳葉位移和槳轂載荷進行分析。在槳葉壽命周期內(0≤t/tF<1),等時間間隔采集60個損傷程度(0≤φ<1)下的仿真記錄,將參考相空間按照坐標分量x1大小劃分成8個等間隔的鄰域,在各鄰域內統(tǒng)計每條記錄的短時狀態(tài)殘差,計算所得每個殘差特征向量ek包含32個列,各殘差特征向量隨使用壽命的變化曲線如圖9所示??梢?,各監(jiān)測信號的殘差特征向量可大致反映槳葉損傷曲線的變化趨勢,但數據依然有較大起伏,難以判斷其中含有幾種損傷發(fā)展模態(tài)。

圖9 殘差特征向量隨時間的變化Fig.9 Curves of residual tracking metric(Scaled)

2.2 槳葉損傷信息的提取

理論上可以利用殘差特征向量ek與損傷變量之間的函數關系φk=H(ek)來估計槳葉損傷。然而,旋翼的實際損傷情況比較復雜,損傷產生的原因可能來自疲勞、腐蝕和照射等不確定性因素,損傷變量的實際維數往往不能事先知道;有些損傷類型不但損傷擴展規(guī)律不確定,甚至沒有事先可直接測量的損傷變量,難以建立損傷與觀測之間映射關系。因此,需要直接應用ek序列分析損傷變量的維數和損傷演化的模態(tài)信息。奇異譜分析方法[13]能夠研究多維時序數據在模態(tài)空間中的能量分布情況,可用于分析旋翼的損傷維數和演化趨勢信息。

(8)

式(8)是一種Rayleigh商的優(yōu)化問題[13],可化為廣義特征值問題

(9)

其中:

Δt為序列的時間間隔;D為微分算子矩陣,大小為(n-1)n;λ為特征值;q為特征向量。

在Matlab軟件中廣義特征值方程(8)可采用奇異譜分解法求解:(Q,λ) =eig(YTY, (DY)TDY)。矩陣Y在模態(tài)坐標系中的投影分量為

φk=Yqk

(10)

其中:qk為第k個廣義模態(tài)向量;φk為第k個模態(tài)投影。

矩陣Y中既包含損傷變量緩慢演化信號,又包含旋翼系統(tǒng)的狀態(tài)變量快速運動的信號。如果殘差矩陣中有k個特征值比其它特征值高出一個數量級,則表明系統(tǒng)存在著k個緩慢的損傷演化模式,這些特征值對應的投影可指示損傷演化過程。

采用奇異譜分解算法對BO-105旋翼各監(jiān)測信號的殘差矩陣Y進行分析,得到的32個廣義特征值如圖10(a)所示,相應的,各殘差矩陣在32個模態(tài)坐標系中的投影分量如圖10(b)所示??梢?,殘差信號的能量按模態(tài)頻率重新分配后,有一個特征值遠遠大于其他特征值,說明殘差矩陣中有一種確定的趨勢分量。與之對應模態(tài)的投影分量較為清晰地反映槳葉損傷擴展的趨勢信息,容易驗證槳葉損傷與其之間有較為顯著的線性關系。將該模態(tài)投影分量作為一種名義損傷變量,可對旋翼的損傷狀態(tài)進行定量跟蹤。

圖10 廣義特征值和模態(tài)坐標投影Fig.10 Eigenvalues and modal projections

2.3 槳葉剩余壽命的預測

為估計槳葉剩余壽命需要建立合理的損傷趨勢預測模型。槳葉剛度損傷過程是一種非線性隨機過程,雙指數平滑方法[11]既能靈活適應回歸變量的非線性動態(tài)特性,又能快速計算多變量回歸預測模型的統(tǒng)計分析參數,是一種優(yōu)選的自適應預測技術。DES方法的基本預測公式為

(11)

當給定失效條件φF后,根據式(10)求得槳葉剩余壽命。采用DES方法對實例旋翼各殘差矩陣的模態(tài)投影分量進行分析,得到槳葉剩余壽命估計結果如圖11所示。在使用時間t≥0.2tF區(qū)間內,采用各特征模態(tài)投影分量對槳葉剩余壽命跟蹤的誤差如表1所示??梢?,這種方法能較為有效地跟蹤槳葉壽命過程。

圖11 應用不同觀測數據對剩余壽命估計Fig.11 Prediction of remaining time to failure using different measurements

表1 槳葉剩余壽命的跟蹤誤差

3 結束語

對結構損傷機理和材料失效理論的研究通常都是在微觀尺度上進行的。這種研究是必要的,但這類研究不能闡明這些損傷過程在整個系統(tǒng)宏觀行為中是如何表現的,也不能提供一種系統(tǒng)級的損傷跟蹤手段。源自混沌信號分析方法的相空間重構技術可以從不同的時空角度觀察系統(tǒng)的損傷行為特性。筆者利用信號重構技術和Volterra預測模型發(fā)展了一種新的損傷觀測殘差生成方法,對損傷觀測殘差序列分析可提取損傷變量的維度和趨勢信息。所提方法能充分利用被監(jiān)測信號的非線性特性在高維相空間中重建系統(tǒng)動力學本質,以不同的時間尺度來觀測系統(tǒng)的演化特性。該方法具有損傷模式的自動識別能力,可用于損傷演化模型或維數信息難以事先確定的場合。

[1] Azzam H, Andrew M J.The use of math-dynamic model to aid the development of integrated health and usage monitoring[J].Journal of Aerospace Engineering,1992,206(1):71-96.

[2] Ganguli R,Chopra I,Haas D J.Detection of helicopter rotor system simulated faults using neural networks[J].Journal of American Helicopter Society,1997,42(2):161-171.

[3] Pawar P M, Ganguli R.Fuzzy logic based health monitoring and residual life prediction for composite rotor[J].Journal of Aircraft,2007,44(3): 981-995.

[4] Pawar P M, Ganguli R. Helicopter rotor health monitoring-a review[J]. Journal of Aerospace Engineering,2007,221(1): 631-647.

[5] Johnson W. Comprehensive analytical model for rotorcraft aerodynamics and dynamics volumes ii:components theory [M].Palo Alto:Johnson Aeronautics,2007:280-283.

[6] Miles T J, Lucas M.Torsional and bending vibration measurements on rotors using laser technology[J]. Journal of Sound Vibration, 1999,266(3):1441-1467.

[7] Pawar P M,Ganguli R. On the effect of progressive damage on composite helicopter rotor system behavior[J]. Composite Structures, 2007,78:410-423.

[8] Chelidze D, Liu M.Dynamical systems approach to fatigue damage identification [J]. Journal of Sound Vibration,2005, 281(3):887-904.

[9] Hilborn R C.Chaos and nonlinear dynamics[M].Oxford:Oxford University Press,2006:375-389.

[10]Cao Liangyue.Practical method for determining the minimum embedding dimension of a scalar time series[J]. Physica D:Nonlinear Phenomena,1997,110(1-2):43-50.

[11]Guo Zhihao, Shahdi M, Shaya S. Multi-objective OLSR for proactive routing in manet with delay,energy, and link lifetime predictions[J].Applied Mathematical Modeling,2011,35:1413-1426.

[12]張家樹.混沌信號的非線性自適應預測技術及其應用研究[D].成都:中國電子科技大學, 2001.

[13]錢征文,程禮,李應紅.利用奇異值分解的信號降噪方法[J].振動、測試與診斷,2011,31(4):459-463.

Qian Zhengwen,Cheng Li,Li Yinghong. Noise reduction method based on singular value decomposition[J].Journal of Vibration,Measurement & Diagnosis, 2011,31(4):459-463.(in Chinese)

[14]Marquardt D. An algorithm for least-squares estimation of nonlinear parameters [J]. SIAM Journal on Applied Mathematics, 1963, 11: 431-441.

歡迎訂閱《振動、測試與診斷》

《振動、測試與診斷》由工業(yè)和信息化部主管,南京航空航天大學和全國高校機械工程測試技術研究會聯(lián)合主辦,是反映振動、動態(tài)測試及故障診斷學科領域的科研成果及其應用情況的技術性刊物。主要刊登國內外以振動測試與故障診斷為中心的動態(tài)測試理論、方法和手段的研究及應用方面的技術文獻,包括實驗測試技術、測試儀器的研制、方法和系統(tǒng)組成、信號分析、數據處理、參數識別與故障診斷以及有關裝置的設計、使用、控制、標定和校準等,不拘泥于行業(yè)和測試項目。

本刊為EICompendex數據庫收錄期刊和中文核心期刊,雙月刊,每逢雙月末出版,每本定價20元,全年120元。歡迎訂閱和投稿,歡迎在本刊刊登各類廣告和科技信息。

編輯部地址:南京市御道街29號 郵政編碼:210016 電話:(025)84893332

傳真:(025)84893332 E-mail:qchen@nuaa.edu.cn 網址:http:∥zdcs.nuaa.edu.cn

*國家自然科學基金資助項目(51375490)

2013-01-07;

2013-02-24

10.16450/j.cnki.issn.1004-6801.2015.01.020

V275.1; TH113.1

朱旭程,男,1973年11月,博士、副教授。主要研究方向為航空機械維修、綜合保障。曾發(fā)表《直升機旋翼健康監(jiān)測系統(tǒng)地面試驗臺方案設計與開發(fā)》(《兵工學報》2008年第29卷第4期)等論文。 E-mail:zhu133@189.cn

猜你喜歡
相空間槳葉維數
探究奇偶旋翼對雷達回波的影響
β-變換中一致丟番圖逼近問題的維數理論
束團相空間分布重建技術在西安200 MeV質子應用裝置的應用
一類齊次Moran集的上盒維數
立式捏合機槳葉結構與槳葉變形量的CFD仿真*
關于齊次Moran集的packing維數結果
非對易空間中的三維諧振子Wigner函數
涉及相變問題Julia集的Hausdorff維數
基于相空間重構的電磁繼電器電性能參數預測研究
直升機槳葉/吸振器系統(tǒng)的組合共振研究
嘉峪关市| 昌图县| 密云县| 西乌珠穆沁旗| 台州市| 灯塔市| 镶黄旗| 龙陵县| 东乌珠穆沁旗| 南投市| 苏尼特右旗| 德令哈市| 濉溪县| 龙游县| 安阳市| 西昌市| 宁蒗| 台北县| 察雅县| 乡城县| 高要市| 中方县| 博爱县| 平阳县| 祁连县| 晋江市| 理塘县| 阿城市| 象州县| 东台市| 华蓥市| 台南县| 和林格尔县| 体育| 安远县| 辉县市| 华宁县| 雷波县| 雷州市| 莱州市| 交口县|