近58 a川西北高原的氣候變化及其生態效應

苑全治， 任 平

（1.四川師范大學 地理與資源科學學院，四川 成都610101； 2.四川師范大學 川西資源環境與可持續發展研究中心，四川 成都610066；3.四川師范大學 西南土地資源評價與監測教育部重點實驗室，四川 成都610066）

川西北高原位于青藏高原東緣，范圍包括四川省甘孜藏族自治州和阿壩藏族羌族自治州，以及涼山彝族自治州的一部分，面積約為2.55×105km2，是我國第一級階梯向第二級階梯的過渡區域.區域內大小河流交織，是長江、黃河上游重要的水源涵養區，也被譽為四川盆地的“水塔”.草地面積占60%以上，牧草產量高、質量好，使川西北高原成為我國五大牧區之一［1］.但近幾十年來，川西北生態環境退化嚴重，高寒草地沙化日益加重，部分區域沙化面積急劇擴張，有逆向演替的趨勢，即由高寒草甸→草原→荒漠草原→荒漠逐漸轉變，生態系統的脆弱性在升高，已經影響到牧區經濟的可持續發展，甚至威脅到長江、黃河源區的生態安全［2-3］；這使得該區域及周邊的草地狀況備受學者關注［4-5］.截至2009年，川西北沙化土地面積已經達到8.22×103km2，占四川省沙化土地面積的89.9%.1994—2009年沙化總面積增加了28.1%，表明近十幾年沙化加劇.沙化土地主要分布于高寒草地區和干旱河谷區，高寒草地沙化土地占全區沙化面積的75.1%.草地沙化會顯著影響生態系統結構和功能的穩定，沙化的擴張會造成草地嚴重退化、土壤肥力降低、多樣性降低.因此，該區域的草原沙化問題亟待研究［6］.

高寒草地的沙化驅動力可以歸納為2方面，即人為因素和自然因素.人為因素主要有過度放牧、草地開墾、毒藥滅鼠、挖藥取薪、修路采礦和牧民定居等，其中長時間的持續過度放牧被認為是導致高寒草地沙化的最主要人為驅動因素［7-8］.連續超載放牧，草地被過度踐踏、過量啃食，植被得不到充分恢復，覆蓋度持續下降，土壤流失加劇，草地出現退化和沙化［9］.導致草地沙化的自然因素主要包括氣候變化、土壤質地和鼠兔泛濫等，其中以氣候變暖為標志的氣候變化是導致高寒草地沙化的主要誘因［10-12］，其他因素加劇了草地的沙化.近幾年，通過封育限牧、退耕還草、招鷹滅鼠等一系列生態恢復手段，已經在一定程度上控制了人為因素的影響，并且隨著人們保護生態環境的意識不斷增強，高寒草地沙化的人為影響將繼續減輕.因此，可以預見，在未來的高寒草地沙化過程中，氣候變化等自然因素的影響將會越來越凸顯.本文以此為出發點，利用川西北高原石渠、德格、色達、甘孜、新龍、巴塘、理塘、稻城、九龍、康定、道孚、小金、馬爾康、松潘、紅原、若爾蓋16個氣象站點數據，研究過去58 a與草地沙化相關的關鍵氣候變化特征，并結合已有的研究討論氣候變化在川西北高寒草地沙化過程中的作用.

1 數據與方法

1.1 數據本文所用氣象數據來自國家氣象信息中心，包括1961—2018年川西北高原16個全國基本氣象站點的最低氣溫、平均氣溫、最高氣溫、風速、降水量、相對濕度和日照時數的日值數據.

1.2 潛在蒸散與氣候干燥度計算潛在蒸散（potential evapotranspiration，ETo）是假設水分供給充足時，一定氣象條件下區域的最大蒸散量.潛在蒸散是實際蒸散的理論上限，能夠表征區域氣候的水分蒸散能力.由于實際蒸散觀測資料的缺乏，人們往往參照潛在蒸散來估算實際蒸散量［13］.

