基于Landsat-8時(shí)間序列影像分析西昆侖山地區(qū)冰川滑移特征

鄭　茜　孫建寶　張　永

1　中國(guó)地震局地質(zhì)研究所地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京市華嚴(yán)里甲1號(hào),100029

基于Landsat-8時(shí)間序列影像分析西昆侖山地區(qū)冰川滑移特征

鄭茜1孫建寶1張永1

1中國(guó)地震局地質(zhì)研究所地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京市華嚴(yán)里甲1號(hào),100029

摘要：利用圖像亞像元互相關(guān)分析方法處理了Landsat-8衛(wèi)星獲取的時(shí)間序列影像數(shù)據(jù)，得到中國(guó)青藏高原西北部地區(qū)西昆侖峰區(qū)冰川勻速滑移的時(shí)空演化過程。利用亞像元影像互相關(guān)技術(shù)對(duì)Landsat-8光學(xué)影像精確配準(zhǔn)，配準(zhǔn)精度達(dá)到0.01像元，即該光學(xué)影像的水平形變監(jiān)測(cè)精度達(dá)到0.15 m。通過對(duì)2013-07～2014-08的15景Landsat-8影像進(jìn)行互相關(guān)和形變時(shí)間序列反演分析，獲得了西昆侖峰區(qū)兩條冰川的滑動(dòng)位移場(chǎng)和速度場(chǎng)。研究表明，該區(qū)域的冰川基本處于勻速滑移狀態(tài)(無明顯加速和減速現(xiàn)象)；同時(shí)也驗(yàn)證了Landsat-8光學(xué)影像在監(jiān)測(cè)較大地表位移和地殼形變事件(如沙丘移動(dòng)、地震、滑坡、火山等)上的應(yīng)用潛力。

關(guān)鍵詞：西昆侖峰區(qū)；Landsat-8；配準(zhǔn)；冰川滑移；亞像元圖像互相關(guān)

依據(jù)中國(guó)冰川編目，西昆侖山共有冰川5 485條，冰川面積8 817.78 km2，約占整個(gè)昆侖山冰川總面積的3/4以上[1]。來自昆侖山冰川的冰川水是塔里木盆地的主要水源，為生活在這一流域的近630萬人口提供飲用水[2-3]。分析研究這些高山冰川的位移場(chǎng)和它們的動(dòng)態(tài)時(shí)空演化過程，對(duì)了解青藏高原地區(qū)的氣候變化、冰川活動(dòng)對(duì)人類活動(dòng)的影響等有重要意義。

目前，常用的形變監(jiān)測(cè)技術(shù)包括差分全球定位系統(tǒng)技術(shù)(DGPS)、干涉合成孔徑雷達(dá)技術(shù) (InSAR)和亞像元光學(xué)影像互相關(guān)技術(shù)等。前兩種方法雖然可以獲得較高的形變觀測(cè)精度，但在監(jiān)測(cè)高山冰川方面具有較大的局限性[4]。在過去的幾十年里，美國(guó)的Landsat-1～7系列衛(wèi)星所獲得的數(shù)據(jù)廣泛用于地球科學(xué)研究和人類日常生活的各個(gè)領(lǐng)域。特別是在冰川研究領(lǐng)域，通過分析數(shù)10 a尺度的Landsat數(shù)據(jù)，可以直觀地測(cè)量冰川的演進(jìn)速度。然而，由于軌道控制方面不夠精準(zhǔn)，應(yīng)用Landsat-1～7系列衛(wèi)星數(shù)據(jù)需要手工選擇地面控制點(diǎn)，以實(shí)現(xiàn)精確的地面定位，形變的量算也只能通過肉眼對(duì)圖像進(jìn)行比較獲得。Landsat-8解決了上述問題，為冰川科學(xué)研究提供了非常有價(jià)值的數(shù)據(jù)源。

在本研究中，我們通過亞像元級(jí)的圖像匹配技術(shù)對(duì)軌道和幀號(hào)為145/35的15景Landsat-8影像數(shù)據(jù)進(jìn)行互相關(guān)計(jì)算，并對(duì)其組成的時(shí)間序列進(jìn)行SVD反演，獲得了西昆侖地區(qū)目標(biāo)冰川的二維位移場(chǎng)及速度場(chǎng)。

