本發(fā)明涉及一種城市面源污染控制領(lǐng)域，尤其是涉及一種基于swmm與matlab數(shù)據(jù)交互的參數(shù)靈敏度優(yōu)化方法。
背景技術(shù)：
為了提高swmm模型的參數(shù)率定效率，國(guó)內(nèi)外眾多科研和工程人員提出各種方法進(jìn)行swmm模型參數(shù)靈敏度分析，包括：以徑流深度、洪峰流量和流量過程線作為目標(biāo)函數(shù)的對(duì)swmm模型中的產(chǎn)匯流以及管網(wǎng)運(yùn)移模塊進(jìn)行參數(shù)靈敏度分析的局部靈敏度分析方法；以徑流體積，峰值流量以及峰值所需時(shí)間作為目標(biāo)函數(shù)，使用局部靈敏度分析方法對(duì)不透水面積比例、區(qū)域?qū)挾?、地面洼蓄深?透水和不透水區(qū)域)、管道曼寧系數(shù)以及透水和不透水區(qū)的曼寧糙率進(jìn)行的靈敏度分析；基于monte-carlo采樣方法和區(qū)域靈敏度分析(rsa)的方法確定swmm模型參數(shù)的靈敏性；采用morris方法分別對(duì)swmm模型的產(chǎn)流模塊和水文水力學(xué)參數(shù)靈敏度進(jìn)行了分析，分別考慮影響徑流深度和峰值流量的參數(shù)的靈敏度；以徑流深度、洪峰流量和洪峰時(shí)間作為輸出參數(shù)的基于montecarlo采樣方法的逐步回歸分析法對(duì)模型的水文水動(dòng)力學(xué)參數(shù)進(jìn)行靈敏度分析。
然而上述這些方法未考慮參數(shù)之間的耦合作用對(duì)模擬結(jié)果的影響加上模型參數(shù)的選定以及輸入?yún)?shù)樣本數(shù)等的限制，降低了靈敏度分析結(jié)果的代表性。因此我們選擇更合適的方法——efast對(duì)swmm模型的輸入?yún)?shù)的模型靈敏度進(jìn)行更加全面系統(tǒng)科學(xué)地分析。該方法是一種基于方差分析的參數(shù)定量靈敏度分析方法，是由sailtelli等綜合了sobol法和傅里葉幅度靈敏度檢驗(yàn)法的優(yōu)點(diǎn)提出的一種全局靈敏度分析方法，具有穩(wěn)健、計(jì)算高效且所需樣本數(shù)較低等優(yōu)點(diǎn)。然而，雖然efast方法相比于其他方法所需樣本較少，根據(jù)研究發(fā)現(xiàn)其進(jìn)行靈敏度分析時(shí)所需的參數(shù)樣本數(shù)量至少時(shí)65*k個(gè)，k為輸入?yún)?shù)個(gè)數(shù)，但是考慮到swmm作為多參數(shù)的復(fù)雜模型，其計(jì)算仍然不是能夠通過人工輕易處理的。需要尋找一種替代方法來提高靈敏度分析的效率、降低分析誤差。
技術(shù)實(shí)現(xiàn)要素：
本發(fā)明的目的就是為了克服上述現(xiàn)有技術(shù)存在的缺陷而提供一種基于swmm與matlab數(shù)據(jù)交互的參數(shù)靈敏度優(yōu)化方法。
本發(fā)明的目的可以通過以下技術(shù)方案來實(shí)現(xiàn)：
一種基于swmm與matlab數(shù)據(jù)交互的參數(shù)靈敏度優(yōu)化方法，包括步驟：
s1：通過matlab模擬efast分析方法，將輸入數(shù)據(jù)樣本傳遞給swmm所需的inp格式文件；
s2：采用swmm的運(yùn)算引擎調(diào)用并提取inp格式文件中的信息進(jìn)行運(yùn)算，并將運(yùn)算結(jié)果分別存儲(chǔ)在rpt格式文件和out格式文件中；
s3：通過matlab調(diào)用并提取rpt格式文件中相應(yīng)的輸出信息并生成相應(yīng)的輸出數(shù)據(jù)樣本；
s4：分別將輸入數(shù)據(jù)樣本和輸出數(shù)據(jù)樣本通過efast靈敏度分析算法分別計(jì)算相應(yīng)參數(shù)的一階靈敏度因子；
s5：根據(jù)獲得的一階靈敏度因子優(yōu)化針對(duì)性優(yōu)化一階靈敏度因子大的參數(shù)對(duì)應(yīng)的條件。
所述步驟s1具體包括步驟：
s11：運(yùn)行matlab，根據(jù)efast分析方法所需的參數(shù)的輸入數(shù)據(jù)樣本；
s12：將生成的輸入數(shù)據(jù)樣本傳遞給swmm所需的inp格式文件。
所述步驟s4具體包括：
s41：構(gòu)建多維函數(shù)模型y＝f(x)，其中x(x1,x2,…,xn)為輸入?yún)?shù)；
s42：設(shè)定搜尋函數(shù)gi，并采用搜尋函數(shù)將多維函數(shù)模型轉(zhuǎn)換為一維函數(shù)模型y＝f(s)：
gi＝(sinωis),i＝1,2,…,n
其中：i為參數(shù)xi的參數(shù)序號(hào)，ωi為參數(shù)xi的整數(shù)頻率，s為實(shí)數(shù)；
s43：對(duì)一維函數(shù)模型進(jìn)行傅里葉變換：
其中：p為傅里葉因子，p∈z，ap,bp為傅里葉振幅，且：
s44：計(jì)算由各參數(shù)引起的方差以及模型總方差：
其中：vi為參數(shù)xi引起的方差，v為模型總方差，λj為傅里葉級(jí)數(shù)的頻譜，j∈z；
s45：根據(jù)獲得的由各參數(shù)引起的方差以及模型總方差，得到各參數(shù)的一階靈敏度因子：
其中：si為參數(shù)xi的一階靈敏度因子。
傅里葉級(jí)數(shù)具體為：
其中：λp為傅里葉級(jí)數(shù)的頻譜。
所述輸入?yún)?shù)至少包括區(qū)域?qū)挾?、區(qū)域坡度、區(qū)域不透水率、透水性區(qū)域的貯水深度、不透水區(qū)域的貯水深度、不透水曲曼寧系數(shù)、透水區(qū)曼寧系數(shù)和不透水區(qū)無洼蓄地面技比例。
與現(xiàn)有技術(shù)相比，本發(fā)明具有以下優(yōu)點(diǎn)：
1)efast方法是一種基于方差分析的參數(shù)定量靈敏度分析方法，是由sailtelli等綜合了sobol法和傅里葉幅度靈敏度檢驗(yàn)法的優(yōu)點(diǎn)提出的一種全局靈敏度分析方法，該方法穩(wěn)健、計(jì)算高效且所需樣本數(shù)較低。
2)efast靈敏度分析方法屬于全局靈敏度分析方法，不同于當(dāng)前常用的局部靈敏度分析方法和morris方法，該方法考慮了swmm模型的復(fù)雜性以及不同參數(shù)的耦合作用對(duì)輸出結(jié)果的影響，即考慮參數(shù)之間的相互作用對(duì)于模型輸出結(jié)果的靈敏性。
3)利用matlab強(qiáng)大的數(shù)據(jù)處理功能，實(shí)現(xiàn)efast算法涉及的大量計(jì)算，快速、準(zhǔn)確地明確swmm模型參數(shù)靈敏度情況，尋找出對(duì)模型模擬影響最為顯著的一些參數(shù)，以便模型率定時(shí)提高效率，增加模型模擬準(zhǔn)確度。
附圖說明
圖1為本發(fā)明的主要方法流程示意圖。
具體實(shí)施方式
下面結(jié)合附圖和具體實(shí)施例對(duì)本發(fā)明進(jìn)行詳細(xì)說明。本實(shí)施例以本發(fā)明技術(shù)方案為前提進(jìn)行實(shí)施，給出了詳細(xì)的實(shí)施方式和具體的操作過程，但本發(fā)明的保護(hù)范圍不限于下述的實(shí)施例。
本發(fā)明目的在于提供一種可以提高用efast靈敏度分析方法對(duì)swmm模型參數(shù)進(jìn)行靈敏性排序的工作效率，降低處理大量數(shù)據(jù)時(shí)的誤差的方法。方法具有高效，準(zhǔn)確，省時(shí)省力的優(yōu)點(diǎn)。
通過數(shù)學(xué)分析軟件matlab編程軟件，編程實(shí)現(xiàn)matlab解決efast算法計(jì)算問題以及matlab與swmm軟件之間的數(shù)據(jù)交互，完成swmm模型參數(shù)的快速靈敏度分析，最終依賴得到的結(jié)果進(jìn)行模型的快速優(yōu)化，有針對(duì)性的控制特定變量達(dá)到優(yōu)化目標(biāo)。
一種基于swmm與matlab數(shù)據(jù)交互的參數(shù)靈敏度優(yōu)化方法，包括步驟：
s1：通過matlab模擬efast分析方法，將輸入數(shù)據(jù)樣本傳遞給swmm所需的inp格式文件，具體包括步驟：
s11：運(yùn)行matlab，根據(jù)efast分析方法所需的參數(shù)的輸入數(shù)據(jù)樣本；
s12：將生成的輸入數(shù)據(jù)樣本傳遞給swmm所需的inp格式文件。
s2：采用swmm的運(yùn)算引擎調(diào)用并提取inp格式文件中的信息進(jìn)行運(yùn)算，并將運(yùn)算結(jié)果分別存儲(chǔ)在rpt格式文件和out格式文件中；
s3：通過matlab調(diào)用并提取rpt格式文件中相應(yīng)的輸出信息并生成相應(yīng)的輸出數(shù)據(jù)樣本；
s4：分別將輸入數(shù)據(jù)樣本和輸出數(shù)據(jù)樣本通過efast靈敏度分析算法分別計(jì)算相應(yīng)參數(shù)的一階靈敏度因子，具體包括步驟：
s41：構(gòu)建多維函數(shù)模型y＝f(x)，其中x(x1,x2,…,xn)為輸入?yún)?shù)，輸入?yún)?shù)至少包括區(qū)域?qū)挾取^(qū)域坡度、區(qū)域不透水率、透水性區(qū)域的貯水深度、不透水區(qū)域的貯水深度、不透水曲曼寧系數(shù)、透水區(qū)曼寧系數(shù)和不透水區(qū)無洼蓄地面技比例，各參數(shù)具有特定變化范圍及分布形式，構(gòu)成一個(gè)多為參數(shù)空間rn；
s42：設(shè)定搜尋函數(shù)gi，并采用搜尋函數(shù)將多維函數(shù)模型轉(zhuǎn)換為一維函數(shù)模型y＝f(s)：
gi＝(sinωis),i＝1,2,…,n
其中：i為參數(shù)xi的參數(shù)序號(hào)，ωi為參數(shù)xi的整數(shù)頻率，s為實(shí)數(shù)；
s43：對(duì)一維函數(shù)模型進(jìn)行傅里葉變換：
其中：p為傅里葉因子，p∈z，ap,bp為傅里葉振幅，且：
s44：計(jì)算由各參數(shù)引起的方差以及模型總方差：
其中：vi為參數(shù)xi引起的方差，v為模型總方差，λj為傅里葉級(jí)數(shù)的頻譜，j∈z；
傅里葉級(jí)數(shù)的頻譜曲線具體為：
其中：λp為傅里葉級(jí)數(shù)。
根據(jù)傅里葉變換的特性，ap＝a-p，bp＝b-p，λp＝λ-p，則由參數(shù)xi的輸入變化所引起的模型結(jié)果的方差vi可以表示為：
其中：z0為z的去心鄰域。則模型的總方差v計(jì)算方法如下：
對(duì)s在區(qū)間[-π,+π]之間等間距取樣，把取樣值通過搜尋函數(shù)轉(zhuǎn)換成為每一個(gè)參數(shù)的取值，輸入模型，經(jīng)過多次運(yùn)行模型后，可分別獲得ap,bp的近似值：
其中：
ns為取樣數(shù)。由ap和bp，以及參數(shù)xi所對(duì)應(yīng)的頻率ωi通過各參數(shù)引起的方差以及模型總方差的計(jì)算式即可獲得各個(gè)參數(shù)所引起的模型輸出的方差vi以及模型輸出的總方差v。由于模型輸出的總方差是由各個(gè)參數(shù)及參數(shù)間相互作用共同得到，可將總方差v分解為如下形式：
v＝∑ivi+∑i≠jvij+∑i≠j≠mvijm+…+v1,2,…,n
其中：vij為參數(shù)xi通過與參數(shù)xj的相互作用貢獻(xiàn)的方差，vijm為參數(shù)xi與xj及xm的耦合總用貢獻(xiàn)的方差，v11...k為xi與x1，x2，...，xk共同作用貢獻(xiàn)的方差。因此，可以分別獲得參數(shù)xi的一階敏感因子si為：
其中：vij為參數(shù)xi通過與參數(shù)xi的相互作用貢獻(xiàn)的方差，vijm為參數(shù)xi與xj及xm的耦合總用貢獻(xiàn)的方差，v1,2,…,n為xi與除xi以外所有參數(shù)共同作用貢獻(xiàn)的方差。因此，可以分別獲得參數(shù)xi的一階敏感因子。
s45：根據(jù)獲得的由各參數(shù)引起的方差以及模型總方差，得到各參數(shù)的一階靈敏度因子：
其中：si為參數(shù)xi的一階靈敏度因子。
一階靈敏度因子反映了參數(shù)對(duì)于模型輸出總方差的獨(dú)立貢獻(xiàn)率。同理，可以獲得參數(shù)的二階以及各階靈敏度因子，計(jì)算如下：
參數(shù)的總靈敏度因子為參數(shù)的各階靈敏度因子之和，計(jì)算方式如下：
sti＝si+sij+sijm+…+s1,2,…,n
參數(shù)的總靈敏度因子代表了參數(shù)同其他所有具有交互作用的參數(shù)之間的相互作用對(duì)于輸出總方差的貢獻(xiàn)之和。
s5：根據(jù)獲得的一階靈敏度因子優(yōu)化針對(duì)性優(yōu)化一階靈敏度因子大的參數(shù)對(duì)應(yīng)的條件，提高模型優(yōu)化效率，從而加快模型優(yōu)化速度。
申請(qǐng)方案的主要技術(shù)難題有兩個(gè)：一是efast算法的matlab實(shí)現(xiàn)問題，二是matlab與swmm軟件之間的數(shù)據(jù)交互的實(shí)現(xiàn)。
根據(jù)前人研究成果，在efast算法的基礎(chǔ)上，通過對(duì)其部分算法進(jìn)行修改或者重新編寫成適合matlab的算法；在matlab與swmm的混合編程中，通過matlab調(diào)用swmm的動(dòng)態(tài)數(shù)據(jù)鏈接數(shù)據(jù)庫(kù)文件實(shí)現(xiàn)swmm的調(diào)用，swmm調(diào)用時(shí)需要提供運(yùn)行所需的數(shù)據(jù)文件格式。
城市雨洪模擬軟件swmm有其獨(dú)立的輸入文件格式*.inp，其中包含了模型運(yùn)行所需的各個(gè)參數(shù)取值以及建模的空間數(shù)據(jù)，模型運(yùn)行時(shí)自動(dòng)生產(chǎn)另外兩種格式的輸出文件*.inp和*.out，分別表示swmm的分析報(bào)告文件和運(yùn)行結(jié)果文件。通過在matlab中實(shí)現(xiàn)efast靈敏度的計(jì)算，輸入swmm模型所需參數(shù)，matlab根據(jù)輸入的數(shù)據(jù)樣本生成傳遞swmm的輸入文件*.inp，模型運(yùn)算生成*.rpt和*.out的結(jié)果文件，通過matlab調(diào)用、提取*.rpt中的輸出信息并生成可供其處理的與輸入數(shù)據(jù)樣本相對(duì)應(yīng)的輸出數(shù)據(jù)樣本，繼而將輸入樣本與對(duì)應(yīng)輸出樣本通過matlab的efast靈敏度分析算法計(jì)算相應(yīng)參數(shù)的靈敏度因子。
swmm有自己獨(dú)立的輸入文件格式*.inp，文件中包含了所有swmm模型運(yùn)行所需的各個(gè)參數(shù)的取值以及建模的空間數(shù)據(jù)，在swmm調(diào)用時(shí)，同時(shí)需要提供另外兩種格式的輸入文件*.rpt和*.out，兩者分別是swmm的分析報(bào)告文件和swmm的運(yùn)行結(jié)果文件，這兩個(gè)文件的內(nèi)容可以在swmm的運(yùn)行時(shí)自動(dòng)生成，而*.inp文件則需要在swmm運(yùn)算前嚴(yán)格按照swmm所需參數(shù)及高程信息生成。為此，首先通過運(yùn)行matlab，根據(jù)efast分析方法所需的特定數(shù)據(jù)樣本生成方x(根據(jù)特定搜尋曲線，searchingcurve)，再將生成的數(shù)據(jù)樣本x分別傳遞給swmm的輸入文件*.inp，每次傳遞都分別用一組新的數(shù)據(jù)樣本代替舊數(shù)據(jù)樣本。之后，swmm的運(yùn)算引擎swmmengine將調(diào)用提取輸入文件中信息進(jìn)行運(yùn)算，并將運(yùn)算結(jié)果分別存儲(chǔ)在*.rpt和*.out文件中，然后通過matlab調(diào)用并提取*.rpt文件中相應(yīng)的輸出信息并生成相應(yīng)的輸出數(shù)據(jù)樣本y，最后，分別將輸入樣本x和輸出樣本y通過efast靈敏度分析特定的算法分別計(jì)算相應(yīng)參數(shù)的一階和總靈敏度因子。上述過程將會(huì)根據(jù)所預(yù)設(shè)的參數(shù)個(gè)數(shù)k和搜尋曲線的取樣個(gè)數(shù)ns確定，一次完整的靈敏度分析將需要進(jìn)行k*ns次swmm與matlab間的數(shù)據(jù)交互，而matlab同swmm的*.inp文件之間的數(shù)據(jù)交互是通過正則表達(dá)式完成的
以下各處一個(gè)實(shí)驗(yàn)案例說明：
基于chj排水系統(tǒng)，考慮不同降雨條件下(不同降雨強(qiáng)度、不同降雨歷時(shí))，swmm模型的不同輸入?yún)?shù)對(duì)系統(tǒng)輸出(徑流深度，洪峰流量，達(dá)到洪峰的時(shí)間以及系統(tǒng)總排水量)。案例分別考慮了降雨入滲模型horton、green-ampt以及cn曲線法下，不同降雨條件下，模型參數(shù)靈敏度的分析結(jié)horton入滲模型下，不透水面比例是影響徑流深度、峰值流量、累計(jì)排水量的首要影響因子；green-ampt入滲模型下，不透水面百分比是影響徑流深度的首要影響因子、不透水面曼寧系數(shù)是影響峰值流量的。為我們明確參數(shù)對(duì)模型結(jié)果影響較大的參數(shù)類型，減少模型率定工作所需時(shí)間，提高模擬結(jié)果的準(zhǔn)確度。