?

植被覆蓋區高光譜遙感影像上蝕變巖與蝕變礦物信息的提取

2015-12-13 05:57徐元進馬洪超孟鵬燕楊明國
地球學報 2015年2期
關鍵詞:變巖礦物光譜

徐元進 ,馬洪超,孟鵬燕,楊明國

1)中國地質大學(武漢)資源學院數學地質遙感地質研究所,湖北武漢 430074;2)武漢大學遙感信息工程學院,湖北武漢 430079

資源勘查中,高光譜遙感已得到了廣泛的應用(Rowan et al.,2000;甘甫平等,2002;闞明哲等,2005;劉圣偉等,2006;Van der Meer,2006a,2012;Chen et al.,2007;Gersman et al.,2008;Bedini,2011;Bishop et al.,2011;高建陽,2011;王潤生等,2011;郭娜等,2012;Pour et al.,2013)。在該應用中,大多數研究和熱液蝕變作用有關,幾乎都是通過識別蝕變礦物去研究資源前景(Van der Meer et al.,2012),使用的參考光譜來自于JPL等光譜庫中礦物光譜或從影像上提取的礦物光譜。但是,在植被覆蓋區植被等會導致像元光譜中一些蝕變礦物信息微弱,難以被有效提取(Bishop et al.,2011)。同時,由于一些礦物(如高嶺石)既能在蝕變作用下形成,也能在非蝕變作用下形成,所以從影像上提取的蝕變礦物信息中,除了有真正找礦意義的蝕變礦物信息外,還包含一些形成于非蝕變作用下無找礦意義的礦物信息。

由于蝕變巖除了包含蝕變礦物的特征外,還包含了其它一些特征(如原巖),因而在遙感影像上,它的信息較其表面蝕變礦物信息更豐富。所以,通過識別信息相對較多、與成礦直接相關的蝕變巖去研究資源前景是很有意義的,這不僅為微弱礦物信息難提取的問題提供一種思路,也可避免一些無找礦意義的礦物信息。本研究選取植被覆蓋的云南普朗斑巖銅礦區作為研究區,從影像上分別提取蝕變巖和蝕變礦物信息,分析它們對找礦的有用性。

在高光譜遙感影像上識別地物,光譜匹配技術是進行定量比較、最直接有效的方法,該技術通過比較參考光譜和像元光譜的形態,表達兩者的相似性(Van der Meer,2006b,2012)。目前,該技術已得到了廣泛應用(Baugh et al.,1998;闞明哲等,2005;Debba et al.,2005;劉圣偉等,2006;Van der Meer,2006b;Kumara et al.,2010;Mountrakis et al.,2011)。SAM(spectral angle mapper)是該技術中應用較多的方法之一(Crósta et al.,1998;Dennison et al.,2004;闞明哲等,2005;Bishop et al.,2011;高建陽,2011;Pour et al.,2013)。本研究將使用SAM,從光譜整體形態出發進行蝕變礦物信息提取。在蝕變巖信息提取中,本研究不僅考慮光譜的整體形態,而且考慮光譜局部形態(本文著重于吸收谷位置形態);識別結果除提供可能擁有蝕變巖的像元外,還提供這些像元光譜的匹配值,便于判斷像元中存在蝕變巖信息的可能性。本文在整體形態匹配結果上比較蝕變礦物和蝕變巖的信息,判斷它們對找礦的有用性;同時,比較局部形態匹配前后的蝕變巖信息,目的是為了論證本研究方法的有效性。

1 方法

1.1 蝕變礦物信息提取思路

本研究使用的高光譜遙感影像數據來自于美國EO-1衛星傳感器Hyperion,為L1R產品。該產品共有242個波段,其中 198個波段(VNIR8-57和SWIR77-224)經過了輻射定標處理,同時由于VNIR56、57與 SWIR77、78重疊,所以只使用了196個波段(426.82—2395.50 nm)。從影像上提取蝕變礦物信息,包含了最小噪聲變換(MNF)、像元純度指數(PPI)、n維可視化與端元選取、端元蝕變礦物識別、SAM方法進行蝕變礦物成圖等過程(Kruse,2012),這些過程是在ENVI軟件上實現的。在此之前,還必須進行大氣校正,本文使用的是 FLAASH模型校正后的結果。