1數(shù)據(jù)源

研究區(qū)為西昆侖峰區(qū),位于青藏高原西北部的新疆自治區(qū)境內(nèi)，是我國(guó)冰川分布最密集的區(qū)域之一(圖1(a))。

圖1　研究區(qū)域Fig.1　The studied area

Landsat-8圖像覆蓋范圍為185 km×185 km的區(qū)域，因此一景圖像可以覆蓋多條冰川。在一次圖像互相關(guān)形變計(jì)算中，可以得到圖像內(nèi)所有冰川的水平位移場(chǎng)，因而具有很高的空間監(jiān)測(cè)效率。同時(shí)，16 d的固定重訪周期和不間斷的數(shù)據(jù)獲取，使得該數(shù)據(jù)的時(shí)間分辨率比其他衛(wèi)星數(shù)據(jù)源具有更大優(yōu)勢(shì)。圖1(b)為2013-09-11獲取的Landsat-8影像；圖1(c)顯示了本研究中的目標(biāo)冰川分布，即圖1(b)中矩形框所示的范圍。本研究將分析該景圖像范圍內(nèi)的形變信息。

我們從USGS網(wǎng)站下載了2013-07～2014-08共15景覆蓋西昆侖峰區(qū)的Landsat-8影像，構(gòu)成形變研究的時(shí)間序列，數(shù)據(jù)信息如表1所示。為了提高估算冰川速度場(chǎng)的精度，我們只提取空間分辨率為15 m的全色波段圖像進(jìn)行圖像互相關(guān)計(jì)算。同時(shí)，選擇云量在30%以下的影像數(shù)據(jù)，以確保圖像可以有效地精確配準(zhǔn)。

表1　研究區(qū) Landsat-8影像數(shù)據(jù)

2數(shù)據(jù)分析方法

從USGS網(wǎng)站下載的Landsat-8影像為使用SRTM DEM數(shù)據(jù)正射糾正過的圖像。只有經(jīng)過正射校正的圖像，排除成像幾何畸變，才能互相進(jìn)行精確配準(zhǔn)并分析形變。

圖像亞像元互相關(guān)分析技術(shù)，最先用于SAR圖像幅值信息的分析[5]，一方面用于計(jì)算SAR斜距向和方位向的位移場(chǎng)，另一方面用于SAR圖像的精確配準(zhǔn)。其基本原理是在兩幅幅值圖像上，按照一定間隔選取一定的配準(zhǔn)網(wǎng)格點(diǎn)，在每個(gè)網(wǎng)格點(diǎn)附近按照一定的窗口大小搜索另一幅圖像的相應(yīng)位置，直到找到窗口內(nèi)互相關(guān)系數(shù)最大的點(diǎn)，認(rèn)為這兩個(gè)圖像互相關(guān)系數(shù)最大的點(diǎn)是地面同一個(gè)散射體目標(biāo)。通過計(jì)算同一個(gè)目標(biāo)在不同時(shí)間獲取的SAR圖像上的位置變化，就能知道該點(diǎn)地面發(fā)生的位移大小。為了實(shí)現(xiàn)高精度測(cè)量，SAR圖像的計(jì)算都是在單視復(fù)圖像上、雷達(dá)坐標(biāo)系下進(jìn)行的。而在光學(xué)圖像的分析中，圖像經(jīng)過了地理編碼和正射校正，所有形變計(jì)算直接在地理坐標(biāo)系下進(jìn)行。采用同樣的方法，通過比較兩幅圖像的地面反射信息，可以實(shí)現(xiàn)亞像元級(jí)的圖像精確匹配和位移場(chǎng)計(jì)算。

使用SAR圖像配準(zhǔn)方法，將兩幅光學(xué)圖像配準(zhǔn)。一般配準(zhǔn)精度(或稱配準(zhǔn)誤差、偏差)在0.001～0.01像元范圍內(nèi)，可以達(dá)到理想的形變監(jiān)測(cè)信噪比。在圖像精確配準(zhǔn)過程中，有必要控制獲取影像的云量。 在配準(zhǔn)之后，把一幅圖像(稱為“輔圖像”或者“從圖像”)重采樣到另一幅圖像(稱為“主圖像”)的幾何坐標(biāo)下[6]。圖像重采樣采用SINC插值方法。將任意兩幅配準(zhǔn)和重采樣后的圖像在COSI-CORR軟件[7]中作互相關(guān)計(jì)算，可得到初步的位移結(jié)果。

