基于任務分解模型的離散數據格網化并行優化

王家潤，謝海峰

(華北計算技術研究所基礎三部，北京 100083)

0 引 言

在國家“核高基”重大專項支持下，自主可控的軟硬件發展迅速。但國產基礎軟硬件技術成熟度不高、技術儲備不足、信息化生態環境不完善、體系化能力弱、整體安全性尚待提高等系列問題，還需要重點解決應用遷移、系統適配、性能優化等問題[1]。文獻[1]針對典型基礎軟硬件，優化參數配置、語句性能、界面渲染、數據管理等；文獻[2]通過綜合運用并發處理、數據緩存、分片處理和分層顯示等技術，對信息的產生、傳輸、處理、查詢及顯示的全過程進行性能優化；文獻[3]研究了嵌入式環境中的國產化GPU JM5400對漢字的顯示及大數據量圖形的顯示優化，但與商用相比，性能差距比較大。

基于CPU多核的OpenMP[4,5]及基于GPU的CUDA[5]的并行優化技術已在各領域中大量應用。文獻[5]對現代處理器特性、串行代碼優化、常用并行編程模型及設計方法等進行了深入分析；文獻[6]研究了OpenMP在Kriging插值中的并行優化；文獻[7]中將CUDA應用在運動目標檢測處理；文獻[8]研究了OpenMP在矢量空間數據拓撲算法中的并行優化；文獻[9]研究了CUDA在電磁可視化中的性能加速；文獻[10]研究了GPU性能優化的主要途徑并提出了一個可視化的性能指導模型。上述研究針對不同的實際問題進行了并行優化，但是處理過程相對復雜、多樣。

本文以離散數據格網化處理中的并行優化研究為例，研究國產應用性能提升技術。借鑒文獻[10]的性能指導模型思想，針對并行編程中的線程任務劃分這一共性難點，建立了一個用于指導任務劃分的線程任務分解通用處理模型。

1 軟硬件結合并行優化策略及處理流程

1.1 軟硬件結合并行優化策略

國產CPU與商用相比雖然還存在差距，但因為普遍采用了多核技術，在一定程度上可彌補主頻不足的問題；國產GPU目前沒有面向桌面的產品，主要使用商用GPU，由于缺乏適配的驅動程序，只能采用開源的驅動程序替代，與商用性能差距較大，目前國產下GPU通用計算還無法實施。因此國產下性能優化應以國產CPU處理器提供的多核處理技術為主，以商用GPU為輔(主要用于對比驗證及技術儲備)。

1.2 軟硬件結合并行優化處理流程

(1)性能分析

通過性能測試工具Intel VTune等找出軟件中的性能瓶頸。

(2)并行處理

并行編程的思想是分治策略，即將大問題劃分為一些小問題，再把這些小問題交給相應的處理單元并行處理。分解后的線程處理任務是否不相關(或相關性很低)，這是并行優化的前提。

(3)對比驗證

采用并行處理的軟件，復雜性及調試難度較高，建議先采用串行編程，再并行優化，串行編程用于驗證并行優化算法的正確性及性能優化對比。

(4)性能調優

根據CPU及GPU處理器硬件的結構特征，對并行優化過程中的細節調優(編譯器選項、數據分塊、內存使用、負載平衡、調度優化、緩存優化等[5,11,12])。

2 離散數據格網化IDW算法及并行優化設計

2.1 離散數據格網化IDW算法

離散數據格網化有廣泛的應用，例如：氣象觀測站的分布是離散不規則的，需要先根據這些觀測站數據來插值出規則網格點上的屬性值(溫度、氣壓等)，基于規則網格，再采用等值線、填色圖等進行可視化等。基于空間相近相似原理的反距離權重IDW插值算法是一種常用的方法。

反距離權重IDW插值算法[13]：

(1)輸入

離散采樣點數據(位置、屬性等)集合；待插值規則網格的范圍、行列數目等。

(2)過程

A：鄰近搜索(Searching)

針對需要構建的規則網格中的每個格點，依據某種搜索規則，從離散采樣點數據集合中搜索與該格點相距較近的若干個離散采樣點，使用這些離散采樣點插值計算該格點的屬性值。目前有多種搜索規則：圓域、四方向、八方向、基于Delaunay三角網一階鄰近點、基于遮蔽程度的搜索等[13-15]，本文采用四方向搜索法。

B：插值計算(Computing)

