基于貝葉斯模式平均與區間二型模糊集的分區灌溉方法

邵東國，鄒亮峰，顧文權，農翕智，王 鶴，王 柏

（1.武漢大學水資源與水電工程科學國家重點實驗室，湖北武漢 430072；2.黑龍江省水利科學研究院，黑龍江哈爾濱 150080）

1 研究背景

松嫩平原是我國重要的糧食產區，受氣候變化與糧食市場的影響，松嫩平原水稻種植面積不斷擴大［1-2］，水田灌溉用水占農業用水80%以上［3］，局部地區出現地下水超采現象，對當地河湖濕地生態系統構成嚴重威脅［4］。因此，合理優化水稻灌溉制度，推行節水控采，緩解地下水超采現象，對保護當地糧食生產與河湖濕地生態系統穩定具有重要意義。

作物灌水量與地下水埋深之間存在一定相關性［5-6］。在干濕交替稻田中，地下水通過毛管上升成為可供水稻根系吸收利用的土壤水，從而影響根區土壤水分分布［7-9］。地下水埋深的變化將改變稻田地表水-土壤水-地下水之間的轉化運動規律，進而影響稻田滲漏與作物需水量［10］。現有水稻灌溉制度優化研究，大多基于作物耗水與田面水層等試驗監測數據，通過多目標遺傳算法與水分生產函數、水量平衡等模型耦合，忽略了稻田地下水埋深對灌水量的影響［11-14］。實際上，隨著東北等地大規模集約化水稻生產的推進和降雨、地下水開采等非均勻分布影響，水稻種植區存在明顯的地下水埋深時空變化［15-17］。因此，需要根據水稻種植區不同地下水埋深對稻田滲漏與水稻作物需水等的影響差異，通過不同區域地下水埋深差異分區來進一步優化水稻灌溉制度，最大限度地提高毛管水的利用率，減少稻田深層滲漏量，實現大規模集約化水稻種植區域的高效節水灌溉。

灌溉制度優化研究中，產量模擬的精度決定了優化結果的可靠性［18］。AquaCrop 作物模型內置地下水模塊［19］，能模擬動態地下水埋深條件，所需參數少且便于多情景模擬，但產量模擬精度較差［20-22］。作物水分生產函數模型根據當地田間試驗數據確定，在產量模擬方面更具針對性，模擬精度更高，但缺少對田間地下水影響的考慮［23］。因此，需要深入研究地下水埋深動態變化條件下，提高作物產量模擬精度的有效方法。

考慮到現有灌溉管理決策常憑經驗或用水計劃，忽略了農業生產中氣候、灌溉等農業活動的隨機性，以及決策者的風險偏好與判斷模糊性，導致灌溉制度優化結果過于理想化，而經驗性灌溉制度又常發生灌溉退水現象，造成灌溉水資源浪費［24-26］。如何綜合考慮農業集約化生產過程中灌溉決策的不確定性，結合決策者主觀風險傾向進行灌溉制度優化決策，是提高灌溉制度優化結果實用性與有效性的重要途徑。

本文以松嫩平原北部和平灌區水稻種植區2017-2019年降雨-地下水-產量監測試驗為基礎，通過分析地下水埋深變化規律，基于空間拓撲關系聚類分析與空間疊加分析確定研究區地下水埋深的分區界限；設定AquaCrop 模型［19，22］中不同區域地下水埋深模塊，運用AquaCrop 和Jensen 模型［23］對平、枯水年各區不同灌溉情景進行產量模擬，并用貝葉斯模式平均方法（BMA）［27-29］融合兩個模型的產量模擬系列，基于區間二型模糊集理論［30-32］從節水、增產與地下水穩定三個屬性對灌溉情景進行模糊綜合評價，結合決策者風險偏好，采用有序加權算術平均算子（OWA）集結各典型年不同分區所有灌溉情景評價值，實現灌溉制度模糊綜合優化決策，以期為水稻規模化種植區提供基于地下水埋深動態的分區節水控采新方法。

2 灌溉區域地下水埋深動態分區

2.1 地下水動態監測試驗及分區方法