2.1圖像配準(zhǔn)

配準(zhǔn)步驟是整個(gè)形變計(jì)算的關(guān)鍵。使用圖像亞像元互相關(guān)的方法，尋找兩幅圖像中共同的地面目標(biāo)，用最小二乘的方法建立主從圖像的多項(xiàng)式變換關(guān)系。配準(zhǔn)和重采樣過程消除了兩幅圖像的投影偏差，因而能夠保證獲得精確的地表形變信息[4]。另外一種配準(zhǔn)方法是在主從圖像之間建立映射關(guān)系，即查找表。查找表記錄了由輔圖像幾何坐標(biāo)到主圖像的映射信息。建立查找表的同時(shí)，考慮地形的影響，可以進(jìn)一步優(yōu)化圖像的配準(zhǔn)過程。在建立了主從圖像的查找表之后，繼續(xù)用前面的配準(zhǔn)多項(xiàng)式來校正查找表，并對(duì)其進(jìn)行更新?；诟潞蟮牟檎冶?，重采樣輔圖像到參考圖像坐標(biāo)系下。經(jīng)過該配準(zhǔn)步驟，主從圖像的配準(zhǔn)精度可以進(jìn)一步提升，配準(zhǔn)偏差更小，因而更有利于形變的監(jiān)測(cè)。在本研究中，我們選用獲取時(shí)間較早的圖像作為主圖像來配準(zhǔn)獲取時(shí)間較晚的另一幅，當(dāng)這幅圖像精確配準(zhǔn)之后，也可將它作為一幅新的主圖像來配準(zhǔn)其他從圖像。以此類推，配準(zhǔn)所有獲得的Landsat-8數(shù)據(jù)?？偟膩碚f，較短時(shí)間間隔內(nèi)的圖像對(duì)，其配準(zhǔn)結(jié)果相對(duì)較好，信噪比較高。

2.2圖像互相關(guān)計(jì)算

地表水平位移場(chǎng)的獲取是通過多時(shí)相正射影像圖的像元互相關(guān)計(jì)算得到的。本研究采用的圖像互相關(guān)方法是在傅里葉域里對(duì)圖像進(jìn)行反復(fù)迭代、無偏估計(jì)的過程[8]。這一過程已經(jīng)在COSI-CORR 軟件中實(shí)現(xiàn)，它將生成兩幅互相關(guān)形變圖像，分別是地表水平位移場(chǎng)的東西向和南北向分量；同時(shí)給出一幅評(píng)估形變監(jiān)測(cè)結(jié)果精度的信噪比圖像(圖2)。圖2給出的形變監(jiān)測(cè)結(jié)果中，EW向形變不明顯，而NS向位移場(chǎng)中可以看到幾條明顯的白色高值區(qū)域，顯示幾條最活躍冰川形成的地面位移場(chǎng)。由于這些冰川主要為近南北走向，與地形梯度方向一致，因而冰川位移也以南北向?yàn)橹?，其NS向的形變比較顯著，這與本次觀測(cè)的結(jié)果完全一致。整個(gè)圖像互相關(guān)的信噪比SNR比較高，圖像顯示為高亮色。

圖2　主輔圖像互相關(guān)結(jié)果Fig.2　Cross-correlation products from a pair of master and slave images

2.3后處理

2.3.1條帶去除

目前，光學(xué)圖像一般通過多傳感器機(jī)械掃描得到，這是為了避免單一光電傳感器在掃描速率上的局限性[9]。例如，Landsat所使用的傳感器系統(tǒng)，每一個(gè)頻譜波段數(shù)據(jù)都由6臺(tái)傳感器同時(shí)掃描得到。因此，一次單相掃描能夠產(chǎn)生圖像上的6行，當(dāng)下一次掃描時(shí)，衛(wèi)星軌道適量前進(jìn)，之前的6臺(tái)傳感器掃描生成下一個(gè)6行圖像(LANDSAT 用戶手冊(cè))。然而，這些傳感器并沒有統(tǒng)一的增益函數(shù)，因此導(dǎo)致互相關(guān)結(jié)果圖上有明顯的平行軌道的條帶，并在垂直軌道方向起伏波動(dòng)。例如，在ASTER和SPOT影像上出現(xiàn)過的類似振蕩偏差或CCD陣列的錯(cuò)誤振幅[10-11]。這種因傳感器偏差導(dǎo)致的條帶誤差，可以通過COSI-CORR軟件的后處理工具Destrip移除。

