平坦沙質(zhì)地表蝕積量計(jì)算模型研究

安志山，張克存，屈建軍，譚立海，張 號(hào)

（中國(guó)科學(xué)院 寒區(qū)旱區(qū)環(huán)境與工程研究所 沙漠與沙漠化重點(diǎn)實(shí)驗(yàn)室，敦煌戈壁荒漠生態(tài)與環(huán)境研究站，甘肅省風(fēng)沙災(zāi)害防治工程技術(shù)研究中心，甘肅 蘭州730000）

風(fēng)是塑造地貌形態(tài)的基本營(yíng)力之一，也是沙粒發(fā)生運(yùn)動(dòng)的動(dòng)力基礎(chǔ)。半個(gè)多世紀(jì)以來(lái)，許多學(xué)者對(duì)氣流和風(fēng)沙活動(dòng)進(jìn)行了研究，提出各種各樣的表示方法，如輸沙勢(shì)、最大可能輸沙量和輸沙率公式等。輸沙勢(shì)是衡量區(qū)域風(fēng)沙活動(dòng)強(qiáng)度及風(fēng)沙地貌演變的重要指標(biāo)，也是目前風(fēng)沙活動(dòng)強(qiáng)度計(jì)算、應(yīng)用最為廣泛的方法，表示潛在的輸沙能力［1－2］。輸沙率是指單位時(shí)間通過(guò)單位床面寬度的沙粒質(zhì)量，國(guó)內(nèi)學(xué)者［3－6］在前人研究基礎(chǔ)上，通過(guò)實(shí)驗(yàn)修正了Bagnold輸沙率公式，提出最大可能輸沙量概念，表示輸沙率的理論極限值，并給出計(jì)算公式，但兩者皆無(wú)法驗(yàn)證且無(wú)法直觀表示區(qū)域風(fēng)沙活動(dòng)強(qiáng)度。蝕積量和蝕積強(qiáng)度是闡明風(fēng)成地貌發(fā)育以及防治風(fēng)沙災(zāi)害的基礎(chǔ)，可以直觀地表示區(qū)域風(fēng)沙活動(dòng)強(qiáng)度，但目前風(fēng)沙蝕積量的計(jì)算較多采用輸沙率間接計(jì)算，誤差較大［7－10］。為此，本文通過(guò)野外觀測(cè)，參考輸沙勢(shì)和最大可能輸沙量公式，得出平坦沙質(zhì)地表蝕積量的計(jì)算模型，希望能夠?yàn)樯衬臄U(kuò)張與收縮提供量化預(yù)測(cè)依據(jù)。

1 觀測(cè)區(qū)與試驗(yàn)設(shè)計(jì)

騰格里沙漠位于阿拉善高原東南部，地理坐標(biāo)為37°54′—39°33′N，103°52′—105°36′E，沙漠面積3.67×104km2，是中國(guó)第4大沙漠。以流動(dòng)沙丘為主，兼有一定比例的固定、半固定沙丘，沙漠內(nèi)部湖泊較多，多數(shù)已經(jīng)干涸并演化為鹽湖。風(fēng)向以NW和EN為主，沙丘類型以格狀沙丘面積最大，由NE—SW走向的主梁和與之垂直的SE—NW的副梁為主，其次為新月形沙丘。沙丘高度變化幅度較小，以3～20m為主，沙丘間距較小，為30～170m［11－13］。

試驗(yàn)在中國(guó)科學(xué)院寒區(qū)旱區(qū)環(huán)境與工程研究所吊坡梁野外觀測(cè)場(chǎng)完成。觀測(cè)場(chǎng)建立在騰格里沙漠東南緣，中衛(wèi)市的西北角（37°32′N，105°02′E），觀測(cè)場(chǎng)內(nèi)經(jīng)推土機(jī)平整，觀測(cè)場(chǎng)規(guī)格為500m×400m，觀測(cè)場(chǎng)中心安裝一座觀測(cè)塔，塔高50m，風(fēng)速風(fēng)向傳感器安裝在距地面12m高度處，傳感器數(shù)據(jù)采集間隔為1min。本文參考前人的實(shí)驗(yàn)結(jié)果，定義沙粒啟動(dòng)的風(fēng)速為5.9m／s。

