粒子群局部優化的反距離權重插值算法

向 峰，李中志*，熊 熙，李斌勇

（1.成都信息工程大學 網絡空間安全學院，成都 610255；2.先進密碼技術與系統安全四川省重點實驗室（成都信息工程大學），成都 610225）

0 引言

目前，在各處布置自動監測設備是獲取環境中污染物或其他研究成分實時數據，并繪制監測指標在整個研究區域的空間分布情況、分析其時空演變過程的重要途徑。由于這些自動化監測設備數量眾多，而且離散地分布在研究區域中，不便于及時維護，部分監測設備會受各種環境因素的影響導致數據上傳不及時或丟失，不利于監測任務的進行。此外，由于還受監測成本的限制，很難通過大面積布置監測設備來獲取空間中任意位置的真實數據，但可以對這些丟失數據或未布置監測的站點處的值進行估計，空間插值算法可利用在空間中已獲取到的真實數據，根據待估計點與當前已知點之間的空間相關關系對未知區域的取值進行估算插值，以滿足對損壞節點處的丟失數據進行補齊、得到某一未知點的估計數據以及更進一步研究的需要。

目前已有多種空間插值算法被提出，其中較為經典的有泰森多邊形法、反距離權重（Inverse Distance Weight，IDW）法、克里金（Kriging）以及樣條法等［1］，這些算法在不同的場景中有各自的優勢［2］，并且在氣象領域的降水［3］、氣溫［4］、空氣質量［5］的分布規律，以及水土資源［6］、環境污染物［7］分布的研究領域應用廣泛。

近年來由于神經網絡較強的非線性擬合能力，許多學者也將神經網絡應用于空間插值中，如，胡廣義等［8］使用反向傳播（Back-Propagation，BP）網絡對降雨量進行插值；陳飛香等［9］使用徑向基函數（Radial Basis Function，RBF）網絡預測土壤鉻含量，相較于簡單的BP 網絡，收斂速度更快，結果也更好；杜景林等［10］使用改進的聚類算法確定RBF 的參數；曾金迪等［11］提出了空間自回歸網絡（Spatial Auto-Regressive Neural Network，SARNN），通過空間距離求解權重。然而上述模型只是簡單地擬合了研究區域某一時刻的空間相關關系。Ma 等［12］提出一種整合時空關系的Geo-LSTM（Geographic Long Short-Term Memory neural network）模型，利用Geo 層提取數據集的空間特征，而LSTM 則從序列數據中學習數據間的時間依賴性，最終使該模型同時具有插值及預測能力，但需要足夠多的歷史數據記錄進行訓練。Zhu 等［13］利用條件生成對抗網絡獲取空間分布數據集中的深層空間特征，通過地形數據訓練生成器，最后使用得到的生成器能夠準確地生成空間數據，但對其他類型數據的生成能力還需更多驗證；此外該模型需要一定數量的高精度訓練樣本數據，且計算量偏大。

盡管上述研究使用神經網絡取得了較好的效果，但還是需要一定的學習與訓練以保證準確度，否則所得結果會不穩定。IDW 算法原理簡單，運算速度較快，易于實現，仍具有一定的實用價值。近年來學者不斷改進IDW 以提高準確性，如參考點選擇方式［14-15］、冪指數［16］等。李慧晴等［17］考慮了地形對IDW 的影響，通過增設與地形相關的權重系數提高了精確度，但通用性不佳。空間數據中存在的空間異質性［18］包括某一處的屬性與鄰近區域不同（局部異質），以及研究區域中有多個屬性互不相同的子區域（空間分層異質）兩種情況。文獻［19］中使 用K 近鄰算法（K Nearest Neighbors，KNN）判斷待插值點的空間屬性，通過降低分層異質參考點的權重以提高最終結果的準確性，但存在誤判的問題。文獻［20］中考慮了要素空間分布的各向異性，得到的參考點權重能更好地反映數據實際的空間分布情況。文獻［21］中使用粒子群優化算法，對參考點個數、冪指數、各向異性的參數進行全局尋優，提出多參數協同優化反距離權重插值（multi-Parameter co-optimization Inverse Distance Weighting，PIDW）算法，具有較好的普適性，但忽略了局部特性，因此在部分場景反而會降低準確度。