2.3.2數(shù)據(jù)濾波和異常值消除

互相關(guān)計(jì)算得到的形變監(jiān)測(cè)結(jié)果必須通過濾波來降低噪聲，但某些合理可信的結(jié)果常會(huì)被一些低相關(guān)系數(shù)的斑塊所掩蓋，必須先從形變場(chǎng)中移除一些過大或過小的異常值。首先計(jì)算位移值的概率密度分布，剔除低信噪比異常值，如圖3所示(藍(lán)色標(biāo)識(shí)的部分認(rèn)為是異常值，予以剔除)。

圖3　位移值概率分布Fig.3　The probability distribution of a displacement field

我們選用基于非局部均值算法的濾波方法[12]對(duì)結(jié)果圖像進(jìn)行濾波處理，其在保證局部細(xì)節(jié)的同時(shí)降低了加性高斯白噪聲。非局部均值算法是通過逐步調(diào)試來實(shí)現(xiàn)數(shù)據(jù)的濾波，具體方法是：對(duì)每個(gè)像元點(diǎn)，將其周圍的小區(qū)域作為一個(gè)中心斑塊，斑塊大小一般為5×5或7×7像元。依據(jù)中心斑塊與其他斑塊的相似度來確定中心斑塊周圍更大范圍的搜索區(qū)域，最終以該搜索區(qū)域內(nèi)所有像元的平均值作為中心像元的濾波值。

3觀測(cè)結(jié)果分析

3.1西昆侖峰區(qū)冰川位移場(chǎng)

在西昆侖山地區(qū)，大部分高山冰川呈南北走向，因此在互相關(guān)結(jié)果中，東西分量的水平位移信號(hào)相對(duì)較弱(圖2)。所以，我們以南北向水平位移結(jié)果為例，研究西昆侖峰區(qū)冰川滑移及其隨時(shí)間的演化。圖4(a)給出了西昆侖峰區(qū)的地表南北向水平位移場(chǎng)，這是由2013-09-27和2013-12-16獲取的兩幅Landsat-8正射影像圖作互相關(guān)得到的，其中一幅覆蓋目標(biāo)冰川的正射影像如圖4(b)所示。該地區(qū)的位移最大值出現(xiàn)在圖4(a)中的灰色矩形框內(nèi)。在80 d內(nèi)，兩條目標(biāo)冰川的滑移速率達(dá)到了1.5 m/d。此外，該地區(qū)內(nèi)其他冰川上稍小的滑移速率也比較明顯。圖上紅色表示冰川向南移動(dòng)，藍(lán)色表示向北移動(dòng)。我們將滑移速率最大的兩條目標(biāo)冰川作為研究對(duì)象，其位移場(chǎng)在圖4(c)中顯示。

圖4　西昆侖峰區(qū)的地表南北向水平位移場(chǎng)Fig.4　Horizontal displacement field of the glaciers in the west Kunlun mountains

圖5　目標(biāo)冰川位移場(chǎng)的時(shí)間序列演化Fig.5　The time-series plot of the displacement fields of the target glaciers

3.2 時(shí)間序列分析

我們處理分析了2013-07～2014-08的15景Landsat-8影像(表1)，并且利用奇異值分解(SVD)算法計(jì)算了目標(biāo)冰川的年地表形變[13]。圖5顯示西昆侖山地區(qū)的冰川在觀測(cè)時(shí)間段內(nèi)，滑移量逐漸累積(南北向)。在384 d內(nèi)，冰川位移隨時(shí)間呈明顯遞增趨勢(shì)。