2.1.1 地下水埋深監測試驗 試驗于2017—2019年在松嫩平原北部（黑龍江省慶安縣）和平灌區水稻灌溉試驗站開展，該區域為平原區，地勢平坦，屬于北溫帶大陸性季風氣候，多年平均降雨量558 mm，多年平均蒸發量764.5 mm，多年平均氣溫2.5 ℃，年內氣溫變幅較大，夏季極端最高氣溫36.7 ℃，冬季極端最低氣溫-44.9 ℃，無霜期128 d。試驗區位于寒地黑土核心區，土壤種類為白漿型水稻土，土壤容重為1.02 g/cm3，土壤基本理化性質為：pH值6.40，全氮15.10 g/kg，全磷15.21 g/kg，全鉀20.09 g/kg，有機質41.5 g/kg。區域內設有7 眼地下水觀測井，根據試驗區地形條件，均勻分布在主要排水口附近，在生育期內采用地下水監測系統逐日觀測地下水埋深（Groundwater depth，GD），具體分布如圖1所示。

圖1 地下水監測井分布

2.1.2 分區方法 采用線性函數歸一化方法，對生育期內地下水埋深監測數據進行歸一化處理。基于K-Means聚類算法對歸一化后的埋深數據進行聚類分析，確定灌溉區域根據地下水埋深分區的節點。K-Means 聚類算法的思路為，使樣本空間所有的非中心點到各自所屬簇的中心點的距離的平方和最小，目標函數F為：

式中：Xnorm為歸一化后數據；X為原始數據；Xmax、Xmin分別為原始數據集的最大值與最小值；C1、C2、…、Ck分別是k個簇的中心點；C（Xi）為Xi點所屬的簇的中心點；d2（Xi，C（Xi））為求兩點的的距離平方。

根據聚類分析結果，基于ArcGIS平臺采用Kriging插值方法進行多次插值，通過分析地理實體的拓撲關系，對插值結果進行空間疊加分析，最終實現區域分區。

2.2 地下水動態分析及分區結果

2.2.1 地下水埋深動態分析 根據2017、2018年地下水監測數據，計算生育期平均地下水埋深（Av?erage groundwater depth，AGD），結果見圖2。從圖2 可以看出，2017年AGD在1.25～2.12 m 之間變化，變幅為0.87 m；2018年AGD在0.65～2.03 m之間變化，變幅為1.38 m；2018年為特豐水年，降雨量較大，因此地下水埋深波動更加劇烈，2018年AGD變幅相比2017年增加了58.6%。

圖2 生育期區域平均地下水埋深與降雨量關系

受降雨、灌溉等因素影響，水稻不同生育期內地下水埋深不同。分蘗期以灌溉為主，降雨量少，GD在1.55～2.03 m 之間；孕穗期、抽穗期降雨豐富且集中，因此GD逐漸減小，2017、2018年GD分別為1.32 ～ 1.76 m 和0.65 ～ 1.49 m；乳熟期降雨、灌溉減少，使得區域GD增加，GD在1.52～1.85 m之間。

2.2.2 地下水埋深分區結果 根據2017年、2018年監測數據對地下水埋深進行頻數分布統計，其結果如圖3所示，GD近似滿足正態分布。調用R語言中kmeans（）函數，設置類別數k=2，分別對2017年、2018年GD數據（一維數據）進行聚類分析，結果如表1所示，將GD監測數據分為兩簇后，簇1 以1.8 m 左右為數據中心，簇2 以1.48 ～1.17 m為數據中心。

圖3 地下水埋深監測數據頻數分布直方圖

表1 地下水埋深數據聚類分析后各簇數據質心結果

土壤水-地下水垂向運動規律研究表明，地下水埋深增加將導致稻田滲漏量增大，作物對地下水的利用率降低［33］。GD>2.8 m后，垂向運動以深層滲漏為主，地下水無法補給作物需水［9］；GD<0.7 m 時，地下水蒸發量較大，易發生土壤鹽漬化，阻礙作物生長及其根系發育［11］；埋深1.5 m介于兩簇質心之間，與聚類分析結果相契合，以此為分區界限可以兼顧深層滲漏與毛管上升對稻田灌水量的影響。基于GD聚類分析結果，結合土壤水-地下水垂向運動規律，最終確定分區灌溉的地下水埋深分區界限為1.5 m。

以埋深1.5 m為界限，結合2017、2018年水稻不同生育期區域地下水埋深的多次Kriging插值結果，對研究區進行初步分區。根據地理實體的拓撲關系，采用空間疊加分析方法對初步分區結果進行相交分析，最終分區結果如圖4所示（Ⅰ區表示GD≤1.5 m的區域，Ⅱ區為GD>1.5 m的區域），為確定不同地下水埋深區域的水稻最優灌溉模式、保證地下水可持續利用奠定了基礎。

圖4 結合聚類分析、多次Kriging插值與相交分析對灌域地下水埋深分區結果

3 基于貝葉斯模式平均方法的分區產量模擬模型

