莫林利,趙秀紹,鄭 偉 ,王 旭
(華東交通大學1.軟件學院; 2.土木建筑學院,江西 南昌330013)
雙線性插值在測量[1-2]、物理[3-4]、巖土[5]、計算機圖像[6-7]等方面有著廣泛的應用,是二維空間中最簡單和最實用的插值計算方法。 許多學者建立解決實際問題的二維插值數學模型,例如王軍討論了基于雙線性內插模型的高程異常值獲取的方法[1],通過內插模型建立高精度、高分辨率的區域高程異常數字模型,從而內插出該區域內任何一點或一小片區域的高程異常值。傳統進行復雜雙線性插值計算時一般利用計算機語言編程實現,但存在著對于變化的實際問題無法變更源程序的不便。 Excel 有著強大的計算功能,在工業建筑、土工試驗方面有著不可比擬的計算優勢[8-12],許多復雜或繁雜的數值計算可以通過Excel 內部函數或公式來實現。
在巖土工程中,許多土工試驗、土力學計算都涉及到二維表查詢并進行雙線性插值計算,例如第四系粉土、粘性土和粉質粘土地基承載力、利用布辛奈斯克求解半無限空間的應力,均是應用二維表格進行查表計算。當要查詢的值正好在表格中時,可輕易的查得對應的承載力或豎向應力。當要查的雙指標中有一個指標不是表中值時,為一維線性插值問題;當要查的雙指標均不在表中時,為雙線性插值。 雙線性插值時,一般要進行三次一維插值才能解決問題,計算工作量大,許多現場工作者采用估值法進行插值,估值因人而異,造成雙線性插值非唯一性。 因此,設計能夠自動進行二維表查詢和線性插值計算的Excel 表格,是快速查表求解地基承載力的有力工具,也能給其它二維表查詢和線性插值提供有價值的參考方法。
以《鐵路橋涵規范》中對于第四紀(Q4)沖洪積粘土求解承載力為例,其相關數據已經錄入Excel 數據表中,如表1所示,其中A 列為孔隙比值,第1 行為液性指數值,其承載力值如表1中B2 到N8 單元格區域所示。

表1 第四紀沖洪積土承載力Excel 數據Tab.1 Soil bearing capacity data for Q4 alluvial soil
某工程地基為Q4沖洪積土,地下水位在地表下2 m 處,水位以下含水率w=32%,土粒比重GS=2.68,土的液限wL=40%,土的塑限wP=18%,地下水位以上土體重度γ=18 kN·m-3,寬度4 m 的條形基礎埋置于地表下5 m 處,按《鐵路橋涵規范》確定該地基承載力。
對于此類問題,首先根據土力學中公式計算e 和IL,計算公式如式(1)和式(2)所示

很顯然,所得孔隙比e 和液性指數IL均不能在表1所示的數據表中直接查出。 許多現場工作者認為插值計算公式(見式(6))計算量大因而一般會利用表1中的值估算一個承載力,而估值因人而異,可能造成承載力求解不夠準確。
設要查表的e3在e1和e2之間,IL3在IL1和IL2之間, 與承載力相關關系如表2所示。
如圖1所示,當液性指數IL1不變時,已知孔隙比e1值對應的承載力為f1,孔隙比e2對應的承載力為f2,求孔隙比e3對應的承載力f5。

表2 孔隙比、液性指數與承載力對應關系Tab.2 Correspondence between void ratio, liquid index and bearing capacity

圖1 二維插值原理圖Fig.1 Principle diagram for two-dimensional interpolation
在圖1(a)中,e2>e3>e1,可知△XVW?△XUY,可得式(3)

可以證明,當e3<e1或e3>e2時,用式(3)計算結果仍然正確,式(3)即為求解線性插值的通用公式。 同理根據圖1中(b)可得式(4)

根據式(3)和式(4)可得e3,IL1時承載力為f5,e3、IL2時承載力為f6。 此時設孔隙比e3不變,IL3在IL1與IL2之間插值,則仍根據線性插值計算方法和表1可得f7,如式(5)

則f7即為e3、IL3時對應的承載力,把式(3)和(4)代入式(5)可得線性插值通式如式(6)所示

當計算e3=0.858, IL3=0.636 對應的承載力f7時,利用表1和式(6)進行二維插值計算,查表1可得查找的最近區域為e1=0.8, e2=0.9, IL1=0.6, IL2=0.7;對應的承載力為f1=230, f2=190, f3=210, f4=180,代入式(6)可得f7=201.688 kPa, f7就是此地基的基本承載力。 然后根據深度和基礎寬度修正可得地基的容許承載力,此處不再詳述。
利用公式(6)進行線性插值時每次均要進行公式推導,計算工作量大且易出錯。此時利用Excel 中TREND函數實現自動插值計算,將大大減少計算工作量且語法明確。 TREND 函數用途是:返回一條線性回歸擬合線的一組縱坐標值(y 值),即找到適合給定的數組known_y′s 和known_x′s 的直線,并返回指定數組new_x′s值在直線上對應的y 值。 TREND 函數語法為