為了確定目標(biāo)冰川精確的滑移速率，我們?cè)趦蓷l冰川上選取6個(gè)參考點(diǎn)(圖4(c))，分析它們的位移(圖6)和速率(圖7)差異。由于圖像處理過程中存在一些不可避免的隨機(jī)噪聲，導(dǎo)致利用奇異值分解算法反演位移場(chǎng)時(shí)間序列的結(jié)果存在部分“負(fù)增長(zhǎng)”值(圖6)，但依然可以明顯看到，在384 d內(nèi)位移場(chǎng)隨時(shí)間近似呈線性遞增，同時(shí)曲線的斜率在不同的參考點(diǎn)上是不同的，說明各點(diǎn)的滑動(dòng)速率并不完全一致。圖7給出了6個(gè)參考點(diǎn)的速率變化?？梢钥闯?，2013-11～2014-08冰川滑移速率保持在0.2～0.6 m/d，而對(duì)于每個(gè)參考點(diǎn)，冰川速度場(chǎng)呈勻速運(yùn)動(dòng)趨勢(shì)。因此，研究區(qū)的這兩條冰川始終保持了極為高速的滑動(dòng)，其能夠被Landsat-8這種圖像分辨率的數(shù)據(jù)檢測(cè)到。與其相比，InSAR形變監(jiān)測(cè)方法，由于時(shí)間去相干效應(yīng)的影響，監(jiān)測(cè)如此高速的地面冰川滑動(dòng)存在很大的局限性，這也反映了Landsat-8數(shù)據(jù)的觀測(cè)優(yōu)勢(shì)。

圖6　位移時(shí)間序列曲線Fig.6　The curves of the displacement time-series

圖7　速度時(shí)間序列曲線Fig.7　The curves of the velocity time-series

4結(jié)語

本研究首次利用Landsat-8光學(xué)圖像形變分析方法，獲得西昆侖峰區(qū)冰川的形變場(chǎng)時(shí)間序列演化特征。基于光學(xué)影像資料提取地表形變，應(yīng)用時(shí)間序列分析的方法、原理，對(duì)Landsat-8光學(xué)影像進(jìn)行形變分析，獲得西昆侖峰區(qū)冰川滑移的地表水平形變及其時(shí)間序列演化過程。結(jié)論如下：

1)光學(xué)影像互相關(guān)技術(shù)監(jiān)測(cè)地表形變，特別是在監(jiān)測(cè)較大形變和較快移動(dòng)的地面目標(biāo)時(shí)具有一定的優(yōu)勢(shì)。

2)基于幅值信息的影像配準(zhǔn)方法、亞像元影像互相關(guān)形變計(jì)算、SVD反演時(shí)間序列——這些方法能夠有效獲得研究區(qū)時(shí)空形變場(chǎng)，提高測(cè)量精度。

3)在384 d觀測(cè)期內(nèi)，西昆侖峰區(qū)冰川位移場(chǎng)隨時(shí)間近似呈線性變化，滑移速率在0.2～0.6 m/d。

4)本研究的成果不僅對(duì)冰川動(dòng)力學(xué)和氣候變化方面的研究有一定意義，而且驗(yàn)證了Landsat-8光學(xué)影像在形變分析中的應(yīng)用潛力。該方法也適用于分析大地震和滑坡等形變場(chǎng)。

致謝：USGS為本研究提供了所有的Landsat-8影像數(shù)據(jù),在此表示感謝!

參考文獻(xiàn)

[1]焦克勤,姚檀棟,李世杰.西昆侖山32 Ka來的冰川與環(huán)境演變[J].冰川凍土,2000，22(3),250-256(Jiao Keqin,Yao Tandong,Li Shijie.Evolution of Glaciers and Environment in the West Kunlun Mountains During the Past 32 Ka[J].Journal of Glaciology and Geocryology,2000,22(3),250-256)

[2]李成秀,楊太保,田洪陣.1990～2011 年西昆侖峰區(qū)冰川變化的遙感監(jiān)測(cè)[J].地理科學(xué)進(jìn)展,2013,32(4)：548-559 (Li Chengxiu,Yang Taibao,Tian Hongzhen.Variation of West Kunlun Mountains Glacier during 1990-2011[J].Progress in Geography,2013,32(4):548-559)

[3]李成秀. 昆侖山冰川和積雪變化的遙感監(jiān)測(cè)[D]. 蘭州:蘭州大學(xué),2014(Li Chengxiu. Remote Sensing Monitoring of Glacier and Snow Cover Changes in the Kunlun Mountain[D]. Lanzhou :Lanzhou University,2014)

