論文2-通過三維數(shù)據(jù)插值進行電阻率的擬合及估計

.z.摘要在實際問題中，我們無法對一個物體進行密集測量，一般等間隔的選取部分點進行測量，這時候就需要采用數(shù)學(xué)方法，根據(jù)已知位置的數(shù)據(jù)計算未知位置的數(shù)據(jù)。本文根據(jù)題目要求建立了兩種合理的差值加密模型——反距離加權(quán)插值模型、克里金插值模型，對文件中給出的坐標網(wǎng)格大小為10m*10m*10m的三維電阻率數(shù)據(jù)進行插值，得到坐標網(wǎng)格大小為1m*1m*1m的數(shù)據(jù)。我們借助多種數(shù)學(xué)軟件，處理大量數(shù)據(jù)并對其進行可視化處理，繪制出顏色圖反映插值效果，根據(jù)插值擬合特點對插值效果給出定量指標。針對問題一，首先對所給出的三維的電阻率數(shù)據(jù)進行分析，根據(jù)不同空間插值方法的優(yōu)勢，結(jié)合已知數(shù)據(jù)的網(wǎng)格化特點，確定插值模型——反距離加權(quán)插值模型、克里金插值模型。反距離加權(quán)插值模型算法簡單易于實現(xiàn)，但對權(quán)重函數(shù)的選擇十分敏感，而網(wǎng)格化的數(shù)據(jù)恰好規(guī)避了這一劣勢?？死锝鸩逯凳窃谧儺惡瘮?shù)理論及結(jié)構(gòu)分析基礎(chǔ)上，進行無偏、最優(yōu)估計的一種方法，不僅考慮了觀測點與待估計點的相對位置，而且考慮了各觀測點之間的相對關(guān)系，插值效果較好。利用兩種方法分別計算出空間*點(45.8,-32.7,68.2)處的電阻率數(shù)值。反距離加權(quán)插值法結(jié)果為Z1=196.18；克里金插值結(jié)果為Z2=194.13。針對問題二：利用問題一中兩種插值方法，分別計算出網(wǎng)格大小為1m*1m*1m的三維電阻率數(shù)據(jù)（由于數(shù)據(jù)量過大，具體數(shù)據(jù)見電子版1、2），分析兩種模型的計算流程可知，克里金算法計算過程較為復(fù)雜，數(shù)據(jù)量較大，計算時間與占用內(nèi)存大，而反距離加權(quán)插值計算過程較為簡單，計算時間明顯短于克里金方法，但計算精度較低。將原網(wǎng)格數(shù)據(jù)及兩種方法加密網(wǎng)格后數(shù)據(jù)利用E*cel計算出平均值與標準差，見下表。平均值E方差S2原始數(shù)據(jù)197.26893713.73324911反距離加權(quán)插值數(shù)據(jù)196.54990515.70823291克里金差值數(shù)據(jù)197.64752115.47241092針對問題三：對插值加密網(wǎng)格后的數(shù)據(jù)，利用顏色圖直觀展示。將電阻率最小值置為純藍色，中間值置為純綠色，最大值置為純紅色，中間數(shù)值采用過渡的顏色，利用matlab對切片Z=0,50分別給出原數(shù)據(jù)，兩種方法加密網(wǎng)格后數(shù)據(jù)的顏色對比圖。另外分別對切片*=82，Y=47,Z=88，給出兩種加密方法得到數(shù)據(jù)的顏色圖。通過同一切面顏色圖的比較，直觀反映出原數(shù)據(jù)、反距離加權(quán)插值加密、與克里金插值加密的插值效果。針對問題四：利用插值加密前后電阻率均值的偏離程度以及電阻率沿*，y，z三個軸向變化趨勢的偏離程度給出定量指標，定量反應(yīng)兩種插值方法的差值效果。關(guān)鍵詞：反距離加權(quán)插值變差函數(shù)加權(quán)系數(shù)克里金插值一、問題重述在實際問題中，由于技術(shù)與測量成本等原因，只能等間隔的選取部分點進行測量。實際中需要采用數(shù)學(xué)方法獲得更多位置數(shù)據(jù)，根據(jù)已知位置的數(shù)據(jù)計算未知位置的數(shù)據(jù)。根據(jù)坐標網(wǎng)格大小為10m*10m*10m的空間三維體電阻率數(shù)據(jù)，通過插值加密后獲得網(wǎng)格大小為1m*1m*1m的三維電阻率數(shù)據(jù)，必須保證插值后的電阻率數(shù)據(jù)極值不變，且坐標位置相同。另外，必須保證插值后的三維成像結(jié)果與插值前的成像結(jié)果形態(tài)基本一致，只是前者像素更高。（1）請給出符合條件的兩種計算方法，并給出相應(yīng)數(shù)學(xué)公式，證明該方法插值后的數(shù)據(jù)極值與位置不變。分別計算出空間*點(45.8,-32.7,68.2)處的電阻率數(shù)值。(2)利用（1）中兩種方法分別計算出網(wǎng)格大小為1m*1m*1m的三維電阻率數(shù)據(jù)，給出計算流程，并對兩種方法的計算量或計算復(fù)雜性進行評估。同時給出原網(wǎng)格數(shù)據(jù)及兩種方法加密網(wǎng)格后數(shù)據(jù)的平均值與標準差。(3)對加密網(wǎng)格后的直觀效果，可采用顏色圖展示。對每一幅需要對比顯示效果的圖，請將最小值置為純藍色（RGB為(0,0,255)），中間值置為純綠色（RGB為(0,255,0)），最大值置為純紅色（RGB為(255,0,0)），中間數(shù)值采用過渡的顏色，可自行設(shè)計。采用這種方法，請對切片Z=0,50分別給出原數(shù)據(jù)，兩種方法加密網(wǎng)格后數(shù)據(jù)的顏色對比圖。另外分別對切片*=82，Y=47,Z=88，請給出兩種加密方法得到數(shù)據(jù)的顏色圖。（4）對不同的插值加密方法的效果，給出定量的指標，分別計算出兩種不同加密方法得到數(shù)據(jù)的指標值，并給出評價。二、模型假設(shè)1、基本假設(shè)：1)假設(shè)data3D_.t*t所給電阻率數(shù)據(jù)背景值真實；2)被測電阻各向同性；3)不考慮環(huán)境因素對電阻率的影響；2、針對反距離加權(quán)插值法：未知值的點受較近控制點的影響比較遠控制點的影響更大；3、針對克里金插值法：本征假設(shè)：設(shè)Z（*）為空間點*處的屬性值即電阻率，Z(*)與Z(*+h)之間的相關(guān)性不依賴于在電阻內(nèi)的特定位置。實際中只需假設(shè)其1,2階矩存在且平穩(wěn)。當(dāng)區(qū)域變化量的協(xié)方差函數(shù)不存在時，只需滿足式（1）、（2），即滿足本征假設(shè)。①區(qū)域化變量Z(*)的增量[Z(*)-Z(*+h)]的數(shù)學(xué)期望為0,即:E[Z(*)-Z(*+h)]=0(*,h)(1)②增量[Z(*)-Z(*+h)]的方差存在且平穩(wěn),即:Var[Z(*)-Z(*+h)]=E[Z(*)-Z(*+h)]-{E[Z(*)-Z(*+h)]}2=E[Z(*)-Z(*+h)]2,*,h(2)三、符號說明空間中待插值點待插值點與其鄰域內(nèi)第i個點之間的距離待插值點P處的電阻率空間點*處的電阻率滯后距實驗變差函數(shù)理論變差函數(shù)塊金常數(shù)基臺值變程原始數(shù)據(jù)與加密后數(shù)據(jù)期望的偏離程度沿*軸向加密數(shù)據(jù)與原始數(shù)據(jù)變化率的偏離程度沿y軸向加密數(shù)據(jù)與原始數(shù)據(jù)變化率的偏離程度沿z軸向加密數(shù)據(jù)與原始數(shù)據(jù)變化率的偏離程度四、模型的分析、建立與求解4.1問題一問題分析由于技術(shù)和測量成本等問題，我們只能獲得有限個電阻率數(shù)據(jù)，而為了獲得更多位置的數(shù)據(jù)，則需要利用空間插值加密的方法建立電阻率空間分布模型，并保證插值前后電阻率數(shù)據(jù)極值與位置不變。反距離加權(quán)插值法設(shè)電阻中待插點為P(*,y,z),P點鄰域內(nèi)有已知散亂點Qi(*i,yi,zi),i=1,2,…,n，利用反距離加權(quán)插值法對P點的屬性值（即電阻率）進行插值，其插值原理是以插值點與樣本點間的距離為權(quán)重進行加權(quán)平均，離插值點越近的樣本點賦予的權(quán)值越大，權(quán)的大小是距離k（k一般取2）次方的倒數(shù)。即：Zp=（3）其中：di為待插點與其鄰域內(nèi)第i個點之間的距離。即：（4）4.1.3本題中克里金的研究對象(電阻率)是區(qū)域化變量,而表征區(qū)域化變量空間相關(guān)性的工具就叫做變差函數(shù)。使用已知點的屬性值所計算的變差函數(shù)叫做實驗變差函數(shù)。數(shù)學(xué)表達式為：（5）其中,表示在*軸向上的實驗變差函數(shù)值;表示在點*處的屬性值;表示在點處的屬性值。（1）計算實驗變差函數(shù)在本征假設(shè)的前提下，區(qū)域化變量Z(*)的增量[Z(*)-Z(*-h)]只依賴于向量h，并不依賴于位置*，被向量h分割的每對數(shù)據(jù)可以看成一對隨機變量的一個不同實現(xiàn)。這時，對于每個滯后距h，我們可以算出的值，其計算公式為：（6）而本題需要計算不同軸向方向的多個實驗變差函數(shù)，所以將公式變形為：（7）（2）理論變差函數(shù)球狀模型理論變差函數(shù)模型是未知的，所以需要從有限的空間已知樣本中去估計，通過不同的滯后距計算不同的，再用理論變差函數(shù)球狀模型來擬合。球狀模型的基本公式為：（8）式中：為塊金常數(shù)，為基臺值，a為變程。接近原點時，變差函數(shù)呈線性形狀，在變程處達到基臺值。其圖形如下：圖4-1-3-1（3）擬合實驗變差函數(shù)假設(shè)對于不同的滯后距已經(jīng)計算出相應(yīng)的實驗變差函數(shù)值，又對每個滯后距參與計算的數(shù)據(jù)對的數(shù)目為Ni。球狀模型表達式：（9）當(dāng)時，令，得到：，則有：（10）目標函數(shù)取最小值時，可求得其中參數(shù)。目標函數(shù)分別對求偏導(dǎo)，得到：（11）將上述方程組整理并寫成矩陣形式，有：（12）解得：（13）（4）求解待插值點函數(shù)值（14）問題求解：（1）由于以上兩種插值加密方法的插值條件保證了插值函數(shù)在插值節(jié)點處函數(shù)值與原數(shù)據(jù)相同，故插值后的電阻率數(shù)據(jù)極值與位置不變。（2）由反距離插值方法得：Z（45.8,-32.7,68.2）=196.18由克里金插值方法得：Z（45.8,-32.7,68.2）=194.134.2問題二：三維電阻率數(shù)據(jù)分別采用問題一中所建立的反距離插值模型與克里金插值模型計算出網(wǎng)格大小為1m*1m*1m的三維電阻率數(shù)據(jù)見1、2。計算流程將反距離插值法與克里金插值法的計算流程分別用圖（4-2-1）、圖（4-2-2）表示。待插值節(jié)點P（*，y，z）待插值節(jié)點P（*，y，z）P點鄰域內(nèi)Qi(*i,yi,zi)圖4-2-1反距離加權(quán)插值法計算實驗變差函數(shù)計算已知數(shù)據(jù)點*i與*0（0,0,0）之間距離di由di找出不同滯后距h，及其個數(shù)N（h）理論變差函數(shù)基本模型加權(quán)最小二乘擬合理論變差函數(shù)利用計算圖4-2-2克里金插值法計算量與計算復(fù)雜性評估（1）在反距離插值算法中：需要根據(jù)目標點到已知數(shù)據(jù)點間的距離的平方的倒數(shù)來進行加權(quán)求出加權(quán)系數(shù)，進而求出待測點的屬性值，隨著數(shù)據(jù)量及擬和精度的提高的增加，計算量以及存儲數(shù)據(jù)所需要的空間會成倍增加，計算數(shù)據(jù)所需時間較長，所需內(nèi)存空間較大。在計算復(fù)雜度方面主要是進行乘法和加法運算，算法的時間復(fù)雜度和空間復(fù)雜度較大。且差值效果一般。（2）在克里金插值算法中：主要是對變差函數(shù)的計算以及求解加權(quán)系數(shù)。在計算變差函數(shù)時，由于所給區(qū)域化變量為三維空間內(nèi)的函數(shù)，在進行計算時，主要由三個參數(shù)：搜索方向，最大搜索半徑，基本滯后距。本文采用等距搜索，由于數(shù)據(jù)量較大搜索次數(shù)較多，采用簡化方法后，以原點為基本點，減少了搜索次數(shù)。在運算中主要解一個矛盾方程組（克里金方程組），求解加權(quán)系數(shù)過程與反距離插值相同，因此克里金插值法德計算量比反距離插值法大得多，在時間和空間復(fù)雜度上均是反距離插值法的數(shù)倍，但插值效果較好。平均值與標準差平均值E方差S2原始數(shù)據(jù)197.26893713.73324911反距離加權(quán)插值數(shù)據(jù)196.54990515.70823291克里金差值數(shù)據(jù)197.64752115.47241092表4-2-34.3問題三問題分析對原始數(shù)據(jù)進行插值加密后，可以得到更加密集的數(shù)據(jù)，利用matlab繪制顏色圖，使加密后繪制圖的像素明顯提高。繪制顏色圖（1）Z=0原數(shù)據(jù)：圖4-3-2-1反距離加權(quán)插值后數(shù)據(jù)：圖4-3-2-2克里金插值后數(shù)據(jù)：圖4-3-2-3（2）Z=50原數(shù)據(jù)：圖4-3-2-4反距離加權(quán)插值后數(shù)據(jù)：圖4-3-2-5克里金插值后數(shù)據(jù)：圖4-3-2-6（3）*=82反距離插值后數(shù)據(jù)：圖4-3-2-7克里金插值后數(shù)據(jù)：圖4-3-2-8（4）Y=47反距離加權(quán)插值后數(shù)據(jù)：圖4-3-2-9克里金插值后數(shù)據(jù)：圖4-3-2-10（5）Z=88反距離加權(quán)插值后數(shù)據(jù)：圖4-3-2-11克里金插值后數(shù)據(jù)：圖4-3-2-12問題四：4.1問題分析：數(shù)據(jù)插值加密效果主要從數(shù)據(jù)均值和數(shù)據(jù)變化趨勢兩方面分析，故給出定量衡量指標與。4.2衡量指標：(1)指標：表征原始數(shù)據(jù)與差值加密后數(shù)據(jù)期望的偏離程度設(shè)Z（*,y,z）為三維空間點（*,y,z）的電阻率值。其中表示原始數(shù)據(jù)的期望，表示差值加密后的數(shù)據(jù)期望。(2)指標S：可分別估計出沿*，y，z三個方向上插值加密數(shù)據(jù)與原數(shù)據(jù)變化率的偏離程度。設(shè)三維空間點（*,y,z）取z為插值*圍內(nèi)任意常數(shù)z0，則在z=z0平面內(nèi)：設(shè)沿*軸的插值節(jié)點為，兩點間的一階差商為，利用差值方法計算點處的函數(shù)值。表示在z=z0平面內(nèi)，沿y軸方向，插值加密數(shù)據(jù)與原數(shù)據(jù)變化率的偏離程度。設(shè)沿y軸的插值節(jié)點為，兩點間的一階差商為，利用差值方法計算點處的函數(shù)值。表示在z=z0平面內(nèi)，沿*軸方向，插值加密數(shù)據(jù)與原數(shù)據(jù)變化率的偏離程度。設(shè)三維空間點（*,y,z）取*為插值*圍內(nèi)任意常數(shù)*0，則在*=*0平面內(nèi)：設(shè)沿z軸的插值節(jié)點為，兩點間的一階差商為，利用差值方法計算點處的函數(shù)值。表示在*=*0平面內(nèi)，沿z軸方向，插值加密數(shù)據(jù)與原數(shù)據(jù)變化率的偏離程度。4.3數(shù)據(jù)分析：平均值分析均值E平面原數(shù)據(jù)反距離加權(quán)插值克里金插值*=0204.4205.18204.4724109266y=0204.5205.78204.2246618829z=0199.0197.22199.0055699984表4-31：整體差值效果分析：反距離加權(quán)插值：2.07克里金插值：=0.79由平均值計算結(jié)果可知克里金插值結(jié)果較好，與原數(shù)據(jù)相差較小，說明整體差值結(jié)果與原數(shù)據(jù)吻合度較高，而反距離插值效果相對較差。2：數(shù)據(jù)整體變化趨勢分析：為了更好地衡量差值效果本模型提出用相鄰點的一階差分與原數(shù)據(jù)的一階差分進行對比，通過原數(shù)據(jù)點相鄰區(qū)間內(nèi)差值點的差商與原數(shù)據(jù)的差商的差值作為偏離程度，并繪出各軸向偏離程度圖線進行對比如下：由圖線可以看出在各軸向上都有偏離程度均有少許震蕩，但都穩(wěn)定在零偏附近，說明克里金插值效果與電阻率的實際分布規(guī)律較為吻合，進一步說明差值效果較好。五、模型評估5.1模型優(yōu)點本文采用的差值模型均為精確差值方法，較全面的運用了已經(jīng)觀測點數(shù)據(jù)，能夠較好的與實際模型相吻合，相對精確地再現(xiàn)了模型的真實值，在整體效果與變化趨勢方面均能與已有數(shù)據(jù)相吻合，通過matlab繪圖所得結(jié)果與原數(shù)據(jù)繪圖結(jié)果基本吻合，具有較高精確度，屬于較精確的差值模型。5.2缺點：在計算量和計算復(fù)雜度方面，需要根據(jù)目標點到已知數(shù)據(jù)點間的距離的平方的倒數(shù)來進行加權(quán)求出加權(quán)系數(shù)，進而求出待測點的屬性值，隨著數(shù)據(jù)量及擬和精度的提高的增加，計算量以及存儲數(shù)據(jù)所需要的空間會成倍增加，計算數(shù)據(jù)所需時間較長，所需內(nèi)存空間較大。在計算復(fù)雜度方面主要是進行乘法和加法運算，算法的時間復(fù)雜度和空間復(fù)雜度較大。5.3模型改進：（1）在計算權(quán)重時，我們采用的是簡化模型，計算精度較低，若采用其他克里金插值如指數(shù)克里金插值，泛克里金插值法，高斯模型，差值效果可能會有所提高。有待進一步測試。（2）在計算變差函數(shù)進行搜索時，考慮到算法復(fù)雜度，采用模型過于簡化，可能會影響計算結(jié)果的精度，下一步可才能考有關(guān)論文，對搜索方法進行進一步改善提高插值的精確度和縮小計算步驟，降低算法復(fù)雜度。（3）在進行算法的推導(dǎo)方面，由于時間倉促，沒有進行大量修改，因此程序運行時間較長，后期可對相關(guān)算法程序進行改進。六、參考文獻[1]翟進乾.克里金插值方法在煤層分布檢測中的應(yīng)用研究.[學(xué)位論文].**：**理工大學(xué).2008[2]蔡元菲.快速克里金算法的研究與實現(xiàn).[學(xué)位論文].**：電子科技大學(xué).2013[3]靳國棟，*衍聰，牛文杰等.距離加權(quán)反比插值與克里金插值法的比較.**大學(xué)學(xué)學(xué)報.2003[4]阮曉青，周義倉等數(shù)學(xué)建模引論..高等教育，2005.7[5]清源計算機工作室編.MATLAB基礎(chǔ)與應(yīng)用..機械工業(yè)3：相關(guān)算法程序1、Z切面程序formatlongl=1;for*2=-100:1:100;m=1;fory2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^2+(y2-y1)^2+(88-z1)^2)^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end>>V=r;pcolor(V);figure(gcf);gca=pcolor(*,y,V);set(gca,'LineStyle','none')title('Z=0數(shù)據(jù)成像')title('Z=88數(shù)據(jù)成像')2、*切面程序l=1;forz2=0:100m=1;fory2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/(sqrt((82-*1)^2+(y2-y1)^2+(z2-z1)^2))^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end>>V=r;[*,y]=meshgrid(1:201,1:101);pcolor(*,y,r);gca=pcolor(*,y,r);set(gca,'linestyle','none');>>title('*=82數(shù)據(jù)成像')3、Y切面程序l=1;forz2=0:100m=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^2+(47-y1)^2+(z2-z1)^2);endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;end4、平均值程序：B=data3D_(:,4);formatlongl=1;forz2=0:1:100;m=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((*2-*1)^4+(1-y1)^4+(z2-z1)^4);endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(l,m)=e/f;m=m+1;endl=l+1;ends1=sum(sum(r))/203015、變化趨勢偏離度：Z軸向B=data3D_(:,4);fori=1:11forj=1:21A(i)=B((10+(i-1)*21)*21+11);endendformatlongm=1;forz2=0:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/10,(y1+110)/10,(z1+10)/10)=1/((1-*1)^2+(1-y1)^2+(z2-z1)^2)^3;endendendfori=1:11;forj=1:21fork=1:21C(j,k,i)=B(k+(j-1)*21+(i-1)*21*21);endendendf=0;fori=1:4851;f=f+D(i);endE=C.*D;e=0;fori=1:4851;e=e+E(i);endr(m)=e/f;m=m+1;endfori=2:11C1(i-1)=(A(i)-A(i-1))/10;endfori=2:101C2(i-1)=(r(i)-r(i-1));endfori=1:10forj=1:10C3((i-1)*10+j)=(C2((i-1)*10+j)-C1(i));endend6、*軸向fori=1:21forj=1:21D1(i)=B((i-1)*21+11)endendformatlongm=1;for*2=-100:1:100forz1=0:10:100for*1=-100:10:100fory1=-100:10:100D((*1+110)/