張永奇,韓美濤,韓曉飛,鄭增記,翟宏光

(1.陜西省地震局，陜西 西安 710068；2.中國(guó)地震局地質(zhì)研究所 地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029)

Kriging法在區(qū)域重力場(chǎng)插值中的適用性研究

張永奇1,2,韓美濤1,韓曉飛1,鄭增記1,翟宏光1

(1.陜西省地震局，陜西 西安 710068；2.中國(guó)地震局地質(zhì)研究所 地震動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029)

由于重力觀測(cè)網(wǎng)點(diǎn)位分布不均勻，密度不夠，在精細(xì)刻畫(huà)區(qū)域重力場(chǎng)變化時(shí)受到很大限制，因此需采用合理的插值算法。文中介紹Kriging插值算法和變異函數(shù)理論模型的基礎(chǔ)上，結(jié)合陜西地區(qū)2011—2012年離散重力變化數(shù)據(jù)，利用最小二乘法得到各個(gè)理論變異函數(shù)模型的擬合參數(shù)，并將實(shí)驗(yàn)變異函數(shù)模型用于Kriging插值算法，同時(shí)綜合考慮交叉驗(yàn)證精度統(tǒng)計(jì)、插值精度統(tǒng)計(jì)結(jié)果，以此來(lái)研究Kriging插值算法的適用性，研究結(jié)果表明：基于最小二乘法獲得的球形模型、指數(shù)模型擬合參數(shù)精度高，高斯模型、對(duì)數(shù)模型精度最差；基于球形模型、指數(shù)模型的實(shí)驗(yàn)變異函數(shù)用于Kriging插值算法得到的插值結(jié)果精度最高、圖像平滑自然、異常區(qū)明顯，是一種適合陜西地區(qū)重力離散數(shù)據(jù)進(jìn)行插值計(jì)算的有效算法。

Kriging法；變異函數(shù)模型；陜西地區(qū)；重力場(chǎng)；插值

重力資料在地震長(zhǎng)期預(yù)測(cè)預(yù)報(bào)中占有重要地位，合理的重力網(wǎng)布設(shè)是獲取真實(shí)可靠資料的基礎(chǔ)。但由于重力網(wǎng)布設(shè)易受各種因素的限制，點(diǎn)位分布往往呈現(xiàn)不均勻、密度不夠等特點(diǎn)，這對(duì)于精細(xì)刻畫(huà)某個(gè)區(qū)域重力場(chǎng)變化來(lái)說(shuō)，顯得力不從心。因此，針對(duì)離散觀測(cè)值采用合理的插值算法進(jìn)行插值計(jì)算就成為一個(gè)重要的解決途徑。目前，重力場(chǎng)插值應(yīng)用較多的算法有最小二乘配置法、樣條函數(shù)法、Kriging法等[1-3]。Kriging法是一種依據(jù)已知點(diǎn)進(jìn)行空間插值的地質(zhì)學(xué)統(tǒng)計(jì)方法，該方法是由南非礦產(chǎn)地理學(xué)家P.G.Krige最先提出，該法綜合考慮空間點(diǎn)的結(jié)構(gòu)性和隨機(jī)性，估值具有最優(yōu)、線性、無(wú)偏等特點(diǎn)，在地學(xué)領(lǐng)域得到廣泛的應(yīng)用，并取得很好的效果[4-9]。目前，基于Kriging法進(jìn)行數(shù)據(jù)插值的研究主要集中在幾種插值算法的比較，研究表明Kriging法具有明顯優(yōu)越性[3,6,9-10]；基于不同變異函數(shù)的Kriging法進(jìn)行重力場(chǎng)插值的研究相對(duì)較少。變異函數(shù)是Kriging法進(jìn)行插值的核心，常用的變異函數(shù)理論模型有球面模型、高斯模型、指數(shù)模型等[12-14]，變異函數(shù)模型選擇是否恰當(dāng)，直接影響插值效果的好壞，甚至影響計(jì)算結(jié)果的可靠性，因此對(duì)不同變異函數(shù)模型的適用性進(jìn)行研究，有利于提高Kriging法插值的精度和可靠性?；谏鲜龇治?，本文在介紹Kriging法及變異函數(shù)模型的基礎(chǔ)上，以陜西區(qū)域離散重力變化數(shù)據(jù)作為研究對(duì)象，對(duì)基于不同變異函數(shù)的Kriging法插值結(jié)果進(jìn)行精度統(tǒng)計(jì)，以此評(píng)價(jià)變異函數(shù)模型的適用性，最終得到適合該區(qū)域重力場(chǎng)插值的可靠算法。