平坦沙質(zhì)地表試驗(yàn)區(qū)位于觀測(cè)塔旁側(cè)，面積規(guī)格為10m×10m，蝕積量測(cè)定采用測(cè)釬法，以正東和正北為正方向設(shè)置觀測(cè)點(diǎn)，每隔1m設(shè)置1根測(cè)釬，供布設(shè)121根測(cè)釬。每根測(cè)釬高度為80cm，初露高度統(tǒng)一為40cm。

2 蝕積量計(jì)算方法

先對(duì)121根測(cè)釬的出露高度按照以下公式進(jìn)行處理：

式中：h——測(cè)量周期內(nèi)蝕積深度（cm）；h1——沙下測(cè)釬深度（cm）；h0——沙下測(cè)釬初始深度（cm）。

運(yùn)用Table Curve 2Dv 5.01軟件，采用最小二乘法對(duì)h進(jìn)行函數(shù)擬合f（x），用決定系數(shù)（R2＞0.98）控制，用 MATLAB對(duì)f（x）進(jìn)行積分，求出風(fēng)沙蝕積量Q（cm3／m2）。當(dāng)Q＞0時(shí)，表示風(fēng)沙堆積；Q＜0，表示風(fēng)蝕；Q＝0，表示蝕積平衡［7］。

3 風(fēng)沙活動(dòng)特征

從圖1可知，觀測(cè)期間內(nèi)，起沙風(fēng)以E和NW風(fēng)為主，與前人［14－15］觀測(cè)結(jié)果一致，分別占起沙風(fēng)頻率的21%和22.7%，合成起沙方向?yàn)?0.6°，合成輸沙方向?yàn)?90.6°。其次起沙風(fēng)頻率較大者為NNW，ESE，WNW 和 ENE，頻率分別為10.4%，7.7%，7.3%和6.8%。其余風(fēng)向頻率則較小，起沙風(fēng)頻率之和僅為24.1%。進(jìn)一步比較偏東方向（ENE，E，ESE）和偏西北方向可知，偏西北方向頻率較偏東方向大4.9%。

圖1 研究區(qū)起沙風(fēng)分布情況

通過(guò)對(duì)測(cè)量期內(nèi)風(fēng)速資料的統(tǒng)計(jì)可知，風(fēng)速大于5.9m／s的起沙風(fēng)頻數(shù)占總頻數(shù)的11.69%。在大于5.9m／s起沙風(fēng)速中，以5.9～6.9m／s的風(fēng)速為主，占全部起沙風(fēng)頻數(shù)的58.12%；其次為7.0～7.9m／s的風(fēng)速，占24.94%；8.0～8.9m／s的風(fēng)速占12.35%；三者之和占95.4%，而大于9.0m／s的風(fēng)速僅占4.6%?？梢姡L(fēng)速越高，頻率越低，且觀測(cè)期間起沙風(fēng)以5.9～9.0m／s風(fēng)速為主。進(jìn)一步分析不同時(shí)段起沙風(fēng)頻率分布圖（圖2）可知，起沙風(fēng)頻率在各個(gè)統(tǒng)計(jì)時(shí)期內(nèi)差別較大，其中7月30日至8月30日內(nèi)最大，占總起沙風(fēng)頻數(shù)的36.68%，風(fēng)沙活動(dòng)強(qiáng)度較大，其次為6月4日至7月6日，占總起沙風(fēng)頻數(shù)的31.86%，風(fēng)沙活動(dòng)強(qiáng)度其次，最小為10月30日至11月15日，占1.09%，其余觀測(cè)期間起沙風(fēng)頻率差別較小，即風(fēng)沙活動(dòng)強(qiáng)度相差較小。

圖2 研究區(qū)不同時(shí)間段起沙風(fēng)頻率分布

在觀測(cè)期內(nèi)，觀測(cè)場(chǎng)平均速度為3.10m／s，小于起沙風(fēng)風(fēng)速平均值（6.60m／s）。且各統(tǒng)計(jì)期間內(nèi)平均風(fēng)速均小于起沙風(fēng)，其中以7月30日至8月30日平均風(fēng)速最大，為4.48m／s；10月31日至11月15日平均風(fēng)速最小，為1.59m／s。平均起沙風(fēng)速最大時(shí)期為8月31日至9月31日，為7.06m／s；平均起沙風(fēng)風(fēng)速最小時(shí)期是10月31日至11月15日，平均起沙風(fēng)為6.05m／s。

4 蝕積過(guò)程

