本發(fā)明涉及一種基于貢獻(xiàn)度的科普勒(Copula)函數(shù)地下水位分析方法，屬于聯(lián)合概率分析技術(shù)領(lǐng)域。

背景技術(shù)：

水作為實(shí)現(xiàn)經(jīng)濟(jì)社會(huì)可持續(xù)發(fā)展的重要戰(zhàn)略性資源,具有不可替代的作用。而地下水資源作為其中重要的組成部分,不僅彌補(bǔ)了地表水資源時(shí)空分布不均勻引起的區(qū)域供水不足,同時(shí)也支撐了生態(tài)環(huán)境的和諧。地下水為國民經(jīng)濟(jì)發(fā)展和人民生活水平提高提供了水源保證,做出了巨大貢獻(xiàn),但近年來，由于無序開采,地下水嚴(yán)重超采,特別是在河北、山西等地，引發(fā)了地下水位持續(xù)下降及其一系列的水環(huán)境問題。為實(shí)施最嚴(yán)格的水資源管理制度,加強(qiáng)對(duì)地下水資源優(yōu)化配置,在及時(shí)掌握地下水動(dòng)態(tài)的前提下,了解其的變化十分必要,地下水預(yù)測(cè)分析工作也是當(dāng)前地下水研究的熱點(diǎn)。

最早對(duì)地下水的分析研究方法采用的是比較簡(jiǎn)單的水均衡方法和水文地質(zhì)比擬法，后來隨著計(jì)算機(jī)的普遍應(yīng)用，數(shù)值模擬法的使用日益廣泛。由于確定性模型對(duì)于復(fù)雜的水文地質(zhì)條件的模擬難以保證精度，隨機(jī)模型方法也應(yīng)運(yùn)而生，主要有回歸分析法、灰色模型和時(shí)間序列模型等。

地下水的補(bǔ)給來源主要是降雨，消耗主要是人類的開采，目前進(jìn)行地下水位分析研究具有局限性。首先，目前的研究并沒有考慮二者對(duì)于地下水位變化的貢獻(xiàn)度問題，僅僅將二者作為影響因素帶入模型進(jìn)行預(yù)測(cè)分析，缺少對(duì)權(quán)重的考慮；其次，并沒有從聯(lián)合分布的角度去探求地下水位變化補(bǔ)給消耗之間的關(guān)系，從概率學(xué)的角度給出可能存在的風(fēng)險(xiǎn)，為決策者提供依據(jù)。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是提出一種基于貢獻(xiàn)度的科普勒函數(shù)地下水位分析方法，以解決目前在進(jìn)行地下水位分析是缺少對(duì)不同因素對(duì)地下水位變化貢獻(xiàn)度的考慮，同時(shí)沒有從聯(lián)合概率分布的角度去研究地下水位和影響地下水位因素的補(bǔ)給消耗之間的概率關(guān)系。

本發(fā)明提出的基于貢獻(xiàn)度的科普勒函數(shù)地下水位分析方法，包括以下步驟：

(1)獲取待分析區(qū)域的相關(guān)資料，包括：從待分析區(qū)域所屬水文部門獲取待分析區(qū)域最近30年中各年平均地下水位資料A_年和最低地下水位資料A_低，從待分析區(qū)域所屬政府規(guī)劃部門獲取待分析區(qū)域最近30年中各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)分別的地下水使用量C₁、C₂、C₃和獲取待分析區(qū)域最近30年中各年地下水總使用量C，從待分析區(qū)域所屬水文部門獲取待分析區(qū)域最近30年中年降雨資料B_年，從待分析區(qū)域所屬政府管理部門獲取待分析區(qū)域最近30年中各年第一產(chǎn)業(yè)增加值D₁、第二產(chǎn)業(yè)增加值D₂和第三產(chǎn)業(yè)增加值D₃；

(2)根據(jù)上述步驟(1)獲取的數(shù)據(jù)，分別計(jì)算第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)的地下水消耗占地下水總消耗的權(quán)重α₁、α₂和α₃如下：