1 數(shù)據(jù)處理理論與模型

1.1 普通Kriging法

假定Z(x)服從二階平穩(wěn)假設(shè)，其在空間位置x1,x2,…,xn有n個(gè)觀測(cè)值Z(xi),i=1,2,…,n,則在未知點(diǎn)x0處，其觀測(cè)值Z(x0)可由已知的n個(gè)觀測(cè)值Z(xi)估計(jì)。

(1)

式中:λi為待定的權(quán)重系數(shù)，為達(dá)到線性無(wú)偏估計(jì)，使估計(jì)方差最小，權(quán)系數(shù)λi可由Kriging方程組決定。

(2)

式中:μ是拉格朗日乘數(shù);γij表示第i,j個(gè)觀測(cè)點(diǎn)間向量所對(duì)應(yīng)的變異函數(shù)值;γi0表示第i個(gè)觀測(cè)點(diǎn)到估值點(diǎn)間向量對(duì)應(yīng)的變異函數(shù)值，這些量可通過(guò)實(shí)驗(yàn)變異函數(shù)擬合理論模型計(jì)算得到。在式(2)約束條件下，求出λj，再回代式(1)，即可根據(jù)已知點(diǎn)觀測(cè)數(shù)據(jù)內(nèi)插出待定點(diǎn)的觀測(cè)估值[14]。

1.2 變異函數(shù)

(3)

式中:γ*(h)為理論變異函數(shù)值γ(h)的估計(jì)值;N為被向量h相隔的實(shí)驗(yàn)數(shù)據(jù)對(duì)的數(shù)目。

按式(3)計(jì)算出空間相距h距離的變異函數(shù)值γ*(h)后，再選擇合適的變異函數(shù)理論模型進(jìn)行擬合，進(jìn)而求出適合研究區(qū)域的變異函數(shù)理論模型。

1.3 變異函數(shù)理論模型

變異曲線圖表示的是一定滯后距h的變異函數(shù)值γ*(h)與該h的對(duì)應(yīng)圖，如圖1所示。圖中C0稱(chēng)為塊金效應(yīng)，表示h很小時(shí)兩點(diǎn)x和x+h之間空間信息的變化；a稱(chēng)為變程，其大小反映了研究對(duì)象中某一區(qū)域變化量的變化程度。CS稱(chēng)為總基臺(tái)值，反映了某區(qū)域化變量在研究范圍內(nèi)變異的強(qiáng)度，它是最大滯后距的可遷移性變異函數(shù)的極限值，當(dāng)h→∞時(shí)，γ(∞)=C(0)=Var[Z(x)]=CS，即當(dāng)h→∞時(shí)，變異函數(shù)值近于先驗(yàn)方差C(0)，當(dāng)無(wú)塊金效應(yīng)C0時(shí)，CS=C，當(dāng)有塊金效應(yīng)時(shí)，CS=C+C0；而C稱(chēng)為基臺(tái)值，它是先驗(yàn)方差與塊金效應(yīng)(常數(shù))之變差，C=CS-C[14]。

圖1 變異曲線圖

變異函數(shù)理論模型分為有基臺(tái)值和無(wú)基臺(tái)值兩大類(lèi)。其中有基臺(tái)類(lèi)模型包括球形模型、高斯模型、指數(shù)模型、純塊金模型、線性模型(有基臺(tái))，有基臺(tái)模型一般認(rèn)為區(qū)域變量滿足二階平穩(wěn)假設(shè)；無(wú)基臺(tái)類(lèi)模型包括冪函數(shù)模型、對(duì)數(shù)模型、線性模型(無(wú)基臺(tái))等，一般認(rèn)為只滿足本征假設(shè)[14]。具體模型如下：

1)球面模型:

(4)

2)高斯模型:

(5)

3)指數(shù)模型:

(6)

4)純塊金效應(yīng)模型:

(7)

