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

?

生物天敵暴發(fā)導(dǎo)致珊瑚礁退化的高分遙感監(jiān)測與分析
——以南海太平島為例

2023-11-13 09:25:56鄭金輝任廣波胡亞斌張飛飛李明杰王瑞富
熱帶地理 2023年10期
關(guān)鍵詞:太平島珊瑚礁黑皮

鄭金輝,任廣波,胡亞斌,張飛飛,馬 毅,李明杰,王瑞富

(1.山東科技大學(xué) 測繪與空間信息學(xué)院,山東 青島 266590;2.自然資源部海洋環(huán)境探測技術(shù)與應(yīng)用重點實驗室,廣州 510300;3.自然資源部第一海洋研究所,山東 青島 266061;4.自然資源部海洋遙測技術(shù)創(chuàng)新中心,山東 青島 266061;5.自然資源部南海發(fā)展研究院,廣州 510310;6.海南南沙珊瑚礁生態(tài)系統(tǒng)國家野外科學(xué)觀測研究站,廣州 510399)

珊瑚礁主要是以造礁石珊瑚石灰質(zhì)骨骼為主干,附著珊瑚藻、仙掌藻、軟體動物殼、有孔蟲等鈣質(zhì)生物的巖石體(余克服 等,2018)。珊瑚礁是海底及整個地球的重要生態(tài)系統(tǒng),具有極高的初級生產(chǎn)力(黃榮永 等,2019),其面積僅占海洋面積的0.25%,但卻是25%的海洋魚類和其他生物的棲息地,具有巨大的生物資源和社會價值。

20 世紀60 年代以來,黑皮海綿、長棘海星、食珊瑚蝸牛(Julianna et al., 2022)、藤壺(Jennie et al., 2016)等珊瑚生物天敵暴發(fā)導(dǎo)致主要珊瑚分布區(qū)大面積受損,由此造成海洋生物多樣性降低、珊瑚礁生態(tài)系統(tǒng)退化等一系列環(huán)境問題(李元超 等,2019)。黑皮海綿和長棘海星是2種最常見的珊瑚生物天敵,其中,黑皮海綿以間接包裹侵入方式破壞珊瑚,1971 年在關(guān)島暴發(fā)后造成該地區(qū)30%~80%珊瑚礁死亡(Brayn et al., 1973),隨后在印度洋-太平洋地區(qū)的多個地點也發(fā)現(xiàn)有黑皮海綿的暴發(fā)現(xiàn)象,造成該地區(qū)20%~90%的珊瑚礁損毀(Diraviya et al., 2018; Masashi et al., 2017; Liao et al., 2007);而長棘海星以直接啃噬方式破壞珊瑚,導(dǎo)致活珊瑚覆蓋率①活珊瑚覆蓋率:指珊瑚生長覆蓋的面積占島礁總體面積的百分比。急劇下降,對珊瑚礁生態(tài)系統(tǒng)結(jié)構(gòu)組成和功能將產(chǎn)生長期影響,如1985—2012年大堡礁一項有關(guān)活珊瑚覆蓋率的研究發(fā)現(xiàn),在長棘海星暴發(fā)后活珊瑚覆蓋面積下降達50%(Glenn et al., 2012);2002—2010 年巴布亞新幾內(nèi)亞的布特萊斯灣(Bootless Bay)的實地調(diào)查發(fā)現(xiàn),在長棘海星暴發(fā)事件后超過55%的活珊瑚死亡,珊瑚覆蓋率從2005年的42.4%下降到2006 年的19.1%(Pratchett et al.,2009)。綜上,珊瑚生物天敵暴發(fā)具有擴張速度快、對珊瑚礁損毀程度高的特點(Mohsen et al., 2012),在活珊瑚覆蓋率大幅下降的區(qū)域,珊瑚砂、珊瑚碎屑等底質(zhì)急劇增多,可通過長時序、高頻次監(jiān)測高分遙感影像中珊瑚礁地貌類型的變化,來評估珊瑚礁退化現(xiàn)象,進而開展珊瑚礁健康狀態(tài)監(jiān)測。

高空間分辨率遙感影像具有長時間、大范圍區(qū)域覆蓋的能力,能精細表征珊瑚礁地貌的空間、結(jié)構(gòu)等信息,在珊瑚礁地貌類型提取和健康狀況監(jiān)測等方面具有應(yīng)用價值(John et al., 2018; Xu et al.,2021)。Zuo (2017)、董 娟(2020)、張 飛 飛(2023)等基于覆蓋南海西沙群島的GF-2、World-View-2 高分辨率遙感影像開展珊瑚礁地貌類型分類,證明了利用遙感影像研究珊瑚礁的可行性,為相關(guān)研究提供了新方法。在此基礎(chǔ)上,可進一步開展珊瑚礁健康狀況研究,如逄今朝等(2021)基于長時序QuickBird、WorldView-2多源遙感影像,結(jié)合機器學(xué)習(xí)開展了中國西沙群島中永樂群島的珊瑚礁底質(zhì)類型分類,并評定了研究區(qū)的珊瑚礁白化等級,研究發(fā)現(xiàn)總體時段內(nèi)14 個島礁中有11 個達到重度白化及以上等級,導(dǎo)致珊瑚礁退化的因素有厄爾尼諾現(xiàn)象、人為活動等;La 等(2022)使用Landsat TM、Landsat ETM+、Landsat 8多源遙感影像進行了Tiworo海峽保護區(qū)珊瑚礁退化研究,結(jié)果發(fā)現(xiàn)活珊瑚覆蓋面積從1994 年的78.30 hm2減少到2019 年的8.01 hm2,珊瑚礁退化速率為2.81 hm2/a,該研究綜合評估了經(jīng)濟、生態(tài)、社會、法律、技術(shù)、制度等因素對珊瑚礁的影響,認為生態(tài)因素為主要致災(zāi)因素。上述研究對珊瑚礁退化因素進行了較為系統(tǒng)的探討,但多為海表面升溫、人類活動影響等綜合因素,對珊瑚生物天敵暴發(fā)導(dǎo)致短期珊瑚礁地貌類型退化探討較少。