對于存在空間異質性的研究區域，使用固定的全局參數無法反映局部情況，會造成一定的局部誤差，而針對區域局部特性設置相關系數的改進策略，雖然能夠提高準確度但也會耗費一定的人力。因此，本文改進PIDW 算法，提出粒子群算法局部優化的反距離權重（Particle swarm Local optimization Inverse Distance Weight，PLIDW）算法，并使用模擬數據與氣象要素數據進行對比驗證。

1 算法概述

1.1 各向異性的反距離權重模型

1.1.1 反距離權重法

反距離權重法假設觀測點對于插值點的影響隨距離的增加而減弱，即成反比關系，表達式如下：

其中：Z*為預測值；n為選取的參照采樣點數目；λi為觀測點對估計插值點的權重系數；Zi為第i個檢測點所測的實際值。權重系數由式（2）給出：

其中：di為待插值點與第i個觀測點的距離；b為距離衰減系數。

1.1.2 各向異性模型

文獻［21］引入“橢圓”距離，可以體現出參考點選取時，不同方向的差異：

對于坐標系中的兩點（xi，yi）、（xj，yj），參數λ調節了二者在Y軸方向上的距離：λ>1 時增加Y軸方向的距離；λ<1 時則會使X軸方向的距離相對增加；λ=1 即退化為歐氏距離。

式（3）只調節了坐標軸方向上的距離，將坐標按式（4）變換，即可反映各向異性：

由式（3）、（4）最終可得到橢圓距離如下：

其中：代表調整后的橢圓距離；Δx與Δy為原數據間的歐氏距離；θ∈（0，2π）為旋轉角度。

1.2 局部參數尋優反距離權重法

考慮到空間異質性的存在，本文通過粒子群優化算法，使用交叉驗證的策略，分別對已知的參考點以式（2）中的參數n、b與式（5）中的參數θ、λ進行尋優并保留局部優化結果。

1.2.1 粒子群優化算法

粒子群優化算法模仿了鳥群的覓食行為，通過信息共享獲取當前群體最優值［22］，并經過迭代得到最終的全局最優值，有全局性能好、簡單易實現的優點［23］。

算法運行時，粒子的位置向量p根據自身經過的最佳位置pbest和種群最佳位置gbest進行更新：

其中：ω為慣性系數；c1和c2為學習因子，反映粒子的學習能力；k代表迭代次數；xid代表當前d維的位置值；pid與pgd分別代表粒子自身第d維的最佳位置與群體最佳位置值，φ1與φ2是［0，1］區間內均勻分布的隨機數。

1.2.2 優化策略

已知研究區域共有N個樣本點｛Z1，Z2，…，ZN｝，對其中任意的Zi利用粒子群優化算法對參數P=（n，b，θ，λ）進行尋優，得到Zi的一組最佳參數Pbesti，如式（7）所示：

其中：代表估計時不包含點Zi本身；Pk為某一組參數向量，對應了一個粒子的坐標；得到的Pbesti即為在點Zi處尋找到的一組最優參數，需要將它進行存儲。

實際使用環境不一定處處都存在各向異性的情況，為了減少多參數尋優帶來的不穩定，本文在對每個站點的參數尋優過程中使用兩個粒子種群：一個種群對P=（n，b，θ，λ）進行完整的尋優，而另一個種群固定參數θ=0 和λ=1，即只對參數n、b進行尋優，最后只保留最優種群的全局最優結果。

1.2.3 算法流程

本文算法主要包含參數尋優過程以及插值使用過程，完成對所有已知樣本點的最佳參數尋優后，得到相應的Pbesti（i=1，2，…，N）。假設空間距離越近，最優參數的取值也越相似，因此對于未知點的取值，算法將使用最近樣本點的最優參數進行近似優化估計。

算法的整體流程如圖1 所示。由于各個參考點的參數尋優過程相互獨立，因此可進行并發優化處理，提高程序的運算速度。此外，由于算法需要反復地計算節點之間的距離，找到最近的參考點，為了提高效率，本文使用K 維樹（KDimensional Tree，KD-Tree）優化查詢過程，可減少選擇最近參考點時的計算時間［24］。