該模型對(duì)純隨機(jī)變量才適用。

5)有基臺(tái)的線性模型:

(8)

6)冪模型:

γ(h)=hθ,0<θ<2.

(9)

當(dāng)θ=1，即為無(wú)基臺(tái)線性模型，γ(h)=ah,h>0，a為常數(shù)。

7)對(duì)數(shù)模型:

γ(h)=log(h).

(10)

由于h→0時(shí)，log(h)→-∞，這與變異函數(shù)性質(zhì)不符合，因此在區(qū)域化變量為離散點(diǎn)時(shí)不適用。

利用式(3)求出各γ(h)，按最小二乘擬合法求出式(4)～式(10)中的參數(shù)C0，C，a，則可確定出樣本空間的實(shí)驗(yàn)變異函數(shù)模型,也就是適合研究區(qū)域的變異函數(shù)模型[11-14]。

2 實(shí)例及結(jié)果分析

2.1 研究區(qū)概況

本文選取陜西地區(qū)2011～2012年離散重力變化數(shù)據(jù)作為研究對(duì)象，重力網(wǎng)范圍為32.42°～36.25°N，106.15°～110.56°E，總共包括221個(gè)點(diǎn)，其中關(guān)中地區(qū)測(cè)點(diǎn)密度最大，平均點(diǎn)間距約10～15 km；陜南地區(qū)測(cè)點(diǎn)密度次之，平均點(diǎn)間距約20～30 km；陜北地區(qū)點(diǎn)位分布最為稀疏，平均點(diǎn)間距達(dá)到40～50 km。陜西地區(qū)地形地貌差異較大，陜北地區(qū)以黃土塬為主，海拔大致為800～1 300 m；關(guān)中地區(qū)為斷陷盆地，平均海拔約520 m；陜南地區(qū)為秦巴山地，海拔基本上為1 500～2 000 m。為了同時(shí)顧及點(diǎn)位分布不均以及地形差異的影響，本文選擇Kriging法對(duì)離散觀測(cè)值進(jìn)行插值，其目的既考慮觀測(cè)值的空間性也考慮觀測(cè)點(diǎn)間的相關(guān)性。通過(guò)對(duì)不同變異函數(shù)模型的參數(shù)進(jìn)行擬合，獲得不同的實(shí)驗(yàn)變異函數(shù)，并將其用于Kriging插值算法，通過(guò)精度對(duì)比分析，得到可靠性和精度都比較滿意的Kriging插值算法。

2.2 研究結(jié)果與分析

首先利用研究區(qū)分布的211個(gè)離散點(diǎn)數(shù)據(jù)，建立適合本研究區(qū)的實(shí)驗(yàn)變異函數(shù)。根據(jù)離散點(diǎn)空間分布特征及其相關(guān)性，獲得實(shí)驗(yàn)變異函數(shù)參數(shù)，具體參數(shù)見(jiàn)表1。由表1知，變異函數(shù)最大滯后距為1.9°，滯后組數(shù)劃分為25組，滯后寬度為0.076°，變異函數(shù)高值、滯后方向、滯后方向公差分別為644、60°、90°。其次求取變異函數(shù)模型參數(shù)，選擇變異函數(shù)理論模型(球形模型、指數(shù)模型、冪函數(shù)模型、塊金效應(yīng)模型等)，并給定初始參數(shù)值，利用最小二乘擬合法對(duì)各模型參數(shù)進(jìn)行迭代處理，設(shè)置迭代次數(shù)為50，經(jīng)過(guò)計(jì)算獲得各組合變異函數(shù)模型(塊金效應(yīng)+其它模型)參數(shù)并繪制擬合圖，如表2、圖2所示。表2中參數(shù)就是實(shí)驗(yàn)變異函數(shù)模型所需的參數(shù)C0，C，a，將參數(shù)代入模型(4)～模型(10)，就可以獲得適合研究區(qū)的變異函數(shù)模型，根據(jù)模型(1)～模型(3)就可以求取研究區(qū)任一點(diǎn)的重力異常值。從圖2直觀可知，球形模型、指數(shù)模型、線性模型、冪函數(shù)模型擬合效果較好，高斯模型、對(duì)數(shù)模型擬合效果較差，研究表明對(duì)數(shù)模型不適合進(jìn)行離散點(diǎn)擬合[15]。