2017 年,南海太平島珊瑚礁受到黑皮海綿入侵,最密集處黑皮海綿的覆蓋度達27.4%(Yang et al., 2018);2021 年,長棘海星持續(xù)暴發(fā)(Konstantin et al., 2022),最高密度可達1 920 個/hm2,遠遠高于15 個/hm2的暴發(fā)閾值,導(dǎo)致活珊瑚覆蓋率從2017 年的33%降低到2021 年的0.9%(Wei et al.,2022)。在短短5年內(nèi)發(fā)生2次不同種類的珊瑚生物天敵暴發(fā)事件,嚴重威脅了太平島的珊瑚礁生態(tài)健康。上述相關(guān)學(xué)者對太平島黑皮海綿、長棘海星2次暴發(fā)事件進行的調(diào)查研究,側(cè)重于單一時間節(jié)點的活珊瑚覆蓋率狀況監(jiān)測,未對持續(xù)的生物天敵暴發(fā)后珊瑚礁地貌變化過程開展詳細分析。而通過開展長時序珊瑚礁地貌類型監(jiān)測,可反映珊瑚礁覆蓋度整體變化,進而評估珊瑚礁生態(tài)系統(tǒng)健康狀況。因此,本文針對太平島2016—2022年黑皮海綿和長棘海星暴發(fā)事件,選取覆蓋太平島的26期Sentinel-2 遙感影像,結(jié)合GF-2(PMS)和Google Earth 平臺的高分辨率遙感影像數(shù)據(jù),開展生物天敵暴發(fā)前后珊瑚礁地貌類型分類實驗,從珊瑚礁地貌的類型變化、面積變化、退化率與恢復(fù)率等角度分析珊瑚礁地貌類型變遷特征,評估生物天敵暴發(fā)事件對太平島珊瑚礁生態(tài)系統(tǒng)退化的影響。以期為太平島海域珊瑚礁生態(tài)資源保護與研究提供理論依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

中國南海太平島位于南沙群島西北部,具有重要的地理位置和戰(zhàn)略意義。本研究區(qū)為太平島及其周邊礁盤(圖1),其地理坐標為10°22′38″ N、114°21′59″ E。太平島是南沙群島中最大的、擁有淡水資源的天然島嶼,出露水面的島嶼由東至西長約1 289.3 m,南北寬約365.7 m;海岸線長度為3 470 m,總面積為360 870 m2;平均海拔3.8 m,東部最高達4.18 m,低潮時可顯露出環(huán)礁寬約30 m的珊瑚帶(夏小明 等,2012);水下礁盤大而連續(xù),平均水溫為28.98℃,其上孕育著豐富的珊瑚礁生態(tài)系統(tǒng)。據(jù)報道,2018 年太平島活珊瑚覆蓋率在60%~80%,包括至少9種大型無脊椎動物、35種藻類和46種珊瑚礁②資料來源:https://www.agriharvest.tw/archives/20987#:~:text=%E5%8D%97%E6%B2%99%E5%A4%AA%E5%B9%B3%E5%B3%B6%E 5%91%A8%E9%82%8A%E6%B5%B7,UCN%E6%98%93%E5%8D%B1%E7%89%A9%E7%A8%AE。有學(xué)者曾對南沙群島55個島礁從自然、經(jīng)濟、航運、政治、軍事、災(zāi)害等6個方面的綜合價值進行研究,得出其綜合價值評估分值排第四位(李弘毅 等,2018)。

圖1 南沙太平島區(qū)位和遙感影像示意Fig.1 Schematic diagram of the location and remote sensing images of the Taiping Island in Nansha

1.2 數(shù)據(jù)與處理

本文所用遙感數(shù)據(jù)有空間分辨率為10 m的Sentinel-2 影 像、4 m GF-2 (PMS) 影 像 和1 m 的Google Earth 平臺的影像數(shù)據(jù),不同年份對應(yīng)影像數(shù)量分布情況見圖2,影像詳細信息見表1。

表1 遙感影像信息Table 1 Remote sensing image information

圖2 不同遙感影像數(shù)量與成像時間統(tǒng)計Fig.2 Statistics of the number and imaging time of different remote sensing images

分類實驗數(shù)據(jù)為Sentinel-2 遙感數(shù)據(jù),影像來源為歐空局哥白尼數(shù)據(jù)中心③https://scihub.copernicus.eu/dhus/#/home。Sentinel-2 影像具有較高的時間分辨率和空間分辨率,Sentinel-2A 與Sentinel-2B兩顆衛(wèi)星互補使得重訪周期為5 d,空間分辨率為10 m,在這種時間分辨率與空間分辨率的配合下,可以刻畫出較精細的珊瑚礁地貌類型,適用于面向小尺度上具有高度異質(zhì)性的珊瑚礁監(jiān)測。Sentinel-2 L1C(Level-1C)級數(shù)據(jù)是經(jīng)過幾何校正的大氣表觀反射率產(chǎn)品,Sentinel-2 L2A(Level-2A)級數(shù)據(jù)是經(jīng)過輻射定標、大氣校正和幾何校正的標準數(shù)據(jù)產(chǎn)品,對L1C 數(shù)據(jù)進行輻射定標和大氣校正處理為L2A 級數(shù)據(jù)后,所得影像數(shù)據(jù)可滿足本研究需求。Sentinel-2 遙感影像的波段2(中心波長為0.490 nm的藍光波段)、3(中心波長為0.560 nm 的綠光波段)可較好穿透水體,表征水下珊瑚礁地貌信息;4(中心波長為0.665 nm的紅光波段)、8(中心波長為0.842 nm的近紅外波段)可有效監(jiān)測出露水面的珊瑚礁地物信息,因此選用這4 個波段開展2016—2022 年珊瑚礁地貌類型長時間序列變化監(jiān)測。

還采用GF-2(PMS)與Google Earth 平臺的影像數(shù)據(jù)作為驗證數(shù)據(jù)的源數(shù)據(jù),其中,GF-2(PMS)衛(wèi)星影像來源于中國海洋衛(wèi)星數(shù)據(jù)服務(wù)系統(tǒng),Google Earth 平臺的影像數(shù)據(jù)來源于https://earth.google.com,二者空間分辨率分別為4和1 m,優(yōu)于Sentinel-2 影像。結(jié)合專家解譯知識和GF-2(PMS)、Google Earth 平臺的高分辨率影像數(shù)據(jù),獲取的珊瑚礁地貌類型專家解譯結(jié)果作為驗證數(shù)據(jù)。