在平坦沙質(zhì)地表觀測(cè)初期，沙地整體變化趨勢(shì)較為緩慢，監(jiān)測(cè)區(qū)內(nèi)測(cè)點(diǎn)高度變化較小，在沙地西部和西北部，出現(xiàn)輕微的剝蝕現(xiàn)象，而沙地局部出現(xiàn)堆積現(xiàn)象，形成小塊沙包。由圖3可見，7月7日至7月15日，沙包分別向東北和西南方向移動(dòng)，中間出現(xiàn)NW—ES走向的輕微風(fēng)蝕凹槽。7月16日至7月29日，風(fēng)蝕凹槽進(jìn)一步發(fā)育，凹槽底部風(fēng)蝕程度進(jìn)一步加劇，兩側(cè)堆積高度增加。7月30日至8月30日，風(fēng)蝕凹槽逐漸發(fā)育成熟，凹槽底部風(fēng)蝕程度和兩側(cè)堆積高度變化較小，而凹槽NW—ES走向明顯，但凹槽風(fēng)蝕程度由西北向東南方向逐漸減弱。8月31日至9月15日，風(fēng)蝕凹槽發(fā)育成熟，風(fēng)蝕等值線無(wú)明顯波動(dòng)，堆積和風(fēng)蝕相對(duì)穩(wěn)定，試驗(yàn)區(qū)域地勢(shì)無(wú)明顯變化。9月16日至9月30日，在風(fēng)蝕凹槽中心產(chǎn)生一個(gè)小型風(fēng)蝕坑，風(fēng)蝕凹槽再次產(chǎn)生變化。10月1日至10月15日，在原小型風(fēng)蝕坑的位置處出現(xiàn)堆積現(xiàn)象，形成一塊相對(duì)凹槽前后地勢(shì)較高的平地。10月16日至11月15日，在西北風(fēng)向影響下，堆積平地繼續(xù)發(fā)育，平地面積增大，高度增加，且向東南風(fēng)向移動(dòng)。

圖3 平坦沙質(zhì)地表蝕積動(dòng)態(tài)變化

5 模型計(jì)算

5.1 基礎(chǔ)模型

假設(shè)研究區(qū)域的風(fēng)沙環(huán)境為完全的開放系統(tǒng)，并與周邊環(huán)境存在風(fēng)沙物質(zhì)交換，則研究區(qū)域的蝕積狀況取決于該研究區(qū)域內(nèi)沙物質(zhì)的輸入和輸出。假定研究區(qū)域在某一時(shí)間段內(nèi)風(fēng)沙物質(zhì)的輸入量為M0，輸出量為Me，則該區(qū)域沙物質(zhì)蝕積狀況可以表示為：

式中：We——特定時(shí)段內(nèi)風(fēng)沙蝕積量（cm3）；M0——特定時(shí)段內(nèi)風(fēng)沙物質(zhì)輸入量（cm3）；Me——特定時(shí)段內(nèi)風(fēng)沙物質(zhì)輸出量（cm3）。

當(dāng)We＞0，物質(zhì)輸入量大于輸出量，為風(fēng)積區(qū)；當(dāng)We＜0，物質(zhì)輸入量小于輸出量，為風(fēng)蝕區(qū)；當(dāng)We＝0，物質(zhì)輸入量等于輸出量，為風(fēng)沙蝕積平衡區(qū)［7］。

5.2 模型推導(dǎo)與計(jì)算

數(shù)學(xué)模型參考Fryberger輸沙勢(shì)公式和凌裕泉的最大可能輸沙量的計(jì)算公式［17］。

假設(shè)：

式中：Q——蝕積量，特定時(shí)間段內(nèi)單位面積蝕積量（cm3／m2）；v——大 于 臨 界 啟 動(dòng)值的 風(fēng) 速 （m／s）；vt——臨界啟動(dòng)風(fēng)速，這里取5.9m／s，高度12m；f——對(duì)應(yīng)時(shí)間段內(nèi)大于臨界啟動(dòng)風(fēng)速的頻數(shù)；t——測(cè)量時(shí)間間隔（min）；a，b——待定常數(shù)。使用最小二乘法進(jìn)行計(jì)算，計(jì)算過(guò)程如下：

求得Q極值時(shí)a和b的值，進(jìn)而確定最優(yōu)值。由以上解法可得：

5.3 模型驗(yàn)證

