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

?

基于線性剪切空間調制快拍成像測偏動態(tài)定標技術的溫度寬域變化研究

2024-05-21 00:00:00蔣思悅賈辰凌張晶潘楊柳唐金鳳江敏馮健萱毛智遠樊東鑫鄧婷王華華

DOI:10.16601/j.cnki.issn2096-7330.2024.01.008"文章編號:2096-7330(2024)01-0060-10

摘"要:空間調制快拍成像測偏技術通過單次測量即可獲得目標圖像的全部斯托克斯參數(shù),在民用、軍事、航空等領域有廣泛應用. 線性剪切空間調制快拍成像測偏技術(SMSIP)可在測量目標的同時實現(xiàn)系統(tǒng)定標,無需預知參考目標,從而能顯著提高測量的速度和精度. 然而,前期的動態(tài)定標技術研究并不適用于溫度寬域變化情況,因為溫度變化會對系統(tǒng)的空間調制相位產生影響,這將導致自定標解調算法失效.該文針對該問題提出一種動態(tài)定標算法,通過引入校準系數(shù)對SMSIP系統(tǒng)的動態(tài)定標原理進行理論推導,并用計算機仿真進行原理論證,對仿真結果進行評估,為線性剪切空間調制快拍成像偏振動態(tài)定標算法的在實際工程化應用提供理論指導.

關鍵詞:線性剪切;快拍成像側偏技術;空間調制;溫度寬域變化;動態(tài)定標算法

中圖分類號:TN256文獻標志碼:A

偏振是光除光譜、光強、相位以外的又一重要物理特性[1].傳統(tǒng)成像技術通常只能采集光強信息或波長信息,難以直接獲得目標的偏振信息.偏振成像技術在傳統(tǒng)成像的基礎上,通過增加偏振維度信息,可直接獲得目標光學輻射的光強信息和偏振信息[2].目標的偏振信息與其表面狀態(tài)、固有理化屬性等息息相關,各種物質表現(xiàn)出獨特的偏振特征.目前偏振成像技術已經得到了廣泛的關注和應用,如在軍事、民用、航空等領域,國內外研究人員進行了深入的探索和實踐[3,4,5].

偏振成像技術根據(jù)調制方式的不同,分為分時型和快照型(也稱為快拍式)兩種[2].分時型偏振成像技術在獲取圖像時使用固定的偏振片與旋轉波片,多次分時采集目標信息,獲得的偏振圖像具有較高的成像精度和準確性[6].但該技術只適用于靜態(tài)成像場景,且需要反復測量多次,延長了系統(tǒng)的數(shù)據(jù)采集時間,這就導致它的應用范圍受到限制.目前的快照型偏振成像技術可劃分為分振幅型、分孔徑型、分焦平面型和空間調制型四種[7].分振幅型的光學系統(tǒng)結構復雜、體積龐大;分孔徑型的可靠性低,光路設計復雜,系統(tǒng)成本較高;分焦平面型功耗低、消光比高,但通常需要在計算目標場景的偏振信息之前,需要預處理來對圖像原始分辨率進行恢復[8].前三種快照型偏振成像技術都需要對不同通道輸出的偏振圖的4個斯托克斯(Stokes)參量信息全部融合在圖像信息中,重構包含干涉信息的圖像,最后完成對目標偏振信息的重建[9].

近年來,國內外的研究學者在空間調制快拍成像側偏系統(tǒng)的核心調制器件上開展了大量研究.Oka等研發(fā)出以雙折射楔形棱鏡為核心調制器件的偏振成像技術[10].張淳民等研究出以Savart偏光鏡為核心調制器件的偏振成像技術[11].曹奇志等研發(fā)了以改進型Savart 偏光鏡為核心調制器件的偏振成像技術[12].在對包含干涉數(shù)據(jù)的目標圖像重構前,傳統(tǒng)方法采用的是參考光定標法,即選擇已知偏振狀態(tài)的光源作為參考光源,一般可以使用線偏振光或圓偏振光,將它們的偏振態(tài)作為參考數(shù)據(jù),再將目標圖像數(shù)據(jù)與參考數(shù)據(jù)進行歸一化即可重構出目標的偏振信息. 但參考光定標法對環(huán)境條件比較敏感,尤其是在室外環(huán)境溫度變化較大的情況下,成像系統(tǒng)的非線性響應、光學器件的畸變等都會影響定標的準確性.