利用ENVI 軟件中的Radiometric Correction、FLAASH Atmospheric Correction 工 具 對 GF-2(PMS)數(shù)據(jù)進行輻射定標、大氣校正,利用Registration 工具以Google Earth 平臺的影像數(shù)據(jù)為基準對GF-2(PMS)數(shù)據(jù)進行幾何校正,配準誤差<0.5 像元。最后對太平島礁體外的區(qū)域進行掩膜,完成影像預(yù)處理。

2 研究方法

選用2016—2022 年覆蓋太平島的26 期Sentinel-2遙感影像數(shù)據(jù)作為實驗分類數(shù)據(jù),基于構(gòu)建的珊瑚礁地貌類型分類體系,采用SVM 分類方法開展珊瑚礁地貌類型分類實驗,結(jié)合11 期GF-2 及Google Earth 平臺的高分辨率影像獲取的專家解譯結(jié)果,進行珊瑚礁地貌類型分類修正后結(jié)果精度評價;在長時間序列珊瑚礁地貌類型分類實驗結(jié)果的基礎(chǔ)上,從珊瑚礁地貌類型變遷、轉(zhuǎn)移面積、珊瑚礁退化率與恢復(fù)率等角度開展黑皮海綿和長棘海星暴發(fā)影響下的珊瑚礁地貌類型特征變化分析。技術(shù)流程如圖3所示。

圖3 技術(shù)流程Fig.3 Technical flowchart

2.1 太平島珊瑚礁地貌類型分類體系構(gòu)建

珊瑚礁地貌類型的分類體系眾多,Sarah 等(2016)從物理指標如沉積物類型和海水動力成因角度,將其劃分為砂質(zhì)前礁、活珊瑚、固結(jié)碎石、藻平面等9 類;程益鋒等(2018)則從物理指標結(jié)合生物指標如生物生長狀況角度將其劃分為潟湖、礁坪、灰沙島、礁脊等6個一級類及潟湖盆、點礁、內(nèi)礁坪、沙灘等10個二級類;董娟等(2020)將活珊瑚覆蓋率納入分類標準,劃分地貌類型為礁坪、向海坡、潟湖、灰沙島等4 個一級類和礁前階地、礁脊、珊瑚叢生區(qū)、珊瑚沉積區(qū)等11個二級類。在以上學(xué)者分類體系基礎(chǔ)上,本文結(jié)合太平島實際情況構(gòu)建太平島珊瑚礁地貌類型分類體系(表2),包括深礁前斜坡、淺礁前斜坡、珊瑚叢生區(qū)、稀疏珊瑚沉積區(qū)、密集珊瑚沉積區(qū)、沙坪、陸地等7 類,不同地貌類型示意見圖4所示。由于黑皮海綿與長棘海星暴發(fā)均會造成珊瑚礁退化(Konstantin et al.,2022),因此通過分析遙感影像中不同珊瑚礁地貌類型的變化,可進行生物天敵暴發(fā)下的珊瑚礁生態(tài)系統(tǒng)退化監(jiān)測研究。

表2 太平島珊瑚礁地貌類型分類解譯標志Table 2 Classification and interpretation signs of coral reef geomorphology types of the Taiping Island

圖4 珊瑚礁地貌類型示意(以Sentinel-2影像為例)Fig.4 Schematic diagram of coral reef landform types(Taking Sentinel-2 image as an example)注:Sentinel-2影像波段組合方式:3(紅)、2(綠)、1(藍)。

結(jié)合董娟(2021)活珊瑚覆蓋率與珊瑚礁地貌類型的關(guān)系,在本文構(gòu)建的珊瑚礁地貌類型分類體系中,確定活珊瑚覆蓋率由高到低的地貌類型順序為:珊瑚礁叢生區(qū)>淺礁前斜坡>密集珊瑚沉積區(qū)>稀疏珊瑚沉積區(qū)>深礁前斜坡>沙坪=陸地。

2.2 基于SVM的珊瑚礁地貌類型分類方法

2.2.1 分類方法 確定研究區(qū)珊瑚礁地貌不同類別后,對預(yù)處理后的遙感影像選取不同類別影像區(qū)域制作樣本,采用ENVI軟件的SVM模塊進行分類實驗。SVM是一種結(jié)合核函數(shù)和優(yōu)化理論的線性分類器,通過訓(xùn)練算法建立模型,將選取的訓(xùn)練樣本表示為空間上的點,根據(jù)映射函數(shù)將不同樣本劃分出盡可能寬的分類平面,測試樣本在映射到相同空間后,便可在訓(xùn)練好模型的基礎(chǔ)上預(yù)測所屬類別,從而對影像提供多種類別分類(Nyan et al., 2021)。其公式可表達為:

式中:{xi,yj} 為訓(xùn)練樣本,其中xi為向量,yi屬于{ - 1, + 1};Ke為核函數(shù),又稱映射函數(shù),可將輸入樣本映射為更高維樣本;α為拉格朗日系數(shù),其作用是將目標函數(shù)與對應(yīng)的限制條件整合為一個函數(shù),αi、αj為對應(yīng)約束條件的系數(shù);C為懲罰系數(shù),其大小影響算法的正則化程度。該公式表示對于給定的一組樣本{xi,yj}和對應(yīng)Ke、C,每個yi={ - 1, + 1}被分為兩種類別之一。

該分類方法主要優(yōu)勢為:1)采用非線性算法將線性不可分樣本轉(zhuǎn)化至高維特征空間,重新構(gòu)建最優(yōu)分割超平面;2)SVM 更適合中小樣本數(shù)據(jù)量分類,能在樣本量較小的情況下更快達到較高精度(Liu et al., 2017)。太平島研究區(qū)珊瑚礁地貌類型的面積分布,決定其分類樣本屬于不平衡樣本,且不同類型區(qū)域在交界處存在相近的顏色、紋理特征,由此采用SVM分類方法,在保證分類精度和分類速度的同時,與專家解譯結(jié)合從而較好完成分類。對Sentinel-2遙感影像進行SVM分類后,利用先驗地學(xué)分類知識,即影像中珊瑚礁地貌類型的空間分層特性、顏色、紋理、位置等解譯標志,對分類結(jié)果中明顯錯分的圖斑進行修正與優(yōu)化。