針對規則網格中的每個格點搜索出的若干個離散采樣點，采用反距離權重IDW插值公式對這若干個離散采樣點屬性數據進行計算，得到該網格格點的屬性值。IDW插值計算公式

其中，權重系數

其中，Z(S0)為S0處的預測值，N為計算過程中使用的鄰近采樣點數量，Zi為第i個采樣點的值，di為插值點S0到各采樣點的距離，2是距離的冪。

(3)輸出

插值計算出的各網格的格點屬性值。根據需要可進一步進行交叉驗證，計算插值誤差等。

由IDW插值算法過程可知，規則格網的格點數目一般比較大、計算密集，每個格點的搜索及插值與鄰近的格點沒有依賴，相互獨立，非常符合并行化處理。下面遵循軟硬件結合的并行優化策略，進行并行優化。

2.2 基于OpenMP的并行優化PIDW算法

OpenMP使用Fork-Join并行編程模型，程序創建一個獨立的串行執行的主線程，執行到OpenMP標識的并行域時，即進入Fork過程：主線程會根據任務的劃分創建一組并行的子線程，然后并行域中的程序代碼在不同的線程中并行執行；Join過程：當主線程創建的并行子線程在并行域中執行完任務，退出并行域時，OpenMP內部同步各子線程，結束各子線程，返回到只有主線程運行狀態。

基于Fork-Join并行編程模型的并行優化PIDW算法：

(1)輸入

待插值網格格點數據、離散采樣點數據、網格范圍及行列數目等。

(2)過程

步驟1 設置OpenMP運行參數：獲取CPU核數，根據核數設置要啟動的子線程數。設置子線程數目與核數相同，可最大程度上利用各核提供的硬件能力，使各核具有較好的負載均衡。

步驟2 進入OpenMP并行域(Fork)：根據子線程數目進行任務劃分，確定每個線程負責的任務。子線程計算網格格點的屬性值(搜索鄰近采樣點，并采用IDW插值公式計算)。各子線程間的計算相對獨立，互不干擾。詳細可參見基于OpenMP的線程任務分解模型部分。

步驟3 退出OpenMP并行域(Join)：OpenMP并行域執行完畢后，即結束并行計算模式，此時網格格點數據插值處理完畢。

(3)輸出

經過插值處理的網格格點數據。

基于OpenMP實現的代碼示例：

//gridData：待插值網格格點數據，待插值網格列

//數、行數：gridNumX，girdNumY

//離散采樣點數據：sampleData

//獲取CPU處理器的核數目

intthreads=omp_get_Num_threads()；

//設置并行啟用threads個子線程，啟用CPU所有//核

omp_set_threads(threads)；

//進入并行域(Fork)，啟用threads個子線程

//各線程相互獨立，并行處理計算任務

#pragmaompparallel//OpenMP并行域標識

