通过对蒲江县2021物流发展现状进行分析，运用数据爬取，使用高德地图API实现批量查询各小区中每一个地址（省、区/市县、乡镇）的经纬度数据，运用Geopy库中的默认WGS-84模型计算出误差在0.5%以内的小区距离，并以矩阵的形式导入Excel文件中，最后运用所获数据结合双约束重力模型，用编程迭代实现了蒲江县2025、2030、2035年的物流分布量预测。2021小区矩阵如下：

在双约束重力模型中物流分布量的计算公式如下：

确定阻抗矩阵

       考虑蒲江与各小区的运输费用难以查询，并与各地区流通的交通工具并不相似，故费用与时间均不能作为阻抗矩阵的依据。故此处，以各小区的距离为阻抗矩阵。因每个小区的区域广阔，故取小区间的几何中心作为参考点，以各小区几何中心点间的距离为各小区之间的距离。

        小区间距离获取的方法如下：首先在高德开发平台上申请用户key，使用高德地图API实现批量查询各小区中每一个地址（省、区/市县、乡镇）的经纬度，整理表格，用分类汇总方法计算出每一个L-OD小区的平均经纬度。

        将小区经纬度数据导入python中，用Geopy库中的默认WGS-84模型计算出误差在0.5%以内的小区距离，并以矩阵的形式导入excel表中。

物流分布量预测

        将2025年、2030年、2035年的预测物流发生量、吸引量以及各小区之间的距离阻抗矩阵数据导入Python中，用编程的方法进行重力模型系数迭代，由于所给数据显示小区5到小区16之间没有物流交换活动，因此进行重力模型系数迭代时不应考虑5-16小区。

        将误差率err设置为0.000001，结果显示2025年预测迭代了7344次，2030年预测迭代了9019次，2035年预测迭代了10095次。将分布量预测结果导入xls文件中。应注意由于小区5到小区16之间没有物流交换，重力模型系数迭代时也未考虑5-16小区，因此在物流分布量预测表中也应人为将5-16小区之间的分布量设置为0，物流分布量迭代实现代码如下：

import copy
import numpy as np
import pandas as pd
import xlwt  
import xlrd

#读取发生量或吸引量
df=pd.read_excel(r"D:\预测2025.xlsx")
dis=pd.read_excel('D:\距离.xlsx',header=None)  #读取距离矩阵
A=df['A']  #读取第一列
P=df['P']  #读取第二列

#获取行列数
row=len(dis)
col=dis.columns.size
cut=4

k2=np.ones(16)  #初始化重力模型系数
k1=np.ones(16)   
err=0.000001  #设定误差率
x=0  #控制循环条件
ite=0  #迭代次数
mat1=np.zeros(16)  
mat2=np.ones(16)  #初始化迭代矩阵

while x==0 :
    for i in range(cut):
        a=0
        for j in range(col):
            if dis[i][j]!=0:
                a=a+k2[j]*A[j]/dis[i][j]
        a=1/a
        k1[i]=a 
    for i in range(cut,row):
        for j in range(cut):
            if dis[i][j]!=0:
                a=a+k2[j]*A[j]/dis[i][j]
        a=1/a
        k1[i]=a 

    for j in range(cut):
        a=0
        for i in range(row):
            if dis[i][j]!=0:
                a+=k1[i]*P[i]/dis[i][j]
        a=1/a
        k2[j]=a 
    for j in range(cut,col):
        for i in range(cut):
            if dis[i][j]!=0:
                a=a+k1[i]*P[i]/dis[i][j]
        a=1/a
        k2[j]=a 
    
    ite=ite+1
    mat2=np.row_stack((mat2,k2))  
    mat1=np.row_stack((mat1,k1))  #更新矩阵
    
    #判断是否继续迭代
    for j in range(col):
        if abs(mat1[ite][j]-mat1[ite-1][j])>=err:
            x1=0
            break
        else:
            x1=1
    for j in range(col):
        if abs(mat2[ite][j]-mat2[ite-1][j])>=err:
            x2=0
            break
        else:
            x2=1

    if x2==1&x1==1 :
        x=1
    else:
        x=0

print("mat1矩阵：",mat1)
print("mat2矩阵：",mat2)
print("ite=",ite)
 
#计算分布量
l=np.zeros((16,16))
for i in range(row):
    for j in range(col):
        l[i][j]=k1[i]*k2[j]*P[i]*A[j]/dis[i][j]
for i in range(cut,row):
    for j in range(cut,col):
        l[i][j]=0
filename =xlwt.Workbook() #创建工作簿
sheet1 = filename.add_sheet(u'sheet1',cell_overwrite_ok=True) #创建sheet
for i in range (0,16):
    for j in range (0,16):
        sheet1.write(i,j,l[i,j])
filename.save('D:\物流分布量2.xls')

代码运行结果如下：

将预测分布量导入excel，以2025年为例，预测结果整理如下：