2.2.2 精度評價方法 主要評價指標包括總體精度(Overall Accuracy, OA)、Kappa 系數(shù),其中OA 評估分類正確的樣本數(shù)與總樣本數(shù)量的比例,OA值越高,表示分類模型的整體性能越好。Kappa系數(shù)是一種綜合評估分類模型分類準確性和隨機分類差異的指標,該值越接近1,表示模型的分類準確性越高;越接近0,表示模型的分類結(jié)果接近隨機分類;越接近-1,表示模型的分類準確性越低。通常,Kappa系數(shù)>0.8 被認為是較好的分類結(jié)果,0.6~0.8 為中等,<0.6 表示分類結(jié)果較差。

2.3 珊瑚礁地貌類型時空變化監(jiān)測方法

對太平島珊瑚礁地貌類型變化監(jiān)測主要方式為統(tǒng)計2016―2022年各地貌類型面積,結(jié)合馬爾科夫轉(zhuǎn)移矩陣(transition matrix by Markov)計算不同地貌類型間相互轉(zhuǎn)移量。轉(zhuǎn)移矩陣各行元素均為非負數(shù),表示不同類型的轉(zhuǎn)移狀態(tài),在一定條件下是相互轉(zhuǎn)移的,其本質(zhì)特征不僅可反映單個地貌類型時空變化情況,還可反映不同類型之間隨時間變化的相互轉(zhuǎn)化量,從而進行更深入的珊瑚礁地貌類型變化定量分析(徐新良 等,2014)。轉(zhuǎn)移矩陣原理如圖5所示。

圖5 轉(zhuǎn)移矩陣原理(a.珊瑚礁地貌類型轉(zhuǎn)移矩陣;b.珊瑚礁地貌類型轉(zhuǎn)移量示意)Fig.5 Schematic diagram of transfer matrix principle (a.Coral reef landform type transfer matrix; b.Schematic diagram of coral reef landform type transfer)

2.4 退化與恢復(fù)定量計算方法

利用轉(zhuǎn)移矩陣分析不同珊瑚礁地貌類型變化情況,再依據(jù)逄今朝等(2021)的方法,將地貌類型由珊瑚礁覆蓋率④珊瑚礁覆蓋率:指有活珊瑚生長覆蓋的地貌類型面積占島礁總體面積的百分比。高轉(zhuǎn)化為低的過程定義為珊瑚礁退化,退化面積為DA(Degradation Area),反之為珊瑚礁恢復(fù),恢復(fù)面積為RA(Recovery Area)。珊瑚礁退化率DR(Degradation Rate)計算公式為:

珊瑚礁恢復(fù)率(Recovery Rate, RR)計算公式為:

式中:TA(Total Area)為珊瑚礁地貌類型總面積。

3 結(jié)果與分析

3.1 分類結(jié)果精度分析與評價

采用專家解譯知識和高分辨率遙感影像相結(jié)合的方式,針對26 期Sentinel-2 遙感影像中7 種不同的珊瑚礁地貌類型開展SVM 分類實驗。為保證樣本數(shù)量,提高分類精度,依據(jù)每種珊瑚礁地貌類型的解譯特征,最終選擇訓(xùn)練樣區(qū)3 640 個,測試樣區(qū)1 560個,每個樣區(qū)包含3~6個像元。

不同的珊瑚礁地貌類型在遙感影像中位置、顏色、紋理均呈現(xiàn)明顯差異,如淺礁前斜坡與深礁前斜坡相鄰,位于內(nèi)環(huán);沙坪在影像中反射率、亮度較高且位于礁盤之上;稀疏珊瑚沉積區(qū)、密集珊瑚沉積區(qū)、珊瑚叢生區(qū)的亮度依次降低,紋理明顯;陸地位于島礁上中心位置,覆有植被,特征明顯,故不同珊瑚礁地貌類型可分性較高。

實驗處理過程見圖6。選用2 種珊瑚礁生物天敵暴發(fā)前后6 期影像(2016-11-20、2017-06-28、2018-08-27、 2019-08-12、 2020-06-22、 2021-07-12)的分類結(jié)果為例,結(jié)果見圖7。由圖7可知,各珊瑚礁地貌類型分類精度存在較為明顯的差異性,深礁前斜坡、陸地受珊瑚礁變化影響最小,且解譯時區(qū)分度較高,故精度穩(wěn)定且較高;而珊瑚叢生區(qū)、沙坪(部分時期)面積占比較小,樣本選擇時像元數(shù)量少,屬于小樣本分類,影響SVM 分類精度。

圖6 太平島珊瑚礁地貌類型分類過程(以2017-06-28所選影像為例)Fig.6 Classification process of coral reef landform types in Taiping Island (taking the image selected on June 28th, 2017 as an example)

圖7 6景影像分類后珊瑚礁地貌類型精度Fig.7 Accuracy of coral reef landform types after classification of 6 images

為保證對太平島黑皮海綿、長棘海星天敵暴發(fā)后珊瑚礁地貌類型影響的后續(xù)研究,需結(jié)合專家解譯知識,對SVM 分類后的珊瑚礁地貌類型分類結(jié)果進行修正,由此可最大程度去除錯分區(qū)域(圖8)。對該修正后珊瑚礁地貌類型分類結(jié)果進行精度評定,得到6 期珊瑚礁地貌類型總體精度和Kappa系數(shù)均值分別為95.30%和0.89,其中最高分類精度和Kappa系數(shù)達96.46%和0.94(表3)。這證明采用SVM結(jié)合專家解譯進行太平島珊瑚礁地貌類型分類具有較高準確性,可用于監(jiān)測珊瑚生物天敵暴發(fā)造成的珊瑚礁退化現(xiàn)象。

表3 太平島珊瑚礁地貌類型精度評定(專家解譯修正后)Table 3 Precision evaluation table for the landform types of the Taiping Island coral reef (after expert interpretation and correction)

圖8 選取影像分類結(jié)果Fig.8 Selected image classification result graph

3.2 黑皮海綿暴發(fā)下珊瑚礁地貌類型演變特征分析

2017年5月,在南海太平島監(jiān)測到黑皮海綿暴發(fā)(Yang et al., 2018),通過對該地珊瑚礁生物多樣性調(diào)查,得出黑皮海綿的過度擴張覆蓋率接近27.4%。相關(guān)研究證明,黑皮海綿暴發(fā)期通常持續(xù)1月至1 年,最長可達3 年(Chow et al., 2022),因此,基于2016-11-20—2019-12-20 的13 幅影像,使用SVM 開展黑皮海綿暴發(fā)下珊瑚礁退化遙感監(jiān)測與分析。需要指出的是,深礁前斜坡珊瑚覆蓋率極低,實驗對比發(fā)現(xiàn),長期的時間序列影像中其位置、形狀變化不明顯,對后續(xù)珊瑚礁地貌類型監(jiān)測與分析影響較小,故只針對其他6種類型進行變化分析。