在云南普朗地區的遙感找礦中,Bishop等(2011)使用高光譜Hyperion數據,使用上述過程進行了該區蝕變礦物信息提取,取得了一些成果??紤]到影像上提取的礦物光譜與JPL等光譜庫中對應礦物的光譜存在一定差異,而與影像上其它區域的該礦物光譜有相同尺度的獲取環境,Bishop等(2011)使用的參考光譜是從影像上提取的礦物端元光譜。本研究將采用 Bishop等的思路進行蝕變礦物信息提取,目的是為了得到和前人相同思路的結果。

圖1 使用蝕變巖光譜進行高光譜遙感影像上蝕變巖識別流程圖Fig.1 Flow chart of the identification of altered rocks in hyperspectral image using the spectra of altered rocks

1.2 蝕變巖信息提取思路

使用不同蝕變類型的野外蝕變巖光譜作為參考光譜,從光譜的整體形態和吸收谷位置光譜形態兩方面,進行高光譜遙感影像上蝕變巖識別,其流程如圖1,作者在VC++下實現了該流程(Xu et al.,2010)。

1.2.1 數據預處理

蝕變巖光譜測試采用美國ASD公司FieldSpec Pro FR野外光譜輻射儀,該儀器能獲取 2151個波段(350.00—2500.00 nm)的光譜。由于野外光譜在1360—1400 nm、1810—1915 nm和2380—2500 nm處被大氣水汽吸收干擾(圖 3),所以光譜匹配前對這些波段數據進行了干擾剔除。

由于蝕變巖光譜的波譜范圍寬于像元光譜,且波段數遠大于像元光譜的波段數,所以光譜匹配前應對蝕變巖光譜進行取值,獲取一致的光譜范圍及波段。取值方法為:在像元光譜的對應波段處獲取蝕變巖光譜的反射率。

1.2.2 光譜匹配

考慮到吸收谷位置匹配的復雜性,本研究首先進行光譜整體形態匹配,獲取那些可能存在蝕變巖的像元,然后在該結果上進行吸收谷位置光譜形態匹配。吸收谷位置的匹配是通過計算各個對應吸收谷位置的相關系數,然后對光譜整體形態匹配的結果進行進一步判定。

(1)光譜整體形態的匹配

通過計算每個像元光譜和所有蝕變巖光譜的皮爾森相關系數(Pearson correlation coefficient),選出最大相關系數的蝕變巖光譜(以下簡稱“最可能光譜”)。如果這個最大相關系數大于一個指定值(比如下文第2部分中的0.80),認為該像元可能存在對應“最可能光譜”的蝕變巖。

(2)吸收谷位置光譜形態的匹配

首先,必須確定吸收谷位置波長。在連續統歸一化后光譜曲線(Van der Meer,2004)上,每個吸收谷中反射率最小處的波長即為吸收谷位置波長(Xu et al.,2010)。

圖2 研究區侵入巖的地質簡圖Fig.2 Simplified geological map of the intrusive rocks in the study area

圖3 六種蝕變巖光譜曲線Fig.3 Curves of six kinds of spectra of altered rocks in Table 1

吸收谷位置光譜形態的匹配,是像元光譜和它的“最可能光譜”,在每個吸收谷波段范圍內,計算皮爾森相關系數。在每個吸收谷位置,存在左右兩端反射率的極大值點,這兩點的波長范圍即為該吸收谷的波長范圍(Xu et al.,2010)。

1.2.3 識別結果成圖

識別結果可通過兩幅圖表達出來。第一幅圖是蝕變巖信息分布圖(如圖4b),第二幅圖是相關系數圖(如圖4c)。由于二維圖像的局限性,第二幅圖表達的是光譜整體形態的匹配系數(如圖4c)。

對于第一幅圖中那些擁有蝕變巖的像元,即使在第二幅圖中整體形態的相關系數較大,如果該像元的大部分吸收谷相關系數太小(比如大于80%吸收谷的相關系數為0.5以下),還是認為該像元存在蝕變巖的可能性較小。在下文第2部分中,我們研究的主要是那些存在蝕變巖可能性較大的像元,如圖4d(即整體形態的相關系數不小于0.80,且所有吸收谷位置的相關系數都不小于0.85)。對于整體形