[4]Berthier E,Vadon H,Baratoux D,et al. Surface Motion of Mountain Glaciers Derived from Satellite Optical Imagery[J].Remote Sensing of Environment,2005,95(1):14-28

[5]劉方堅(jiān). 星載雷達(dá)影像配準(zhǔn)方法的研究與實(shí)現(xiàn)[D]. 北京:中國(guó)地質(zhì)大學(xué)(北京),2006(Liu Fangjian. The Registration Algorihtm Research of Spaceborne InSAR Images[D]. Beijing :China University of Geosciences(Beijing),2006)

[6]Massonnet D,Feigl K. Radar Interferometry and Its Application to Changes in the Earth’s Surface[J].Reviews of Geophysics,1998,36(4): 441-500

[7]Leprince S,Berthier E,Ayoub F,et al. Monitoring Earth Surface Dynamics with Optical Imagery[J].EOS Transactions,2008,89(1):1-12

[8]Scherler D,Leprince S,Strecker M R. Glacier-Surface Velocities in Alpine Terrain from Optical Satellite Imagery -Accuracy Improvement and Quality Assessment [J].Remote Sensing of Environment,2008,112(10):3 806-3 819

[9]Horn B K P, Woodham R J. Dest-Riping Landsat MSS Images by Histogram Modification[J].Computer Graphics and Image Processing,1979(10):69-83

[10] Leprince S,Muse P,Avouac J P. In-Flight CCD Distortion Calibration for Pushbroom Satellites Based on Subpixel Correlation[J]. IEEE Transactions on Geoscience and Remote Sensing,2008,46(9):2 675-2 683

[11] Teshima Y,Iwasaki A. Correction of Attitude Fluctuation of Terra Correction of Attitude Fluctuation of Terra with Parallax Observation[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(1):222-227

[12] Ayoub F,Leprince S,Keene L.User’s Guide to COSI-CORR[Z].Co-Registration of Optically Sensed Images and Correlation,2009

[13] Berardino P,Fornaro G,Lanari R,et al. A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2 375-2 383

Foundation support:National Natural Science Foundation of China,No. 41374040.

About the first author:ZHENG Qian,postgraduate,majors in time-series analysis of SAR interferometry and crustal deformation application,E-mail: 359300096@qq.com.

Fast and Uniformly Slipping Western-Kunlun Glaciers from Time-series Deformation Analysis Using Periodically Captured Landsat-8 Imagery

ZHENGQian1SUNJianbao1ZHANGYong1

1State Key Laboratory of Earthquake Dynamics,Institute of Geology,CEA,A1 Huayanli,Beijing 100029,China

Abstract:We processed time-series data,using a sub-pixel cross-correlation analysis method, to obtain the temporal evolution of the surface displacements of the west Kunlun mountain glaciers in the north-western Tibetan plateau. This work co-registers periodically acquired Landsat-8 images using the sub-pixel cross-correlation technique at an accuracy of 0.01 pixels, implying a horizontal displacements detecting precision of the optical image 0.15 m. By cross-correlation analysis of the Landsat-8 images acquired from July 2013 to August 2014 and time-series inversion, this study constructs the displacement and velocity fields of two glaciers in the west Kunlun mountains. It demonstrates that the glaciers in the study area are slipping fast and uniformly (no obvious acceleration or deceleration). Based on Landsat-8 images, the method of this study can potentially be used for measuring ground deformation caused by dunes shifting, earthquakes, landslides, or volcanoes etc.

Key words:west Kunlun mountains; Landsat-8; co-registration; glacial motion; sub-pixel cross-correlation

收稿日期：2016-03-16

第一作者簡(jiǎn)介：鄭茜，碩士生，研究方向?yàn)镮nSAR時(shí)間序列分析與地殼形變,E-mail: 359300096@qq.com;zhengqianabc@sina.com。

DOI：10.14075/j.jgg.2016.07.010

文章編號(hào)：1671-5942(2016)07-0604-05

中圖分類號(hào)：P228

文獻(xiàn)標(biāo)識(shí)碼：A

項(xiàng)目來源：國(guó)家自然科學(xué)基金(41374040)。