由圖9 可知,在2017-06-28 影像中珊瑚礁覆蓋率相對較高的地貌類型受災(zāi)嚴重,稀疏珊瑚沉積區(qū)面積的變化趨勢與密集珊瑚沉積區(qū)并不完全一致,但在之后下一時段的2017-09-16統(tǒng)計中的面積也經(jīng)歷了明顯減少,相比均值(38.26 hm2) 下降了18.10 hm2;密集珊瑚沉積區(qū)面積在2017-06-28影像中顯示有明顯下降,相較于均值(57.64 hm2)下降了23.48 hm2;而沙坪則在2017-06-28 有明顯上升,指示沙坪面積有異常增加,相較均值(8.29 hm2)增加了19.93 hm2;珊瑚叢生區(qū)活珊瑚覆蓋率最高,在2017-06-28 該地貌類型的面積有較明顯的下降,2017-11-20—2018-02-18 珊瑚叢生區(qū)有一增加的峰值隨后回落至均值(1.06 hm2);因陸地周圍圍繞沙坪,故其面積受沙坪影響略有變化;淺礁前斜坡面積無較大變化。

圖9 2016—2019年太平島珊瑚礁地貌類型面積變化Fig.9 Statistics of changes in the area of the Taiping Island coral reef landform types from 2016 to 2019

與2016-11-20初始狀態(tài)相比,2017-6-28后太平島各珊瑚礁地貌類型變化達到峰值:珊瑚叢生區(qū)面積降低64.10%,密集珊瑚沉積區(qū)降低34.43%,稀疏珊瑚沉積區(qū)降低13.71%,即珊瑚礁覆蓋率高的地貌類型都受到較大破壞;而無活珊瑚覆蓋的沙坪面積有激增情況,增幅達3 327.59%。陸地、淺礁前斜坡變化較小,為-1.6%、-0.39%。因此以影像2017-06-28為黑皮海綿暴發(fā)節(jié)點⑤不等于實際暴發(fā)時間,影像獲取時間代表的是這一時間點的狀況信息,而具體一段時間的整體情況及變化趨勢須通過多時相以及多源數(shù)據(jù)才能完整掌握。。

根據(jù)影像分類結(jié)果,得到珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣(表4、5,對應(yīng)轉(zhuǎn)移矩陣變遷情況請見附圖-a)。

表4 2016-11-20—2017-06-28太平島珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣Table 4 Transition matrix of changes in coral reef geomorphic types on the Taiping Island during November 20th, 2016 to June 28th, 2017

由表4 可知,起始監(jiān)測時間2016-11-20 與時間節(jié)點2017-06-28相比,變化較大的珊瑚礁地貌類型集中于稀疏珊瑚沉積區(qū)、沙坪、密集珊瑚沉積區(qū)。密集珊瑚沉積區(qū)向沙坪轉(zhuǎn)移量為12.66 hm2,稀疏珊瑚沉積區(qū)向沙坪轉(zhuǎn)移量12.97 hm2,而沙坪向密集/稀疏珊瑚沉積區(qū)轉(zhuǎn)移量合計不足1 hm2。除此之外,面積占比較小的珊瑚叢生區(qū)在轉(zhuǎn)移為沙坪量為0.51 hm2,也證實了該時段呈珊瑚覆蓋率高的珊瑚礁地貌類型向低珊瑚覆蓋率高的珊瑚礁地貌類型轉(zhuǎn)移的趨勢。由表5 可知,2017-06-28 至監(jiān)測時間段末2019-12-20,珊瑚礁覆蓋率有所提高,沙坪向其他珊瑚礁地貌類型轉(zhuǎn)移總量為28.22 hm2,其他轉(zhuǎn)為沙坪總量則為10.51 hm2,沙坪整體呈現(xiàn)減少趨勢,珊瑚礁在黑皮海綿暴發(fā)后有逐年恢復(fù)趨勢。

表5 2017-06-28—2019-12-20太平島珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣Table 5 Transition matrix of changes in coral reef geomorphic types on the Taiping Island during June 28th, 2017 to December 20th, 2019

通過統(tǒng)計2016—2019年黑皮海綿暴發(fā)后的珊瑚礁地貌類型變化增減情況(表6)發(fā)現(xiàn),研究區(qū)珊瑚礁地貌類型在黑皮海綿暴發(fā)前后的3—10 月變化顯著,2017-03-20影像中已有密集珊瑚沉積區(qū)和稀疏珊瑚沉積區(qū)劇烈變化。這與Siti 等(2021)的研究相吻合:黑皮海綿擴展速率在暖季(4—10 月)較快,涼季(11至次年3月)較慢,即黑皮海綿在不同月份和不同季節(jié)的擴展速率可能存在差異。

表6 黑皮海綿暴發(fā)研究時段太平島珊瑚礁地貌類型增減統(tǒng)計Table 6 Changes in coral reef geomorphic types on the Taiping Island during the period of Terpios hoshinota sponge outbreak study hm2

以2016-11-20中各珊瑚礁地貌類型面積為基準(見表6)進行分析,可知密集珊瑚沉積區(qū)的減少比例最高,為72.92%;珊瑚礁叢生區(qū)次之,為55.24%;稀疏珊瑚沉積區(qū)第三,為44.99%。

3.3 長棘海星暴發(fā)下珊瑚礁地貌類型演變特征分析

由Wei等(2022)的2021年4月實地考察記錄顯示,太平島上6個站點、共12個樣帶(樣帶長度為50 m)發(fā)現(xiàn)長棘海星分布的平均密度為630 個/hm2,最高可達1 920個/hm2,大規(guī)模暴發(fā)的長棘海星使得珊瑚礁遭到嚴重破壞,該現(xiàn)象引發(fā)的珊瑚礁地貌類型變化在遙感影像上體現(xiàn)更為直觀。在本文實驗中選用2020-02-13—2022-06-07 共13 幅影像作為研究時段二。