草地沙化與表層土壤水分的下降有密切聯系，而蒸散的變化直接影響土壤水分及草地植被的生長.因此，本文將研究川西北高原ETo的變化趨勢，以分析其對草地沙化可能的影響.因為長時間序列的潛在蒸散觀測資料很難獲取，所以學者通常采用數學模型進行模擬計算.Penman-Monteith模型（簡稱為P-M模型）由于考慮了蒸散中的能量平衡和水汽傳輸的空氣動力學結構［14］，能夠反映各氣候要素的綜合影響，具有明確的物理意義，因此得到很多學者的驗證和推薦.聯合國糧農組織（FAO）根據假想的參考面又進一步修改了P-M模型（即FAO56 P-M模型），改進后的模型能夠適應不同氣候類型區域的ETo計算，目前已經得到廣泛應用.Yin等［15］校正了FAO56 P-M模型中計算太陽凈輻射Rn的經驗參數，使其能夠更好地適應中國的ETo計算.計算方法如下：

其中，ETo為日潛在蒸散量（mm），Δ為飽和水汽壓曲線斜率（kPa/℃），T為日平均氣溫（℃），Rn為太陽 凈輻射（MJ/（m2·d）），γ為干濕常數（kPa/℃），G為土壤熱通量（MJ/（m2·d）），σ為Stefan-Boltzmann常數（4.903×10-9MJ/（K4·m2·d）），U2為2 m高處的風速（m/s），ea為實際水汽壓（kPa），es為平均飽和水汽壓（kPa），N為可照時數（h），n為實際日照時數（h），Rso為晴天太陽總輻射（MJ/（m2·d）），Tmin，k為日最低氣溫絕對溫標（K），Tmax，k為日最高氣溫絕對溫標（K）.ETo計算公式中其他參量的計算方法參考文獻［16］.

氣候的干濕狀況對于植被分布和生產力有重大影響，本文將研究川西北高原的氣候干濕狀況變化，以討論其對草地沙化的作用.計算氣候干濕狀況的方法常用的是道庫恰耶夫和維索茨基在20世紀初提出的降水和蒸發之比［16］，即基于大氣水分收支的2個分量：降水量h和潛在蒸散量ETo.本研究采用干燥度指數作為衡量干濕狀況的指標，該指數已在很多研究中得到了驗證：

其中，K為干燥度指數，h為降水量（mm），ETo為潛在蒸散量（mm）.

1.3 氣候傾向率計算氣候傾向率s是氣候因子或特征值的每10 a線性變化率，是分析氣候變化趨勢的重要方法.氣候數據往往是離散的時間序列數據，因此在計算氣候傾向率時一般需要對離散的氣候數據進行線性擬合，即用一元線性方程表示因變量（氣候數據）與自變量（時間）之間的函數關系，所得直線斜率（slope）的10倍即為氣候傾向率.直線擬合通常采用最小二乘法（least square method），尋找一條直線擬合氣象數據，使擬合值與實際值之間的偏差平方和最小，即假設回歸方程為

根據測試點（X1，Y1），（X2，Y2），…，（Xn，Yn），利用極值法可推導出直線的斜率a和截距b：

任意時間序列的離散氣候數據，都可以用最小二乘法進行直線擬合，但并非所有氣候因子都有顯著的線性變化趨勢.因此，需要對擬合的直線進行顯著性檢驗.方程參數估計F檢驗的P值代表氣象數據線性趨勢的顯著程度.P值越小越顯著；反之，亦然.當P值高于假設的顯著性水平0.05時，統計不顯著；當P值≤0.01時，方程為高度顯著，文中記為**；當0.01＜P≤0.05時，方程統計顯著，記為*.除了分析川西北高原內16個氣象站點的氣候因子變化趨勢外，為了避免單個氣象站點氣候數據誤差帶來的不確定性影響，本文將對高原整體氣候變化特征進行分析.雖然區域內氣象站點較多且在空間上的分布相對均勻，但由于地形起伏，各個氣象站點所能夠代表的區域范圍不盡相同.因此，把16個氣象站點的氣象數據相加或者平均后來分析變化趨勢，并不能較好地反映川西北高原的整體氣候變化特征.本文將利用空間插值的方法，首先把離散的站點數據轉換成空間上連續的柵格數據，然后把各個柵格的氣候數據進行平均并分析其變化，以代表川西北高原的整體氣候變化趨勢.采用薄盤光滑樣條插值法（thin plate smoothing splines），以經度和緯度作為樣條函數自變量，對氣象站點觀測數據進行空間插值（其中對氣溫、潛在蒸散和風速的插值用高程作為協變量），數據空間分辨率為0.1°.對插值結果的精度進行檢驗，結果顯示溫度、潛在蒸散、降水和風速的實際值和擬合值之間的擬合優度R2分別為0.95、0.93、0.90和0.81，其中風速的擬合度稍差，這可能與風速更容易受到小尺度地理要素差異的影響有關系，但也能夠滿足本研究的需求.氣候干燥度柵格數據由潛在蒸散和降水量根據（3）式計算得到.