圖1 本文算法流程Fig.1 Flowchart of proposed algorithm

2 仿真實驗與結果分析

為了評估算法的效果，本文首先使用仿真環境進行對比實驗。將本文算法與經典的IDW 算法、Kriging、PIDW、SARNN 進行對比。

2.1 評價方法

本文通過計算預測值與真實值之間的平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）、均方根誤差（Root Mean Square Error，RMSE）以及平均絕對誤差（Mean Absolute Error，MAE）并統計誤差的最大值（MAXimum，MAX），對算法的效果進行評估。其中MAPE、RMSE 與MAE的數學表達式如式（8）～（10）所示，數值越小越準確。

2.2 仿真數據

本文對曲面進行插值測試，曲面的函數方程如下：

插值實驗區間為（0，1），曲面的圖像如圖2 所示。

圖2 測試曲面函數圖像Fig.2 Image of test surface function

2.3 實驗方法

實驗隨機生成50 組坐標點作為參考樣本點，另外生成100 組坐標點進行效果評估，選擇較少數量的參考點更能反映出不同算法間的差異。部分算法的參數設置與實現如下。

1）IDW。經過多次實驗，設置參考最近7 個站點，距離衰減系數為2。

2）SARNN。輸入層傳遞函數選擇Softmax，隱含層激活函數選擇ReLU（Rectified Linear Unit），損失函數選擇RMSE，優化器選擇隨機梯度下降（Stochastic Gradient Descent，SGD），輸入輸出層神經元數為n，隱含層數為2，神經元數目均為2n，n為參考的樣本點數目。使用Dropout 與Early stop防止神經網絡過度擬合［25］，神經網絡由PyTorch 實現。

3）PIDW。粒子群算法需要尋優的參數為：n、b、θ、λ。經過多次實驗測試，本文參照文獻［15，17，19-21］的應用情況及建議，選取n∈[4，10]，該范圍能夠應對多數場景；b∈[1，4]；θ∈[0，2π]，無需更改；λ∈[0.01，100]，該參數范圍可以根據實際情況調整。

4）PLIDW。參數n、b、θ、λ的設置與PIDW 一致；另一種群固定設置θ=0 和λ=1，其他也與PIDW 的參數保持一致。

2.4 仿真實驗結果及討論

仿真實驗結果如表1 所示。預測值（Predictive Value，PV）與真實值（Actual Value，AV）的對比如圖3 所示，數據點越貼近軸線PV=AV 說明越準確。

表1 仿真實驗結果 單位：%Tab.1 Results of simulation experiments unit：%

如表1 所示，PLIDW 在RMSE、MAE 上表現最佳，相較于PIDW 與IDW，PLIDW 的四項指標都更低，準確度更高。SARNN 的MAPE 最低，這與它的MAX 最小有關，如圖3（c）所示，在AV（0，0.2）的區間內，SARNN 的預測點偏離PV=AV 軸線的程度最低，但總體偏差更大，這可以從MAE 值得到反映，它的MAE 為7.572%，高于其他對比算法。Kriging 也出現了相同問題，總體上也出現偏差，MAE 值僅次于SARNN。由此可見，相較于對比算法，PLIDW 具有明顯優勢。

圖3 各算法模擬數據散點圖Fig.3 Simulated data scatter plot of different algorithms

3 氣象要素數據實驗與結果分析

為了得到更加可靠的結論，本文以四川省2021 年12 月某日整點的氣溫為例，在真實場景中對各算法進行比較，所有的數據來源于中國氣象數據網，數據集中共有145 個自動化監測站點的監測數據。

本文首先以不放回的形式從所有的站點中隨機抽取20個站點用作效果評估，再從剩余的站點中抽取站點作為參考點。第一組實驗抽取20 個，后續實驗以20 的倍數遞增抽取，以此對數個算法進行多次對比實驗，相關參數的設置與2.3節保持一致。

為了與SARNN 進行對比，需要將數據進行歸一化，歸一化方式如式（12）所示：

3.1 評價指標