上式中，為上述各年地下水總使用量C的多年平均值，分別為各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)分別的地下水使用量C₁、C₂、C₃的多年平均值；

通過下式，計(jì)算得到地下水消耗指數(shù)W：

W_i＝D_1i×α₁+D_2i×α₂+D_3i×α₃ i＝1,2,3…..30

上式中，下標(biāo)i代表年份，W_i表示各年的地下水消耗指數(shù)，D_1i、D_2i、D_3i分別為各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)、第三產(chǎn)業(yè)的增加值，α₁、α₂、α₃分別為第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)地下水消耗占地下水總消耗的權(quán)重；

(3)根據(jù)上述步驟(1)獲取的數(shù)據(jù)，畫出年最低地下水位A_低和時(shí)間的線性關(guān)系圖，橫坐標(biāo)為年份，縱坐標(biāo)為各年最低地下水位值，從最低地下水位A_低和時(shí)間的線性關(guān)系圖中得到地下水位變化的拐點(diǎn)年份；

(4)根據(jù)步驟(1)獲取的數(shù)據(jù)，分別畫出平均地下水位資料A_年、年降雨資料B_年和地下水消耗指數(shù)W與年份之間的線性關(guān)系圖，橫坐標(biāo)為年份，縱坐標(biāo)分別為各年A_年、B_年和W，根據(jù)上述步驟(3)確定的地下水位拐點(diǎn)年份，將上述三個(gè)不同線性關(guān)系圖均劃分為地下水位拐點(diǎn)前和后兩部分，分別添加各部分的趨勢(shì)線，假定在拐點(diǎn)前后，地下水位A_年與年份之間線性關(guān)系斜率分別為S_Ab和S_Aa，降雨B_年與年份之間的線性關(guān)系斜率分別為S_Bb和S_Ba，地下水消耗指數(shù)W與年份之間的線性關(guān)系斜率分別為S_Wb和S_Wa，則降水量對(duì)地下水位變化的貢獻(xiàn)度C_B為：

C_B＝((S_Bb+S_Ba)/2)/((S_Ab+S_Aa)/2)

地下水消耗指數(shù)W對(duì)地下水位變化的貢獻(xiàn)度C_W為：

C_W＝-((S_Wb+S_Wa)/2)/((S_Ab+S_Aa)/2)；

(5)根據(jù)上述步驟(1)-步驟(4)，得到按照年份先后排列的三類數(shù)據(jù)：第一類為A_年系列，第二類為不同年份降雨量B_年與C_b的乘積B‘_年系列，第三類為不同年份地下水消耗指數(shù)W與C_g的乘積W'系列；

(6)根據(jù)上述步驟(5)得到的三類數(shù)據(jù)系列A_年、B‘_年、W'，分別計(jì)算每個(gè)系列的邊緣概率分布，選取在水文頻率分析中的多種分布線型分別進(jìn)行曲線擬合，然后對(duì)每一種分布線型進(jìn)行檢驗(yàn)，得到每一類數(shù)據(jù)的最優(yōu)擬合分布線型，比較分布線型的理論值與實(shí)際值之間的均方根R_RMSE，以均方根R_RMSE最小為原則選擇擬合度較好的分布函數(shù)：

(7)采用阿基米德型科普勒家族中的Clayton Copula和Gumbel Copula函數(shù)分別構(gòu)建三維聯(lián)合概率分布：

Frank Copula形式為：

Clayton Copula形式為：

式中，θ為copula函數(shù)的參數(shù)，u₁、u₂、u₃為邊際分布函數(shù)；

(8)選用離差平方和最小準(zhǔn)則法對(duì)Copula函數(shù)的擬合優(yōu)度進(jìn)行評(píng)價(jià)，選取離差平方和最小準(zhǔn)則法最小的copula作為優(yōu)選科普勒函數(shù)，離差平方和最小準(zhǔn)則法的計(jì)算公式為：

上式中，p_i分別為經(jīng)驗(yàn)頻率和理論頻率；i為樣本序號(hào)；