表1 實(shí)驗(yàn)變異函數(shù)參數(shù)

表2 不同變異函數(shù)模型擬合參數(shù)

備注：線性模型斜率slope=225;冪函數(shù)模型因變量θ=0.692

圖2 不同模型的變異函數(shù)擬合圖

為了進(jìn)一步比較不同變異函數(shù)模型的適用性，采用交叉驗(yàn)證的方法進(jìn)行分析。將221個(gè)離散點(diǎn)數(shù)據(jù)全部作為檢核點(diǎn)，計(jì)算擬合點(diǎn)值與離散數(shù)據(jù)點(diǎn)值的殘差精度以及擬合點(diǎn)自身精度，精度統(tǒng)計(jì)如表3所示。分析表3發(fā)現(xiàn)，從擬合點(diǎn)的均方根來(lái)分析，精度從高到低順序?yàn)椋焊咚鼓Ｐ?指數(shù)模型>球形模型=冪函數(shù)模型>線性模型>對(duì)數(shù)模型；但從擬合點(diǎn)值與離散數(shù)據(jù)點(diǎn)值的殘差均方根來(lái)分析，精度從高到低順序?yàn)椋呵蛐文Ｐ?冪函數(shù)模型>線性模型>指數(shù)模型>高斯模型>對(duì)數(shù)模型；從標(biāo)準(zhǔn)差來(lái)分析，也可以得到和均方根大致相同的結(jié)論，因此綜合4個(gè)精度指標(biāo)考慮，球形模型、冪函數(shù)模型擬合精度較高，更加適合研究區(qū)離散數(shù)據(jù)Kriging法插值，而高斯模型、對(duì)數(shù)模型擬合精度最差。

表3 不同變異函數(shù)模型的交差驗(yàn)證精度統(tǒng)計(jì)

利用上述各個(gè)實(shí)驗(yàn)變異函數(shù)模型進(jìn)行Kriging法插值，總共獲得8 700個(gè)節(jié)點(diǎn)的插值數(shù)據(jù)，緯度方向插值節(jié)點(diǎn)為100行，經(jīng)度方向插值節(jié)點(diǎn)為87列。根據(jù)插值結(jié)果繪制插值等值線圖，超出插值邊界最大值的數(shù)據(jù)采用白化值，如圖3所示，圖3中(a)～(f)分別代表球形模型、高斯模型、指數(shù)模型、線性模型、冪函數(shù)模型、對(duì)數(shù)模型的插值結(jié)果。同時(shí)統(tǒng)計(jì)基于不同變異函數(shù)模型對(duì)應(yīng)的Kriging法插值精度，具體統(tǒng)計(jì)指標(biāo)有插值花費(fèi)時(shí)間、插值最大值、最小值、均值、標(biāo)準(zhǔn)差，具體見(jiàn)表4。從圖3直觀看出，(a)、(c)～(e)結(jié)果較好，很好地體現(xiàn)了重力場(chǎng)異常變化的細(xì)節(jié)特征，是對(duì)離散觀測(cè)值的有效插值和擬合，而(b)、(f)插值效果最差，(b)反映的是研究區(qū)的趨勢(shì)變化，其中的細(xì)節(jié)體現(xiàn)不完整，而(f)完全不適合離散點(diǎn)插值，對(duì)重力場(chǎng)異常變化反映不明顯，而且量值遠(yuǎn)遠(yuǎn)大于采用其它模型進(jìn)行插值的結(jié)果。從細(xì)節(jié)上分析，基于球形模型和指數(shù)模型的插值結(jié)果，“牛眼”較少，曲線比較光滑，由南向北過(guò)渡比較自然；而基于線性模型和冪函數(shù)模型的插值結(jié)果反映細(xì)節(jié)比較清晰，形成很多局部異常，因此基于這兩種模型得到的插值結(jié)果最好是進(jìn)行高通濾波，得到更加真實(shí)的異常變化圖。