由于本文選用的評價指標較多，為了更加直觀地對各個算法進行比較，將所有不同參考點數量的實驗結果的四個評估指標按照式（13）計算綜合評分：

式（13）綜合考慮了4 個評估指標，得分越高者綜合表現更佳，本文以此衡量算法的準確度。

3.2 實驗結果及討論

將各個算法插值的結果通過式（8）～（10）計算誤差，并統計誤差的最大值IMAX。將所有實驗的結果按照式（13）計算，得到的相關數據如表2 所示。PLIDW 的各項指標表現相較于對比算法更平穩，各個評測指標處于最優或偏上的水平。在6 組不同參考點數的測試中，PLIDW 有4 組的綜合表現優于其他對比算法，雖然隨著參考點數目的增加，SARNN 因為神經網絡的擬合能力使優勢得以顯現，但收斂也因此減慢。而PIDW 由于未考慮區域中存在的局部特性，在少數場景中出現弱于經典IDW 的情況。總體上看，PLIDW 算法的表現明顯優于IDW 與PIDW 算法，相較于PIDW 算法，PLIDW 算法的準確度提高了4.18%～40.08%。

表2 氣象要素數據的實驗結果Tab.2 Experimental results of meteorological element data

綜上所述，本文提出的PLIDW 算法考慮了局部參數選取對精確度的影響，并使用兩個粒子種群進行參數尋優，緩解多參數尋優的不穩定性對結果產生的影響，提高了算法的適應能力，在多數的實驗對比場景中占據優勢，不需要其他的輔助數據，易于實現和使用。

4 運行效率分析

由于比傳統算法多了粒子群參數尋優的過程，因此增加了本文算法的時間復雜度。為在精度與耗時之間進行平衡，首先需對本文算法的時間復雜度進行分析。假設參考點總數為N，選擇n個樣本點參與計算，尋優區間為［a，b］，a≤n≤b，迭代次數為K，種群規模為M。算法首先遍歷N個參考點，分別找到其中距離最近的b個點，復雜度為O（Nbf（N）），f（x）表示某種尋找最近點的算法復雜度。設對一個粒子進行評估耗時為O′(1)，可知對一個參考點的優化耗時為O′(KM)，整個優化過程的時間復雜度為O′(KMN)，最后還需尋找一次最近鄰點再計算，與IDW 算法的時間復雜度相同，為O′(nf(N))。最終可得整體的時間復雜度為：O(Nbf(N)) +O′(KMN) +O′(nf(N))，實際應用場景中KM的量級與N較為接近，對計算時間影響較大，可針對它進行調整以減少耗時。

因為對K、M的取值作手動觀察調整較為困難，本文根據文獻［26］的建議，固定設置種群規模M為30，以在第3 章的數據中隨機抽取120 個采樣點為例，逐漸提高K的取值并計算準確度（式（13））得到的結果如圖4 所示。

圖4 準確度結果Fig.4 Results of accuracy

如圖4 可以看出，迭代次數增加到一定程度時，算法的準確度很難再有明顯提升，計算耗時反而成倍增加，因此迭代20 次左右即可。雖然本文算法的時間復雜度高于傳統算法，但由于實際使用時參考點的總數通常在幾十到數百之間，因此仍舊可以接受，在使用時逐漸增加迭代次數進行嘗試，并盡量使量級不超過待尋優的參考點的數目N，可減少計算耗時。

5 結語

本文提出的PLIDW 算法提高了IDW 算法的準確度，相較于全局尋優的PIDW 算法，主要具有以下的特點：

1）考慮了研究區域內的局部特性，相較于使用不變的全局參數，能夠減小局部誤差。

2）算法同時還考慮了對各向同性模型的參數優化，與各向同性模型尋優結果的適應度進行比較，保存最佳的模型參數，減少尋優難度。

3）在多參數優化的過程中，各站點之間相互獨立，易于設計并發程序以提升運行效率。

由于加入了粒子群算法，本文算法的復雜度有所提高，探索并使用更簡單快速的優化算法可進一步提高運行效率。此外，本文算法對未知點插值參數的選取方式還可進一步改進，本文直接根據最近距離進行選擇，而綜合衡量周邊多個節點的情況進行推斷可能會更加準確。