本發(fā)明涉及涉及一種針對大規(guī)模問題的空間曲面建模的方法，可用于大數(shù)據(jù)背景下數(shù)字高程模型的生成、氣象、生物量、土壤的空間分布模擬等領(lǐng)域，也可以視為曲面柵格逼近的一種方法，可用于大規(guī)模的物理、化學(xué)、醫(yī)學(xué)等方面的曲面造型。
背景技術(shù)：
：為了解決了半個(gè)世紀(jì)以來困擾曲面建模的誤差問題和多尺度問題，我們以微分幾何原理和優(yōu)化控制論為理論基礎(chǔ)，建立了一個(gè)以全局性近似數(shù)據(jù)(包括遙感數(shù)據(jù)和全球模型粗分辨率模擬數(shù)據(jù))為驅(qū)動(dòng)場、以局地高精度數(shù)據(jù)(包括監(jiān)測網(wǎng)數(shù)據(jù)和調(diào)查采樣數(shù)據(jù))為優(yōu)化控制條件的高精度曲面建模(HighAccuracySurfaceModelling，簡寫為HASM)方法，并在20多年大量應(yīng)用研究基礎(chǔ)上，提煉形成了地球表層建?；径?。HASM最后可轉(zhuǎn)化為一個(gè)由地面采樣約束的等式約束最小二乘問題(Equalityconstraintleastsquaresproblem:LSE)，目的是為了在保證采樣點(diǎn)處模擬值等于采樣值的條件下，保持整體模擬誤差最小。充分利用采樣信息，也是保證迭代趨近于最佳模擬效果的有效手段。盡管HASM方法比傳統(tǒng)的插值方法在模擬精度上有了很大的改善，但研究表明，當(dāng)前HASM只能處理小規(guī)模的問題，上述求解過程是將LSE問題轉(zhuǎn)化為法方程組求解。一方面，轉(zhuǎn)化后的方程組系數(shù)矩陣條件數(shù)非常大，造成了問題的敏感性及病態(tài)性，這使得求解方程組的迭代方法收斂速度非常慢，或者不收斂，并由此造成會(huì)得不到原求解問題的近似解；另一方面，當(dāng)前采用的預(yù)處理共軛梯度法求解HASM方程組需要將方程組系數(shù)矩陣全部的一次性輸入內(nèi)存，隨著求解問題規(guī)模的增加，其存儲(chǔ)空間會(huì)明顯增長，對于大規(guī)模問題當(dāng)前求解策略已不再使用。當(dāng)前HASM對于全國1km分辨率數(shù)據(jù)的模擬，往往采用分片處理策略，這增加了分區(qū)邊界處的誤差，從而降低了模擬結(jié)果的整體精度。技術(shù)實(shí)現(xiàn)要素：本發(fā)明針對現(xiàn)有的HASM模型在求解大規(guī)模問題中的缺陷，提供了一種改進(jìn)的HASM，所需存儲(chǔ)空間小，能夠在保證精度的情況下，解決了大規(guī)模問題。本發(fā)明解決上述技術(shù)問題的技術(shù)方案如下：一種基于大數(shù)據(jù)的高精度曲面建模方法，包括以下步驟：步驟1：創(chuàng)建各采樣點(diǎn)的地理坐標(biāo)信息和待測變量采樣值，所述地理坐標(biāo)信息中包括采樣點(diǎn)的經(jīng)度信息和采樣點(diǎn)的緯度信息；步驟2：將待測區(qū)域空間離散化為網(wǎng)格點(diǎn)形式，得到網(wǎng)格點(diǎn)離散值，根據(jù)地理坐標(biāo)信息和待測變量采樣值建立采樣方程，所述采樣方程用于判斷采樣點(diǎn)是否在網(wǎng)格點(diǎn)上，其中，如果采樣點(diǎn)在所述網(wǎng)格點(diǎn)上，則該網(wǎng)格點(diǎn)的值即為待測變量采樣值；如果采樣點(diǎn)在網(wǎng)格內(nèi)，則將距該采樣點(diǎn)最近的網(wǎng)格點(diǎn)上利用泰勒展開得到該網(wǎng)格點(diǎn)上的近似采樣值；步驟3：根據(jù)網(wǎng)格點(diǎn)離散值計(jì)算待測區(qū)域每個(gè)網(wǎng)格點(diǎn)的第一類基本量E、F、G和第二類基本量L、M、N，其中所述第一基本量用于表示模擬曲面上曲線的長度、模擬曲面的面積和模擬曲面上曲線的曲率，所述第二基本量用于表示模擬曲面的局部彎曲變化程度，將用第一基本量和第二基本量表示的曲面的偏微分方程組進(jìn)行高階差分離散，獲得離散方程組對應(yīng)的代數(shù)系統(tǒng)，將所述代數(shù)系統(tǒng)與所述采樣方程組合成等式約束最小二乘問題；步驟4：將等式約束最小二乘為題轉(zhuǎn)化為求解截?cái)嗄繕?biāo)函數(shù)極小值問題；步驟5：把求解截?cái)嗄繕?biāo)函數(shù)極小值問題轉(zhuǎn)化為求解對稱不定方程組，所述對稱不定方程組為高精度曲面建模HASM方程組；具體的，步驟6：隨機(jī)選取HASM方程組的迭代初值；步驟7：將HASM方程組中的系數(shù)矩陣進(jìn)行行分塊，并對系數(shù)矩陣分解后的塊矩陣進(jìn)行存儲(chǔ)；步驟8：對迭代初值采用塊行投影迭代法，求解HASM方程組，并判斷求解結(jié)果是否收斂；步驟9：當(dāng)求解結(jié)果不收斂時(shí)，對求解結(jié)果重新采用塊行投影迭代法，求解HASM方程組，并判斷求解結(jié)果是否收斂，如果收斂，執(zhí)行步驟10，否則，重新執(zhí)行步驟9；步驟10：當(dāng)HASM方程組的解收斂時(shí)，進(jìn)一步判斷HASM方程組的解是否滿足高斯科達(dá)齊方程組，若不滿足，則執(zhí)行步驟6；若滿足，則根據(jù)HASM方程組的解輸出關(guān)于待測變量的高精度模擬曲面模型。本發(fā)明旨在解決空間曲面建模領(lǐng)域中的大數(shù)據(jù)問題，具有如下優(yōu)點(diǎn)：1、完善并發(fā)展了已有的模型，改進(jìn)后的模型具有較低的存儲(chǔ)需求及較少的計(jì)算時(shí)間；2、引入截?cái)嘞到y(tǒng)，給出了求解約束最小二乘問題的優(yōu)化求解方法，結(jié)果更穩(wěn)定；3、給出了大型方程組系數(shù)矩陣特定的分塊方式，保證了矩陣各行之間的線性獨(dú)立性，避免了信息冗余，從而給出了行投影算法的加速算法即塊行投影迭代法，提高了方法的收斂速度；4、改進(jìn)后的方法可求解大規(guī)模問題，計(jì)算時(shí)間與計(jì)算網(wǎng)格數(shù)成線性關(guān)系，所需存儲(chǔ)空間遠(yuǎn)遠(yuǎn)小于HASM方程組系數(shù)矩陣中非零元素個(gè)數(shù)；在保證精度的情況下，解決了HASM在求解大規(guī)模問題中的缺陷。在上述技術(shù)方案的基礎(chǔ)上，本發(fā)明還可以做如下改進(jìn)。進(jìn)一步，所述曲面的偏微分方程組為：其中，E＝1+fx2，F(xiàn)＝fx·fy，x為空間離散化的網(wǎng)格點(diǎn)的橫坐標(biāo)，y為空間離散化的網(wǎng)格點(diǎn)的縱坐標(biāo)，fx為函數(shù)f對x的一階偏導(dǎo)數(shù)，fxx為函數(shù)f對x的二階偏導(dǎo)數(shù)，fy為函數(shù)f對y的一階偏導(dǎo)數(shù)，fyy為函數(shù)f對y的二階偏導(dǎo)數(shù)Ex、Fx、Gx分別為E、F、G對x的一階偏導(dǎo)數(shù)，Ey、Fy、Gy分別為E、F、G對y的一階偏導(dǎo)數(shù)，為第二類克里斯托弗爾符號。采用上述進(jìn)一步方案的有益效果是，使得HASM建立在完整的微分幾何學(xué)理論基礎(chǔ)之上，保證了HASM的穩(wěn)定性，提高了HASM的模擬精度。進(jìn)一步，所述步驟3具體包括：步驟3.1：根據(jù)第一類基本量E、F、G和第二類基本量L、M、N的有限差分逼近計(jì)算待測區(qū)域每個(gè)網(wǎng)格點(diǎn)的第一類基本量和第二類基本量的值；步驟3.2：結(jié)合第一類基本量的有限差分逼近、第二類基本量的有限差分逼近和第二類克里斯托弗爾符號的有限差分得到離散方程組；步驟3.3：將所述離散方程組轉(zhuǎn)換成矩陣形式的代數(shù)系統(tǒng)；步驟3.4：將所述代數(shù)系統(tǒng)與采樣方程組合成等式約束最小二乘問題。采用上述進(jìn)一步方案的有益效果是，HASM最后轉(zhuǎn)化為一個(gè)由地面采樣約束的等式約束最小二乘問題，目的是為了在保證采樣點(diǎn)處模擬值等于采樣值的條件下，保持整體模擬誤差最小。充分利用采樣信息，也是保證迭代趨近于最佳模擬效果的有效手段。進(jìn)一步，所述步驟7中將HASM方程組中的系數(shù)矩陣進(jìn)行行分塊的具體方法為：設(shè)系數(shù)矩陣A的每個(gè)子塊Ai最多有u行，且Ai∈Ru×n,u＜＜n均為滿秩，即每塊的行是獨(dú)立的，假設(shè)正在生成一個(gè)子塊Ai，該子塊已包含矩陣A中的j行(j＜u)，且已知AiT的QR分解為AiT＝QjRj，塊AiT的行相關(guān)性的估計(jì)其中rhh為位于Rj的(h,h)位置的元素，a1為矩陣A中的一行，令根據(jù)AiT的QR分解，采用QR分解更新方法來得到的QR分解，然后計(jì)算的行相關(guān)性估計(jì)，如果該估計(jì)仍小于κ，則a1加入塊Ai，其中，κ為預(yù)先設(shè)定的一個(gè)正數(shù)。進(jìn)一步，所述QR分解采用FrancisQR分解策略。采用上述進(jìn)一步方案的有益效果是，對于涉及到的QR分解，采用FrancisQR分解策略，可將其計(jì)算花費(fèi)由O(u3)降低為O(u2),u＜＜n，從而提高計(jì)算效率。對矩陣分塊可加速行投影迭代法的收斂，減少HASM方程組求解計(jì)算時(shí)間。進(jìn)一步，所述步驟7中對系數(shù)矩陣分解后的塊矩陣進(jìn)行存儲(chǔ)的方式為行壓縮存儲(chǔ)方式CSR。采用上述進(jìn)一步方案的有益效果是，CSR格式是最常見和靈活的壓縮格式，它將矩陣的非零元素進(jìn)行按行壓縮存儲(chǔ)，并用專用數(shù)組來記錄非零元素原有的位置，壓縮效率高，壓縮過程便于理解，由此顯著降低了內(nèi)存需求。進(jìn)一步，所述步驟8和步驟9具體包括：假設(shè)當(dāng)前迭代點(diǎn)為，判斷xk是否收斂，若不收斂，則找出離xk直交距離最遠(yuǎn)的塊，選取xk在所述塊上的投影作為下一個(gè)迭代點(diǎn)。采用上述進(jìn)一步方案的有益效果是，加快收斂速度。附圖說明圖1為本發(fā)明實(shí)施例提供的一種基于大數(shù)據(jù)的高精度曲面建模方法步驟流程圖；圖2為標(biāo)準(zhǔn)數(shù)學(xué)曲面及采樣點(diǎn)分布圖；圖3為HASM、FRK、IDW和FRS模擬曲面與真實(shí)曲面的誤差圖。具體實(shí)施方式以下結(jié)合附圖對本發(fā)明的原理和特征進(jìn)行描述，所舉實(shí)例只用于解釋本發(fā)明，并非用于限定本發(fā)明的范圍。如圖1所示，本發(fā)明實(shí)施例提供的一種基于大數(shù)據(jù)的高精度曲面建模方法，包括以下步驟：步驟1：創(chuàng)建各采樣點(diǎn)的地理坐標(biāo)信息和待測變量采樣值，所述地理坐標(biāo)信息中包括采樣點(diǎn)的經(jīng)度信息和采樣點(diǎn)的緯度信息；步驟2：將待測區(qū)域空間離散化為網(wǎng)格點(diǎn)形式，得到網(wǎng)格點(diǎn)離散值，根據(jù)地理坐標(biāo)信息和待測變量采樣值建立采樣方程，所述采樣方程用于判斷采樣點(diǎn)是否在網(wǎng)格點(diǎn)上，其中，如果采樣點(diǎn)在所述網(wǎng)格點(diǎn)上，則該網(wǎng)格點(diǎn)的值即為待測變量采樣值；如果采樣點(diǎn)在網(wǎng)格內(nèi)，則將距該采樣點(diǎn)最近的網(wǎng)格點(diǎn)上利用泰勒展開得到該網(wǎng)格點(diǎn)上的近似采樣值；步驟3：根據(jù)網(wǎng)格點(diǎn)離散值計(jì)算待測區(qū)域每個(gè)網(wǎng)格點(diǎn)的第一類基本量E、F、G和第二類基本量L、M、N，其中所述第一基本量用于表示模擬曲面上曲線的長度、模擬曲面的面積和模擬曲面上曲線的曲率，所述第二基本量用于表示模擬曲面的局部彎曲變化程度，將用第一基本量和第二基本量表示的曲面的偏微分方程組進(jìn)行高階差分離散，獲得離散方程組對應(yīng)的代數(shù)系統(tǒng)，將所述代數(shù)系統(tǒng)與所述采樣方程組合成等式約束最小二乘問題；所述步驟3具體包括：步驟3.1：根據(jù)第一類基本量E、F、G和第二類基本量L、M、N的有限差分逼近計(jì)算待測區(qū)域每個(gè)網(wǎng)格點(diǎn)的第一類基本量和第二類基本量的值；步驟3.2：結(jié)合第一類基本量的有限差分逼近、第二類基本量的有限差分逼近和第二類克里斯托弗爾符號的有限差分得到離散方程組；步驟3.3：將所述離散方程組轉(zhuǎn)換成矩陣形式的代數(shù)系統(tǒng)；步驟3.4：將所述代數(shù)系統(tǒng)與采樣方程組合成等式約束最小二乘問題。具體的，根據(jù)曲面論基本定理：設(shè)曲面的第一類和第二類基本量E、F、G、L、M和N滿足對稱性，E、F、G正定，E、F、G、L、M和N滿足曲面的偏微分方程組，即Gauss方程組，則全微分方程組在f(x,y)＝f(x0,y0)(x＝x0,y＝y(tǒng)0)初始條件下，存在著唯一的解z＝f(x,y)。Gauss方程組的表達(dá)式為，x為空間離散化的網(wǎng)格點(diǎn)的橫坐標(biāo)，y為空間離散化的網(wǎng)格點(diǎn)的縱坐標(biāo)，fx為函數(shù)f對x的一階偏導(dǎo)數(shù)，fxx為函數(shù)f對x的二階偏導(dǎo)數(shù)，fy為函數(shù)f對y的一階偏導(dǎo)數(shù)，fyy為函數(shù)f對y的二階偏導(dǎo)數(shù)Ex、Fx、Gx分別為E、F、G對x的一階偏導(dǎo)數(shù)，Ey、Fy、Gy分別為E、F、G對y的一階偏導(dǎo)數(shù)，為第二類克里斯托弗爾符號。若{(xi,yj)}是計(jì)算域Ω的正交剖分、[0,Lx]×[0,Ly]為無量綱準(zhǔn)化計(jì)算域、max{Lx,Ly}＝1、Lx為研究區(qū)域水平長度的歸一化后的取整數(shù)值。Ly為研究區(qū)域垂直方向長度的歸一化后的取整數(shù)值。為計(jì)算步長、{(xi,yj)|0≤i≤I+1,0≤j≤J+1}為標(biāo)準(zhǔn)化計(jì)算域的柵格，則第一類基本量的有限差分逼近為，其中，Ei,j，F(xiàn)i,j，Gi,j、分別為E、F、G在網(wǎng)格點(diǎn)(i,j)上的值，fi,j為模擬曲面在網(wǎng)格點(diǎn)(i,j)上的值；第二類基本量的有限差分逼近為，其中，Li,j，Mi,j，Ni,j分別為L、M、N在網(wǎng)格點(diǎn)(i,j)上的值；第二類克里斯托弗爾符號的有限差分表達(dá)為，其中，為第二類克里斯托弗爾符號在網(wǎng)格點(diǎn)(i,j)上的值；高斯方程組的有限差分形式為,(1.1.2)的矩陣形式可寫為：其中，其中，為上述(1.1.2)式右端項(xiàng)采用有限差分格式離散后的值。結(jié)合采樣信息的有效約束控制，上述約束最小二乘問題(1.1.3)可表達(dá)為HASM所求解的等式約束的最小二乘問題，Sx＝k為采樣點(diǎn)滿足的采樣方程，其中S和k分別為采樣矩陣和采樣向量；如果是z＝f(x,y)在第p采樣點(diǎn)(xi,yj)的值，則sp,(i-1)×J+j＝1，HASM最后轉(zhuǎn)化為一個(gè)由地面采樣約束的等式約束最小二乘問題(Equalityconstraintleastsquaresproblem:LSE)步驟4：將等式約束最小二乘為題轉(zhuǎn)化為求解截?cái)嗄繕?biāo)函數(shù)極小值問題；具體的，通過研究發(fā)現(xiàn)，引入一參數(shù)λ，LSE問題可首先通過Lagrange乘子法將其轉(zhuǎn)化為如下求截?cái)嗄繕?biāo)函數(shù)極小值問題：其中步驟5：把求解截?cái)嗄繕?biāo)函數(shù)極小值問題轉(zhuǎn)化為求解對稱不定方程組，所述對稱不定方程組為高精度曲面建模HASM方程組；具體的，對求導(dǎo)，另其導(dǎo)數(shù)為0，可得結(jié)合及Sx＝k可得如下對稱不定系統(tǒng)：其中，I為對角線元素全為1，其他元素為0的單位矩陣。為矩陣的轉(zhuǎn)置矩陣。ST為S的轉(zhuǎn)置矩陣。λ為采樣點(diǎn)權(quán)重構(gòu)成的權(quán)重矩陣，在HASM中取值為1～2之間。步驟6：隨機(jī)選取HASM方程組的迭代初值；步驟7：將HASM方程組中的系數(shù)矩陣進(jìn)行行分塊，并對系數(shù)矩陣分解后的塊矩陣進(jìn)行存儲(chǔ)；具體的，對于系數(shù)矩陣的分塊方式，本研究給出了A的一分塊算法，應(yīng)用該算法每塊最多有u行，且每塊Ai∈Ru×n,u＜＜n均為滿秩，即每塊的行是獨(dú)立的。假設(shè)正在生成一個(gè)塊Ai，該塊已包含矩陣A中的j行(j＜u)，且已知AiT的QR分解為AiT＝QjRj，令rhh為位于Rj的(h,h)位置的元素。塊AiT的行相關(guān)性的估計(jì)其中對于矩陣A中的一行a1是否加入該塊Ai的準(zhǔn)則為a1沒有被分配給任何其他的塊，且a1加入后Ai的行相關(guān)性估計(jì)βi不超過κ，κ為預(yù)先設(shè)定的一個(gè)正數(shù)。為了決定一個(gè)新的行a1是否加入塊Ai，令且根據(jù)已有的AiT的QR分解來得到的QR分解，對此采用QR分解更新方法來實(shí)現(xiàn)。然后計(jì)算的行相關(guān)性估計(jì)，如果該估計(jì)仍小于κ，則a1加入塊Ai。該塊分解算法描述如下，該算法基于C++實(shí)現(xiàn)。令Γ＝1,2,…,m,Γs＝{l:A中的a1列已分配給某塊}初始化：令i＝1,Γs＝φ，whileΓs≠Γ,do令Γc＝Γ\Γs，j＝1,Γi＝φ選取l∈Γc，令A(yù)i＝[a1],R1＝[||a1||2,0]T,計(jì)算Q1＝In-2vvT令dmax＝||a1||2，dmin＝||a1||2Γc←Γc\{l}，Γs←Γs∪{l}，Γi←Γi∪{l}while(j＜u且Γc≠φ)do選取l∈Γc令v＝[||a||2,0]T∈Rn-j計(jì)算tmax←max{dmax,||a||2},tmin←min{dmin,||a||2},iftmin≠0,且則更新Ai←[Ai,a1]Γc←Γc\{l}，Γs←Γs∪{l}，Γi←Γi∪{l}j＝j(luò)+1dmax←tmax,dmin←tminelseΓc←Γc\{l}endendi←i+1end為了提高計(jì)算效率，對于涉及到的QR分解，我們采用了FrancisQR分解策略，并將其計(jì)算花費(fèi)由O(u3)降低為O(u2),u＜＜n.上述算法中，涉及到稀疏矩陣塊Ai間的乘法及稀疏矩陣與向量乘的運(yùn)算，向量內(nèi)積及向量加法運(yùn)算。稀疏矩陣通常存儲(chǔ)為壓縮格式。隨著應(yīng)用場景及計(jì)算平臺(tái)的不同，矩陣通常被壓縮為不同的存儲(chǔ)格式。壓縮格式的選擇要綜合考慮矩陣的稀疏特點(diǎn)及計(jì)算平臺(tái)。本研究對矩陣存儲(chǔ)采用行壓縮存儲(chǔ)方式(CSR)，CSR格式是最常見和靈活的壓縮格式，它將矩陣的非零元素進(jìn)行按行壓縮存儲(chǔ)，并用專用數(shù)組來記錄非零元素原有的位置，壓縮效率高，壓縮過程便于理解。CSR格式并于被移植到不同平臺(tái)上。許多新的壓縮格式也都基于CSR格式修改而來。這種方法存儲(chǔ)n階矩陣A時(shí)，假設(shè)A中共有l(wèi)個(gè)非零元素，則需要用一個(gè)l維向量x按先行后列的順序依次存放A中的非零元素，用一個(gè)l維向量xJ按同樣的順序依次存放A中的這些非零元素列號，同時(shí)還需要引入一個(gè)n+1維整型向量x(R)，指明A中第i行中第一個(gè)非零元素被存儲(chǔ)在x中的位置，此研究中不定方程組(4)的系數(shù)矩陣為對稱矩陣，因此實(shí)際中只需存儲(chǔ)上三角部分非零元素。對于Huge型，對應(yīng)的采用了塊行壓縮存儲(chǔ)方式。基于上述格式實(shí)現(xiàn)的矩陣與向量乘法y＝Av可寫為：y(x(R)(i))＝x(i)×v(xJ(i))end基于CSR行壓縮存儲(chǔ)方式，可分析得上述算法的計(jì)算花費(fèi)最大為O(n)，其中n為計(jì)算網(wǎng)格數(shù)。步驟8：對迭代初值采用塊行投影迭代法，求解HASM方程組，并判斷求解結(jié)果是否收斂；步驟9：當(dāng)求解結(jié)果不收斂時(shí)，對求解結(jié)果重新采用塊行投影迭代法，求解HASM方程組，并判斷求解結(jié)果是否收斂，如果收斂，執(zhí)行步驟10，否則，重新執(zhí)行步驟9；具體的，對于上述對稱不定系統(tǒng)的求解方法一般為直接法及迭代法。直接法涉及到矩陣的分解等對大規(guī)模問題并不適合。常用的迭代法為Krylov子空間迭代法，比較預(yù)處理共軛梯度法PCG、基于最小二乘的QR分解方法LSQR,及廣義極小殘差化方法GMRES等。由于此類迭代法需要將方程組系數(shù)矩陣全部一次性輸入內(nèi)存，隨著計(jì)算規(guī)模的增加，此類方法已不再使用。行投影迭代法可避免此類問題，行投影迭代法應(yīng)用廣泛，算法簡單，經(jīng)典的行投影迭代法如Cimmio算法、Kaczmar算法等。定義Hi為集合：Rn為n維實(shí)數(shù)向量構(gòu)成的空間，i＝1,…,p。因此任何即x*為p個(gè)子空間Hi的交集，為線性系統(tǒng)ATx＝b的解。設(shè)所求(1.1.5)中的對稱不定系統(tǒng)為Ax＝b，根據(jù)計(jì)算區(qū)域的網(wǎng)格刨分將矩陣A分解為塊行形式：并定義Pk(AiT)為AiT值域上的投影算子?？傮w思路是先將xk投影到p個(gè)超平面Algorithm:塊Cimmino行投影法:選取x(0),令k＝0,重復(fù)上述步驟直到收斂，begindoinparalleli＝1,…pδi(k)＝Ai+bi-PA(AiT)x(k)＝Ai+(bi-Aix(k))endparallelsetk＝k+1end此方法可避免不用每次都將矩陣元素全部輸入，每次只需輸入矩陣的一塊行，保證了大規(guī)模問題的求解可能性。但缺點(diǎn)是收斂速度慢，因此在實(shí)際中較少應(yīng)用。為此本研究基于Kaczmar算法，給出了一個(gè)新的行投影塊迭代算法。通過選取離當(dāng)前迭代點(diǎn)距離最遠(yuǎn)的塊來進(jìn)行直交投影，并將投影作為下一個(gè)迭代點(diǎn)，從而加速了收斂速度。假設(shè)當(dāng)前迭代點(diǎn)為xk，首先計(jì)算xk在所有Hi上的直交投影Pi(xk)，對于i＝1,…,p，然后比較所有的||Pi(xk)-xk||2并求出其中最大者，即找出離xk直交距離最遠(yuǎn)的塊Hi，假設(shè)||Pj(xk)-xk||2最大，即選取xk在Hj上的投影Pj(xk)作為下一個(gè)迭代點(diǎn)，即xk+1＝Pj(xk)，繼續(xù)上面的過程，得到一序列{xk}，不斷迭代直至收斂。該算法描述如下：初始化：對矩陣A進(jìn)行分塊，給定初始點(diǎn)x0∈Rn，0≤ε＜1，令k＝0步驟10：當(dāng)HASM方程組的解收斂時(shí)，進(jìn)一步判斷HASM方程組的解是否滿足高斯科達(dá)齊方程組，若不滿足，則執(zhí)行步驟6；若滿足，則根據(jù)HASM方程組的解輸出關(guān)于待測變量的高精度模擬曲面模型。。下面介紹本發(fā)明實(shí)施例的驗(yàn)證，1、數(shù)值試驗(yàn)首先以標(biāo)準(zhǔn)數(shù)學(xué)曲面為例，研究算法的有效性。比較改進(jìn)后的HASM方法與可適應(yīng)大規(guī)模求解問題的Kriging方法(FRK)、可適應(yīng)大規(guī)模求解問題的Spline方法(FRS)及反距離權(quán)重法IDW的模擬性能差異。定義域?yàn)閇0,1]×[0,1]。如圖1所示，在標(biāo)準(zhǔn)數(shù)學(xué)曲面上隨機(jī)選取1/100的采樣點(diǎn)，所有實(shí)驗(yàn)在64位機(jī)上運(yùn)行，HASM中停機(jī)準(zhǔn)則為Gauss-Codazii方程組。將計(jì)算區(qū)域刨分為1000×1000.u＝100,κ＝5.對計(jì)算誤差的衡量，采用平均絕對誤差RMSE，其計(jì)算公式如下：其中fi，為第ith個(gè)真實(shí)值與模擬值，N為采樣點(diǎn)個(gè)數(shù)，計(jì)算結(jié)果如表1所示：表1不同方法計(jì)算效率分析比較方法HASMG-IDWFRKFRSRMSE1.3076×10-41.66×10-22.20×10-35.78×10-4計(jì)算時(shí)間(s)54.06s20.11s40.00s24.79s結(jié)果表明，HASM計(jì)算精度最高，其次是FRS方法。HASM計(jì)算精度分別為IDW,FRK及FRS的127倍，17倍及4倍。但HASM的計(jì)算速度低于其他方法。這是因?yàn)镠ASM計(jì)算過程中需要對偏微分方程組有限差分離散，需要計(jì)算方程組右端項(xiàng)，及求解方程組等，計(jì)算過程更加復(fù)雜。圖2為根據(jù)(a)HASM,(b)FRK,(c)G-IDW,和(d)FRS模擬曲面與真實(shí)曲面的誤差圖，其中，可以看出IDW模擬誤差最大，最大誤差出現(xiàn)在邊界及峰值處。最大偏差為0.2139.FRK也產(chǎn)生了較大誤差，最大誤差為0.02737.FRS的最大偏差出現(xiàn)在邊界處為-0.0123.誤差最小的是HASM方法，可以看出誤差分布相對均勻。盡管HASM的計(jì)算時(shí)間多于其他方法，但其計(jì)算時(shí)間與計(jì)算網(wǎng)格成線性比例關(guān)系。表2給出了不同計(jì)算規(guī)模下HASM的計(jì)算時(shí)間。計(jì)算時(shí)間與計(jì)算網(wǎng)格的關(guān)系式可表達(dá)為：t＝1.1099×10-5n+26.7447,R2＝0.91表2.不同計(jì)算規(guī)模下HASM的計(jì)算時(shí)間計(jì)算網(wǎng)格數(shù)1million4million9million16million25million計(jì)算時(shí)間(s)54.0688.21107.54155.19339.152、實(shí)際案例以中國為研究區(qū)，采用上述方法計(jì)算過去50年1961-2010年多年平均降水的分布，空間分辨率為1km，對應(yīng)計(jì)算網(wǎng)格數(shù)為19,606,916。從HASM、FRK、IDW和FRS模擬的中國范圍多年平均降水分布結(jié)果可以看出，所有方法都能較好的反映出降水的分布趨勢。HASM及FRK方法模擬的曲面相對其他方法較為更加光滑。IDW及FRS方法出現(xiàn)了很多牛眼現(xiàn)象，且FRS方法在邊界處產(chǎn)生了較大幅度的振蕩。隨機(jī)選擇15％的作為驗(yàn)證點(diǎn)驗(yàn)證不同方法。此過程重復(fù)了10次，并計(jì)算這10次的平均RMSE。結(jié)果如表3所示?？煽吹紿ASM計(jì)算精度是最高的，其模擬精度分別我FRK,IDW及FRS的1.4,1.5及1.6倍。FRS的模擬誤差最大。IDW計(jì)算時(shí)間最少，其次是HASM.表3.實(shí)際案例中不同方法的比較方法HASMFRKIDWFRSRMSE(mm)90.87123.51133.44147.10計(jì)算時(shí)間(s)16623763172數(shù)值試驗(yàn)與實(shí)際案例，表明改進(jìn)后的HASM方法在保證精度的同時(shí)，使得其計(jì)算時(shí)間與計(jì)算網(wǎng)格數(shù)成線性關(guān)系，在存儲(chǔ)上，由于在求解約束最小二乘問題時(shí)，采用了塊行投影算法，只需存儲(chǔ)矩陣的一塊行，并采用行壓縮存儲(chǔ)方式對其進(jìn)行存儲(chǔ)，保證了求解問題的規(guī)模。以上所述僅為本發(fā)明的較佳實(shí)施例，并不用以限制本發(fā)明，凡在本發(fā)明的精神和原則之內(nèi)，所作的任何修改、等同替換、改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。當(dāng)前第1頁1 2 3