由圖10可知,珊瑚礁覆蓋率高的地貌類型面積在2021 年3 月開始有不同程度的變化。圖10-c 所示,稀疏珊瑚沉積區(qū)面積在2021-03-29 有一谷值,面積減少了7.43 hm2,相當于稀疏珊瑚沉積區(qū)在該時段面積均值(41.19 hm2)的36.08%,在2021-04-28有一暫時性的增加,隨后呈下降趨勢;而密集珊瑚沉積區(qū)面積在2021-03-29出現(xiàn)異常升高,但在之后的2021-04-21急劇下降,下降了13.56 hm2,表明該珊瑚礁地貌類型在受到強擾動后面積有所減少,到同年9月份縮減至最低,較均值(50.55 hm2)減少近1/4;沙坪面積在2021-03-29—07-21 呈遞增趨勢,增加了17.16 hm2,接近原面積的3倍;珊瑚叢生區(qū)的面積變化有所滯后,于2021-07-21 達最低,但同年9月25日有一異常升高,與稀疏珊瑚叢生區(qū)的變化一致,這一異??赡苁情L棘海星暴發(fā)后種群聚集,導(dǎo)致影像中產(chǎn)生與珊瑚礁地貌類型相似特征。進一步分析可發(fā)現(xiàn),與2020-02-13初始狀態(tài)相比,2021-06-28后,珊瑚礁覆蓋率較高的密集珊瑚沉積區(qū)下降量最大,為22.46%,但降幅較黑皮海綿暴發(fā)后小;沙坪面積增加最多,達到70.52%。

圖10 2020—2022年太平島珊瑚礁地貌類型面積變化統(tǒng)計Fig.10 Statistics of changes in the area of the Taiping Island coral reef landform types from 2020 to 2022

根據(jù)影像分類結(jié)果,得到珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣(表7、8,對應(yīng)轉(zhuǎn)移矩陣變遷情況圖請見附圖-b)。

表7 2020-02-13—2021-06-27太平島珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣Table 7 Transition matrix of changes in coral reef geomorphic types on the Taiping Island during February 13th, 2020 to June 27th, 2021

表8 2021-06-27—2022-06-07太平島珊瑚礁地貌類型變遷轉(zhuǎn)移矩陣Table 8 Transition matrix of changes in coral reef geomorphic types on the Taiping Island during June 27th, 2021 to June 7th, 2022

由表7、8可知,在長棘海星暴發(fā)影響下,面積最大變化集中于密集珊瑚沉積區(qū)、稀疏珊瑚沉積區(qū)、沙坪。長棘海星暴發(fā)前期(2020-02-13—2021-06-27),密集珊瑚沉積區(qū)轉(zhuǎn)移為稀疏珊瑚沉積區(qū)占比較大,轉(zhuǎn)移量為20.42 hm2,稀疏珊瑚沉積區(qū)向密集珊瑚沉積區(qū)轉(zhuǎn)移量為12.11 hm2;同時稀疏珊瑚沉積區(qū)和密集珊瑚沉積區(qū)向沙坪轉(zhuǎn)移量在長棘海星暴發(fā)后有所減少。長棘海星暴發(fā)后期(2021-06-27—2022-06-07)珊瑚礁覆蓋率高的地貌類型轉(zhuǎn)移為珊瑚礁覆蓋率低的類型面積均未超過10 hm2,且珊瑚礁覆蓋率高的地貌類型面積開始有所增加,如最高增加量為由稀疏珊瑚沉積區(qū)轉(zhuǎn)移為密集珊瑚沉積區(qū)的18.23 hm2。

以2020-02-13中各地貌類型面積為基準(表9)分析發(fā)現(xiàn),珊瑚叢生區(qū)的減少量占原本面積比例最高,達59.17%,稀疏珊瑚沉積區(qū)減少量次之,達34.91%,密集珊瑚沉積區(qū)減少量第三,為23.97%。這與長棘海星的食性高度吻合,珊瑚礁覆蓋率高的地貌類型成為受災(zāi)最為嚴重的區(qū)域。

表9 長棘海星暴發(fā)研究時段太平島各珊瑚礁地貌類型增減統(tǒng)計Table 9 Changes in each coral reef geomorphic type on the Taiping Island during the period of crown-of-thorns starfish outbreak study hm2

綜上可知,黑皮海綿與長棘海星暴發(fā)前后太平島珊瑚礁地貌類型的變化過程中,與活珊瑚覆蓋率呈正相關(guān)的密集珊瑚沉積區(qū)、稀疏珊瑚沉積區(qū)、珊瑚叢生區(qū)在實地調(diào)查報道時間點前后均有異常,以減少為主,黑皮海綿對密集珊瑚沉積區(qū)的破壞程度達72.92%,該珊瑚礁地貌類型退化面積約為37.99 hm2;而長棘海星對珊瑚叢生區(qū)破壞程度達59.17%,該珊瑚礁地貌類型退化面積約為0.71 hm2。這說明2種珊瑚生物天敵暴發(fā)后對太平島的珊瑚礁地貌類型影響顯著,生態(tài)健康狀況有明顯退化趨勢。

4 討論

4.1 珊瑚生物天敵暴發(fā)致不同珊瑚礁地貌類型變化光譜對比

選取太平島Sentinel-2影像中4個包含不同珊瑚礁地貌類型的區(qū)域,進行典型區(qū)域光譜變化對比分析。由于多光譜影像中每個波段光譜標準差可表現(xiàn)像元反射率在空間中的變化,進而反映珊瑚礁健康與退化變化規(guī)律,即珊瑚退化后可見光反射率升高,光譜標準差先波動變化后波形趨于平滑(陳啟東 等,2015),該現(xiàn)象可指示珊瑚礁地貌類型發(fā)生不同變化。由典型區(qū)域影像及反射率標準差變化可知,典型區(qū)Ⅰ~Ⅳ在黑皮海綿暴發(fā)期間有3次較明顯的光譜標準差波動升高,參考線①、②、③分別對應(yīng)2017-06-28、2018-02-18、2019-02-18;長棘海星暴發(fā)期間有2次較明顯的光譜標準差波動升高,參考線④、⑤分別對應(yīng)2021-03-29、2021-09-25(圖11)。

圖11 典型區(qū)域影像及反射率標準差變化Fig.11 Typical regional images and changes in reflectance standard deviation