2 結果與分析

把形成的1961—2018年川西北高原各年的氣溫、降水、風速、潛在蒸散和氣候干燥度柵格數據進行平均形成散點圖，并在此基礎上做58 a平均線和線性趨勢線.把川西北高原各個氣象站點1961—2018年的氣候數據進行平均，并做線性擬合，所得平均值和氣候傾向率分別標記在氣象站點旁的橫線上面和下面.

2.1 氣溫變化特征如圖1所示，從川西北高原氣溫的總體變化情況看，雖然各年平均氣溫圍繞趨勢線上下波動，但總體呈明顯的上升趨勢，趨勢線的擬合優度R2達到了0.58.過去58 a川西北高原氣候傾向率s為0.25℃/Da.在1986年以前，區域總體處于溫度比較低的階段；1987—2000年，氣溫圍繞平均值上下波動；2001年以后，川西北高原整體進入相對溫暖的時期.

圖1 川西北高原氣溫整體變化趨勢Fig.1 Overall trend of air temperature of Northwestern Sichuan Plateau

從川西北高原氣溫的總體分布來看，呈北低南高的趨勢，即越深入青藏高原氣溫越低（圖2）.氣溫的變化趨勢顯示：除小金站外，其他氣象站點的氣溫都有顯著的增暖趨勢，但增暖幅度差異較大.氣候傾向率最大的站點是理塘站，達到0.43℃/Da.區域中部地區的增暖幅度比南北兩側要小.

圖2 各站點平均氣溫及氣候傾向率Fig.2 Mean air temperature and climatic tendency of each station

2.2 降水變化特征由圖3可知，川西北高原多年平均降水量h為700.68 mm，1961—2018年的氣候傾向率為8.3 mm/Da，趨勢線斜率小，與平均線幾乎重合（圖3）.但趨勢線的統計不顯著，直線的R2也僅為0.07，說明川西北高原降水變化的線性趨勢不明顯，更多的表現出非線性特征.1998年降水量最大，達825.66 mm，比平均值高19.4%，這一年也是厄爾尼諾年，但并非所有厄爾尼諾年都會給川西北高原帶來大量降水，如1986—1987年、2009—2010年等，這可能是因為厄爾尼諾主要影響太平洋季風區的氣候.而受高原地形阻擋，東南夏季風所帶來的太平洋水汽對川西北高原降水影響有限，隨西南季風沿著縱向嶺谷區北上的印度洋水汽卻會對該區域降水產生較大影響［17］.

圖3 川西北高原降水量整體變化趨勢Fig.3 Overall trend of precipitation of Northwestern Sichuan Plateau

區域內降水量由東南向西北遞減，東南的九龍和康定站的年降水量較高，均大于800 mm，年降水量最小的是海拔最高的石渠站（圖4）.從變化趨勢來看，若爾蓋站的降水量在這58 a里平均減少，其余站點降水量則平均增加，其中馬爾康、新龍、小金、道孚、康定等站的降水量呈顯著增長趨勢，但顯著性不同.

圖4 各站點平均降水量及氣候傾向率Fig.4 Mean precipitation and climatic tendency of each station

2.3 風速變化特征1961—2018年川西北高原的年平均風速v為2.7 m/s，并呈現出降低趨勢，但氣候傾向率較小，為-0.03 m/（s·Da），但下降趨勢不顯著，R2也只有0.02.從年平均風速的波動情況來看，大體分成4個階段：第一階段從1961—1968年，年平均風速變化不大，且比58 a平均值偏低，為低風速階段；第二階段從1969—1988年，這段時間的年平均風速變化同樣不大，但比58 a平均值偏高，屬高風速階段；第三階段從1989—2002年，年平均風速呈線性下降趨勢，為風速下降階段；第四階段從2003—2018年，年平均風速有一定的波動，但總體為低風速階段（圖5）.