本文在線性剪切空間調制快拍成像動態(tài)定標技術[2]的基礎上,針對溫度寬域變化情況提出了一種動態(tài)定標算法,該方法直接通過校準系數(shù)對空間相位實時測量來動態(tài)定標.本文第1節(jié)是系統(tǒng)的光路設計,系統(tǒng)的核心器件為兩塊MSP.第2節(jié)是溫度動態(tài)定標原理,在一定的溫度變化范圍內,改進后的算法能準確的進行目標重建.第3節(jié)是數(shù)值仿真,使用Matlab軟件,給定輸入Stokes圖像,模擬Stokes圖像經過系統(tǒng)之后產生的變換,再由動態(tài)定標算法重建Stokes圖像,并對重建Stokes圖像做了相應的系數(shù)分析.第4節(jié)是結論,指出了本文的創(chuàng)新點以及動態(tài)定標算法的應用范圍.

1"SMSIP的光路設計

SMSIP的光學設計如圖1所示.改進型薩瓦偏光鏡MSP1和MSP2的主截面與xoz平面平行,半波片HWP的快軸方向與x軸的夾角為22.5°,分析器A的偏振方向與x軸的夾角為45°.從目標發(fā)出的光經過準直后,平行于紙面入射,并通過濾光片得到準單色光. 接下來,入射光經過MSP1,被分成兩束線偏振光,分別位于xoz平面和yoz平面內且彼此垂直.這兩束光經過半波片HWP后,偏振方向將旋轉45°,然后垂直入射進入MSP2,并分成四束線偏振光.最后,這四束光經過分析器A和成像鏡L后干涉,將目標的圖像和干涉條紋疊加在CCD上.

2"溫度寬域變化下SMSIP系統(tǒng)動態(tài)定標原理

輸入光束的偏振狀態(tài)用Stokes矢量來描述:

Sin=[S0,S1,S2,S3]T,JY,2 (1)

其中S0是總光強,S1和S2分別代表兩個方向上的線偏振光,S3代表圓偏振光[13].入射光經過SMSIP系統(tǒng)的傳輸方程可表示為

Sout=MA(45°)M2MM(22.5°)M1Sin,JY,2(2)

MA(45°)表示分析器A的Muller矩陣、M1和M2表示兩個改進型薩瓦偏光鏡的Muller矩陣、MH(22.5°)表示HWP的Muller矩陣[14].

將式(2)展開計算后可得到焦平面上光強公式:

I=S0,out=0.5S0+0.5S1cos(φ)+0.25S2[cos(φ2+φ1)-cos(φ2-φ1)]+

0.25S3[sin(φ2+φ1)+sin(φ2-φ1)],JY,2 (3)

其中φ1=2π×2Ω×X,φ2=2π×2Ω×2X,Ω=Δ/λf是空間載波頻率,λ是入射光波的中心波長,f是成像透鏡 L2的焦距,Δ是MSP1的單板橫向剪切量,φ1和φ2分別是 MSP1和 MSP2的空間調制相位,改進型薩瓦偏光鏡的厚度比決定φ1和φ2的比值[15].

對光強I進行傅里葉變換后可將Stokes矢量分離到不同的頻譜通道C[16]i:

Ci=F(I)·WiJY,2 (4)

其中F為傅里葉變換符號,i為通道號,Wi為窗函數(shù).變換后可得到含7個峰值的干涉圖.

通過對通道(峰)進行反傅里葉變換,將各個通道頻譜信息表示為:

Fi=F-1(Ci).JY,2(5)

當MSP1和MSP2的厚度按1∶2比例配置時,各個通道頻譜信息與Stokes矢量分量之間的關系為