如圖2所示,當TREND 函數所用參數數組僅有兩組數據x1、x2和y1、y2時,則其返回的線性擬合的直線一定是由m(x1,y1)和n(x2,y2)兩點所連的直線,且其擬合的相關系數R2=1。 當所求第3 點C 點的橫坐標為x3時, 則函數返回值即為o 點的縱坐標y3, 即o 點是由m 點和n 點線性插入的。 實現上述計算的公式為“=TREND({y1,y2},{x1,x2},x3)”,在Excel 中大括號代表數組。
在具體應用于線性插值時,known_y′s、known_x′s 兩個參數必須是數組或Excel 連續單元格引用。 因此必須將要線性插值的原始數據布置到相鄰的單元格。 根據表1可得,e3=0.858 的位置在0.8~0.9 (表中列表值) 之間,IL3=0.636 在0.6~0.7 之間, 則對應的f1=230 kPa, f2=190 kPa, f3=210 kPa,f4=180 kPa,數據如圖3中R5~T7 區域所示。
計算e3=0.858, IL1=0.6 時的承載力f5的方法為:在S8 單元格中輸入公式(7)所示的計算公式,計算結果為206.8 kPa。 公式(7)如下所示


圖2 Trend 函數應用于線性插值原理Fig.2 Principle of TREND function applied to the linear interpolation
式(7)中,S6:S7 為f1和f2組成的承載力引用,為參數known_y′s 的具體形式;R6:R7 為e1和e2組成的孔隙比引用, 為參數known_x′s 的具體形式;R2 為e3,新的要插值的孔隙比,為參數new_x′s 的具體形式。
計算e3=0.858,IL2=0.7 時的承載力f6的方法為在T8 單元格輸入公式(8)所示的計算公式,計算結果為192.6 kPa。 公式(8)如下所示

計算e3=0.858,IL3=0.636 時的承載力f7的方法為在T9 單元格輸入式(9)所示的計算公式,計算結果f7=201.69 kPa。 公式(9)如下所示

式(8)和(9)各參數意義與式(7)類似。圖3所示就是二維線性插值的結果,可見進行二維插值需要調用三次TREND 函數。

圖3 使用TREND 函數實現二維插值算法圖Fig.3 Application of TREND function to achieve twodimensional interpolation
在上面的表格已經實現了二維線性插值,但在圖3中需要的f1, f2, f3, f4是手動查出來并輸入的。若e3,IL3的值變更時, f1, f2, f3, f4又要重新查找,如何查找e3和IL3對應的最近的f1, f2, f3, f4是待解決的問題。
Excel 查表時,在孔隙比和液性指數均為有序排列的情況下,就是確定指定孔隙比和液性指數在表中對應的位置區間。 在Excel 中,查找數據在某個單元格區域的位置可以使用MATCH 函數,其使用語法為

該函數返回符合特定值特定順序的項在數組中的相對位置,三個參數的含義分別為:lookup_value 是需要在數據表(lookup_array)中查找的值,lookup_array 包含有所要查找數值的連續的單元格區域,match_type 為1代表升序排列,0 代表任意順序排列,-1 代表降序排列。
查e3=0.858 時,首先要查e3在0.5~1.1 中的位置,首先確定小于e3=0.858 的最大值為0.8,對應數組[0.5,0.6,0.7,0.8,0.9,1.0,1.1]中的位置就是4,其查找方法可以在T2 單元格中用式(10)來計算

同理查得小于IL3=0.636 在數組[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2]的位置計算方法為在T3 單元格中用式(11)來計算,計算結果為7。 公式(11)如下所示

T2 和T3 單元格計算結果顯示4 和7,即0.8(小于0.858 的最大值)和0.6(小于0.636 的最大值)所在的位置,在承載力區域B2:N8 的相對位置為第4 行和第7 列,其查得的承載力值為f1, 5 行7 列對應的為f2,4行8 列對應的為f3,5 行8 列對應的為f4。
根據上面的查找方法,可得第4 行第7 列承載力值為f1, 5 行7 列對應的為f2,4 行8 列對應的為f3,5行8 列對應的為f4。 則可以通過表1查知f1, f2, f3, f4的具體數據分別為230,190,210,180。 根據相對位置獲得f1~f4的承載力數據,需要用到INDEX 函數,其語法為

該函數可返回指定單元格數組(array)指定行(row_num)指定列(column_num)的數據。 具體應用時,可以在S6,S7,T6,T7 單元格中分別輸入以下的式(12)~式(15),即獲得插值區域