{

//獲取當前子線程索引，索引范圍：0～threads-1

intthreadIdx=omp_get_thread_num();

//設置該子線程的循環計算任務：不越界

//(i

//(i+=threads)，不重復且保證所有的線程計算

//覆蓋所有的網格格點計算

for(inti=threadIdx;i

for(intj=0;j

//搜索與插值計算

SearchCompute(i,j,gridData,

gridNumX,gridNumY，sampleData);

}//退出并行域,結束threads個子線程，返回主線

//程(Join)，此時網格gridData的各格點數據已

//經過插值處理，可用等值線跟蹤算法等進一步處//理，可視化顯示等。

2.3 基于OpenMP的線程任務分解模型

OpenMP總共分配有N個子線程(在整個并行處理域內，OpenMP會為每個子線程分配唯一的索引threadIdx，索引號范圍：0～N-1)，需要處理columnNum列、rowNum行的二維網格格點的計算任務(網格格點從左下角(0,0)開始，沿X、Y方向編號)。

在一般情況下，子線程數目遠遠小于網格格點數目，因此，每個子線程需要負責多個格點的計算。基本劃分思想是將網格格點按列(或行)分解到N個組中，每個組由一個子線程負責計算，重點是建立線程索引與網格格點索引的映射關系，下面通過二重循環來完成每個子線程的任務分解：

intthreadIdx=omp_get_thread_num()；

for(inti=threadIdx;i

for(intj=0;j

GridNodeCompute(…);

即針對并行域中的當前線程threadIdx，由該線程負責處理網格中從第i=threadIdx列(即將線程索引threadIdx映射到網格的第i列)開始到columnNum-1列范圍內的各列格點的所有計算(j:0～rowNum)，其中列i的遞增間隔為線程總數N，原因是已分配有N個線程，前0～N-1列已分配到這N個子線程中并行處理，所以前N列無需重復計算，這是線程并行任務劃分的關鍵。i=i+N表示沿X方向繼續搜尋需要該線程計算的列，保證計算不遺漏。i

線程任務分解模型如圖1所示，共使用4個子線程：0～3，將網格中的前4列直接映射到4個線程(圖中左側細線矩形范圍內的4列網格格點)，其中線程0負責網格的列0與列4格點計算(參見圓符號標記)，線程1、2與線程0處理類似，線程3只處理網格列3(參見矩形符號標記)。

圖1 基于OpenMP的線程任務分解模型

線程任務分解模型的計算均衡(負載平衡)：一般情況下，首先每個線程參與計算，對應使用了CPU的各核，對計算任務進行了基于全部核的分解；其次，通過任務分解模型示意圖可看出，各線程之間最多相差網格的一列或一行格點，各線程的計算量相差不太大，整體來說，計算比較均衡，即線程任務分解模型具有較好的負載均衡。

2.4 基于CUDA的并行優化PIDW算法

根據CUDA架構的編程模型(HostDevice)[5]，基于CUDA的并行計算主要分為以下3個過程：將待處理數據由Host端(CPU)讀入Device端(GPU)顯存中；在Device端(GPU)中，進行并行計算：采用多線程并行處理數據，并將處理結果存儲在顯存中；把Device端(GPU)中的處理過的數據讀回Host端(CPU)的內存中。CUDA內核函數運行在GPU上，啟動后CUDA中的每一個線程都將會同時并行地執行內核函數中的代碼。GPU是由多個流多處理器構成，以塊(Block)為基本調度單元，因此對于流處理器較多的GPU，一次可以處理的塊(Block)更多，從而運算速度更快。

使用CUDA編程主要注意：①線程的設置。獲取GPU硬件參數：流處理器及線程塊最大線程個數等參數及處理問題規模等。②線程的索引標識。每個線程都有線程索引threadIdx、所在的塊索引blockIdx、所在塊的維數blockDim、所在的線程網格維數gridDim等CUDA內置的信息(圖2)，通過這些信息可計算出線程在整個線程網格Grid中的索引位置，如在線程網格的X索引方向的索引位置

threadIdx.x+blockIdx.x*blockDim.x

基于CUDA架構編程模型的并行優化PIDW算法：

(1)輸入

待插值網格格點數據、離散采樣點數據、網格范圍及行列數目等。

(2)過程

步驟1 獲取GPU設備的硬件參數：GPU的流處理器(SP)個數、線程塊的最大線程數量等。

步驟2 根據GPU設備參數及網格格點數據維數(columnNum，rowNum)，設置CUDA中線程網格的劃分gridDim；線程塊的劃分blockDim，這樣在CUDA中分配的總線程數：gridDim.x*gridDim.y*；blockDim.x*blockDim.y。

因為不同的GPU硬件的限制及網格數據規模等的限制，網格及塊劃分要根據實際環境調整。gridDim及blockDim的劃分有多種方式，本文采用二維形式的gridDim及blockDim為例。一種常用的劃分如下：gridDim(ceil(columnNum/16),ceil(rowNum/16));blockDim(16,16)，其中ceil取上限整數。

步驟3 在顯存(GPU)和內存(CPU)中為待插值網格格點屬性數據開辟存儲空間，在顯存中為離散采樣點數據開辟存儲空間，將離散采樣點的數據由內存拷貝到顯存。

步驟4 根據CUDA中設置的線程數和任務量進行任務劃分，指定每個線程負責的任務。各線程計算網格格點的屬性值(搜索鄰近采樣點，并采用IDW插值公式計算)。各線程間的計算相對獨立，互不干擾。詳細參見基于CUDA的線程任務分解模型。

步驟5 啟動并行計算模式(即CPU端啟動核函數)，在GPU端由GPU中的多個線程并行計算網格格點的屬性值。

步驟6 結束并行計算模式(即GPU端核函數結束，返回CPU端)，將網格格點數據由顯存(GPU)拷貝到內存(CPU)[11,12]。

(3)輸出

輸出已插值計算的網格格點數據。

CUDA部分實現示例代碼：

//gridData：待插值網格格點數據，待插值網格行//數、列數：gridNumX，girdNumY

//離散采樣點數據：sampleData

//核函數

GridingParalleComputing<<>>()

//核函數類似OpenMP中并行域概念

{

//設置為二維線程網格，當前線程索引CUDA內部

//采用二維(threadIdx.x,threadIdx.y)標識，

//所在的線程塊由(blockIdx.x,blockIdx.y)標

//識

intI=threadIdx.x+blockIdx.x*blockDim.x;

intJ=threadIdx.y+blockIdx.y*blockDim.y;

//該線程的任務分解：不越界(i

//j

//i+=gridDim.x*blockDim.x及

//j+=gridDim.y*blockDim.y,類似OpenMP中的

//計算任務處理

for(inti=I;

i

i+=gridDim.x*blockDim.x)

for(intj=J;

j

j+=gridDim.y*blockDim.y)

//搜索與計算

SearchCompute(i,j,gridData,

gridNumX,gridNumY,

sampleData);

}

2.5 基于CUDA的線程任務分解模型

網格格點數據維數(columnNum，rowNum)，采用二維形式的gridDim及blockDim，CUDA總共分配有threads=gridDim.x*gridDim.y*blockDim.x*blockDim.y個線程。基本劃分思想是將網格格點分解到threads個組中，每個組由一個線程負責計算，重點是建立線程索引與網格格點索引的映射關系。在每個塊內，CUDA會為每個線程分配唯一的索引threadIdx(threadIdx.x,threadIdx.y)，因為采用二維形式的blockDim，所以threadIdx也取二維形式。需要處理columnNum列、rowNum行的網格格點的計算任務(網格格點從左下角(0,0)開始編號)，通過二重循環來完成線程的任務分解：

intcol=threadIdx.x+

blockIdx.x*blockDim.x;

introw=threadIdx.y+

blockIdx.y*blockDim.y;

for(inti=col;

i

i+=gridDim.x*blockDim.x)

for(intj=threadIdx.y;

j

j+=gridDim.y*blockDim.y)

GridNodeCompute(…);

每個線程相對于所在的線程塊，有唯一索引，則在整個線程網格中的索引為：(col,row)=(threadIdx.x+blo-ckIdx.x*blockDim.x，threadIdx.y+blockIdx.y*blockDim.y);

即針對塊中的線程threadIdx，由該線程負責處理對應的網格數據中的col列、row行索引格點的計算(如果col列、row行格點存在)。從第i=col列、第j=row行開始到columnNum-1列、rowNum-1行范圍內的各格點的所有計算，其中列i的遞增間隔為gridDim.x*blockDim.x，其中行j的遞增間隔為gridDim.y*blockDim.y，原因是已分配有threads個線程，其中：x方向間隔gridDim.x*blockDim.x；y方向間隔gridDim.y*blockDim.y。網格數據的前0～gridDim.x*blockDim.x-1列及0～gridDim.y*blockDim.y-1行已映射到threads個線程中并行處理，所以無需重復計算，這是線程并行劃分的關鍵。i+=gridDim.x*blockDim.x與j+=gridDim.y*blockDim.y表示繼續沿x、y兩個方向搜索該線程負責計算的格點，主要保證計算不遺漏。i

線程任務分解模型如圖3所示，其中gridDim(3,2),blockDim(2,2),共分配線程個數：24個，直接映射到網格格點(圖中左下方細線矩形范圍內的24個網格格點)，線程網格索引(0,0)(參見菱形符號標記)，所在的塊索引block(0,0)及該塊中線程索引thread(0,0),分配的計算任務是計算網格格點(0,0)、(6，0)、(0,4)、(6,4)；線程網格索引(5,2)(參見三角符號標記)，所在的塊索引block(2,1)及該塊中線程索引thread(1,0)，分配的計算任務是計算網格格點(5,2)、(5，6)。

圖3 基于CUDA的線程任務分解模型

線程任務分解模型的計算均衡(負載平衡)：一般情況下，首先GPU分配的所有線程參與計算，對計算任務進行了基于全部線程的分解；其次，通過任務分解模型示意圖可看出，各線程處理的格點在二維網格上均勻間隔分布，整體來說，計算比較均衡，即線程任務分解模型具有較好的負載均衡。

3 算法實驗

采用傳統的串行算法(CPU單核)、并行算法(CPU多核OpenMP)、并行算法(GPU多核CUDA)3種方式，統計總消耗時間、加速比(串行算法時間/并行算法時間之比)等。

實驗數據選自全球氣象觀測709個站點，地理范圍：經度-180°～180°，維度-90°～90°，數據網格以360×180為基準(即經緯度按1°劃分)，采用網格系數表達網格疏密程度，例如網格系數為2，表示360×2×180×2。

3.1 實驗一、OpenMP并行優化及分析

(1)實驗環境：國產麒麟Kylin操作系統(32位)，CPU：申威(SW410)、龍芯(Loongson3A2000)、英特爾(Intel I5)，均為四核。

(2)實驗結果：統計消耗時間及加速比，并繪制成曲線圖，結果參見表1、表2、圖4～圖7。

表1 網格系數2時多核并行加速性能測試統計

表2 網格系數4時多核并行加速性能測試統計

圖4 網格系數2時CPU消耗時間對比

圖5 網格系數2時CPU加速比曲線

圖6 網格系數4時CPU消耗時間對比

圖7 網格系數4時CPU加速比曲線

(3)實驗分析：選取3種CPU及兩種網格劃分，對IDW算法并行實驗，從表格及曲線圖可看出：

加速比：經過并行加速后，加速比達到3.75～5.614。龍芯(最大加速比3.82)及申威CPU(最大加速比3.75)接近4倍，這二者基本相當，英特爾CPU超過4倍(最大加速比5.614)，加速比最高。

總消耗時間：經過并行加速后，消耗時間大量減少。龍芯與申威二者基本相當，英特爾消耗最少：6.883。

基本結論：采用OpenMP并行優化IDW算法是有效的，整體加速比基本符合與CPU的核數相當這一經驗值，也說明了OpenMP加速受CPU核心數的限制。從加速比可看出，國產與商用CPU存在差距。

3.2 實驗二、CUDA并行優化及分析

(1)實驗環境：CPU：Intel I2，雙核；GPU：NVIDIA Quadro 240；操作系統：windows XP 32位；CUDA Cores：16個；CUDA v6.5。

(2)實驗結果：統計消耗時間及加速比，并繪制成曲線圖，結果參見表3、圖8、圖9。

(3)實驗分析：通過采用不同疏密網格劃分，在Windows操作系統下對IDW算法并行優化實驗，從表格及曲線圖可看出：

加速比：加速比在1.9～5.9之間，性能提升顯著。隨著網格數據量的增加，加速比更高，可看出GPU比較適合大規模的數據并行計算。

表3 GPU并行加速性能測試統計

圖8 GPU加速前后消耗時間對比

圖9 GPU加速比曲線

總消耗時間：經過GPU并行加速后，優化后消耗時間比優化前消耗的時間大量減少。針對網格系數4，比較表2與表3消耗時間(Intel I5：34.352， GPU：20.187)，可看出CUDA比OpenMP更有效。

基本結論：采用CUDA并行優化IDW算法更加有效。

3.3 結論及應用效果

基于OpenMP與CUDA的并行優化，具有較好的性能提升。將本文并行優化技術應用到國產下某型號氣象業務系統項目中，該應用系統性能得到較大的改善，效果如圖10所示。

圖10 局部區域等值線顯示效果

4 結束語

針對傳統的IDW算法采用并行優化處理后，實現了較好的性能提升，體現了CPUGPU多核的硬件加速能力，驗證了軟硬件結合并行優化思路的正確性，線程任務分解模型也極大地減輕了并行編程任務劃分的難度，為國產軟件的性能提升提供了一定的指導。下一步可研究：①目前CUDA僅支持NVIDIA顯卡，支持多種顯卡的OpenCL通用計算技術；②基于Hadoop MapReduce、Spark等分布式大規模數據的并行計算技術；③針對OpenMP及CUDA的性能優化。

參考文獻：

[1]ZHANG Xiaoqing,GONG Bo,TIAN Liyun,et al.Study on the performance optimization of application running on the platform built with Chinese owned technology[J].Software,2015,36(2):5-9(in Chinese).[張曉清,龔波,田麗韞,等.國產自主可控應用性能優化研究[J].軟件，2015,36(2):5-9.]

[2]XIE Zhouyu,ZANG Fei.Performance optimization technology for information system based on domestic software and hardware[J].Command Information System and Technology,2014,5(3):59-63(in Chinese).[謝宙宇,臧飛.基于國產軟硬件的信息系統性能優化技術[J].指揮信息系統與技術,2014,5(3):59-63.]

[3]FU He,XIE Yongfang.From design of a cockpit display system based on homemade GPU JM5400[J].Computer Engineering & Science,2016,38(10):2083-2090(in Chinese).[符鶴,謝永芳.基于國產化圖形芯片JM5400的座艙顯示系統設計[J].計算機工程與科學,2016,38(10):2083-2090.]

[4]OpenMPI[EB/OL].[2017-04-07].http://www.open-mpi.org.

[5]LIU Wenzhi.Parallel computing and performance optimization[M].Beijing:China Machine Press,2015:2-174(in Chinese).[劉文志.并行算法設計與性能優化[M].北京:機械工業出版社,2015:2-174.]

[6]CHEN Huan,XIE Jian.Kriging interpolation algorithm based on OpenMP[J].Computer Science,2012,39(6A):392-394(in Chinese).[陳歡,謝健.基于OpenMP的Kriging插值算法研究[J].計算機科學,2012,39(6A):392-394.]

[7]LING Bin.DENG Yan.YU Shibo.Processing for accelerated sparse PCNN moving target detection algorithm with CUDA[J].Computer Engineering and Design,2016,37(12):3301-3305(in Chinese).[凌濱,鄧艷,于士博.CUDA并行加速的稀疏PCNN運動目標檢測算法[J].計算機工程與設計,2016,37(12):3301-3305.]

[8]GU Yuhang,ZHAO Wei,LI Li,et al.Design and realization of parallel topology algorithm of vector spatial data based on OpenMP[J].Engineering of Surveying and Mapping,2015,24(11):22-27(in Chinese).[谷宇航,趙偉,李力,等.基于OpenMP的矢量空間數據并行拓撲算法設計與實現[J].測繪工程,2015,24(11):22-27.]

[9]WU Lingda,HAO Liyun,FENG Xiaomeng,et al.Combining isosurface rendering and volume rendering for method of electromagnetic environment visualization[J].Journal of Beijing University of Aeronautics and Astronautics,2017,43(5):887-893(in Chinese).[吳玲達,郝利云,馮曉萌,等.結合等值面繪制與體繪制的電磁環境可視化方法[J].北京航空航天大學學報,2017,43(5):887-893.]

[10]JIA Haipeng.Research of parallel optimization technologies on GPU computing platforms[D].Qingdao:Ocean University of China,2012:1-105(in Chinese).[賈海鵬.面向GPU計算平臺的若干并行優化關鍵技術研究[D].青島:中國海洋大學,2012:1-105.]

[11]DAI Chen,CHEN Peng,YANG Donglei,et al.On parallel programming and optimization for multi-core[J].Computer Application and Software,2013,30(12):198-202(in Chinese).[戴晨,陳鵬,楊冬蕾,等.面向多核的并行編程和優化研究[J].計算機應用與軟件,2013,30(12):198-202.]

[12]SHAO Tianjiao.Research on GPU-based parallel computation performance optimization[D].Changchun:Jilin University,2014:1-48(in Chinese).[邵天驕.基于GPU并行計算的性能優化研究[D].長春:吉林大學,2014:1-48.]

[13]ZHOU Hongwen,DU Yuxing,REN Gaofeng,et al.Analysis of ore reserves based on Arc-Engine and IDW[J].Journal of Wuhan University of Technology,2014,36(2):115-119(in Chinese).[周洪文,杜玉星,任高峰,等.基于Arc-Engine與IDW法的礦石儲量分析研究[J].武漢理工大學學報,2014,36(2):115-119.]

[14]DUAN Ping,SHENG Yehua,LI Jia,et al.Adaptive IDW interpolation method and its application in the temperature field[J].Geographical Research,2014,33(8):1417-1426(in Chinese).[段平,盛業華,李佳,等.自適應的IDW插值方法及其在氣溫場中的應用[J].地理研究,2014,33(8):1417-1426.]

[15]LI Zhengquan,WU Yaoxiang.Inverse distance weighted interpolation involving position shading[J].Acta Geodaetica et Cartographica Sinica,2015,44(1):91-98 (in Chinese).[李正泉,吳堯祥.顧及方向遮蔽性的反距離權重插值法[J].測繪學報,2015,44(1):91-98.]