翻阅论文资料后，气溶胶光学厚度AOD遥感反演算法总结如下：
基于GF-4卫星影像蓝光波段以及MOD09A1地表反射率数据，采用深蓝算法，利用6S模型建立查找表，反演气溶胶光学厚度。
核心原理如下公式：
式中：
是表观反射率，
是气溶胶反射率， 
是瑞利散射， 
是地表反射率
S是大气半球反照率，
是大气下行透过率，
是大气上行透过率
是卫星影像辐射定标后的结果，采用MOD09A1地表反射率产品，、、S、、由6S大气模型计算得到。
6S模型输出结果：
详细步骤如下：
1、对GF-4卫星影像进行预处理，包括辐射定标、正射校正、图像配准、镶嵌、裁剪等步骤，计算得到研究区范围内蓝光波段的表观反射率；
2、分离提取云、水体区域，云体暂时还没有较好的自动化提取算法，水体可用NDWI阈值分割；
3、处理GF-4同时期的MOD09A1数据，获取研究区内的地表反射率数据；
4、将影像转投影为等经纬度坐标，根据影像左上角坐标，计算各像元位置的经纬度；
5、根据各像元经纬度、成像日期时间，计算各像元成像几何条件，主要包括太阳天顶角、太阳方位角、卫星天顶角、卫星方位角；
相关计算公式参考如下(查询了很多网站资料，众说纷纭，公式也不尽相同，还没完全搞清楚计算原理，暂且按照下面的参考资料)：
https://www.it610.com/article/1292398285364994048.html
https://www.pveducation.org/zh-hans/pvcdrom/%E5%A4%AA%E9%98%B3%E6%97%B6
6、逐像元获取成像几何条件执行6S模型，效率太慢，又因一定范围内的像元几何条件可以近似看做相同，故每100*100像元块（GF-4影像空间分辨率50m，即5km*5km）执行一次6S模型，几何条件选取100*100像元的均值，气溶胶厚度依次设为0.0005、0.25、0.5、0.75、1.0、1.25、1.5、1.75、2.0，则每个像元块执行9次6S模型，利用得到气溶胶反射率、瑞利散射、大气半球反照率、大气下行透过率、大气上行透过率代入前面的公式，计算得到表观反射率，9个气溶胶值，则对应9个表观反射率；
7、所有像元块计算完成后，得到9个气溶胶厚度对应的9幅表观反射率数据，逐像素遍历，建立表观反射率与气溶胶厚度的线性关系，再以GF-4影像计算所得表观反射率为自变量代入拟合关系，得到对应的气溶胶光学厚度，逐像素遍历建立回归关系速度慢，可利用多进程机制提高运行速度。
8、将云、水覆盖区域赋为无效值；
9、最终，建立Aeronet站点观测数据与反演结果的线性关系，评价反演精度。
反演结果如下：
注：GF系列卫星幅宽较大，需计算各个像素的几何条件(太阳天顶角、太阳方位角、卫星天顶角、卫星方位角)，不能用中心像素的值代替，而5中的参考方法不适用于计算极轨卫星的卫星天顶角、卫星方位角，故暂时选取GF-4静轨卫星数据反演AOD。
参考文献：
王中挺，厉青，王桥，李莘莘，陈良富，周春艳，张丽娟，徐拥军．2012．利用深蓝算法从HJ-1数据反演陆地气溶胶. 遥感学报
马小雨，陈正华，宿鑫，于会泳，贾丹丹，姚焕玫.2020.GF-4增强型地表反射率库支持法的气溶胶光学厚度反演. 遥感学报
田信鹏，孙林，刘强，李秀红. 2018. 北京地区Landsat8 OLI高空间分辨率气溶胶光学厚度反演. 遥感学报