(9)根據(jù)選定的科普勒函數(shù)計(jì)算風(fēng)險(xiǎn)概率，公式為：

上式中，假定x₁為地下水位系列A_年，x₂為表征地下水消耗的W'系列，x₃為降雨量系列B‘_年系列，得到降雨量等于或者小于設(shè)定概率條件下，地下水位和地下水消耗W'相等的發(fā)生的概率，根據(jù)氣象部門對(duì)降雨量的預(yù)測(cè)，確定一個(gè)降雨發(fā)生的概率，得到經(jīng)濟(jì)增長(zhǎng)一定的值條件下相應(yīng)地下水位可能的值，為制定經(jīng)濟(jì)發(fā)展政策；或假設(shè)x₃為表征地下水消耗的W'系列，x₁為地下水位系列A_年，x₂為降雨量系列B‘_年系列，得到當(dāng)產(chǎn)業(yè)發(fā)展達(dá)到某一個(gè)程度時(shí)，相應(yīng)地下水位和降雨量之間同時(shí)達(dá)到設(shè)定的概率，根據(jù)未來政府的規(guī)劃，確定一個(gè)地下水消耗的W'的值，得到當(dāng)降雨達(dá)到設(shè)定概率值時(shí)，地下水位可能的值，為地下水保護(hù)政策提供依據(jù)。

上述基于貢獻(xiàn)度的科普勒函數(shù)地下水位分析方法的步驟5中，水文頻率分析中的多種分布線型分別為：P－Ⅲ型分布曲線、對(duì)數(shù)正態(tài)分布曲線和廣義極值分布曲線，P－Ⅲ型分布概率密度公式為：

其中，α為形狀參數(shù)，β為尺度參數(shù)，a₀為位置參數(shù)。

廣義極值分布曲線概率密度公式為：

其中，α為尺度參數(shù)，k為形狀參數(shù)，u為位置參數(shù)。

對(duì)數(shù)正態(tài)分布曲線概率密度公式為：

其中μ_y、δ_y分別為x序列取自然對(duì)數(shù)后形成的序列的均值和標(biāo)準(zhǔn)差。

本發(fā)明提出的基于貢獻(xiàn)度的科普勒函數(shù)地下水位分析方法，其優(yōu)點(diǎn)如下：

(1)本發(fā)明選取降水作為地下水補(bǔ)給因素，通過不同產(chǎn)業(yè)的地下水消耗量設(shè)置權(quán)重計(jì)算地下水消耗指數(shù)作為地下水消耗因素，通過線性斜率分析確定補(bǔ)給和消耗對(duì)地下水的貢獻(xiàn)率，直觀的反映了影響地下水位的不同因素之間的關(guān)系；

(2)通過建立三者的聯(lián)合概率分布，從概率學(xué)的角度分析了地下水位變化和消耗補(bǔ)給之間的存在的概率關(guān)系，消除了模型方法無法刻畫地下水動(dòng)態(tài)變化的缺陷，通過耦合氣候變化模型對(duì)降雨的預(yù)測(cè)或區(qū)域經(jīng)濟(jì)發(fā)展規(guī)劃，分析補(bǔ)給消耗對(duì)地下水水位的影響，為決策者提供依據(jù)。

具體實(shí)施方式

本發(fā)明提出的基于貢獻(xiàn)度的科普勒函數(shù)地下水位分析方法，包括以下步驟：

(1)獲取待分析區(qū)域的相關(guān)資料，包括：從待分析區(qū)域所屬水文部門獲取待分析區(qū)域最近30年中各年平均地下水位資料A_年和最低地下水位資料A_低，從待分析區(qū)域所屬政府規(guī)劃部門獲取待分析區(qū)域最近30年中各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)分別的地下水使用量C₁、C₂、C₃和獲取待分析區(qū)域最近30年中各年地下水總使用量C，從待分析區(qū)域所屬水文部門獲取待分析區(qū)域最近30年中年降雨資料B_年，從從待分析區(qū)域所屬政府管理部門獲取待分析區(qū)域最近30年中各年第一產(chǎn)業(yè)增加值D₁、第二產(chǎn)業(yè)增加值D₂和第三產(chǎn)業(yè)增加值D₃；