{F0=0.5·S0

F1=-0.125·(S2+jS3)exp[j(φ2-φ1)]

F2=0.25·S1exp(jφ2)

F3=0.125·(S2-jS3)exp[j(φ2+φ1)].JY,2(6)

用校準系數(shù)Ki來表示式(6)中的幅度項和相位項:

Ki={K0=2

K1=-8·exp(-j[φ2-φ1])

K2=4·exp(-jφ2)

K3=8·exp(-j[φ2+φ1]).JY,2(7)

將(6)與(7)相乘得到

AKS·=[F0·K0

F2·K2-Re(F1·K1)

Im(F1·K1)]或AKS·=[F0·K0F2·K2-Re(F3·K3)

Im(F3·K3)].JY,2(8)

為了精確地重建Stokes矢量,必須精確計算校準系數(shù)Ki,校準系數(shù)可以看作初始校準系數(shù)Kr,i與動態(tài)校準系數(shù)Kd,i的乘積.Kr,i可在系統(tǒng)的初始狀態(tài)下求得.Kd,i在測量過程中隨改進型薩瓦偏光鏡的空間相位變化實時測量.校準系數(shù)表示為:Ki=Kr,i·Kd,i.JY,2(9)

初始校準系數(shù)Kr,i可通過測量一組已知偏振態(tài)的參考光束得到.

S(0°)=[1"1"0"0]T,JY,2(10)

S(45°)=[1"0"1"0],JY,2(11)

其中S(0°)為0°參考光束的斯托克斯矢量,S(45°)為45°參考光束的斯托克斯矢量.

0°參考光束的測量結果用于校準C2通道,45°參考光束的測量結果用于校準C1和C3通道,初始校準系數(shù)Kr,i如下:

{Kr,0=2,

Kr,1=SX(|Kr,0·F(45°)0|F(45°)1SX),

Kr,2=SX(|Kr,0·F(0°)0|F(0°)2SX),

Kr,3=SX(|Kr,0·F(45°)0|F(45°)3SX),JY,2(12)

i為通道號,F(xiàn)(0°)i為系統(tǒng)測量0°參考光束時各通道包含的頻譜信息,F(xiàn)(45°)i為系統(tǒng)測量45°參考光束時各通道包含的頻譜信息.

動態(tài)校準系數(shù)Kd,i在這里主要討論由溫度寬域變化引起MSP雙折射率和厚度變化:

Kd,i(T)(Δn,d),JY,2(13)

其中Δn表示MSP雙折射率,d表示MSP厚度.

溫度對MSP的主要作用是影響SMSIP系統(tǒng)中MSP的空間相位.MSP空間相位可表示為

φi=SX(4πλfSX)·X·di·Δn,JY,2(14)

其中φi(i=1,2)為MSP空間相位,λ是入射光波的中心波長,X是空間位置,di(i=1,2)為MSP厚度.

MSP的空間相位延遲量與MSP雙折射率和MSP厚度有關,溫度變化對MSP的這兩個方面都有影響.當溫度發(fā)生微小變化時,MSP空間相位也相應發(fā)生變化:

Δφi=SX(ΔKG-*1/4φiTSX)·ΔT=

SX(4πλfSX)·X·(SX(diTSX)·Δn+di·SX(ΔKG-*1/4nTSX))·ΔT=

[SX(4πλfSX)·X·(di·Δn·SX(1diSX)·SX(diTSX))+SX(4πλfSX)·X·(di·Δn·SX(1ΔnSX)·SX(ΔKG-*1/4nTSX))]·ΔT=

SX(4πλfSX)·X·di·Δn·(SX(1diSX)·SX(diTSX)+SX(1ΔnSX)·SX(ΔKG-*1/4nTSX))·ΔT=

φi·(SX(1diSX)·SX(diTSX)+SX(1ΔnSX)·SX(ΔKG-*1/4nTSX))·ΔT=r·φi·ΔT,JY,2(15)

r=SX(1diSX)·SX(diTSX)+SX(1ΔnSX)·SX(ΔKG-*1/4nTSX),JY,2(16)

其中Δφi(i=1,2)為MSP相位延遲量,SX(ΔKG-*1/4nTSX)是MSP的雙折射率溫度系數(shù),SX(1diSX)·SX(diTSX)(i=1,2)是MSP在垂直于光軸方向上的線性熱膨脹系數(shù),r是MSP的材料系數(shù)[17].

文獻[18]給出了求解Δφ2的方法.將公式(6)中F2平方,得到

F22=SX(116SX)|S1|2e2jφ2.JY,2(17)

再將F1與F3相乘,得到

F1·F3=-SX(164SX)|S2+jS3|2·e2jφ2.JY,2(18)

公式(17)和(18)得到的相位角與φ2有關,因此可以通過計算公式(17)或公式(18)的相位角得到φ2.但是,兩個公式中|S1|2或|S2+S3|2的值可能為0,因此可通過計算兩個公式的加權和來獲得φ2.

F22-4·F1·F3=SX(116SX)(S21+S22+S23)·e2jφ2JY,2(19)

用angle函數(shù)對公式(19)取相位角可得到2φ2的值,angle函數(shù)提取相位角后,可得到區(qū)間在[-π,π)之間的相位.

angle(F22-4·F1·F3)=2φ2.JY,2(20)

MSP2空間相位延遲量Δφ2的值與Kr,i的相位角有關,根據(jù)公式(12)可計算出K2r,2和Kr,1Kr,3的加權和并取出相位角.

angle(4·K2r,2-Kr,1·Kr,3)*=2φ2+2Δφ2.JY,2(21)

由(20)和(21)便可求得Δφ2的值:

2Δφ2=angleSX(F22-4F1F3(4K2r,2-Kr,1Kr,3)*SX).JY,2(22)

angle函數(shù)取出相位角的范圍為(-π,π),Δφ2的取值范圍為(-π/2,π/2).

因此,由式(15)~(22)可知,SMSIP系統(tǒng)的溫度變化量為

ΔT=SX(Δφ2r·φ1SX)=|SX(π/2r·φ2SX)|.JY,2(23)

由公式(23)可得SMSIP系統(tǒng)的溫度變化范圍為

T0+SX(SX(-π2SX)8πΩXrSX)<T<T0+SX(SX(π2SX)8πΩXrSX).JY,2(24)

SMSIP系統(tǒng)中初始溫度T0為23℃,r的取值為-51×10-6/℃[19],Ω=Δ/λf,Ω的值為13.3703,X的取值為像元尺寸與最大光程差像元空間位置的乘積,X=5.31μm×512,經計算后可得寬域溫度范圍為-10.71℃≤T≤56.71℃.

動態(tài)校準系數(shù)Kd,i可表示為

{Kd,1=exp(-0.5j×r×ΔT×φ2),

Kd,2=exp(-j×r×ΔT×φ2),

Kd,3=exp(-1.5j×r×ΔT×φ2).JY,2(25)

將Kd,i和Kr,i代入公式(8)就能重建Stokes矢量.

在實際工程應用中,SMSIP系統(tǒng)的工作溫度可能高于56.71℃.然而,隨著溫度的增加,MSP2空間相位延遲量Δφ2也進一步增大,當Δφ2大于動態(tài)定標算法能夠校準的最大值時,重建Stokes矢量會出現(xiàn)相位混亂,導致重建Stokes矢量失真.為了擴展動態(tài)定標算法適用的溫度范圍,我們需要進行相位修正.用polyfit函數(shù)對Δφ2與空間位置X進行一階線性擬合,將擬合得到的曲線與原來的Δφ2進行比較,兩條曲線之間有一個大小為n×2π的偏移量(n為任意整數(shù)).

n=argmin2{|aX-(2Δφ2+n·2π)|}.JY,2(26)

將n值代入Kd,i中可擴展SMSIP系統(tǒng)的應用范圍,溫度變化范圍可拓展至-47℃≤T≤ 93℃.修正后的Kd,i如下所示:

{Kd,1=exp[-0.5j×(r×ΔT×φ2+n×π)],

Kd,2=exp[-j×(r×ΔT×φ2+n×π)],

Kd,3=exp[-1.5j×(r×ΔT×φ2+n×π)].JY,2(27)

3"計算機仿真

為了驗證溫度寬域變化時線性剪切空間調制快拍成像動態(tài)定標算法的可行性,我們在MATLAB平臺上進行了仿真實驗,仿真參數(shù)設置為:光源的中心波長為6328nm;CCD 分辨率為1024pixel×1024pixel,像元尺寸為531μm×531μm;改進型薩瓦偏光鏡的折射率分別為no=166527和ne=148956,MSP1和MSP2的尺寸分別為25mm×25mm×16mm 和25mm×25mm×32mm;透鏡 L1和L2的焦距為f1=f2=105mm;空間載頻設置為133703;矩陣點數(shù)N設置為1024.

動態(tài)定標算法相位修正前的溫度適用范圍為-1071℃≤T≤5671℃,當環(huán)境溫度超過5671℃時,重建Stokes矢量相位發(fā)生反轉.

3"相位修正前,當溫度超過56.71℃,重建Stokes矢量出現(xiàn)誤差

當環(huán)境溫度為40℃時,Stokes矢量輸入圖像如圖4所示,輸入S1為peak函數(shù),S2為peak’函數(shù),S3為兩個正弦函數(shù)的乘積,S0為S12+S22+S32的平方根;圖5為模擬干涉圖;圖6為與圖5對應的頻譜圖,可以清晰地看到七個通道(峰值);圖7為重建Stokes目標圖像,重建圖像與輸入圖像基本一致.

當環(huán)境溫度為90℃時,2Stokes矢量輸入圖像不變,重建Stokes目標圖像如圖8所示,重建Stokes目標圖像與Stokes矢量輸入圖像之間的差圖像如圖9所示.與Stokes矢量輸入圖像相比,重建Stokes目標圖像存在明顯的問題.在觀察圖像時,我們可以注意到以下幾個方面的不足:首先,重建S1目標圖像左右兩側出現(xiàn)相位混亂,部分屬于圖像左邊的相位被解調到了圖像右邊.這種相位的錯位會導致圖像失真,影響重建S1目標圖像的準確性和可視化效果. 其次,重建S2目標圖像和重建S3目標圖像左右兩側均有部分相位重建效果不理想,導致圖像中的特征信息丟失,圖像的細節(jié)和結構無法得到準確表達.

為了提高重建Stokes目標圖像的質量,使用含相位修正的動態(tài)定標算法來重建Stokes目標圖像.通過提高動態(tài)定標算法的相位重建精確度來優(yōu)化相位恢復過程,及時糾正相位混亂,提高了重建圖像的精度,有助于重建Stokes目標圖像更好地保留輸入圖像的特征和細節(jié).當環(huán)境溫度為90℃時,相位修正后的重建Stokes目標圖像如圖10所示,重建Stokes目標圖像與Stokes矢量輸入圖像之間的差圖像如圖11所示,相位修正后重建Stokes目標圖像準確的還原了Stokes矢量輸入圖像的形狀、邊緣、紋理等特征.

結構相似度(SSIM)、相關系數(shù)(CORR)、均方誤差(MSE)是常用的統(tǒng)計分析方法,用于判斷兩幅圖像之間的相似性、相關性以及均方誤差.SSIM是一種用于衡量兩個信號或圖像之間結構相似性的評估系數(shù),考慮了亮度、對比度和結構等方面的信息,通過計算三個主要的分量來評估兩個圖像之間的相似程度:亮度相似度(Luminance Similarity)、對比度相似度(Contrast Similarity)和結構相似度(Structure Similarity),三個分量最終綜合得到一個0到1之間的值,值越接近1表示兩個圖像結構越相似.CORR是衡量兩個圖像之間線性相關程度的評估系數(shù),其取值范圍為-1到1,相關系數(shù)接近1時,表示兩個圖像呈正相關. MSE的計算方法是將兩幅圖像對應像素點之間的差值平方后求和并除以像素點的總數(shù),這樣得到的結果就是MSE,MSE越小,則說明兩幅圖像相似度越高;反之,MSE越大,則說明兩幅圖像相似度越低.

由表1和表2對比可知,+1當環(huán)境溫度為90℃時,經過相位修正后的重建Stokes目標圖像的SSIM參數(shù)均在0.987以上,CORR參數(shù)均在0.994以上,MSE參數(shù)均小于3.02×10-5.重建Stokes目標圖像保持了與原始圖像相似的結構特征,并且在亮度、對比度和紋理等方面也能夠準確地復原.此外,重建Stokes目標圖像還具有較高的線性相關度,與原始圖像之間存在著較強的相關性.結果表明,在環(huán)境溫度為90℃時,重建Stokes目標圖像能夠以很高的精度還原輸入圖像的各種特征,圖像重建效果良好.

4"結論

本文在線性剪切空間調制快拍成像定標技術的基礎上,探討了在溫度寬域變化時,動態(tài)定標算法的可行性. 對動態(tài)定標算法進行了理論分析,并通過計算機仿真重建Stokes圖像. 仿真結果表明SMSIP系統(tǒng)的應用溫度范圍為-47℃~93℃.動態(tài)定標算法拓寬了線性剪切空間調制快拍成像定標技術的適用溫度范圍,推動了空間調制快拍成像技術的工程化應用進程.

參考文獻:

[1]"高毅,于津強,張笑東,等.基于偏振自適應融合圖像的水下物證探測方法[J].紅外技術,2023,45(9):962-968.

[2]"曹奇志,唐金鳳,潘楊柳,等.線性剪切空間調制快拍成像動態(tài)定標技術[J].物理學報,2022,71(15):150-156.

[3]"柯劍寒,岳鈞百,程雪岷,等.海洋生物水下原位監(jiān)測技術及其在偏振維度的信息拓展[J].水下無人系統(tǒng)學報,2023,31(4):614-623.

[4]"陳雪嬌.基于超材料的近紅外偏振芯片研制[D]. 長春:長春理工大學, 2022.

[5]"邱鵬,王國聰,張曉明,等.2.16 m望遠鏡偏振光度計[J].紅外與激光工程,2023,52(7):220-233.

[6]"李鑫權.基于PSIM機載大視場寬譜段偏振光譜成像系統(tǒng)光學設計[D].長春:中國科學院大學(中國科學院長春光學精密機械與物理研究所),2023.

[7]"張紫楊,葉松,熊偉,等.基于空間調制偏振成像的Savart棱鏡參數(shù)標定[J].光學學報,2022,42(21):60-68.

[8]"趙永剛,龍錦寧,孫春生.分焦平面型偏振相機成像誤差分析[J].艦船科學技術,2023,45(2):159-162.

[9]"郝平波.空間調制全偏振成像帶寬影響分析及仿真研究[D].桂林:桂林電子科技大學,2021.

[10]Oka K, Kaneko T. Compact complete imaging polarimeter using birefringent wedge prisms[J]. Optics Express, 2003,11(13):1510-1519.

[11]張淳民,楊建峰,原新晶,等.偏振干涉成像光譜技術研究進展[J].光電子·激光,2000(4):444-448.

[12]曹奇志,張淳民.基于改進薩瓦偏光鏡的快拍成像偏振光譜儀[C].陜西省光學學會.2012年西部光子學學術會議論文摘要集,2012:16.

[13]馬佳琪.基于手性二維鈣鈦礦的偏振分辨光電探測器的研究[D].武漢:華中科技大學,2022.

[14]曹奇志,張晶,Edward DeHoog,等.空間調制穩(wěn)態(tài)微型快拍成像測偏技術研究[J].物理學報,2016,65(5):31-42.

[15]Pan Yangliu, Zhang Jing, Jiang Min, et al.Demodulation of polarization information for a spatially modulated snapshot imaging polarimeter based on the coherent demodulation theory[J].Applied Optics,2022,61(21):6349-6355.

[16]Chrysler B D, Oka K, Otani Y, et al. Dynamic calibration for enhancing the stability of a channeled spectropolarimeter[J]. Applied Optics, 2020,59(30):9424-9433.

[17]李景鎮(zhèn).光學手冊[M].西安:陜西科學技術出版社,1986:1304.

[18]Chrysler B D, Otani Y, Hagen N. Dynamic calibration of a channeled spectropolarimeter for extended temperature stability[C].Polarization Science and Remote Sensing IX. SPIE, 2019:156-164.

[19]Kimura T, Saruwatari M. Temperature compensation of birefringent optical filters[J]. Proceedings of the IEEE, 1971,59(8):1273-1274.

[責任編輯:彭喻振]

收稿日期:2023-11-10

*基金項目:國家自然科學基金(11964021;11664004);中國科學院西安光學精密機械研究所瞬態(tài)光學與光子技術國家重點實驗室開放課題(SKLST202220;SKLST202003)

第一作者簡介:蔣思悅(1999—),女,江西樟樹人,碩士研究生,研究方向:快拍偏振成像技術。

通信作者簡介:賈辰凌(1988—),男,講師,研究方向:高光譜偏振成像.Email: xxjsjcl@qq.com;

張晶(1979—),女,教授,研究方向:快拍偏振成像技術.Email:jingzhang2011@163.com。

宁津县| 克拉玛依市| 库尔勒市| 沂水县| 鹤壁市| 凯里市| 青川县| 黄石市| 虞城县| 武功县| 丽江市| 浦东新区| 定远县| 广西| 上高县| 堆龙德庆县| 丰县| 丽水市| 神农架林区| 裕民县| 阳原县| 锦屏县| 洪雅县| 改则县| 岳阳市| 会昌县| 宿松县| 嫩江县| 宜黄县| 马公市| 横山县| 嘉祥县| 邮箱| 乐平市| 常州市| 师宗县| 汶上县| 玛纳斯县| 建瓯市| 来宾市| 南充市|