圖5 川西北高原風速整體變化趨勢Fig.5 Overall trend of wind speed of Northwestern Sichuan Plateau

由圖6可知，川西北高原康定站年平均風速最大，達2.94 m/s；馬爾康站最小，為1.11 m/s.從變化趨勢來看，風速平均增加的站點有紅原、松潘、甘孜和巴塘，但氣候傾向率都不大，甘孜站的線性增加趨勢不顯著.其余站點的年平均風速平均降低，且大部分呈顯著線性降低趨勢，其中石渠站的下降傾向率最大，達到-0.20 m/（s·Da）.總體來看，在風速變化顯著的站點中，風速總體呈下降趨勢.

圖6 各站點平均風速及氣候傾向率Fig.6 Mean precipitation and climatic tendency of each station

2.4 潛在蒸散量變化特征近58 a來，川西北高原的年潛在蒸散量平均為647.60 mm，氣候傾向率為4.4 mm/Da，但線性趨勢不顯著，R2只有0.05，如圖7所示.從ETo的整體波動情況來看，大概經歷了1961—1973年波動上升、1974—2000年波動下降和2001—2018年較平穩波動3個階段.值得注意的是，第一個階段的上升趨勢較強，而下一個階段的波動下降較緩。

圖7 川西北高原ETo整體變化趨勢Fig.7 Overall trend of ETo of Northwestern Sichuan Plateau

從年平均ETo空間分布來看（圖8），總體呈南高北低趨勢，其中ETo最大的是小金站，達990.07 mm，最低的是松潘站，為642.68 mm.松潘、紅原、若爾蓋、巴塘等站ETo呈顯著增加的線性趨勢，而道孚、德格、新龍、稻城和小金等站ETo則顯著降低.

圖8 各站點平均ETo及氣候傾向率Fig.8 Mean ETo and climatic tendency of each station

2.5 氣候干燥度變化特征如圖9所示，近58 a川西北高原的年平均氣候干燥度為0.93，年平均降水量大于ETo，整體氣候濕潤.氣候傾向率為-0.007，說明氣候在平均變濕，但線性擬合沒有通過顯著性檢驗，且R2只有0.012，說明氣候總體干濕狀況的線性變化趨勢不明顯，非線性特征更加突出.

圖9 川西北高原氣候干燥度整體變化趨勢Fig.9 Overall trend of climate dryness of Northwestern Sichuan Plateau

由圖10可知，區域氣候干濕狀況的空間分布趨勢大體上是西干東濕.在所有氣象站點中，氣候最為干燥的是巴塘站，干燥度K為2.10，康定站最濕潤.按照干濕區劃方法，K≤1.0的紅原、康定、馬爾康、松潘和九龍站屬于濕潤氣候；K介于1.0～1.5的若爾蓋、理塘、色達、新龍、德格、石渠、道孚、甘孜和稻城站屬半濕潤氣候；K＞1.5的巴塘和小金站為半干旱氣候.

圖10 各站點平均干燥度及氣候傾向率Fig.10 Mean climate dryness and climatic tendency of each station

3 討論

本文在上述內容中研究了川西北高原各站點

和整體的氣候變化趨勢，目的是為了討論上述關鍵氣候變化特征的環境效應，尤其是氣候變化對高寒草地沙化的影響.氣溫、降水、風速、ETo雖然都是影響川西北高原植被的關鍵氣候因子，但它們之間并不是獨立的，相互之間存在著密切的聯系.因此，首先需要對各個因子的相關關系進行分析，找出影響高寒草地的關鍵氣候因子，然后結合現有研究和遙感數據討論這些關鍵因子的變化對草地沙化的影響作用.

3.1 氣候因子之間的相關關系各因子之間相關關系的大小可通過Pearson系數r來衡量：

其中，xi和yi分別代表2個氣候因子中的第i個數據.若r＞0時，2個因子可能是正相關；而r＜0時，則可能是負相關.r絕對值越接近1，表示2個因子相關性越強；越接近0，則說明相關性越弱.

由表1可見，反映區域熱量條件的潛在蒸散ETo與風速和氣溫的相關性較顯著，而能夠表征區域干濕狀況的干燥度指數則與（3）式中的降水量和ETo相關關系均較為顯著.表中還顯示風速與干燥度指數有一定正相關關系，降水量與ETo存在一定的負相關關系；這可能是因為降水時太陽輻射較弱，降低了蒸散量.