(2)根據(jù)上述步驟(1)獲取的數(shù)據(jù)，分別計(jì)算第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)的地下水消耗占地下水總消耗的權(quán)重α₁、α₂和α₃如下：

上式中，為上述各年地下水總使用量C的多年平均值，分別為各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)分別的地下水使用量C₁、C₂、C₃的多年平均值；

通過下式，計(jì)算得到地下水消耗指數(shù)W：

W_i＝D_1i×α₁+D_2i×α₂+D_3i×α₃ i＝1,2,3…..30

上式中，下標(biāo)i代表年份，W_i表示各年的地下水消耗指數(shù)，D_1i、D_2i、D_3i分別為各年第一產(chǎn)業(yè)、第二產(chǎn)業(yè)、第三產(chǎn)業(yè)的增加值，α₁、α₂、α₃分別為第一產(chǎn)業(yè)、第二產(chǎn)業(yè)和第三產(chǎn)業(yè)地下水消耗占地下水總消耗的權(quán)重；

(3)根據(jù)上述步驟(1)獲取的數(shù)據(jù)，畫出年最低地下水位A_低和時(shí)間的線性關(guān)系圖，橫坐標(biāo)為年份，縱坐標(biāo)為各年最低地下水位值，從最低地下水位A_低和時(shí)間的線性關(guān)系圖中得到地下水位變化的拐點(diǎn)年份；由于近年來地下水位開采日趨嚴(yán)重，添加趨勢(shì)線判斷拐點(diǎn)所在年份,拐點(diǎn)所在年份前后地下水位變化趨勢(shì)有明顯的差異。

(4)根據(jù)步驟(1)獲取的數(shù)據(jù)，分別畫出平均地下水位資料A_年、年降雨資料B_年和地下水消耗指數(shù)W與年份之間的線性關(guān)系圖，橫坐標(biāo)為年份，縱坐標(biāo)分別為各年A_年、B_年和W，根據(jù)上述步驟(3)確定的地下水位拐點(diǎn)年份，將上述三個(gè)不同線性關(guān)系圖均劃分為地下水位拐點(diǎn)前和后兩部分，分別添加各部分的趨勢(shì)線，假定在拐點(diǎn)前后，地下水位A_年與年份之間線性關(guān)系斜率分別為S_Ab和S_Aa，降雨B_年與年份之間的線性關(guān)系斜率分別為S_Bb和S_Ba，地下水消耗指數(shù)W與年份之間的線性關(guān)系斜率分別為S_Wb和S_Wa，則降水量對(duì)地下水位變化的貢獻(xiàn)度C_B為：

C_B＝((S_Bb+S_Ba)/2)/((S_Ab+S_Aa)/2)

地下水消耗指數(shù)W對(duì)地下水位變化的貢獻(xiàn)度C_W為：

C_W＝-((S_Wb+S_Wa)/2)/((S_Ab+S_Aa)/2)；

(5)根據(jù)上述步驟(1)-步驟(4)，得到按照年份先后排列的三類數(shù)據(jù)：第一類為A_年系列，第二類為不同年份降雨量B_年與C_b的乘積B‘_年系列，第三類為不同年份地下水消耗指數(shù)W與C_g的乘積W'系列；

(6)根據(jù)上述步驟(5)得到的三類數(shù)據(jù)系列A_年、B‘_年、W'，分別計(jì)算每個(gè)系列的邊緣概率分布，選取在水文頻率分析中的多種分布線型分別進(jìn)行曲線擬合，然后對(duì)每一種分布線型進(jìn)行檢驗(yàn)，得到每一類數(shù)據(jù)的最優(yōu)擬合分布線型，曲線擬合即用連續(xù)曲線近似地刻畫或比擬平面上離散點(diǎn)組所表示的坐標(biāo)之間的函數(shù)關(guān)系的一種數(shù)據(jù)處理方法，是用解析表達(dá)式逼近離散數(shù)據(jù)的一種方法。主要步驟是首先畫出系列的散點(diǎn)圖，然后選取適當(dāng)?shù)姆植季€型，最后調(diào)整分布線型相關(guān)參數(shù)。比較分布線型的理論值與實(shí)際值之間的均方根R_RMSE，以均方根R_RMSE最小為原則選擇擬合度較好的分布函數(shù)：