態的相關系數不小于 0.80,吸收谷位置處相關系數不全是較大的(比如部分相關系數低于 0.85)的像元或區域,我們對其地面物體及其野外環境進行調查分析,檢驗方法的適用性。

表1 六種蝕變巖光譜信息Table 1 Information concerning the six spectra of altered rocks

圖4 高光譜遙感影像及其識別結果圖Fig.4 Hyperspectral RS image of the study area and diagrams of its recognition results

表2 圖4b或4d中四個圈定點及實地檢驗結果表Table 2 Identification results and field survey results of the four marked targets in Fig.4b or 4d

2 應用與分析

2.1 研究區

研究區位于中國云南省中甸縣普朗斑巖銅礦 區,其 四 角 坐 標 為 :該區屬于原始森林覆蓋區,地物類型主要是巖石、植被和土壤,地表蝕變巖露頭的原巖大部分是侵入巖(主要是石英閃長玢巖)。圖2為研究區侵入巖的地質簡圖,侵入巖的外圍區域是以碳酸鹽等為主的沉積巖區。

2.2 蝕變巖信息提取及結果分析

2003年9月,在研究區進行了野外光譜采集工作。最終選取6種蝕變巖光譜,其信息如表1,光譜曲線如圖3,表1中光譜文件名與圖 3中光譜曲線名對應。高光譜影像數據獲取時間為2003年9月,圖4a為研究區高光譜遙感影像,圖4b、c、d是識別結果圖。圖4b和4d中蝕變巖的色標號與表1中蝕變巖光譜的編號對應。

由圖4b、c、d和圖2可知:圖4b和4d中,蝕變巖信息主要集中在侵入巖及其周圍區域;圖4c中,可能存在蝕變巖的像元的相關系數都超過0.80,這說明這些像元的光譜和圖3中對應的“最佳匹配光譜”有較高相關性。由此說明:使用蝕變巖光譜進行影像中蝕變巖信息提取是可行的。

我們進行了實地檢驗,識別結果圖中數十種蝕變巖圈定點被驗證。如表 2,該表中列舉了 4個典型圈定點的識別結果與野外實地結果。由表2可知:除了部分表面礦物、覆蓋物與實地情況稍微有些不同外,識別結果與實地地物基本一致。

為了論證吸收谷位置識別的重要性,圖4b和4d中蝕變巖區被對比研究。在圖4b和4d中的侵入巖及周圍區域,蝕變巖信息分布區域大致相同,只是存在信息多少的不同。在沉積巖區,圖4b中兩區(A1、A2),主要為表面有褐鐵礦、孔雀石的蝕變石英閃長玢巖(表 1中編號為 3);但是圖4d中,這兩區(A1、A2)無蝕變巖信息。野外調查表明:這兩區沒有發現蝕變巖,但很容易見到“鐵銹”現象,可能是沉積巖中含鐵物質風化氧化所致,因而圖4b的這兩區的信息是假蝕變巖信息。這種結果可能是因為這兩區中含有鐵氧化物信息的像元光譜和表面有褐鐵礦的石英閃長玢巖的光譜有較高相關性,它們中存在共同的信息(即鐵氧化物信息)。因此,為保證識別的精確性,避免一些假的蝕變巖信息,光譜匹配除從整體形態上匹配外,還應從局部吸收谷形態上進行匹配。

2.3 蝕變礦物信息提取及結果分析

由于 Hyperion高光譜數據的短波紅外范圍(SWIR)的信噪比低于可見光近紅外范圍(VNIR),為了避免噪聲的影響,Bishop等(2011)用VNIR(400—1300 nm)的數據提取鐵礦物信息,用SWIR(2000—2400 nm)的數據提取含羥基的礦物信息。本研究采用相同方法提取這兩類礦物信息。圖5a和5c分別是從影像中提取的粘土礦物和鐵礦物的端元光譜曲線,這些端元光譜對應礦物的識別思路為:使用JPL等光譜庫中光譜,對端元光譜進行匹配,挑選出與每個端元光譜匹配較好的前 10位光譜,通過對比分析,最終決定該端元光譜對應的是哪種礦物。圖5b和5d分別是圖5a和5c的端元光譜曲線對應的JPL光譜庫中粘土礦物光譜曲線和鐵礦物的光譜曲線。