由表4知，幾種模型花費(fèi)的時(shí)間差別不大，對(duì)數(shù)模型用時(shí)最多，精度最差。從插值的最大、最小值來(lái)看都沒(méi)有超過(guò)原始離散點(diǎn)的最大值72.7、最小值-52.7，表明插值數(shù)據(jù)可靠性較高，只進(jìn)行內(nèi)插計(jì)算沒(méi)有外推計(jì)算。從標(biāo)準(zhǔn)差來(lái)看高斯模型精度最高，但這是一種“虛假”的高精度，并不適合研究區(qū)的重力異常插值；其余5種模型中球形模型和指數(shù)模型精度較高，與圖3結(jié)果表現(xiàn)一致，是真實(shí)反映插值效果的指標(biāo)；冪函數(shù)模型、線性模型精度稍差，應(yīng)該是局部區(qū)域插值精度差導(dǎo)致的。

表4 不同變異函數(shù)模型的插值精度統(tǒng)計(jì)

圖3 基于不同變異函數(shù)模型的kriging重力場(chǎng)插值圖(ugal)

3 結(jié) 論

本文針對(duì)陜西地區(qū)2011～2012年211個(gè)離散重力變化值，基于不同變異函數(shù)模型構(gòu)建Kriging插值法進(jìn)行模型適用性研究，得出如下結(jié)論：

1)對(duì)于本次試驗(yàn)數(shù)據(jù)，最大滯后距為1.9°，滯后組數(shù)分為25組，滯后方向?yàn)?0°；經(jīng)反復(fù)測(cè)試變異函數(shù)模型選擇純塊金模型和其它模型的組合，擬合效果較好，單獨(dú)模型擬合效果不佳；且利用最小二乘擬合迭代計(jì)算，迭代次數(shù)設(shè)置為50次，基本上可以獲得最優(yōu)解。

2)利用最小二乘擬合獲得的各種模型的參數(shù)差異較大，同時(shí)發(fā)現(xiàn)數(shù)據(jù)具有一定程度的各向異性，量值皆為2，各向異性方向變化范圍為15.61°～16.33°。變異函數(shù)模型參數(shù)精度統(tǒng)計(jì)使用交叉驗(yàn)證方法，交叉驗(yàn)證時(shí)采用全部離散點(diǎn)數(shù)據(jù)，其優(yōu)點(diǎn)是避免由于選擇部分離散點(diǎn)而帶入人為的統(tǒng)計(jì)偏差。交叉驗(yàn)證結(jié)果顯示球形模型、冪函數(shù)模型精度最高，均方根和標(biāo)準(zhǔn)誤差均最小，說(shuō)明這兩種模型建立參數(shù)選取合理，適合用于陜西地區(qū)離散重力變化網(wǎng)格化插值。

3)利用基于不同變異函數(shù)模型的Kriging法進(jìn)行插值時(shí)，考慮插值曲線的光滑性，整個(gè)研究區(qū)域劃分為8 700個(gè)節(jié)點(diǎn)，插值得到各個(gè)節(jié)點(diǎn)的重力變化值。為了減小邊緣效應(yīng)的影響，超出插值最大值的區(qū)域采用白化值。由圖3、表4可以得出：采用高斯模型、對(duì)數(shù)模型的Kriging法進(jìn)行插值效果最差，不適合陜西地區(qū)離散重力變化插值；插值效果最好的是采用球形模型、指數(shù)模型，精度最高、誤差最小、圖像顯示曲線光滑自然，且沒(méi)有出現(xiàn)插值外推現(xiàn)象。

[1] 祝意青，胡斌，李輝，等.用雙三次樣條函數(shù)模擬青藏高原東北緣重力場(chǎng)動(dòng)態(tài)圖像[J].中國(guó)地震，2005，21(2)：165-171.

[2] 祝意青，陳斌，張希，等.景泰5.9級(jí)地震前后的重力變化研究[J].中國(guó)地震，2001，17(4)：356-363.

[3] 李振海，汪海洪.重力數(shù)據(jù)網(wǎng)格化方法比較[J].大地測(cè)量與地球動(dòng)力學(xué)，2010，30(1)：140-144.

[4] 陳小斌.中國(guó)陸地現(xiàn)今水平形變狀況及其驅(qū)動(dòng)機(jī)制[J].中國(guó)科學(xué)(D輯：地球科學(xué))，2007，37(8)：1056-1064.