例如：根據(jù)上述步驟(5)得到的三類數(shù)據(jù)系列A_年、B‘_年、W'，分別計(jì)算每個(gè)系列的邊緣概率分布，選取在水文頻率分析中的多種分布線型分別進(jìn)行曲線擬合，曲線擬合即用連續(xù)曲線近似地刻畫或比擬平面上離散點(diǎn)組所表示的坐標(biāo)之間的函數(shù)關(guān)系的一種數(shù)據(jù)處理方法，是用解析表達(dá)式逼近離散數(shù)據(jù)的一種方法。主要步驟是首先畫出系列的散點(diǎn)圖，然后選取適當(dāng)?shù)姆植季€型，最后調(diào)整分布線型相關(guān)參數(shù)；通過檢驗(yàn)不同的分布線型理論值與實(shí)際值之間的均方根，即R_RMSE,以R_RMSE最小為原則選擇擬合度較好的分布函數(shù)：

上式中,F(x)為理論分布函數(shù)，F(xiàn)₀為經(jīng)驗(yàn)分布函數(shù)，以A_年系列為例，該系列中第一列是時(shí)間(x)，第二列是地下水位值A(chǔ)_年，的含義就是通過算選取的分布函數(shù)F(x)，從x₁開始計(jì)算每一個(gè)理論值F(X_i)和實(shí)際值的均方根，一直計(jì)算到最后一個(gè)x_n；

(7)采用阿基米德型科普勒家族中的Clayton Copula和Gumbel Copula函數(shù)分別構(gòu)建三維聯(lián)合概率分布：

Frank Copula形式為：

Clayton Copula形式為：

式中，θ為copula函數(shù)的參數(shù)，u₁、u₂、u₃為邊際分布函數(shù)；

(8)選用離差平方和最小準(zhǔn)則法(OLS)對(duì)Copula函數(shù)的擬合優(yōu)度進(jìn)行評(píng)價(jià)，選取離差平方和最小準(zhǔn)則法最小的copula作為優(yōu)選科普勒函數(shù)，離差平方和最小準(zhǔn)則法的計(jì)算公式為：

上式中，p_i分別為經(jīng)驗(yàn)頻率和理論頻率；i為樣本序號(hào)；

(9)根據(jù)選定的科普勒函數(shù)計(jì)算風(fēng)險(xiǎn)概率，公式為：

上式中，假定x₁為地下水位系列A_年，x₂為表征地下水消耗的W'系列，x₃為降雨量系列B‘_年系列，得到降雨量等于或者小于設(shè)定概率條件下，地下水位和地下水消耗W'相等的發(fā)生概率，根據(jù)氣象部門對(duì)降雨量的預(yù)測(cè)，確定一個(gè)降雨發(fā)生的概率，得到經(jīng)濟(jì)增長(zhǎng)一定的值條件下相應(yīng)地下水位可能的值，為制定經(jīng)濟(jì)發(fā)展政策；或假設(shè)x₃為表征地下水消耗的W'系列，x₁為地下水位系列A_年，x₂為降雨量系列B‘_年系列，得到當(dāng)產(chǎn)業(yè)發(fā)展達(dá)到某一個(gè)程度時(shí)，相應(yīng)地下水位和降雨量之間同時(shí)達(dá)到設(shè)定概率，根據(jù)未來政府的規(guī)劃，確定一個(gè)地下水消耗的W'的值，得到當(dāng)降雨達(dá)到設(shè)定概率值時(shí)，地下水位可能的值，為地下水保護(hù)政策提供依據(jù)。