式(12)~式(15)實際獲得的插值區域為H5:I6 的引用,在S6,S7,T6,T7 顯示的數據即為所查詢的f1=230,f2=190, f3=210, f4=180,單位均kPa。 則再用前面所述的TREND 函數進行線性插值就可自動完成雙指標的二維自動查表和自動插值計算。
插值計算時,當輸入的e 小于0.5 時,則無法查找到相應的值,計算結果顯示“#N/A”,即在函數中沒有可用數值時,產生錯誤值“#N/A”,因此對可能遇到的邊界問題要特殊處理。 處理方法為當e 小于0.5 時,按0.5 的位置計算,當e 大于1.1 時,按1.1 前一行位置計算小值,需要按1.0~1.1 之間插值。 具體實現方法為:把T2 和T3 單元中的公式修正成式(16)和式(17)

在式(16)中,當R2 單元格的值小于0.5 時,顯然在數據列表中查不到此值,則強制其相對位置為1;當R2大于等于1.1 時,相對位置為1.0 所對應的相對位置,此時強制為COUNT(A2 : A8)-1 即為1.0 對應的位置,其中COUNT(A2 : A8)為孔隙比數據表的總行數。 當R2 在0.5~1.1 之間時,則按MATCH(R2 : A2: A8,1)求解其相對位置。 式(17)解釋與式(16)類似。
以上所有設計結果如圖3所示,經過測試其計算結果與式(6)相同,說明文中的Excel 方法是可靠的。
論文結合巖土工程中的承載力計算的實際工程問題推導了二維線性插值的計算公式,但計算公式復雜且計算容易出錯,不宜在工程中采用。 本文提出的Excel 函數法語法簡單,可以實現自動查表和自動二維線性插值。 經過分析可得出以下結論:
1) TREND 函數參數為(y1,y2)和(x1, x2)兩組數據時,其擬合的直線就是(x1,y1)和(x2,y2)兩點的聯線,且相關系數為1,因此可以用作線性插值計算。TREND 函數使用三次分別計算出e3, IL1對應值f5;e3,IL2對應值f6;e3,IL3對應值f7,即可完成二維表格的二維線性插值。
2) MATCH 函數為Excel 中查找e3(或IL3)在一維數組中小于e3數據表中最大數位置的函數,當所查找值不在給定的數組范圍時,要進行邊界特殊處理。當所查值小于表中最小值時,則位置設為1,當所查找值大于表中最大值,則位置為數組成員數減1。
3) MATCH 函數和INDEX 函數配合可以查找e 和IL在承載力表中最近的插值區域,再配合TREND 函數可以實現自動查表和自動插值,經過與推導公式對比證明文中提出的Excel 方法簡單可靠,設計的計算表格可以作為承載力計算的實用工作。
[1] 王軍.基于雙線性插值模型的高程異常值的獲取[J].云南民族大學學報:自然科學版,2013,22(1): 75-78.
[2] 王勝兵,戴明強,黃登斌.基于雙線性插值擬合的山形曲面面積計算[J].兵工自動化,2012,31(3):42-43.
[3] 馬德堂,楊春雨,李緒宣.基于雙線性插值的三維橫向各向同性介質初至波射線追蹤[J].地球物理學進展,2014,29(3):1201-1205.
[4] 郭曉文,張有會,王志巍,等.參考曲率的雙線性插值濾波[J].計算機工程與應用,2012,48(35):161-165.
[5] 趙元海,史慶濤.土壓力在位移影響下的近似求解方法[J].華東交通大學學報,2013,30(6):49-54.
[6] 趙煌,彭勇.雙線性插值算法的優化及其應用[J].電視技術,2012,36(17):30-32.
[7] 吳桂萍,吳巍,王成,等.基于雙線性插值的魚眼圖像校正方法[J].計算機應用與軟件,2012,29(2):122-124.
[8] 李健,陳國強.Excel 在土工試驗資料整理中的應用[J].電力勘測設計,2011,18(4):8-10.
[9] 張永成,王洪輝,譚桂花.基于Excel VBA 的圖解粒度參數計算[J].成都理工大學學報:自然科學版,2010,37(6):650-653.
[10] 曲蒙,李新勇,唐永寧,等.基于EXCEL 的礦槽料倉計算[J].工業建筑,2012,42(s1):152-154.
[11] 喬玉娥,鄭世棋,翟玉衛,等.Excel 曲線擬合功能在測量不確定度評定中的應用[J].上海計量測試,2013,41(5):25-29.
[12] 路云龍,張轉梅,李文鈺.利用Excel 求解插值問題[J].長春大學學報,2012,22(10):1228-1229.