在起始時間2016-11-20影像中,典型區(qū)Ⅰ主要地貌類型為密集珊瑚沉積區(qū)和稀疏珊瑚沉積區(qū),典型區(qū)Ⅱ包含較多地貌類型:陸地、稀疏珊瑚沉積區(qū)、密集珊瑚沉積區(qū)、珊瑚叢生區(qū)、淺礁前斜坡,典型區(qū)Ⅲ包含大部分密集珊瑚沉積區(qū)和少部分稀疏珊瑚沉積區(qū)、淺礁前斜坡,典型區(qū)Ⅳ全部為密集珊瑚沉積區(qū)。

在黑皮海綿暴發(fā)期,典型區(qū)Ⅰ、Ⅲ在參考線①處有光譜標準差波動升高現(xiàn)象,與起始時間2016-11-20影像中個典型區(qū)的主要地貌類型相比,2個典型區(qū)的稀疏珊瑚沉積區(qū)有所擴大;在參考線②處典型區(qū)Ⅰ、Ⅳ有明顯光譜標準差波動升高現(xiàn)象,典型區(qū)Ⅰ主要地貌類型轉(zhuǎn)移為稀疏珊瑚沉積區(qū)和沙坪,典型區(qū)Ⅳ轉(zhuǎn)移為密集珊瑚沉積區(qū)與稀疏珊瑚沉積區(qū)各占一半,由此說明,參考線①、②對應(yīng)的時間節(jié)點有珊瑚礁退化現(xiàn)象;在參考線③處典型區(qū)Ⅲ、Ⅳ有明顯光譜標準差波動升高現(xiàn)象,典型區(qū)Ⅲ在該時間節(jié)點影像中密集珊瑚沉積區(qū)占主體,典型區(qū)Ⅳ則除密集珊瑚沉積區(qū)與稀疏珊瑚沉積區(qū)外,還新增了沙坪,說明該時間節(jié)點不同典型區(qū)珊瑚礁恢復(fù)與退化同時發(fā)生。

在長棘海星暴發(fā)期,典型區(qū)Ⅰ、Ⅱ、Ⅲ在參考線④處均有光譜標準差波動變化,典型區(qū)Ⅰ轉(zhuǎn)移為起始類型(位置有所變化),表明珊瑚礁有所恢復(fù),典型區(qū)Ⅱ多了沙坪,典型區(qū)Ⅲ為稀疏珊瑚沉積區(qū)占主體,表明珊瑚礁有所退化;典型區(qū)Ⅱ、Ⅲ在參考線⑤處均達到峰值,典型區(qū)Ⅱ沙坪有所擴大,典型區(qū)Ⅲ不含有珊瑚叢生區(qū)、沙坪地貌類型,表明珊瑚礁既有退化也有恢復(fù)現(xiàn)象。

以上地貌類型轉(zhuǎn)移指示參考線①、②、④處的時間段珊瑚礁有退化現(xiàn)象,參考線③、⑤處對應(yīng)珊瑚礁既有退化也有恢復(fù)現(xiàn)象。由此證明黑皮海綿與長棘海星暴發(fā)前期珊瑚礁以退化為主,而后期珊瑚礁既有退化也有恢復(fù)趨勢。

4.2 珊瑚礁退化與恢復(fù)分析

基于太平島SVM 分類實驗珊瑚礁地貌類型特征演變結(jié)果,進一步分析南海太平島珊瑚礁地貌類型總體退化與恢復(fù)率(圖12)。

圖12 珊瑚生物天敵暴發(fā)期間太平島珊瑚礁退化率與恢復(fù)率統(tǒng)計Fig.12 Statistical graph of coral reef degradation and recovery rates on the Taiping Island during the outbreak period of coral biological predators

由圖12 可知,2016—2022 年太平島珊瑚礁總體退化率高于恢復(fù)率,平均珊瑚礁退化率為16.18%,平均珊瑚礁恢復(fù)率為11.94%;其中2016—2019年黑皮海綿暴發(fā)期間,太平島的珊瑚礁退化率最高達到23.76%,在峰值之后珊瑚礁退化率有所下降,而珊瑚礁恢復(fù)率由較低水平升至22.48%。這說明黑皮海綿在暴發(fā)后接近2年時間內(nèi),對珊瑚的影響有所減弱,或是珊瑚的適應(yīng)性有所增強,部分珊瑚從退化中恢復(fù)。在2020—2021年長棘海星暴發(fā)前期的珊瑚礁退化率為最高為17.32%,珊瑚礁恢復(fù)率與黑皮海綿暴發(fā)前期相近;在2021—2022長棘海星暴發(fā)后期珊瑚礁退化率最高達21.25%,珊瑚礁恢復(fù)率有所上升但不及黑皮海綿暴發(fā)后期。黑皮海綿暴發(fā)事件后,長棘海星暴發(fā),此時的珊瑚礁雖有所恢復(fù)但其群落結(jié)構(gòu)發(fā)生較大變化,供給長棘海星的珊瑚礁食物已大幅減少。在2016—2022年,太平島珊瑚礁總體退化率高于恢復(fù)率,其中2017 年3—6月珊瑚礁退化率最高,為23.88%;在2017 年6—9月珊瑚礁恢復(fù)率最高,為18.03%。

綜上,可以發(fā)現(xiàn)不同生物天敵暴發(fā)對珊瑚的影響存在差異。黑皮海綿暴發(fā)后珊瑚礁退化率較高,對珊瑚礁地貌類型的影響程度較大,而長棘海星暴發(fā)后珊瑚礁退化率和恢復(fù)率均較低,對珊瑚礁地貌類型的影響程度相對較小。因此,珊瑚在面對不同暴發(fā)天敵種類時表現(xiàn)出不同的適應(yīng)性和反應(yīng)機制,這些差異可能與珊瑚的生物學(xué)特性、暴發(fā)天敵的生態(tài)位、以及珊瑚和暴發(fā)天敵之間的相互作用等因素有關(guān)。

由圖13可知,太平島沙坪經(jīng)歷珊瑚生物天敵暴發(fā)事件后面積顯著增加,由總面積的0.4%增至4.5%;稀疏珊瑚沉積區(qū)面積有較明顯減少,由23.6% 降至17.3%,密集珊瑚沉積區(qū)有微小增加,而珊瑚叢生區(qū)、陸地和淺礁前斜坡變化較小。

圖13 2016—2022年太平島珊瑚礁地貌類型占比統(tǒng)計Fig.13 Statistical analysis of coral reef geomorphic type proportions on the Taiping Island from 2016 to 2022