5.3.1 蝕積量與輸沙勢(shì)比較 輸沙勢(shì)是衡量區(qū)域風(fēng)沙活動(dòng)的重要標(biāo)志，應(yīng)用較廣［14－15，17－19］，而蝕積量是表明區(qū)域內(nèi)風(fēng)沙物質(zhì)在時(shí)段內(nèi)的輸入與輸出［7］。從三者趨勢(shì)圖（圖4a）中可以看出，野外實(shí)測(cè)值、模型預(yù)測(cè)值和輸沙勢(shì)三者之間變化趨勢(shì)存在同期波動(dòng)特征，但三者波動(dòng)幅度不同。三者同時(shí)在8月30日到達(dá)峰值，分別為：1 7057.8，1 8471.2cm3／m2和184.20VU；在11月15日取得最小值429.75，1 011.5cm3／m2和1.76VU。其中，在10月16日至10月30日，由于接近臨界啟動(dòng)風(fēng)速的頻數(shù)較多，野外實(shí)測(cè)值在10月30日出現(xiàn)較大波動(dòng)，而依據(jù)氣象數(shù)據(jù)進(jìn)行計(jì)算的模型預(yù)測(cè)值和輸沙勢(shì)波動(dòng)較小，與野外實(shí)測(cè)值存在一定差距。圖4b進(jìn)一步說(shuō)明相同時(shí)段內(nèi)平坦沙質(zhì)地表野外實(shí)測(cè)值、模型預(yù)測(cè)量和輸沙勢(shì)三者之間的關(guān)系，三者存在明顯的線性關(guān)系。比較輸沙勢(shì)、模型預(yù)測(cè)值與野外實(shí)測(cè)值的回歸方程（圖4b，其中x表示實(shí)測(cè)蝕積量，y1和y2分別表示預(yù)測(cè)值和輸沙勢(shì)），可知模型預(yù)測(cè)值回歸方程決定系數(shù)較大，擬合度更高。由此可以判定，模型預(yù)測(cè)值能夠更好地表示平坦沙質(zhì)地表蝕積量，其值更加貼近野外實(shí)測(cè)值。

圖4 蝕積量實(shí)測(cè)值、預(yù)測(cè)值和輸沙勢(shì)趨勢(shì)與回歸分析

5.3.2 蝕積量與最大可能輸沙量比較 最大可能輸沙量是指當(dāng)沙源充足，沙面平坦裸露時(shí)，氣流充分作用于流沙表面所具有對(duì)沙物質(zhì)的最大可能搬運(yùn)能力，指輸沙量的理論極限值。其計(jì)算公式為：

式中：Q——最大可能輸沙量（m3／m）；f——對(duì)應(yīng)時(shí)間段內(nèi)大于臨界啟動(dòng)風(fēng)速的頻數(shù)；t——測(cè)量時(shí)間間隔（min）［6］。

由最大可能輸沙量趨勢(shì)圖（圖5）可知：最大可能輸沙量、模型預(yù)測(cè)值和野外實(shí)測(cè)值之間關(guān)系特征與輸沙勢(shì)相同，即存在同期波動(dòng)特征，三者同時(shí)在8月30日達(dá) 到 峰 值，分 別 為 1 690.9m3／m，18 471.2 和17 057.8cm3／m，但最大可能輸沙量與模型預(yù)測(cè)值之間的波動(dòng)趨勢(shì)更相近，主要原因是兩者均依據(jù)氣象數(shù)據(jù)進(jìn)行計(jì)算。進(jìn)一步分析三者回歸方程可知，模型預(yù)測(cè)值與野外實(shí)測(cè)值之間的回歸方程決定系數(shù)較最大可能輸沙量回歸方程大，可見相比最大可能輸沙量，模型預(yù)測(cè)值能夠更加準(zhǔn)確地表示平坦沙質(zhì)地表蝕積量。

圖5 蝕積量實(shí)測(cè)值、預(yù)測(cè)值和最大可能輸沙量趨勢(shì)與回歸分析

6 結(jié)論