本研究使用圖5a和5c中提取的礦物光譜作為參考光譜,進行蝕變礦物信息提取。當所有光譜角的值設置為0.04時,圖6a、b分別是提取的粘土礦物和鐵礦物的信息分布圖。將圖6a、b中信息與圖4a、d中信息進行對比,可知:圖6a、b中信息明顯多于圖4a、d中信息,但是在圖6a、b中四處蝕變巖圈定點(R2、R4、R8和 R9)幾乎無蝕變信息。當增大所有光譜角的值(比如0.1)時,圖6a、b中出現大量新信息,但是這四處圈定點仍無信息。由此說明:遙感影像上,這些圈定點的蝕變礦物信息極其微弱,以致無法被有效提取出;這些圈定點的蝕變巖信息能被有效提取出,可能是因為在影像上這些蝕變巖包含了比它們表面的蝕變礦物更多的信息(如原巖等)。

由圖4b、圖4d、圖6和圖2,可知:

(1)相比較圖6a,圖6b中蝕變信息的位置和輪廓更接近圖4b(或4d);圖4b和圖6b中蝕變信息都主要存在于侵入巖及其周圍。這種現象可能是因為在侵入巖及其周圍的這些鐵礦物主要是熱液蝕變作用的產物或其風化氧化后產物,而本研究中蝕變巖也是熱液蝕變作用的產物。

(2)圖6b中A1、A2區主要是針鐵礦和赤鐵礦,這與 2.2部分中在這兩區驗證的“鐵銹”現象是接近的。但是,這兩區位于沉積巖區,沒有發現蝕變巖,因而這兩區“鐵銹”可能是含鐵物質風化氧化所致。第2.2部分中研究已表明:圖4b中A1、A2區是表面有褐鐵礦的蝕變石英閃長玢巖,雖然該區沒有發現蝕變巖,但卻指示該區有鐵的氧化物;通過增加吸收谷位置的識別,圖4b中A1、A2區假的蝕變巖信息可以避免(如圖4d中A1、A2區)。

圖5 影像上提取的礦物光譜曲線(a,c)與JPL光譜庫中礦物光譜曲線的對比(b,d)Fig.5 Comparison of extracted mineral signatures from Hyperion analysis in comparison with JPL library spectra

圖6 蝕變礦物信息圖Fig.6 Altered mineral maps

(3)在圖6a中,硬石膏是信息最多的礦物,其次是高嶺石。除部分信息存在于侵入巖及其周圍區域外,更多的硬石膏和高嶺石信息存在于沉積巖區。這些存在于沉積巖區的信息可能來自于沉積巖,因而可能是假的蝕變信息。

上述分析表明:遙感影像上,通過蝕變巖光譜獲取的蝕變信息比通過蝕變礦物光譜獲取的蝕變信息更可靠,原因可能是這些蝕變巖僅形成于熱液蝕變作用下,是與研究區斑巖銅礦的形成直接相關,而使用的一些“蝕變礦物”既可在蝕變作用下形成,也可在非蝕變作用下形成。

3 結論

(1)在蝕變巖識別圖上,4個典型蝕變巖點的野外驗證表明:除了部分表面礦物、覆蓋物與實地情況稍微有些不同外,識別結果與實地地物基本一致。但是,通過蝕變礦物光譜獲取的蝕變信息分布圖中,這 4點處無任何信息。由此說明:這些點的蝕變巖信息較其表面蝕變礦物信息更豐富,因而能被有效識別出。

(2)在蝕變巖信息圖中,通過比較吸收谷位置識別前后的結果,發現:只進行光譜整體形態識別的結果中一些區域(如圖4b中A1、A2)存在假的蝕變巖信息。因此,為保證識別的精確性,避免一些假的蝕變巖信息,光譜匹配除從整體形態上匹配外,還應從局部吸收谷形態上進行匹配。

(3)通過比較提取的蝕變巖、粘土礦物、鐵礦物的信息圖,發現:鐵礦物的信息分布與蝕變巖信息分布較接近,主要存在于侵入巖及其周圍,可能是因為鐵礦物主要是熱液蝕變作用的產物或其風化氧化后產物,而蝕變巖也是熱液蝕變作用的產物;在粘土信息分布圖中,硬石膏和高嶺石是主要的礦物,但是大量信息分布在沉積巖區。這些說明:通過蝕變巖光譜獲取的蝕變信息比通過蝕變礦物光譜獲取的蝕變信息更可靠。