上述方法的步驟5中，水文頻率分析中的多種分布線型分別為：P－Ⅲ型分布曲線(P3)、對(duì)數(shù)正態(tài)分布曲線(LN2)和廣義極值分布曲線(GEV)，P－Ⅲ型分布概率密度公式為：

其中，α為形狀參數(shù)，β為尺度參數(shù)，a₀為位置參數(shù)。

廣義極值分布曲線概率密度公式為：

其中，α為尺度參數(shù)，k為形狀參數(shù)，u為位置參數(shù)。

對(duì)數(shù)正態(tài)分布曲線概率密度公式為：

其中μ_y、δ_y分別為x序列取自然對(duì)數(shù)后形成的序列的均值和標(biāo)準(zhǔn)差。

應(yīng)用本發(fā)明的方法，可以根據(jù)氣候模型對(duì)未來降水的預(yù)測(cè)水平，假定降雨量達(dá)到某一頻率時(shí)，研究地下水位和地下水消耗指數(shù)之間風(fēng)險(xiǎn)關(guān)系，即經(jīng)濟(jì)發(fā)展達(dá)到不同的水平時(shí)地下水位可能達(dá)到某一值的概率，為進(jìn)行區(qū)域經(jīng)濟(jì)發(fā)展提供參考；或根據(jù)經(jīng)濟(jì)發(fā)展預(yù)測(cè)，假定地下水消耗指數(shù)達(dá)到某一值的前提下，研究不同的降水條件下地下水位出現(xiàn)的可能，為進(jìn)行地下水保護(hù)提供依據(jù)?？梢约俣˙‘_年等于或小于某一頻率時(shí)，不同區(qū)域經(jīng)濟(jì)發(fā)展W'和地下水位A_年同時(shí)出現(xiàn)的風(fēng)險(xiǎn)概率。

本發(fā)明方法通過確定不同產(chǎn)業(yè)地下水消耗，得到表征消耗權(quán)重的地下水消耗指數(shù)；通過對(duì)降雨和地下水消耗指數(shù)對(duì)地下水位的貢獻(xiàn)度分析，得到補(bǔ)給和消耗對(duì)地下水位的影響程度；通過構(gòu)建三維copula函數(shù)，從概率上分析地下水位變化與消耗補(bǔ)給之間的動(dòng)態(tài)關(guān)系，可以耦合氣候變化模式的降水預(yù)測(cè)或區(qū)域社會(huì)經(jīng)濟(jì)的發(fā)展規(guī)劃，為決策者提供依據(jù)。

本發(fā)明可以應(yīng)用在地下水位分析預(yù)報(bào)上，地下水位的補(bǔ)給來源主要是降雨補(bǔ)給，消耗主要是經(jīng)濟(jì)活動(dòng)帶來的地下水開采，本發(fā)明可以提供給地下水位分析有效的參考信息。具體來說，根據(jù)對(duì)未來降水的預(yù)測(cè)水平，假定降雨量為某一值的時(shí)候，研究地下水位和地下水消耗指數(shù)之間風(fēng)險(xiǎn)關(guān)系，即經(jīng)濟(jì)發(fā)展達(dá)到不同的水平時(shí)地下水位可能達(dá)到某一值的概率，為進(jìn)行區(qū)域經(jīng)濟(jì)發(fā)展提供參考；或根據(jù)經(jīng)濟(jì)發(fā)展預(yù)測(cè)，假定地下水消耗指數(shù)達(dá)到某一值的前提下，研究不同的降水條件下地下水位出現(xiàn)的可能，為進(jìn)行地下水保護(hù)提供依據(jù)。

上面結(jié)合實(shí)施例對(duì)本發(fā)明進(jìn)行了示例性的描述，顯然本發(fā)明的實(shí)現(xiàn)并不受上述方式的限制，只要采用了本發(fā)明的方法構(gòu)思和技術(shù)方案進(jìn)行的各種改進(jìn)，或未經(jīng)改進(jìn)將本發(fā)明的構(gòu)思和技術(shù)方案直接應(yīng)用于其它場(chǎng)合的，均在本發(fā)明的保護(hù)范圍內(nèi)。