表1 各氣候因子之間的Pearson相關系數Tab.1 Pearson correlation coefficient between climatic factors

在上述5個指標中，潛在蒸散ETo與干燥度指數是2個綜合指標，也是常用來反映一個地區水分和熱量綜合水平的指標.從2個指數的總趨勢來看，ETo呈略微下降趨勢，而干燥度指數是輕微上升趨勢，但二者的線性變化都不顯著，都是圍繞平均值上下波動，而且存在一定的階段性.在20世紀70年代以前，二者的變化規律相似，在波動中上升.70年代以后到2000年左右，干燥度指數的波動更加強烈，而ETo的變化則相對和緩.2005年以來，干燥度指數有更加明顯的下降趨勢，ETo沒有很顯著的變化趨勢.

3.2 氣候因子與地表NDVI的相關分析歸一化植被指數（Normalized Difference Vegetation Index，NDVI）是根據植被對近紅外波段強烈反射和對紅波段強烈吸收的特點，基于植被對這2個波段的反射率所構建的植被指數.NDVI與植被生物量、覆蓋度、葉面積指數、凈初級生產力等生態參量有良好的相關關系，且易于獲取，因此被廣泛應用于各類植被信息的遙感監測.本文利用長時間序列的AVHRR-NDVI（8×8 km）數據，把各個站點的氣候因子與NDVI進行相關分析，結果見表2.

表2 各氣候因子與NDVI之間的相關關系Tab.2 Correlation between climate factors and NDVI

表2中顯示各個站點所代表的區域植被類型不同，這是因為高寒草甸雖然是區域典型植被類型，但并非唯一植被類型，在高原面以外的高山深谷地形發育了復雜的植被垂直帯譜，且大多屬于濕潤海洋性帶譜，包括高山上垂帶譜和峽谷下垂帶譜，植被類型有常綠闊葉林、落葉闊葉林、常綠針葉林、灌叢、草甸等.位于高原面上的氣象站點包括石渠、色達、紅原、若爾蓋、理塘和稻城，這些站點所代表的區域植被均以高寒草地為主，而其他站點所處區域則植被類型復雜.從表2也可以看出，植被類型復雜的區域NDVI與單個氣候因子的相關關系大多很微弱，這也證明了生態系統越復雜，對氣候變化的敏感性和脆弱性就越低.高原面上的草地植被NDVI雖然與部分氣候因子有一定相關關系，但相關系數均不高，不確定性較大，如降水與若爾蓋草地NDVI相關關系為正，而與道孚、理塘、九龍等站的草地NDVI相關關系卻為負.溫度與草地NDVI相關關系大多為正，說明溫度的升高對高寒草地的生長有有利的一面，但ETo與色達站NDVI的相關關系為負說明溫度升高導致的蒸散增加對草地植被是不利的.單個氣候因子與NDVI之間相關關系的不確定性較大，一方面是因為遙感所獲取的NDVI包括了人類活動的影響，另一方面是由于氣候變化對植被的影響復雜，有負面的影響，也有正面的影響，二者相互抵消后的結果很難通過統計的方法來反映.

3.3 氣候變化的地表覆被效應分析本文的研究表明，川西北高原溫度呈顯著上升趨勢，而溫度的持續上升會導致地表蒸發量加大、多年凍土退化和土壤凍結時間縮短［18-19］，這些都會引起土壤表層含水量減少［8］，草地植被生長因受水分限制而發生退化，植被覆蓋度下降，進而加速了土壤的水蝕和風蝕，最終造成草地的沙化，沙化后相對干燥的土壤和稀疏的植被為鼠類所喜歡的環境，鼠類大量繁殖并啃食草地植被，又加劇了高寒草地的退化與沙化.溫度升高與高寒草地退化、沙化之間的關系可以總結為圖11.

圖11 氣溫升高所主導的高寒草地、沙化過程Fig.11 Desertification process of alpine grassland dominated by the increase of temperature