3.1 產量模擬模型為深入揭示灌區地下水埋深動態對水稻產量的影響，利用AquaCrop模型地下水模塊［19］模擬不同地下水埋深分區的作物灌水與產量關系，但產量模擬精度次于水分生產函數模型。因此，為提高產量模擬精度，以AquaCrop模型中各區地下水動態埋深為控制條件，將該模型模擬所得水稻蒸發量結果作為Jensen 水分生產函數模型［23］的輸入，得到不同地下水埋深區域的產量模擬值，通過貝葉斯模式平均（BMA）方法［29］將AquaCrop與Jensen模型產量模擬結果融合，得到精度更高的整合產量系列。

3.1.1 AquaCrop作物模型 AquaCrop模型是基于水分驅動的日尺度作物生長模擬模型，主要根據實際蒸騰量與歸一化水分生產率計算生物量，利用收獲指數和生物量的乘積確定最終產量。計算公式為：

式中：B為生物量，t/hm2；Tr為實際蒸騰量，mm；WP*為歸一化水分生產率；HI為收獲指數，%；Y為產量，t/hm2。

模型運行需要輸入：氣象參數、作物參數、田間管理參數、土壤參數［19］。AquaCrop模型地下水模塊中設定不同地下水埋深調節土壤水分脅迫程度。

3.1.2 水分生產函數模型 選取Jensen水分生產函數模型描述水稻不同生育期水分-產量關系，可表示為：

式中：Ya、Ym分別為作物全生育期內的實際產量和最大產量，kg/hm2；ETa、ETm分別為全生育期作物實際騰發量和最大騰發量，mm；λi為i生育階段缺水對作物產量影響的敏感性指數，即水分敏感指數；n為生育階段數。

將其轉化為多元線性方程，采用多元線性回歸分析方法求解各系數，得到適合研究區的Jensen模型［24］為：

式中下標1、2、3、4、5、6分別代表分蘗前中后期、拔節期、抽穗期與乳熟期。

3.1.3 貝葉斯模式平均（BMA）方法 BMA方法［28］是一種集合不同模型模擬值，得到更可靠的綜合模擬值的統計后處理方法。假設S為產量模擬量，R=［X，Y］為輸入數據（其中X代表各模型模擬產量，Y代表實測產量），是K個模型模擬的集合，基于貝葉斯全概率公式可得模擬量S的概率密度函數為：

BMA方法采用后驗概率作為權重，對各模型模擬值進行加權平均，最終輸出綜合模擬結果，精度較高的模型權重值更大。若單個模型模擬值與實測值均符合正態分布，可用BMA方法的校正公式：

3.2 灌溉情景設定根據慶安氣象站歷史數據，采用優化適線法選取Pearson-Ⅲ型曲線進行擬合。選取1992年（P=50%）為典型平水年，2007年（P=85%）為典型枯水年，各典型年份生育期累計降雨量分別為426和274 mm。

根據灌區地下水埋深分區結果（Ⅰ區、Ⅱ區GD數值）設定AquaCrop地下水模塊動態GD情景。通過增加或減少灌水次數、加大或降低灌溉定額的方式，基于解空間動態縮減策略［34］共組合生成576種灌溉模擬情景，組合方式見表2。

表2 灌溉情景設定

采用AquaCrop 模型對各典型年不同分區所有灌溉情景進行蒸發量-產量模擬，同時基于模擬所得蒸發量計算各灌溉情景下Jensen模型產量值。

4 基于區間二型模糊集的灌溉制度風險決策模型

4.1 基于IT2FS 的模糊決策模型在灌溉制度優選中，受管理者主觀因素影響，常存在一定模糊性。為提高優選結果的合理性，本文基于區間二型模糊集（IT2FS）［32］建立了灌溉制度模糊決策模型，主要采用效用函數（U（A））度量IT2FS的優劣性，以區間二型模糊集熵（E（A））度量IT2FS自身的不確定性［31］。同時引入風險偏好因子（θ）刻畫不同決策者的風險態度，各屬性權重由E（A）和θ共同確定。采用有序加權算術平均算子（OWA）集結各方案屬性權重與效用函數得到綜合效用值。

對于有n個方案集xi（0

（1）對各方案的所有屬性進行模糊評價，得到原始決策矩陣A=（aij）n×m（aij均是IT2FS，表示決策者對某方案的評價結果）對其進行規范化得到規范化決策矩陣D（dij）n×m。

（2）計算D矩陣的效用函數U（dij）n×m和區間二型模糊集熵E（dij）n×m。

（3）計算不同決策風險偏好時各屬性權重。

（4）根據屬性權重與效用函數，利用有序加權平均算子集結不同風險偏好下各方案綜合效用值U：