[5] 朱守彪，蔡永恩，石耀霖,等.青藏高原及鄰區(qū)現(xiàn)今地應(yīng)變率場(chǎng)的計(jì)算及其結(jié)果的地球動(dòng)力學(xué)意義[J].地球物理學(xué)報(bào)，2005，48(5)：1053-1061.

[6] 曾懷恩，黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[J].測(cè)繪工程，2007，16(5)：5-8.

[7] 劉曉剛，吳曉平，王寶軍，等.衛(wèi)星重力梯度數(shù)據(jù)的網(wǎng)格化方法研究[J].大地測(cè)量與地球動(dòng)力學(xué)，2010，30(6)：60-65.

[8] 劉曉霞，江在森，武艷強(qiáng)，等.Kriging方法在GPS速度場(chǎng)網(wǎng)格化和應(yīng)變率場(chǎng)計(jì)算中的適用性[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2014，39(4)：457-461.

[9] 劉兆平，楊進(jìn)，武煒，等.地球物理數(shù)據(jù)網(wǎng)格化方法的選取[J].物探與化探，2010，34(1)：93-97.

[10] 姚道榮，鐘波，汪海洪，等.最小二乘配置與普通Kriging法的比較[J].大地測(cè)量與地球動(dòng)力學(xué)，2008，28(3)：77-82.

[11] 熊俊楠，馬洪濱.變異函數(shù)的自動(dòng)擬合研究[J].測(cè)繪信息與工程，2008，33(1)：27-29.

[12] 李章林，王平，李冬梅，等.實(shí)驗(yàn)變差函數(shù)計(jì)算方法的研究與運(yùn)用[J].國(guó)土資源信息化，2008(2)：10-14.

[13] 郭泉河，李秀海.不同變異函數(shù)的泛Kriging法的GPS高程擬合結(jié)果[J].黑龍江工程學(xué)院學(xué)報(bào)(自然科學(xué)版)，2011， 25(4)：26-28.

[14] 沈云中，陶本藻.實(shí)用測(cè)量數(shù)據(jù)處理方法[M].北京：測(cè)繪出版社，2012：194-214.

[15] 侯景儒，尹鎮(zhèn)南，李維明，等.實(shí)用地質(zhì)統(tǒng)計(jì)學(xué)[M].北京：地質(zhì)出版社，1998：31-68.

Research on the applicability of Kriging methodin regional gravity field interpolation

ZHANG Yongqi1,2, HAN Meitao1, HAN Xiaofei1, ZHENG Zengji1，ZHAI Hongguang1

(1. Shaanxi Provincial Earthquake Administration, Xi’an 710068,China;2. State Key Laboratory of Earthquake Dynamics，Institute of Geology, China Earthquake Administration，Beijing 100029，China)

Because the gravity observation network distribution is not uniform and the point density is not enough, it is very limited to describe the variation of regional gravity field, which needs to use the reasonable interpolation algorithm. Based on the detailed introduction of the Kriging interpolation algorithm and the variation function theory, this paper unifies the Shaanxi area 2011～2012 discrete variation of gravity data, obtains the fitting parameters of theoretical variation function model by least square method, and uses the experimental variation model for Kriging interpolation algorithm. Considering the cross validation accuracy statistics, the interpolation precision results, the application to study Kriging interpolation algorithm, the result shows that: the fitting parameter accuracy of spherical model and exponential model are high based on the least square method, and the accuracy of Gauss model and logarithmic model is the worst; the experimental variation function based on the exponential model and spherical model used in the Kriging interpolation algorithm, can get the highest accuracy of interpolation, meanwhile the image is smooth and the abnormal area is obvious, so it is a suitable interpolation calculation algorithm for gravity discrete data of Shaanxi area.

Kriging method；variation function model；Shaanxi area；gravity field；interpolation

2016-10-28

2016年度中國(guó)地震局“三結(jié)合”課題(162704)；陜西省起航與創(chuàng)新基金資助項(xiàng)目(201514)

張永奇(1985-)，男，工程師,博士研究生.

著錄：張永奇,韓美濤,韓曉飛，等.Kriging法在區(qū)域重力場(chǎng)插值中的適用性研究[J].測(cè)繪工程，2018，27(2)：1-6.

10.19349/j.cnki.issn1006-7949.2018.02.001

P207

A

1006-7949(2018)02-0001-06

張德福]