甘甫平,王潤生,楊蘇明.2002.西藏Hyperion數據蝕變礦物識別初步研究[J].國土資源遙感,(4):44-51.

高建陽.2011.Hyperion高光譜數據在福建鐘騰銅鉬礦區的應用研究[J].國土資源遙感,(1):87-90.

郭娜,郭科,張婷婷,劉廷晗,胡斌,汪重午.2012.基于短波紅外勘查技術的西藏甲瑪銅多金屬礦熱液蝕變礦物分布模型研究[J].地球學報,33(4):641-653.

闞明哲,田慶久,張宗貴.2005.新疆哈密三種典型蝕變礦物的HyMap高光譜遙感信息提取[J].國土資源遙感,(1):37-40.

劉圣偉,甘甫平,閆柏琨,楊蘇明,王潤生,王青華,唐攀科.2006.成像光譜技術在典型蝕變礦物識別和填圖中的應用[J].中國地質,33(1):178-186.

王潤生,熊盛青,聶洪峰,梁樹能,齊澤榮,楊金中,閆柏琨,趙福岳,范景輝,童立強,林鍵,甘甫平,陳微,楊蘇明,張瑞江,葛大慶,張曉坤,張振華,王品清,郭小方,李麗.2011.遙感地質勘查技術與應用研究[J].地質學報,85(11):1700-1743.

BAUGH W M,KRUSE F A,JR W W A.1998.Quantitative Geochemical Mapping of Ammonium Minerals in the Southern Ceder Mountains,Nevada,Using the Airborne Visible/Infrared Imaging Spectrometer(AVIRIS)[J].Remote Sensing of Environment,65(3):292-308.

BEDINI E.2011.Mineral mapping in the Kap Simpson complex,central East Greenland,using HyMap and ASTER remote sensing data[J].Advances in Space Research,47(1):60-73.

BISHOP C A,LIU Jian-guo,MASON P J.2011.Hyperspectral remote sensing for mineral exploration in Pulang,Yunnan Province,China[J].International Journal of Remote Sensing,32:2409-2426.

CHEN Xiao-feng,WARNER T A,CAMPAGNA D J.2007.Integrating visible,near-infrared and short wave infrared hyperspectral and multispectral thermal imagery for geological mapping at Cuprite,Nevada[J].Remote Sensing of Environment,110:344-356.

CRO′STA A P,SABINE C,TARANIK J V.1998.Hydrothermal Alteration Mapping at Bodie,California,Using AVIRIS Hyperspectral Data[J].Remote Sensing of Environment,65:309-319.

DEBBA P,VAN RUITENBEEK F J A,VAN DER MEER F D,CARRANZA E J M,STEIN A.2005.Optimal field sampling for targeting minerals using hyperspectral data[J].Remote Sensing of Environment,99:373-386.

DENNISON P E,HALLIGAN K Q,ROBERTS D A.2004.A comparison of error metrics and constraints for multiple endmember spectral mixture analysis and spectral angle mapper[J].Remote Sensing of Environment,93:359-367.

GAN Fu-ping,WANG Run-sheng,YANG Su-ming.2002.Study on the alteration minerals identification using Hyperion data [J].Remote Sensing for Land and Resources,(4):44-51(in Chinese with English abstract).

GAO Jian-yang.2011.The application of the Hyperion hyperspectral image to the Zhongteng Cu-Mo deposit in Pinghe County of Fujing province[J].Remote Sensing for Land and Resources,(1):87-90(in Chinese with English abstract).

GERSMAN R,BEN-DOR E,BEYTH M,AVIGAD D,ABRAHA M,KIBREAB A.2008.Mapping of hydrothermally altered rocks by the EO-1 Hyperion sensor,Northern Danakil Depression,Eritrea[J].International Journal of Remote Sensing,29:3911-3936.

GUO Na,GUO Ke,ZHANG Ting-ting,LIU Ting-han,HU Bin,WANG Chong-wu.2012.Hydrothermal Alteration Distribution Model of the Jiama (Gyama) Copper-Polymetallic Deposit Based on Shortwave Technique[J].Acta Geoscientica Sinica,33(4):641-653(in Chinese with English abstract).