式中：bj為U（dij）第j行元素中第j大的元素；ωj為第j個屬性的權重。

4.2 決策模型屬性與標準設定根據產量模擬結果，從節水效益（C1）、增產效果（C2）和地下水埋深的穩定性（C3）三個屬性對各典型年不同地下水埋深區域所有灌溉情景進行評估。評估結果用語言進行標度，各語言標度對應的區間二型模糊集如表3所示。

表3 語言標度及其對應的區間二型模糊集

5 結果與分析

5.1 產量模擬模型率定驗證結果根據2017、2018年田間實驗數據［23，33］對AquaCrop模型、BMA法參數進行率定、求解。選取體積差DV［35-36］、納什系數NSE和決定系數R2評價模型模擬精度。

5.1.1 AquaCrop模型率定驗證結果 選取土壤含水率指標檢驗AquaCrop模型適用性，以2017年控灌（T1）、淺曬淺灌（T2）、淹灌（T3）處理數據率定模型，調整作物參數，率定后部分參數見表4。以2018年為驗證期，模擬得到土壤含水率與實測數據比較見如圖5，精度評價結果見表5。對比土壤含水率實測值和模擬值，實測點均落在模擬含水率曲線附近，模擬值正負體積差DV不超過3.8%，R2與NSE均在0.8以上，說明土壤含水率模擬值與實測值之間一致性強，AquaCrop模型適用性高。

表4 率定后的模型輸入參數

表5 含水率模擬評價指標

圖5 土壤含水率實測值與模擬值對比

5.1.2 BMA 參數估計與產量模擬精度結果 基于2017、2018年實測產量系列［23］、AquaCrop［19］與Jensen模型［23］模擬產量系列，采用期望最大化（EM）［27］算法并結合R語言BMS包求解BMA方法中各模型所占權重。結果顯示，AquaCrop模型所占權重為0.402，Jensen模型權重為0.598，表明Jensen模型產量模擬精度更高。

表6為貝葉斯模式平均（BMA）方法［29］及其2個模型產量模擬精度評價結果。從表6可知，AquaC?rop模型高估產量（DV<0），Jensen模型低估產量（DV>0），采用BMA方法可以調和兩個產量模擬系列；BMA方法的確定性系數（R2）均在95%以上，大于任一單個模型的R2；BMA方法的納什效率系數NSE最大，表明經過BMA方法融合后的產量模擬精度比單一模型的模擬精度高。

表6 BMA方法與單一模型產量模擬精度評價結果

2017、2018年實測產量與模擬產量對比圖6。從圖6可以發現，BAM方法產量模擬效果最優，產量的相關系數超過0.95，AquaCrop 模型效果最差。因此，在動態地下水埋深條件下，采用BMA方法［29］融合AquaCrop、Jensen模型的模擬產量能夠提高和平灌區產量模擬精度。

圖6 模擬產量與實測產量對比

5.2 分區灌溉模擬結果各典型年不同地下水埋深區域所有灌溉情景產量模擬結果見圖7（Ⅰ區表示GD≤1.5 m 的區域，Ⅱ區表示GD>1.5 m 的區域）。采用Mann-Kendall趨勢檢驗法，對各年型不同區域所有灌溉情景的產量進行趨勢檢驗，結果見表7。結果顯示，不同年型各區域產量趨勢檢驗值均小于-2.58（顯著性水平0.01），隨著灌水量增加，產量總體呈顯著下降趨勢，但枯水年產量變化趨勢與總體變化趨勢存在差異。

圖7 Ⅰ區、Ⅱ區不同灌溉情景產量-灌水量關系

表7 產量Mann-Kendall趨勢檢驗值

篩選相同灌水量下各典型年不同區域產量最大的灌溉情景，擬合產量-灌水量關系曲線，結果如圖8、表8所示。結果表明，水稻產量-灌水量響應規律呈非線性關系。灌水量較小時，土壤含水率較低，作物受到一定程度干旱脅迫，阻礙干物質累積與產量形成，此時增加灌水能夠降低水分脅迫程度，促進產量形成；枯水年水分脅迫更嚴重，產量形成受到抑制，增加灌水能顯著降低水分脅迫對產量形成的抑制作用。灌水量較大時，土壤常處于飽和狀態，土壤含氧量和透氣性低，導致水稻根系發育受阻，進而影響產量形成；灌水量越大，越容易產生深層滲漏，造成無效灌水。因此，隨灌水量增加，水分脅迫逐漸消除，產量為增長趨勢，而后影響根系呼吸進而降低產量，但產量整體仍然呈現隨著灌水量增加而降低的趨勢。