通過(guò)野外布設(shè)試驗(yàn)進(jìn)行觀測(cè)，并在輸沙勢(shì)和最大可能輸沙量計(jì)算公式的基礎(chǔ)上修改參數(shù)，提出了平坦沙質(zhì)地表蝕積量的計(jì)算模型，彌補(bǔ)兩者無(wú)法直觀預(yù)測(cè)風(fēng)沙活動(dòng)的不足，同時(shí)也避免依據(jù)輸沙率進(jìn)行預(yù)測(cè)蝕積量而產(chǎn)生的較大誤差。通過(guò)比較輸沙勢(shì)、最大可能輸沙量、模型預(yù)測(cè)值和野外實(shí)測(cè)蝕積強(qiáng)度之間趨勢(shì)圖和回歸方程，得出模型預(yù)測(cè)值回歸方程更顯著，能夠更準(zhǔn)確表示平坦沙質(zhì)地表蝕積量，更貼近野外實(shí)測(cè)值，更直觀表示區(qū)域風(fēng)沙活動(dòng)狀況。研究結(jié)果可為沙漠?dāng)U張和收縮提供理論依據(jù)。

［1］ Fryberger S G，Dean G.Dune forms and wind regime［C］∥McKee E D.A Study of Global Sand Seas.US Geological Survey Professional Pape，1979，1052：137-169.

［2］ Lettau K，Lettau H.Experimental and micrometeorologicalfield studies of dune migration［C］∥Lettau K，Lettau H.Exporing the World’s Driest Climate.University of Wisconsin-Madison，IES Report，1978，101：110-147.

［3］ 李振山，倪晉仁.挾沙氣流輸沙率研究［J］.泥沙研究，2001（1）：1-10.

［4］ Dong Zhibao，Liu Xiaoping，Wang Hongtao.Aeolian sand transport：A wind tunnel model［J］.Sedimentary Geology，2003，161（1）：71-83.

［5］ 凌裕泉.輸沙量水平分布的非均一性［J］.實(shí)驗(yàn)力學(xué)：1994，9（4）：352-356

［6］ 凌裕泉.最大可能輸沙量的工程計(jì)算［J］.中國(guó)沙漠，1997，17（4）：362-368.

［7］ 劉連友.區(qū)域蝕積量和蝕積強(qiáng)度初步研究：以晉陜蒙接壤區(qū)為例［J］.地理學(xué)報(bào)，1999，54（1）：59-68.

［8］ Hu Xiao，Liu Lianyou，Li Shunjiang.Estimation of sand transportation rate for fixed and semi-fixed dunes using meteorological wind data［J］.Pedosphere，2009，19（1）：129-136.

［9］ Liu Lianyou，Skidmore E，Hasi E.Dune sand transport as influenced by wind directions，speed and frequencies in the Ordos Plateau［J］.Geomorphology，2005，67（3）：283-297.

［10］ Han Zhiwen，Dong Zhibao，Wang Tao.Observations of several characteristics of aeolian sand movement in the Taklimakan Desert［J］.Science in China（D）：Earth Sciences，2004，47（1）：86-96.

［11］ 閆滿存，董光榮，李華梅，等.騰格里沙漠東南緣沙漠演化的初步研究［J］.中國(guó)沙漠，1998，18（2）：111-117.

［12］ 凌裕泉，金炯，甄計(jì)國(guó).騰格里沙漠東南緣的風(fēng)沙活動(dòng)規(guī)律［M］∥騰格里沙漠沙坡頭地區(qū)流沙治理研究.銀川：寧夏人民出版社，1988：309-316.

［13］ 哈斯，董光榮，王貴勇.騰格里沙漠東南緣格狀沙丘的形態(tài)：動(dòng)力學(xué)研究［J］.中國(guó)科學(xué)：D輯，1999（29）：466-471.

［14］ 張正偲，董治寶，趙愛國(guó)，等.沙漠地區(qū)風(fēng)沙活動(dòng)特征［J］.干旱區(qū)研究，2007，24（4）：550-555.

［15］ 張克存，屈建軍，俎瑞平，等.騰格里沙漠東南緣輸沙勢(shì)與最大可能輸沙量之比較［J］.干旱區(qū)地理，2008，31（5）：643-648.

［16］ 黃富祥，牛海山，王明星，等.毛烏素沙地植被覆蓋率與風(fēng)蝕輸沙率定量關(guān)系［J］.地理學(xué)報(bào)，2001，56（6）：700-710.

［17］ 李紅軍，何清，楊青.近40年新疆輸沙勢(shì)的分析［J］.中國(guó)沙漠，2004，24（6）：706-710.

［18］ Fryrear D W.Wind erosion：field measurement and analysis［J］.Transactions of the ASAE，1991，34（1）：155-160.

［19］ Bagnold R A.風(fēng)沙和荒漠沙丘物理學(xué)：中譯本［M］.北京：科學(xué)出版社，1959：168-171.