4.3 其他潛在因素對2016―2022 年太平島珊瑚礁退化影響分析

珊瑚礁生態(tài)系統(tǒng)健康的其他潛在因素主要有海表面溫 度 (Surface Sea Temperature,SST)、人類活動、臺風(fēng)及風(fēng)暴潮等。

4.3.1 SST 影響 溫度是導(dǎo)致珊瑚礁退化的重要因素之一,同時溫度升高對黑皮海綿大量繁殖和長棘海星幼蟲孵化也有促進作用,易引起季節(jié)性生物種群暴發(fā)。根據(jù)2016―2022 年的太平島SST 數(shù)據(jù)統(tǒng)計表明,SST 具有周期變化規(guī)律,溫度最低值位于2―3 月,5―6 月達到最大值(圖14)。造礁珊瑚最適宜生長溫度為25~29℃,36℃為最高極限溫度,在黑皮海綿與長棘海星暴發(fā)的2個關(guān)鍵時間段(2017年4月―2019年2月、2020年7月―2022年4月)處于珊瑚礁生長的適宜溫度,并且超過29℃的天數(shù)熱累積效應(yīng)未達到珊瑚礁熱閾值(高于全球長期最熱月海溫1℃),因此在研究時間段內(nèi),溫度不是導(dǎo)致太平島珊瑚礁地貌類型退化的主要因素。

圖14 2016—2022年中國南海太平島SST年均變化Fig.14 Annual average change of SST on the Taiping Island in the South China Sea from 2016 to 2022

4.3.2 人為活動影響 太平島在2016―2022 年人跡罕至,島上建筑設(shè)施較為穩(wěn)定,無人為擴張(見圖9-f、10-f);島上有駐扎人員,但沒有以養(yǎng)殖為主的漁民,因而大范圍破壞珊瑚礁群的人類捕撈活動極少。故人為影響也不是造成珊瑚礁在某一時間點發(fā)生急劇退化的主要因素。

4.3.3 臺風(fēng)及風(fēng)暴潮影響 中央氣象臺·臺風(fēng)網(wǎng)⑥http://typhoon.nmc.cn記錄了2016—2022年經(jīng)過中國海域的全部臺風(fēng)路徑信息。當臺風(fēng)經(jīng)過珊瑚島礁附近時,其能量通過海流、海浪傳播到珊瑚島礁區(qū)域,受到地形的影響,會產(chǎn)生較大消耗,從而減少對珊瑚礁的沖擊,對珊瑚礁地貌類型影響較小;同時珊瑚礁生態(tài)系統(tǒng)本身具有一定的抗風(fēng)浪的能力,故在本文認為臺風(fēng)也不是造成2016—2022年太平島珊瑚礁急劇退化現(xiàn)象的主要因素。

綜上,SST、人為活動和臺風(fēng)及風(fēng)暴潮災(zāi)害不是2016—2022 年影響太平島珊瑚礁受損退化主要因素。

5 結(jié)論

以南海太平島為研究區(qū),基于2016-11-20—2022-06-07 的26 期Sentinel-2 遙感影像,利用SVM分類方法,開展了7 種珊瑚礁地貌類型信息提取。在此基礎(chǔ)上,分析珊瑚礁地貌類型變遷特征,評估生物天敵暴發(fā)事件對珊瑚礁生態(tài)系統(tǒng)退化的影響。結(jié)果表明:1)結(jié)合專家解譯知識和SVM分類算法獲取的珊瑚礁地貌類型分類結(jié)果,最高總體精度和Kappa 系數(shù)分別為96.46%和0.94。2)2017 年4 月—2019年2月是黑皮海綿暴發(fā)期,至2017年6月黑皮海綿對太平島珊瑚礁地貌類型影響程度達到最大,造成密集珊瑚沉積區(qū)、珊瑚礁叢生區(qū)和稀疏珊瑚沉積區(qū)分別減少72.92%、55.24%、44.99%。3)2020 年7 月—2022 年4 月是長棘海星暴發(fā)期,至2021年6月長棘海星對太平島珊瑚礁地貌類型影響程度最大,造成珊瑚叢生區(qū)、稀疏珊瑚沉積區(qū)、密集珊瑚沉積區(qū)分別減少59.17%、34.91%、23.97%。4)在2016—2022 年太平島珊瑚生物天敵暴發(fā)事件期間,平均珊瑚礁退化率為16.18%,最高達23.88%,發(fā)生在2017 年3—6 月;平均珊瑚礁恢復(fù)率為11.94%,最高達18.03%,發(fā)生在2017 年6—9月??傮w珊瑚礁退化率高于恢復(fù)率。

在生物天敵暴發(fā)期間,珊瑚礁地貌類型處于退化和恢復(fù)動態(tài)過程,整體上呈現(xiàn)退化趨勢。本文對太平島珊瑚生物天敵暴發(fā)造成的珊瑚礁地貌類型變化總體情況進行量化分析,表明遙感手段結(jié)合珊瑚礁地貌類型分類方法可作為監(jiān)測珊瑚礁突發(fā)事件的有效手段,從而為太平島珊瑚礁健康狀況研究提供科學(xué)依據(jù)。

猜你喜歡
太平島珊瑚礁黑皮
終于等到你!ATOLL(珊瑚礁)ST200流媒體播放機、SDA200流媒體播放/功放一體機
珊瑚礁世界的魚兒
跟蹤導(dǎo)練(三)3
雪天
闖禍的車
黑皮系列
臺旅游網(wǎng)推“太平島賞月”
六成臺灣人認為蔡英文應(yīng)登太平島
硨磲采挖對珊瑚礁生態(tài)系統(tǒng)的破壞——以西沙北礁為例
都是話多惹的禍
泰兴市| 麻江县| 巩留县| 虎林市| 绥芬河市| 高雄市| 海口市| 文成县| 安多县| 金湖县| 仁化县| 中卫市| 黑龙江省| 青州市| 南充市| 涿鹿县| 孝昌县| 五峰| 灌云县| 吉林省| 海安县| 长治市| 台安县| 平和县| 固镇县| 乌苏市| 泸州市| 麻城市| 睢宁县| 黑山县| 凤庆县| 连云港市| 大理市| 常州市| 禄丰县| 龙井市| 桐梓县| 罗源县| 静乐县| 沙雅县| 九寨沟县|