表8 各年型不同區域產量-灌水量擬合曲線方程

圖8 產量與灌水量擬合曲線

隨著灌水量增加，Ⅰ區最早出現產量降低現象，Ⅱ區較晚出現。原因是地下水埋深小時，地下水可通過毛管水上升供給水稻需水，此時土壤水分垂向運動為毛管水上升過程，增加灌水量會更早影響根系呼吸作用，最終導致產量降低；埋深較大時，滲漏增強，因此產量降低現象出現較晚。綜上所述，水稻產量與灌水量響應規律表現出高度非線性，且總體下降趨勢明顯，并隨地下水埋深增大表現出滯后現象，能夠為優化水稻分區灌溉制度提供理論依據。

5.3 灌溉制度風險決策結果根據式（11）分別計算各年型不同區域不同決策風險偏好下節水效益（C1）、增產效果（C2）和地下水穩定性（C3）三個屬性的權重值，基于期望效用函數理論，引入絕對風險回避系數A（W）衡量風險回避程度，A（W）>0為風險規避型；A（W）=0為風險中性型；A（W）<0為風險偏好型，計算結果如圖9所示。結果顯示，屬性權重未知時，不同風險態度會影響各屬性權重分布。隨著風險態度逐漸轉變為偏好型，地下水埋深較淺區域（Ⅰ區）C1屬性權重上升，C2屬性權重下降，C3屬性權重保持穩定；埋深較深區域（Ⅱ區）C1屬性權重降低，C2屬性權重略微降低，C3權重上升。

圖9 各典型年不同區域內不同風險偏好下各屬性權重分布

根據各年型各分區不同風險態度的屬性權重，結合效用函數，采用有序加權算術平均算子（OWA）對灌溉情景綜合效用值進行集結，選取綜合效用值最大的灌溉情景為最優灌溉制度，結果見表9。

表9 各年型不同地下水埋深區域的優化灌溉制度

優化后，相比傳統雨養種植，平水年每增加100 mm 灌水，水稻產量提高1100 kg/hm2以上，枯水年每增加100 mm 灌水，產量提高1000 kg/hm2以上，Ⅰ區水稻增產效果更顯著。通過灌溉措施可以緩解降雨時空分布不均導致的水稻水分虧缺，因此松嫩平原水稻分區灌溉模式具有較大灌水增產潛力，優化后的水稻灌溉制度能實現增產目標。相比傳統淹灌［37］，優化后的灌溉制度考慮了稻田地下水的補給作用與田間深層滲漏的影響，平水年平均節水22%以上，枯水年平均節水30%以上，在淺埋區（Ⅰ區）節水效果更顯著，可節水37.6%以上。相較于灌區現行灌溉制度［37］，平水年可節水15%左右，枯水年可達20%。松嫩平原當前種植結構下［1］，采用優化后的水稻灌溉制度至少可節水8.3億m3，減少地下水開采量12%以上。因此，考慮稻田地下水埋深進行分區優化后的水稻灌溉制度，在穩產的同時具有良好的節水控采效果，能有效緩解地下水超采現象，為當地河湖濕地生態系統安全提供保障。

6 結論

（1）通過2017—2019年松嫩平原北部和平灌區降雨-地下水-產量監測試驗與統計分析，揭示了灌區地下水時空動態變化規律，提出了稻田動態地下水埋深條件下的水稻灌溉制度分區優化方法，采用K-Means算法與空間疊加分析方法，確定了和平灌區稻田地下水埋深的分區界限為1.5 m。

（2）采用貝葉斯模式平均（BMA）方法融合AquaCrop、Jensen模型產量模擬系列，提高了動態地下水埋深條件下的產量模擬精度，BMA 產量模擬值與實測值的相關系數超過0.95。通過模擬不同水量分配情景，揭示了各典型年分區灌溉的水稻產量-灌水量差異性響應規律與灌溉增產潛力，為分區灌溉制度優化提供理論基礎。

（3）為提高灌溉制度優化結果的合理性，構建了基于區間二型模糊集理論的灌溉制度模糊風險決策模型，引入風險偏好因子刻畫決策者的風險態度，提出了適合不同決策者的灌域分區灌溉模式，該模式能夠減少稻田灌水量15%～37.6%，減少地下水開采量12%以上。

（4）由于地下水埋深時空差異性大，點尺度歷史地下水埋深數據存在局限性，地下水埋深分區結果的有效性受到挑戰，因此構建區域地下水埋深實時模擬、預測模型，提高分區結果的有效性，獲得更合理的水稻分區灌溉制度，是未來需要完善的重要內容。