KAN Ming-zhe,TIAN Qing-jiu,ZHANG Zong-gui.2005.The extraction of HyMap hyperspectral remote sensing information from three typical altered minerals in Hami area,Xingjing[J].Remote Sensing for Land and Resources,(1):37-40(in Chinese with English abstract).

KRUSE F A.2012.Mapping surface mineralogy using imaging spectrometry[J].Geomorpholgy,137:41-56.

KUMARA A S,KEERTHI V,MANJUNATH A S,VAN DER WERFF H,VAN DER MEER F.2010.Hyperspectral image classification by a variable interval spectral average and spectral curve matching combined algorithm[J].International Journal of Applied Earth Observation and Geoinformation,12:261-269.

LIU Sheng-wei,GAN Fu-ping,YAN Bo-kun,YANG Su-ming,WANG Run-sheng,WANG Qing-hua,TANG Pan-ke.2006.Application of the imaging spectroscopic technique in mineral identification and mapping[J].Geology in China,33(1):178-186(in Chinese with English abstract).

MOUNTRAKIS G,IM J,OGOLE C.2011.Support vector machines in remote sensing:A review[J].ISPRS Journal of Photogrammetry and Remote Sensing,66:247-259.

POUR A B,HASHIM M,VAN G J.2013.Detection of hydrothermal alteration zones in a tropical region using satellite remote sensing data:Bau goldfield,Sarawak,Malaysia[J].Ore Geology Reviews,54:181-196.

ROWAN L C,CROWLEY J K,SCHMIDT R G,AGER C M,MARS J C.2000.Mapping hydrothermally altered rocks by analyzing hyperspectral image (AVIRIS) data of forested areas in the Southeastern United States[J].Journal of Geochemical Exploration,68(3):145-166.

VAN DER MEER F,VAN DER WERFF H,VAN RUITENBEEK F,HECKER C,BAKKER W,NOOMEN M,VAN DER MEIJDE M,CARRANZA E.2012.Multi- and hyperspectral geologic remote sensing:A review[J].International Journal of Applied Earth Observation and Geoinformation,14(1):112-128.

VAN DER MEER F.2004.Analysis of spectral absorption features in hyperspectral imagery[J].International Journal of Applied Earth Observation and Geoinformation,5(1):55-68.

VAN DER MEER F.2006a.Indicator kriging applied to absorption band analysis in hyperspectral imagery:a case study from the Rodalquilar epithermal gold mining area,SE Spain[J].International Journal of Applied Earth Observation and Geoinformation,8:61-72.

VAN DER MEER F.2006b.The effectiveness of spectral similarity measures for the analysis of hyperspectral imagery[J].International Journal of Applied Earth Obser-vation and Geoinformation,8:3-17.

WANG Run-sheng,XIONG Sheng-qing,NIE Hong-feng,LIANG Shu-neng,QI Ze-rong,YANG Jin-zhong,YAN Bai-kun,ZHAO Fu-yue,FAN Jing-hui,TONG Li-qiang,LIN Jian,GAN Fu-ping,CHEN Wei,YANG Su-ming,ZHANG Rui-jiang,GE Da-qing,ZHANG Xiao-kun,ZHANG Zhen-hua,WANG Pin-qing,GUO Xiao-fang,LI Li.2011.Remote Sensing Technology and Its Application in Geological Exploration[J].Acta Geologica Sinica,85(11):1700-1743(in Chinese with English abstract).

XU Yuan-jin,ZHANG Zhen-fei,HU Guang-dao.2010.Ground Object Identification Based on Absorption-band Position Using EO-1 Hyperion Data[J].Photonirvachak-Journal of the Indian Society of Remote Sensing,38:299-308.

猜你喜歡
變巖礦物光譜
基于三維Saab變換的高光譜圖像壓縮方法
煤泥水中煤與不同礦物相互作用的模擬研究
我國首列106節重載列車抵達濟礦物流
楊房溝水電站左岸拱肩槽f27斷層蝕變巖特性及工程實踐
金廠峪蝕變巖制備礦物釉的微觀特征研究
某水電站蝕變巖遇水強度軟化三軸壓縮試驗研究
基于NAIRS和PCA-SVM算法快速鑒別4種含鐵礦物藥
江西相山礦田主要鈾礦化類型及其地球化學特征對比研究
星載近紅外高光譜CO2遙感進展
苦味酸與牛血清蛋白相互作用的光譜研究
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合