降水的增多會在一定程度上緩解溫度升高對高寒草地的不利影響［20］，但根據前文研究，川西北高原降水變化波動性大，線性變化不顯著，即降水不會持續性增加，這對于高寒草地退化和沙化的緩解作用有限；且降水量的陡增反而會加劇對土壤的水蝕強度，陡減又會加劇土壤的風蝕，進而促進土壤的沙化.降水變化的區域差異較大，紅原和若爾蓋站的降水在平均減少，從近年來的研究來看，紅原和若爾蓋草地的沙化問題也更為嚴峻［21-22］.川西北高原在近58 a的風速有降低的趨勢，這有利于降低風蝕強度，但風速的變化同樣存在區域差異，石渠站風速降低明顯，紅原站風速卻有線性增加趨勢.由于風速變化率較小，在草地沙化過程中主要起加劇作用，而并非主導氣候因子.此外，極端氣候事件也會導致小尺度上的草地退化與沙化，如雪崩、冰崩、暴風雪使草地被大雪覆蓋［23］；極端強降水使土壤侵蝕加劇［24］.

氣候變暖對草地生態系統的影響是多方面的，有上述的不利影響，同時也會有正面效應.變暖使生長季節氣溫和土壤溫度增加［25］，有利于植物生長.植物葉片的凈光合速率隨溫度的升高而增加［26］，增溫使冰川退縮、雪線上升，植被可覆蓋的海拔高度更高，草場綠化面積擴大［27］.綜合以上分析，顯著的增暖趨勢可能是川西北高寒草地退化、沙化的主導氣候因子，同時降水的減少和風速的增加可能加速了部分區域（如紅原、若爾蓋地區）的草地沙化過程，但由于氣候變化與地表覆被過程之間的關系非常復雜，同時存在正負兩方面的生態效應，且人為擾動在草地沙化中的作用不能忽略.因此，通過統計方法很難了解氣候變化在高寒草地沙化中的作用.氣候-植被模型能夠模擬無人類擾動的氣候情景下植被的動態變化過程［28-30］，可以成為未來草地沙化與氣候變化關系研究的重要工具.

4 結論

川西北高原位于青藏高原東部，高寒草甸是該區域高寒、濕潤氣候下的典型植被，但近年來草地沙化問題突出，已經影響到長江、黃河源區的生態安全和區域畜牧業的可持續發展，以過度放牧、墾荒為主的人為因素和以氣候變化為主的自然因素被認為是導致沙化的主要驅動力.本文著眼于高寒草地的氣候驅動力，研究了與草地沙化關系密切的氣溫、降水、風速、潛在蒸散和氣候干燥度的變化趨勢，并探討了氣候變化對草地沙化的可能影響，得出的結論如下.

1）1961—2018年川西北高原溫度呈顯著上升趨勢，總體的氣候傾向率為0.25℃/Da，溫度增幅最大的是理塘站，達到0.43℃/Da.降水量呈現平均增多的趨勢，但線性趨勢不顯著，整體氣候傾向率較小，為8.3 mm/Da.而且降水量變化的區域差異較大，康定站增幅最大且線性趨勢顯著，為31.26 mm/Da.風速總體呈下降趨勢，但變化率較小，氣候傾向率為-0.03 m/（s·Da）.

2）近58 a來，川西北高原潛在蒸散呈上升趨勢，氣候傾向率為4.4 mm/Da，期間大概經歷了1961—1973年波動上升、1974—2000年波動下降和2001—2018年較平穩波動3個階段.ETo變化的區域差異也較大，上升幅度最大的是若爾蓋、松潘和巴塘站.氣候干濕狀況在平均變濕潤，線性趨勢不顯著.德格、馬爾康、道孚和康定站的氣候變濕潤趨勢有較高的顯著性.

3）ETo和氣候干燥度2個綜合氣候因子與氣溫、降水、風速等單因子有密切聯系，ETo與風速和氣溫的相關系數較大，氣候干燥度與降水的相關系數最大.但通過與地表植被NDVI的相關性研究發現，無論綜合氣候因子，還是單氣候因子與NDVI的變化大多無顯著相關性，這一方面由于遙感獲取的NDVI包括了人類活動的影響，另一方面氣候變化對植被的影響過程復雜，很難通過統計的方法去解釋二者之間的關系.因此，未來的研究中需要借助氣候-植被模型，模擬無人類擾動的氣候情境下，川西北高寒草地的變化過程，以更加深入的了解氣候變化與草地沙化之間的關系.

致謝四川師范大學川西資源環境與可持續發展研究中心科技項目（KY202002）對本文給予了資助，謹致謝意.